/[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

revision 1556 by gross, Mon May 12 00:54:58 2008 UTC revision 1570 by gross, Sat May 24 21:31:04 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 106  err_t Paso_Solver_PCG( Line 124  err_t Paso_Solver_PCG(
124      tol = *resid;      tol = *resid;
125      Performance_startMonitor(pp,PERFORMANCE_SOLVER);      Performance_startMonitor(pp,PERFORMANCE_SOLVER);
126      /* initialize data */      /* initialize data */
127      #pragma omp parallel      #pragma omp parallel (i0, istart, iend, ipp)
128      {      {
129             #pragma omp for private(i0) schedule(static)         #ifdef USE_DYNAMIC_SCHEDULING
130             for (i0=0;i0<n;i0++) {             #pragma omp for schedule(dynamic, 1)
131                rs[i0]=r[i0];             for (ipp=0; ipp < n_chunks; ++ipp) {
132                x2[i0]=x[i0];                istart=chunk_size*ipp;
133                p[i0]=0;                iend=MIN(istart+chunk_size,n);
134                v[i0]=0;         #else
135               #pragma omp for schedule(static)
136               for (ipp=0; ipp <np; ++ipp) {
137                   istart=len*ipp+MIN(ipp,rest);
138                   iend=len*(ipp+1)+MIN(ipp+1,rest);
139           #endif
140               #pragma ivdep
141               for (i0=istart;i0<iend;i0++) {
142                     rs[i0]=r[i0];
143                     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;      num_iter=0;
154      /* start of iteration */      /* start of iteration */
155      while (!(convergeFlag || maxIterFlag || breakFlag)) {      while (!(convergeFlag || maxIterFlag || breakFlag)) {
156             ++(num_iter);             ++(num_iter);
        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 132  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 parallel 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 (i0, istart, iend, ipp, ss)
166               {
167                      #ifdef USE_DYNAMIC_SCHEDULING
168                          #pragma omp for schedule(dynamic, 1)
169                          for (ipp=0; ipp < n_chunks; ++ipp) {
170                             istart=chunk_size*ipp;
171                             iend=MIN(istart+chunk_size,n);
172                      #else
173                          #pragma omp for schedule(static)
174                          for (ipp=0; ipp <np; ++ipp) {
175                              istart=len*ipp+MIN(ipp,rest);
176                              iend=len*(ipp+1)+MIN(ipp+1,rest);
177                      #endif
178                      ss=0;
179                      #pragma ivdep
180                      for (i0=istart;i0<iend;i0++) ss+=v[i0]*r[i0];
181                      #pragma critical
182                      sum_1+=ss;
183                      #ifdef USE_DYNAMIC_SCHEDULING
184                        }
185                      #else
186                        }
187                      #endif
188               }
189             #ifdef PASO_MPI             #ifdef PASO_MPI
190              /* 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:
191                 OMP master participates in an MPI reduction to get global sum_1 */                 OMP master participates in an MPI reduction to get global sum_1 */
192              loc_sum[0] = sum_1;              loc_saum[0] = sum_1;
193              MPI_Allreduce(loc_sum, &sum_1, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);              MPI_Allreduce(loc_sum, &sum_1, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
194             #endif             #endif
195             tau_old=tau;             tau_old=tau;
196             tau=sum_1;             tau=sum_1;
197             /* p=v+beta*p */             /* p=v+beta*p */
198             if (num_iter==1) {             #pragma omp parallel (i0, istart, iend, ipp,beta)
199                 #pragma omp parallel for private(i0)  schedule(static)             {
200                 for (i0=0;i0<n;i0++) p[i0]=v[i0];                    #ifdef USE_DYNAMIC_SCHEDULING
201             } else {                        #pragma omp for schedule(dynamic, 1)
202                 beta=tau/tau_old;                        for (ipp=0; ipp < n_chunks; ++ipp) {
203                 #pragma omp parallel for private(i0)  schedule(static)                           istart=chunk_size*ipp;
204                 for (i0=0;i0<n;i0++) p[i0]=v[i0]+beta*p[i0];                           iend=MIN(istart+chunk_size,n);
205                      #else
206                          #pragma omp for schedule(static)
207                          for (ipp=0; ipp <np; ++ipp) {
208                              istart=len*ipp+MIN(ipp,rest);
209                              iend=len*(ipp+1)+MIN(ipp+1,rest);
210                      #endif
211                      if (num_iter==1) {
212                          #pragma ivdep
213                          for (i0=istart;i0<iend;i0++) p[i0]=v[i0];
214                      } else {
215                          beta=tau/tau_old;
216                          #pragma ivdep
217                          for (i0=istart;i0<iend;i0++) p[i0]=v[i0]+beta*p[i0];
218                      }
219                      #ifdef USE_DYNAMIC_SCHEDULING
220                        }
221                      #else
222                        }
223                      #endif
224             }             }
225             /* v=A*p */             /* v=A*p */
226             Performance_stopMonitor(pp,PERFORMANCE_SOLVER);             Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
# Line 157  err_t Paso_Solver_PCG( Line 228  err_t Paso_Solver_PCG(
228         Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, p,ZERO,v);         Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, p,ZERO,v);
229             Performance_stopMonitor(pp,PERFORMANCE_MVM);             Performance_stopMonitor(pp,PERFORMANCE_MVM);
230             Performance_startMonitor(pp,PERFORMANCE_SOLVER);             Performance_startMonitor(pp,PERFORMANCE_SOLVER);
231    
232             /* delta=p*v */             /* delta=p*v */
233             #pragma omp parallel for private(i0) reduction(+:sum_2) schedule(static)         sum_2 = 0;
234             for (i0=0;i0<n;i0++) sum_2+=v[i0]*p[i0];             #pragma omp parallel (i0, istart, iend, ipp,ss)
235               {
236                      #ifdef USE_DYNAMIC_SCHEDULING
237                          #pragma omp for schedule(dynamic, 1)
238                          for (ipp=0; ipp < n_chunks; ++ipp) {
239                             istart=chunk_size*ipp;
240                             iend=MIN(istart+chunk_size,n);
241                      #else
242                          #pragma omp for schedule(static)
243                          for (ipp=0; ipp <np; ++ipp) {
244                              istart=len*ipp+MIN(ipp,rest);
245                              iend=len*(ipp+1)+MIN(ipp+1,rest);
246                      #endif
247                      ss=0;
248                      #pragma ivdep
249                      for (i0=istart;i0<iend;i0++) ss+=v[i0]*p[i0];
250                      #pragma critical
251                      sum_2+=ss;
252                      #ifdef USE_DYNAMIC_SCHEDULING
253                        }
254                      #else
255                        }
256                      #endif
257               }
258             #ifdef PASO_MPI             #ifdef PASO_MPI
259            loc_sum[0] = sum_2;            loc_sum[0] = sum_2;
260            MPI_Allreduce(loc_sum, &sum_2, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);            MPI_Allreduce(loc_sum, &sum_2, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
# Line 168  err_t Paso_Solver_PCG( Line 263  err_t Paso_Solver_PCG(
263    
264        
265             if (! (breakFlag = (ABS(delta) <= TOLERANCE_FOR_SCALARS))) {             if (! (breakFlag = (ABS(delta) <= TOLERANCE_FOR_SCALARS))) {
                alpha=tau/delta;  
266                 /* smoother */                 /* smoother */
267                 #pragma omp parallel             sum_3 = 0;
268               sum_4 = 0;
269                   #pragma omp parallel (i0, istart, iend, ipp,d, ss, ss1, alpha)
270                 {                 {
271                    #pragma omp for private(i0) schedule(static)                    #ifdef USE_DYNAMIC_SCHEDULING
272                    for (i0=0;i0<n;i0++) r[i0]-=alpha*v[i0];                        #pragma omp for schedule(dynamic, 1)
273                    #pragma omp for private(i0,d) reduction(+:sum_3,sum_4) schedule(static)                        for (ipp=0; ipp < n_chunks; ++ipp) {
274                    for (i0=0;i0<n;i0++) {                           istart=chunk_size*ipp;
275                             iend=MIN(istart+chunk_size,n);
276                      #else
277                          #pragma omp for schedule(static)
278                          for (ipp=0; ipp <np; ++ipp) {
279                              istart=len*ipp+MIN(ipp,rest);
280                              iend=len*(ipp+1)+MIN(ipp+1,rest);
281                      #endif
282                      ss=0;
283                      ss1=0;
284                      alpha=tau/delta;
285                      #pragma ivdep
286                      for (i0=istart;i0<iend;i0++) {
287                            r[i0]-=alpha*v[i0];
288                          d=r[i0]-rs[i0];                          d=r[i0]-rs[i0];
289                          sum_3+=d*d;                          ss+=d*d;
290                          sum_4+=d*rs[i0];                          ss1+=d*rs[i0];
291                    }                    }
292                      #pragma critical
293                      {
294                          sum_3+=ss;
295                          sum_4+=ss1;
296                      }
297                      #ifdef USE_DYNAMIC_SCHEDULING
298                        }
299                      #else
300                        }
301                      #endif
302                 }                 }
303                 #ifdef PASO_MPI                 #ifdef PASO_MPI
304                 loc_sum[0] = sum_3;                 loc_sum[0] = sum_3;
# Line 188  err_t Paso_Solver_PCG( Line 307  err_t Paso_Solver_PCG(
307                 sum_3=sum[0];                 sum_3=sum[0];
308                 sum_4=sum[1];                 sum_4=sum[1];
309                  #endif                  #endif
310                  gamma_1= ( (ABS(sum_3)<= ZERO) ? 0 : -sum_4/sum_3) ;              sum_5 = 0;
311                  gamma_2= ONE-gamma_1;                  #pragma omp parallel (i0, istart, iend, ipp, ss, gamma_1,gamma_2)
                 #pragma omp parallel  
312                  {                  {
313                      #pragma omp for private(i0,x2_tmp,x_tmp,rs_tmp) schedule(static)                    gamma_1= ( (ABS(sum_3)<= ZERO) ? 0 : -sum_4/sum_3) ;
314                      for (i0=0;i0<n;++i0) {                    gamma_2= ONE-gamma_1;
315                      #ifdef USE_DYNAMIC_SCHEDULING
316                          #pragma omp for schedule(dynamic, 1)
317                          for (ipp=0; ipp < n_chunks; ++ipp) {
318                             istart=chunk_size*ipp;
319                             iend=MIN(istart+chunk_size,n);
320                      #else
321                          #pragma omp for schedule(static)
322                          for (ipp=0; ipp <np; ++ipp) {
323                              istart=len*ipp+MIN(ipp,rest);
324                              iend=len*(ipp+1)+MIN(ipp+1,rest);
325                      #endif
326                      ss=0;
327                      #pragma ivdep
328                      for (i0=istart;i0<iend;i0++) {
329                        rs[i0]=gamma_2*rs[i0]+gamma_1*r[i0];                        rs[i0]=gamma_2*rs[i0]+gamma_1*r[i0];
330                        x2[i0]+=alpha*p[i0];                        x2[i0]+=alpha*p[i0];
331                        x[i0]=gamma_2*x[i0]+gamma_1*x2[i0];                        x[i0]=gamma_2*x[i0]+gamma_1*x2[i0];
332                          ss+=rs[i0]*rs[i0];
333                      }
334                      #pragma omp critical
335                      sum_5+=ss;
336                      #ifdef USE_DYNAMIC_SCHEDULING
337                        }
338                      #else
339                      }                      }
340                      #pragma omp for private(i0) reduction(+:sum_5) schedule(static)                    #endif
                     for (i0=0;i0<n;++i0) sum_5+=rs[i0]*rs[i0];  
341                  }                  }
342                  #ifdef PASO_MPI                  #ifdef PASO_MPI
343                 loc_sum[0] = sum_5;                 loc_sum[0] = sum_5;

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

  ViewVC Help
Powered by ViewVC 1.1.26