/[escript]/trunk/paso/src/PCG.c
ViewVC logotype

Diff of /trunk/paso/src/PCG.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 415 by gross, Wed Jan 4 05:37:33 2006 UTC revision 495 by gross, Mon Feb 6 06:32:06 2006 UTC
# Line 1  Line 1 
1  /* $Id$ */  /* $Id$ */
2    
3    
4  /* PCG iterations */  /* PCG iterations */
5    
6  #include "SystemMatrix.h"  #include "SystemMatrix.h"
# Line 62  err_t Paso_Solver_PCG( Line 63  err_t Paso_Solver_PCG(
63    dim_t n = A->num_cols * A-> col_block_size;    dim_t n = A->num_cols * A-> col_block_size;
64    double *resid = tolerance, *rs=NULL, *p=NULL, *v=NULL, *x2=NULL ;    double *resid = tolerance, *rs=NULL, *p=NULL, *v=NULL, *x2=NULL ;
65    double tau_old,tau,beta,delta,gamma_1,gamma_2,alpha,sum_1,sum_2,sum_3,sum_4,sum_5,tol;    double tau_old,tau,beta,delta,gamma_1,gamma_2,alpha,sum_1,sum_2,sum_3,sum_4,sum_5,tol;
66    double d,norm_of_residual,norm_of_residual_global;    double norm_of_residual,norm_of_residual_global;
67      register double r_tmp,d,rs_tmp,x2_tmp,x_tmp;
68    
69  /*                                                                 */  /*                                                                 */
70  /*-----------------------------------------------------------------*/  /*-----------------------------------------------------------------*/
# Line 93  err_t Paso_Solver_PCG( Line 95  err_t Paso_Solver_PCG(
95         for (i0=0;i0<n;i0++) {         for (i0=0;i0<n;i0++) {
96            rs[i0]=r[i0];            rs[i0]=r[i0];
97            x2[i0]=x[i0];            x2[i0]=x[i0];
98           }
99           #pragma omp for private(i0) schedule(static)
100           for (i0=0;i0<n;i0++) {
101            p[i0]=0;            p[i0]=0;
102            v[i0]=0;            v[i0]=0;
103         }         }
# Line 136  err_t Paso_Solver_PCG( Line 141  err_t Paso_Solver_PCG(
141             if (! (breakFlag = (ABS(delta) <= TOLERANCE_FOR_SCALARS))) {             if (! (breakFlag = (ABS(delta) <= TOLERANCE_FOR_SCALARS))) {
142                 alpha=tau/delta;                 alpha=tau/delta;
143                 /* smoother */                 /* smoother */
144                   #pragma omp for private(i0) schedule(static)
145                   for (i0=0;i0<n;i0++) r[i0]-=alpha*v[i0];
146                 #pragma omp for private(i0,d) reduction(+:sum_3,sum_4) schedule(static)                 #pragma omp for private(i0,d) reduction(+:sum_3,sum_4) schedule(static)
147                 for (i0=0;i0<n;i0++) {                 for (i0=0;i0<n;i0++) {
                      r[i0]-=alpha*v[i0];  
148                       d=r[i0]-rs[i0];                       d=r[i0]-rs[i0];
149                       sum_3=sum_3+d*d;                       sum_3+=d*d;
150                       sum_4=sum_4+d*rs[i0];                       sum_4+=d*rs[i0];
151                  }                  }
152                  gamma_1= ( (ABS(sum_3)<= ZERO) ? 0 : -sum_4/sum_3) ;                  gamma_1= ( (ABS(sum_3)<= ZERO) ? 0 : -sum_4/sum_3) ;
153                  gamma_2= ONE-gamma_1;                  gamma_2= ONE-gamma_1;
154                  #pragma omp for private(i0) reduction(+:sum_5) schedule(static)                  #pragma omp for private(i0,x2_tmp,x_tmp,rs_tmp) schedule(static)
155                  for (i0=0;i0<n;i0++) {                  for (i0=0;i0<n;++i0) {
156                      rs[i0]=gamma_2*rs[i0]+gamma_1*r[i0];
157                    x2[i0]+=alpha*p[i0];                    x2[i0]+=alpha*p[i0];
158                    x[i0]=gamma_2*x[i0]+gamma_1*x2[i0];                    x[i0]=gamma_2*x[i0]+gamma_1*x2[i0];
159                    rs[i0]=gamma_2*rs[i0]+gamma_1*r[i0];                  }
160                    sum_5+=rs[i0]*rs[i0];                  #pragma omp for private(i0) reduction(+:sum_5) schedule(static)
161                }                  for (i0=0;i0<n;++i0) sum_5+=rs[i0]*rs[i0];
162                norm_of_residual=sqrt(sum_5);                  norm_of_residual=sqrt(sum_5);
163                convergeFlag = norm_of_residual <= tol;                  convergeFlag = norm_of_residual <= tol;
164                maxIterFlag = num_iter == maxit;                  maxIterFlag = num_iter == maxit;
165                breakFlag = (ABS(tau) <= TOLERANCE_FOR_SCALARS);                  breakFlag = (ABS(tau) <= TOLERANCE_FOR_SCALARS);
166            }             }
167         }         }
168         /* end of iteration */         /* end of iteration */
169         #pragma omp master         #pragma omp master

Legend:
Removed from v.415  
changed lines
  Added in v.495

  ViewVC Help
Powered by ViewVC 1.1.26