/[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 1388 by trankine, Fri Jan 11 07:45:58 2008 UTC revision 1556 by gross, Mon May 12 00:54:58 2008 UTC
# Line 104  err_t Paso_Solver_PCG( Line 104  err_t Paso_Solver_PCG(
104    } else {    } else {
105      maxit = *iter;      maxit = *iter;
106      tol = *resid;      tol = *resid;
107      #pragma omp parallel firstprivate(maxit,tol,convergeFlag,maxIterFlag,breakFlag) \      Performance_startMonitor(pp,PERFORMANCE_SOLVER);
108                                             private(tau_old,tau,beta,delta,gamma_1,gamma_2,alpha,norm_of_residual,num_iter)      /* initialize data */
109        #pragma omp parallel
110      {      {
111         Performance_startMonitor(pp,PERFORMANCE_SOLVER);             #pragma omp for private(i0) schedule(static)
112         /* initialize data */             for (i0=0;i0<n;i0++) {
113         #pragma omp for private(i0) schedule(static)                rs[i0]=r[i0];
114         for (i0=0;i0<n;i0++) {                x2[i0]=x[i0];
115            rs[i0]=r[i0];                p[i0]=0;
116            x2[i0]=x[i0];                v[i0]=0;
117         }             }
118         #pragma omp for private(i0) schedule(static)      }
119         for (i0=0;i0<n;i0++) {      num_iter=0;
120            p[i0]=0;      /* start of iteration */
121            v[i0]=0;      while (!(convergeFlag || maxIterFlag || breakFlag)) {
        }  
        num_iter=0;  
        /* start of iteration */  
        while (!(convergeFlag || maxIterFlag || breakFlag)) {  
122             ++(num_iter);             ++(num_iter);
123             #pragma omp barrier         sum_1 = 0;
124             #pragma omp master         sum_2 = 0;
125             {         sum_3 = 0;
126             sum_1 = 0;         sum_4 = 0;
127             sum_2 = 0;         sum_5 = 0;
            sum_3 = 0;  
            sum_4 = 0;  
            sum_5 = 0;  
            }  
128             /* v=prec(r)  */             /* v=prec(r)  */
129             Performance_stopMonitor(pp,PERFORMANCE_SOLVER);             Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
130             Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);             Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);
# Line 139  err_t Paso_Solver_PCG( Line 132  err_t Paso_Solver_PCG(
132             Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);             Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);
133             Performance_startMonitor(pp,PERFORMANCE_SOLVER);             Performance_startMonitor(pp,PERFORMANCE_SOLVER);
134             /* tau=v*r    */             /* tau=v*r    */
135             #pragma omp for private(i0) reduction(+:sum_1) schedule(static)             #pragma omp parallel for private(i0) reduction(+:sum_1) schedule(static)
136             for (i0=0;i0<n;i0++) sum_1+=v[i0]*r[i0]; /* Limit to local values of v[] and r[] */             for (i0=0;i0<n;i0++) sum_1+=v[i0]*r[i0]; /* Limit to local values of v[] and r[] */
137             #ifdef PASO_MPI             #ifdef PASO_MPI
138              /* In case we have many MPI processes, each of which may have several OMP threads:              /* In case we have many MPI processes, each of which may have several OMP threads:
139                 OMP master participates in an MPI reduction to get global sum_1 */                 OMP master participates in an MPI reduction to get global sum_1 */
140                  #pragma omp master              loc_sum[0] = sum_1;
141              {              MPI_Allreduce(loc_sum, &sum_1, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
               loc_sum[0] = sum_1;  
               MPI_Allreduce(loc_sum, &sum_1, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);  
             }  
142             #endif             #endif
143             tau_old=tau;             tau_old=tau;
144             tau=sum_1;             tau=sum_1;
145             /* p=v+beta*p */             /* p=v+beta*p */
146             if (num_iter==1) {             if (num_iter==1) {
147                 #pragma omp for private(i0)  schedule(static)                 #pragma omp parallel for private(i0)  schedule(static)
148                 for (i0=0;i0<n;i0++) p[i0]=v[i0];                 for (i0=0;i0<n;i0++) p[i0]=v[i0];
149             } else {             } else {
150                 beta=tau/tau_old;                 beta=tau/tau_old;
151                 #pragma omp for private(i0)  schedule(static)                 #pragma omp parallel for private(i0)  schedule(static)
152                 for (i0=0;i0<n;i0++) p[i0]=v[i0]+beta*p[i0];                 for (i0=0;i0<n;i0++) p[i0]=v[i0]+beta*p[i0];
153             }             }
154             /* v=A*p */             /* v=A*p */
155             Performance_stopMonitor(pp,PERFORMANCE_SOLVER);             Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
156             Performance_startMonitor(pp,PERFORMANCE_MVM);             Performance_startMonitor(pp,PERFORMANCE_MVM);
157         Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, p,ZERO,v);         Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, p,ZERO,v);
        Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, p,ZERO,v);  
158             Performance_stopMonitor(pp,PERFORMANCE_MVM);             Performance_stopMonitor(pp,PERFORMANCE_MVM);
159             Performance_startMonitor(pp,PERFORMANCE_SOLVER);             Performance_startMonitor(pp,PERFORMANCE_SOLVER);
160             /* delta=p*v */             /* delta=p*v */
161             #pragma omp for private(i0) reduction(+:sum_2) schedule(static)             #pragma omp parallel for private(i0) reduction(+:sum_2) schedule(static)
162             for (i0=0;i0<n;i0++) sum_2+=v[i0]*p[i0];             for (i0=0;i0<n;i0++) sum_2+=v[i0]*p[i0];
163             #ifdef PASO_MPI             #ifdef PASO_MPI
164                 #pragma omp master            loc_sum[0] = sum_2;
165             {            MPI_Allreduce(loc_sum, &sum_2, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
              loc_sum[0] = sum_2;  
              MPI_Allreduce(loc_sum, &sum_2, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);  
            }  
166             #endif             #endif
167             delta=sum_2;             delta=sum_2;
168    
# Line 184  err_t Paso_Solver_PCG( Line 170  err_t Paso_Solver_PCG(
170             if (! (breakFlag = (ABS(delta) <= TOLERANCE_FOR_SCALARS))) {             if (! (breakFlag = (ABS(delta) <= TOLERANCE_FOR_SCALARS))) {
171                 alpha=tau/delta;                 alpha=tau/delta;
172                 /* smoother */                 /* smoother */
173                 #pragma omp for private(i0) schedule(static)                 #pragma omp parallel
174                 for (i0=0;i0<n;i0++) r[i0]-=alpha*v[i0];                 {
175                 #pragma omp for private(i0,d) reduction(+:sum_3,sum_4) schedule(static)                    #pragma omp for private(i0) schedule(static)
176                 for (i0=0;i0<n;i0++) {                    for (i0=0;i0<n;i0++) r[i0]-=alpha*v[i0];
177                       d=r[i0]-rs[i0];                    #pragma omp for private(i0,d) reduction(+:sum_3,sum_4) schedule(static)
178                       sum_3+=d*d;                    for (i0=0;i0<n;i0++) {
179                       sum_4+=d*rs[i0];                          d=r[i0]-rs[i0];
180                            sum_3+=d*d;
181                            sum_4+=d*rs[i0];
182                      }
183                 }                 }
184                 #ifdef PASO_MPI                 #ifdef PASO_MPI
185                     #pragma omp master                 loc_sum[0] = sum_3;
186                 {                 loc_sum[1] = sum_4;
187                   loc_sum[0] = sum_3;                 MPI_Allreduce(loc_sum, sum, 2, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
188                   loc_sum[1] = sum_4;                 sum_3=sum[0];
189                   MPI_Allreduce(loc_sum, sum, 2, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);                 sum_4=sum[1];
                  sum_3=sum[0];  
                  sum_4=sum[1];  
                }  
190                  #endif                  #endif
191                  gamma_1= ( (ABS(sum_3)<= ZERO) ? 0 : -sum_4/sum_3) ;                  gamma_1= ( (ABS(sum_3)<= ZERO) ? 0 : -sum_4/sum_3) ;
192                  gamma_2= ONE-gamma_1;                  gamma_2= ONE-gamma_1;
193                  #pragma omp for private(i0,x2_tmp,x_tmp,rs_tmp) schedule(static)                  #pragma omp parallel
194                  for (i0=0;i0<n;++i0) {                  {
195                    rs[i0]=gamma_2*rs[i0]+gamma_1*r[i0];                      #pragma omp for private(i0,x2_tmp,x_tmp,rs_tmp) schedule(static)
196                    x2[i0]+=alpha*p[i0];                      for (i0=0;i0<n;++i0) {
197                    x[i0]=gamma_2*x[i0]+gamma_1*x2[i0];                        rs[i0]=gamma_2*rs[i0]+gamma_1*r[i0];
198                          x2[i0]+=alpha*p[i0];
199                          x[i0]=gamma_2*x[i0]+gamma_1*x2[i0];
200                        }
201                        #pragma omp for private(i0) reduction(+:sum_5) schedule(static)
202                        for (i0=0;i0<n;++i0) sum_5+=rs[i0]*rs[i0];
203                  }                  }
                 #pragma omp for private(i0) reduction(+:sum_5) schedule(static)  
                 for (i0=0;i0<n;++i0) sum_5+=rs[i0]*rs[i0];  
204                  #ifdef PASO_MPI                  #ifdef PASO_MPI
205                     #pragma omp master                 loc_sum[0] = sum_5;
206                 {                 MPI_Allreduce(loc_sum, &sum_5, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
                   loc_sum[0] = sum_5;  
                   MPI_Allreduce(loc_sum, &sum_5, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);  
                }  
207                  #endif                  #endif
208                  norm_of_residual=sqrt(sum_5);                  norm_of_residual=sqrt(sum_5);
209                  convergeFlag = norm_of_residual <= tol;                  convergeFlag = norm_of_residual <= tol;
210                  maxIterFlag = num_iter == maxit;                  maxIterFlag = num_iter == maxit;
211                  breakFlag = (ABS(tau) <= TOLERANCE_FOR_SCALARS);                  breakFlag = (ABS(tau) <= TOLERANCE_FOR_SCALARS);
212             }             }
213         }      }
214         /* end of iteration */      /* end of iteration */
215         #pragma omp master      num_iter_global=num_iter;
216         {      norm_of_residual_global=norm_of_residual;
217             num_iter_global=num_iter;      if (maxIterFlag) {
218             norm_of_residual_global=norm_of_residual;           status = SOLVER_MAXITER_REACHED;
219             if (maxIterFlag) {      } else if (breakFlag) {
220                 status = SOLVER_MAXITER_REACHED;           status = SOLVER_BREAKDOWN;
221             } else if (breakFlag) {      }
222                 status = SOLVER_BREAKDOWN;      Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
            }  
        }  
        Performance_stopMonitor(pp,PERFORMANCE_SOLVER);  
     }  /* end of parallel region */  
223      TMPMEMFREE(rs);      TMPMEMFREE(rs);
224      TMPMEMFREE(x2);      TMPMEMFREE(x2);
225      TMPMEMFREE(v);      TMPMEMFREE(v);

Legend:
Removed from v.1388  
changed lines
  Added in v.1556

  ViewVC Help
Powered by ViewVC 1.1.26