/[escript]/branches/arrexp_2137_win/paso/src/PCG.c
ViewVC logotype

Diff of /branches/arrexp_2137_win/paso/src/PCG.c

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

temp/paso/src/PCG.c revision 1387 by trankine, Fri Jan 11 07:45:26 2008 UTC trunk/paso/src/PCG.c revision 1571 by gross, Sat May 24 22:28:33 2008 UTC
# Line 63  Line 63 
63  *  ==============================================================  *  ==============================================================
64  */  */
65    
66    /* #define PASO_DYNAMIC_SCHEDULING_MVM */
67    
68    #if defined PASO_DYNAMIC_SCHEDULING_MVM && defined __OPENMP
69    #define USE_DYNAMIC_SCHEDULING
70    #endif
71    
72  err_t Paso_Solver_PCG(  err_t Paso_Solver_PCG(
73      Paso_SystemMatrix * A,      Paso_SystemMatrix * A,
74      double * r,      double * r,
# Line 71  err_t Paso_Solver_PCG( Line 77  err_t Paso_Solver_PCG(
77      double * tolerance,      double * tolerance,
78      Paso_Performance* pp) {      Paso_Performance* pp) {
79    
   
80    /* Local variables */    /* Local variables */
81    dim_t num_iter=0,maxit,num_iter_global;    dim_t num_iter=0,maxit,num_iter_global, chunk_size=-1, len,rest, n_chunks, np, ipp;
82    dim_t i0;    register double ss,ss1;
83      dim_t i0, istart, iend;
84    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;
85    err_t status = SOLVER_NO_ERROR;    err_t status = SOLVER_NO_ERROR;
86    dim_t n = Paso_SystemMatrix_getTotalNumRows(A);    dim_t n = Paso_SystemMatrix_getTotalNumRows(A);
# Line 82  err_t Paso_Solver_PCG( Line 88  err_t Paso_Solver_PCG(
88    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;
89    double norm_of_residual,norm_of_residual_global, loc_sum[2], sum[2];    double norm_of_residual,norm_of_residual_global, loc_sum[2], sum[2];
90    register double r_tmp,d,rs_tmp,x2_tmp,x_tmp;    register double r_tmp,d,rs_tmp,x2_tmp,x_tmp;
91      char* chksz_chr;
92      np=omp_get_max_threads();
93    
94    #ifdef USE_DYNAMIC_SCHEDULING
95        chksz_chr=getenv("PASO_CHUNK_SIZE_PCG");
96        if (chksz_chr!=NULL) sscanf(chksz_chr, "%d",&chunk_size);
97        chunk_size=MIN(MAX(1,chunk_size),n/np);
98        n_chunks=n/chunk_size;
99        if (n_chunks*chunk_size<n) n_chunks+=1;
100    #else
101        len=n/np;
102        rest=n-len*np;
103    #endif
104  /*                                                                 */  /*                                                                 */
105  /*-----------------------------------------------------------------*/  /*-----------------------------------------------------------------*/
106  /*                                                                 */  /*                                                                 */
# Line 104  err_t Paso_Solver_PCG( Line 122  err_t Paso_Solver_PCG(
122    } else {    } else {
123      maxit = *iter;      maxit = *iter;
124      tol = *resid;      tol = *resid;
125      #pragma omp parallel firstprivate(maxit,tol,convergeFlag,maxIterFlag,breakFlag) \      Performance_startMonitor(pp,PERFORMANCE_SOLVER);
126                                             private(tau_old,tau,beta,delta,gamma_1,gamma_2,alpha,norm_of_residual,num_iter)      /* initialize data */
127        #pragma omp parallel private(i0, istart, iend, ipp)
128      {      {
129         Performance_startMonitor(pp,PERFORMANCE_SOLVER);         #ifdef USE_DYNAMIC_SCHEDULING
130         /* initialize data */             #pragma omp for schedule(dynamic, 1)
131         #pragma omp for private(i0) schedule(static)             for (ipp=0; ipp < n_chunks; ++ipp) {
132         for (i0=0;i0<n;i0++) {                istart=chunk_size*ipp;
133            rs[i0]=r[i0];                iend=MIN(istart+chunk_size,n);
134            x2[i0]=x[i0];         #else
135         }             #pragma omp for schedule(static)
136         #pragma omp for private(i0) schedule(static)             for (ipp=0; ipp <np; ++ipp) {
137         for (i0=0;i0<n;i0++) {                 istart=len*ipp+MIN(ipp,rest);
138            p[i0]=0;                 iend=len*(ipp+1)+MIN(ipp+1,rest);
139            v[i0]=0;         #endif
140         }             #pragma ivdep
141         num_iter=0;             for (i0=istart;i0<iend;i0++) {
142         /* start of iteration */                   rs[i0]=r[i0];
143         while (!(convergeFlag || maxIterFlag || breakFlag)) {                   x2[i0]=x[i0];
144                     p[i0]=0;
145                     v[i0]=0;
146               }
147           #ifdef USE_DYNAMIC_SCHEDULING
148             }
149           #else
150             }
151           #endif
152        }
153        num_iter=0;
154        /* start of iteration */
155        while (!(convergeFlag || maxIterFlag || breakFlag)) {
156             ++(num_iter);             ++(num_iter);
            #pragma omp barrier  
            #pragma omp master  
            {  
            sum_1 = 0;  
            sum_2 = 0;  
            sum_3 = 0;  
            sum_4 = 0;  
            sum_5 = 0;  
            }  
157             /* v=prec(r)  */             /* v=prec(r)  */
158             Performance_stopMonitor(pp,PERFORMANCE_SOLVER);             Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
159             Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);             Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);
# Line 139  err_t Paso_Solver_PCG( Line 161  err_t Paso_Solver_PCG(
161             Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);             Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);
162             Performance_startMonitor(pp,PERFORMANCE_SOLVER);             Performance_startMonitor(pp,PERFORMANCE_SOLVER);
163             /* tau=v*r    */             /* tau=v*r    */
164             #pragma omp for private(i0) reduction(+:sum_1) schedule(static)         sum_1 = 0;
165             for (i0=0;i0<n;i0++) sum_1+=v[i0]*r[i0]; /* Limit to local values of v[] and r[] */             #pragma omp parallel private(i0, istart, iend, ipp, ss)
166               {
167                      ss=0;
168                      #ifdef USE_DYNAMIC_SCHEDULING
169                          #pragma omp for schedule(dynamic, 1)
170                          for (ipp=0; ipp < n_chunks; ++ipp) {
171                             istart=chunk_size*ipp;
172                             iend=MIN(istart+chunk_size,n);
173                      #else
174                          #pragma omp for schedule(static)
175                          for (ipp=0; ipp <np; ++ipp) {
176                              istart=len*ipp+MIN(ipp,rest);
177                              iend=len*(ipp+1)+MIN(ipp+1,rest);
178                      #endif
179                      #pragma ivdep
180                      for (i0=istart;i0<iend;i0++) ss+=v[i0]*r[i0];
181                      #ifdef USE_DYNAMIC_SCHEDULING
182                        }
183                      #else
184                        }
185                      #endif
186                      #pragma critical
187                      {
188                        sum_1+=ss;
189                      }
190               }
191             #ifdef PASO_MPI             #ifdef PASO_MPI
192              /* 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:
193                 OMP master participates in an MPI reduction to get global sum_1 */                 OMP master participates in an MPI reduction to get global sum_1 */
194                  #pragma omp master              loc_saum[0] = sum_1;
195              {              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);  
             }  
196             #endif             #endif
197             tau_old=tau;             tau_old=tau;
198             tau=sum_1;             tau=sum_1;
199             /* p=v+beta*p */             /* p=v+beta*p */
200             if (num_iter==1) {             #pragma omp parallel private(i0, istart, iend, ipp,beta)
201                 #pragma omp for private(i0)  schedule(static)             {
202                 for (i0=0;i0<n;i0++) p[i0]=v[i0];                    #ifdef USE_DYNAMIC_SCHEDULING
203             } else {                        #pragma omp for schedule(dynamic, 1)
204                 beta=tau/tau_old;                        for (ipp=0; ipp < n_chunks; ++ipp) {
205                 #pragma omp for private(i0)  schedule(static)                           istart=chunk_size*ipp;
206                 for (i0=0;i0<n;i0++) p[i0]=v[i0]+beta*p[i0];                           iend=MIN(istart+chunk_size,n);
207                      #else
208                          #pragma omp for schedule(static)
209                          for (ipp=0; ipp <np; ++ipp) {
210                              istart=len*ipp+MIN(ipp,rest);
211                              iend=len*(ipp+1)+MIN(ipp+1,rest);
212                      #endif
213                      if (num_iter==1) {
214                          #pragma ivdep
215                          for (i0=istart;i0<iend;i0++) p[i0]=v[i0];
216                      } else {
217                          beta=tau/tau_old;
218                          #pragma ivdep
219                          for (i0=istart;i0<iend;i0++) p[i0]=v[i0]+beta*p[i0];
220                      }
221                      #ifdef USE_DYNAMIC_SCHEDULING
222                        }
223                      #else
224                        }
225                      #endif
226             }             }
227             /* v=A*p */             /* v=A*p */
228             Performance_stopMonitor(pp,PERFORMANCE_SOLVER);             Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
229             Performance_startMonitor(pp,PERFORMANCE_MVM);             Performance_startMonitor(pp,PERFORMANCE_MVM);
230         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);  
231             Performance_stopMonitor(pp,PERFORMANCE_MVM);             Performance_stopMonitor(pp,PERFORMANCE_MVM);
232             Performance_startMonitor(pp,PERFORMANCE_SOLVER);             Performance_startMonitor(pp,PERFORMANCE_SOLVER);
233    
234             /* delta=p*v */             /* delta=p*v */
235             #pragma omp for private(i0) reduction(+:sum_2) schedule(static)         sum_2 = 0;
236             for (i0=0;i0<n;i0++) sum_2+=v[i0]*p[i0];             #pragma omp parallel private(i0, istart, iend, ipp,ss)
237               {
238                      ss=0;
239                      #ifdef USE_DYNAMIC_SCHEDULING
240                          #pragma omp for schedule(dynamic, 1)
241                          for (ipp=0; ipp < n_chunks; ++ipp) {
242                             istart=chunk_size*ipp;
243                             iend=MIN(istart+chunk_size,n);
244                      #else
245                          #pragma omp for schedule(static)
246                          for (ipp=0; ipp <np; ++ipp) {
247                              istart=len*ipp+MIN(ipp,rest);
248                              iend=len*(ipp+1)+MIN(ipp+1,rest);
249                      #endif
250                      #pragma ivdep
251                      for (i0=istart;i0<iend;i0++) ss+=v[i0]*p[i0];
252                      #ifdef USE_DYNAMIC_SCHEDULING
253                        }
254                      #else
255                        }
256                      #endif
257                      #pragma critical
258                      {
259                         sum_2+=ss;
260                      }
261               }
262             #ifdef PASO_MPI             #ifdef PASO_MPI
263                 #pragma omp master            loc_sum[0] = sum_2;
264             {            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);  
            }  
265             #endif             #endif
266             delta=sum_2;             delta=sum_2;
267    
268        
269             if (! (breakFlag = (ABS(delta) <= TOLERANCE_FOR_SCALARS))) {             if (! (breakFlag = (ABS(delta) <= TOLERANCE_FOR_SCALARS))) {
                alpha=tau/delta;  
270                 /* smoother */                 /* smoother */
271                 #pragma omp for private(i0) schedule(static)             sum_3 = 0;
272                 for (i0=0;i0<n;i0++) r[i0]-=alpha*v[i0];             sum_4 = 0;
273                 #pragma omp for private(i0,d) reduction(+:sum_3,sum_4) schedule(static)                 #pragma omp parallel private(i0, istart, iend, ipp,d, ss, ss1, alpha)
274                 for (i0=0;i0<n;i0++) {                 {
275                       d=r[i0]-rs[i0];                    ss=0;
276                       sum_3+=d*d;                    ss1=0;
277                       sum_4+=d*rs[i0];                    alpha=tau/delta;
278                      #ifdef USE_DYNAMIC_SCHEDULING
279                          #pragma omp for schedule(dynamic, 1)
280                          for (ipp=0; ipp < n_chunks; ++ipp) {
281                             istart=chunk_size*ipp;
282                             iend=MIN(istart+chunk_size,n);
283                      #else
284                          #pragma omp for schedule(static)
285                          for (ipp=0; ipp <np; ++ipp) {
286                              istart=len*ipp+MIN(ipp,rest);
287                              iend=len*(ipp+1)+MIN(ipp+1,rest);
288                      #endif
289                      #pragma ivdep
290                      for (i0=istart;i0<iend;i0++) {
291                            r[i0]-=alpha*v[i0];
292                            d=r[i0]-rs[i0];
293                            ss+=d*d;
294                            ss1+=d*rs[i0];
295                      }
296                      #ifdef USE_DYNAMIC_SCHEDULING
297                        }
298                      #else
299                        }
300                      #endif
301                      #pragma critical
302                      {
303                          sum_3+=ss;
304                          sum_4+=ss1;
305                      }
306                 }                 }
307                 #ifdef PASO_MPI                 #ifdef PASO_MPI
308                     #pragma omp master                 loc_sum[0] = sum_3;
309                 {                 loc_sum[1] = sum_4;
310                   loc_sum[0] = sum_3;                 MPI_Allreduce(loc_sum, sum, 2, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
311                   loc_sum[1] = sum_4;                 sum_3=sum[0];
312                   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];  
                }  
313                  #endif                  #endif
314                  gamma_1= ( (ABS(sum_3)<= ZERO) ? 0 : -sum_4/sum_3) ;              sum_5 = 0;
315                  gamma_2= ONE-gamma_1;                  #pragma omp parallel private(i0, istart, iend, ipp, ss, gamma_1,gamma_2)
316                  #pragma omp for private(i0,x2_tmp,x_tmp,rs_tmp) schedule(static)                  {
317                  for (i0=0;i0<n;++i0) {                    gamma_1= ( (ABS(sum_3)<= ZERO) ? 0 : -sum_4/sum_3) ;
318                    rs[i0]=gamma_2*rs[i0]+gamma_1*r[i0];                    gamma_2= ONE-gamma_1;
319                    x2[i0]+=alpha*p[i0];                    ss=0;
320                    x[i0]=gamma_2*x[i0]+gamma_1*x2[i0];                    #ifdef USE_DYNAMIC_SCHEDULING
321                          #pragma omp for schedule(dynamic, 1)
322                          for (ipp=0; ipp < n_chunks; ++ipp) {
323                             istart=chunk_size*ipp;
324                             iend=MIN(istart+chunk_size,n);
325                      #else
326                          #pragma omp for schedule(static)
327                          for (ipp=0; ipp <np; ++ipp) {
328                              istart=len*ipp+MIN(ipp,rest);
329                              iend=len*(ipp+1)+MIN(ipp+1,rest);
330                      #endif
331                      #pragma ivdep
332                      for (i0=istart;i0<iend;i0++) {
333                          rs[i0]=gamma_2*rs[i0]+gamma_1*r[i0];
334                          x2[i0]+=alpha*p[i0];
335                          x[i0]=gamma_2*x[i0]+gamma_1*x2[i0];
336                          ss+=rs[i0]*rs[i0];
337                      }
338                      #ifdef USE_DYNAMIC_SCHEDULING
339                        }
340                      #else
341                        }
342                      #endif
343                      #pragma omp critical
344                      {
345                         sum_5+=ss;
346                      }
347                  }                  }
                 #pragma omp for private(i0) reduction(+:sum_5) schedule(static)  
                 for (i0=0;i0<n;++i0) sum_5+=rs[i0]*rs[i0];  
348                  #ifdef PASO_MPI                  #ifdef PASO_MPI
349                     #pragma omp master                 loc_sum[0] = sum_5;
350                 {                 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);  
                }  
351                  #endif                  #endif
352                  norm_of_residual=sqrt(sum_5);                  norm_of_residual=sqrt(sum_5);
353                  convergeFlag = norm_of_residual <= tol;                  convergeFlag = norm_of_residual <= tol;
354                  maxIterFlag = num_iter == maxit;                  maxIterFlag = num_iter == maxit;
355                  breakFlag = (ABS(tau) <= TOLERANCE_FOR_SCALARS);                  breakFlag = (ABS(tau) <= TOLERANCE_FOR_SCALARS);
356             }             }
357         }      }
358         /* end of iteration */      /* end of iteration */
359         #pragma omp master      num_iter_global=num_iter;
360         {      norm_of_residual_global=norm_of_residual;
361             num_iter_global=num_iter;      if (maxIterFlag) {
362             norm_of_residual_global=norm_of_residual;           status = SOLVER_MAXITER_REACHED;
363             if (maxIterFlag) {      } else if (breakFlag) {
364                 status = SOLVER_MAXITER_REACHED;           status = SOLVER_BREAKDOWN;
365             } else if (breakFlag) {      }
366                 status = SOLVER_BREAKDOWN;      Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
            }  
        }  
        Performance_stopMonitor(pp,PERFORMANCE_SOLVER);  
     }  /* end of parallel region */  
367      TMPMEMFREE(rs);      TMPMEMFREE(rs);
368      TMPMEMFREE(x2);      TMPMEMFREE(x2);
369      TMPMEMFREE(v);      TMPMEMFREE(v);

Legend:
Removed from v.1387  
changed lines
  Added in v.1571

  ViewVC Help
Powered by ViewVC 1.1.26