/[escript]/branches/doubleplusgood/paso/src/PCG.cpp
ViewVC logotype

Diff of /branches/doubleplusgood/paso/src/PCG.cpp

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

temp_trunk_copy/paso/src/PCG.c revision 1384 by phornby, Fri Jan 11 02:29:38 2008 UTC trunk/paso/src/PCG.c revision 1811 by ksteube, Thu Sep 25 23:11:13 2008 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2008 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15  /* PCG iterations */  /* PCG iterations */
16    
# Line 63  Line 62 
62  *  ==============================================================  *  ==============================================================
63  */  */
64    
65    /* #define PASO_DYNAMIC_SCHEDULING_MVM */
66    
67    #if defined PASO_DYNAMIC_SCHEDULING_MVM && defined __OPENMP
68    #define USE_DYNAMIC_SCHEDULING
69    #endif
70    
71  err_t Paso_Solver_PCG(  err_t Paso_Solver_PCG(
72      Paso_SystemMatrix * A,      Paso_SystemMatrix * A,
73      double * r,      double * r,
# Line 71  err_t Paso_Solver_PCG( Line 76  err_t Paso_Solver_PCG(
76      double * tolerance,      double * tolerance,
77      Paso_Performance* pp) {      Paso_Performance* pp) {
78    
   
79    /* Local variables */    /* Local variables */
80    dim_t num_iter=0,maxit,num_iter_global;    dim_t num_iter=0,maxit,num_iter_global, chunk_size=-1, len,rest, np, ipp;
81    dim_t i0;    register double ss,ss1;
82      dim_t i0, istart, iend;
83    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;
84    err_t status = SOLVER_NO_ERROR;    err_t status = SOLVER_NO_ERROR;
85    dim_t n = Paso_SystemMatrix_getTotalNumRows(A);    dim_t n = Paso_SystemMatrix_getTotalNumRows(A);
86    double *resid = tolerance, *rs=NULL, *p=NULL, *v=NULL, *x2=NULL ;    double *resid = tolerance, *rs=NULL, *p=NULL, *v=NULL, *x2=NULL ;
87    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;
88    double norm_of_residual,norm_of_residual_global, loc_sum[2], sum[2];  #ifdef PASO_MPI
89    register double r_tmp,d,rs_tmp,x2_tmp,x_tmp;    double loc_sum[2], sum[2];
90    #endif
91      double norm_of_residual,norm_of_residual_global;
92      register double d;
93    
94      /* Should not be any executable code before this ifdef */
95    
96    #ifdef USE_DYNAMIC_SCHEDULING
97    
98        /* Watch out for these declarations (see above) */
99        char* chksz_chr;
100        dim_t n_chunks;
101    
102        chksz_chr=getenv("PASO_CHUNK_SIZE_PCG");
103        if (chksz_chr!=NULL) sscanf(chksz_chr, "%d",&chunk_size);
104        np=omp_get_max_threads();
105        chunk_size=MIN(MAX(1,chunk_size),n/np);
106        n_chunks=n/chunk_size;
107        if (n_chunks*chunk_size<n) n_chunks+=1;
108    #else
109        np=omp_get_max_threads();
110        len=n/np;
111        rest=n-len*np;
112    #endif
113  /*                                                                 */  /*                                                                 */
114  /*-----------------------------------------------------------------*/  /*-----------------------------------------------------------------*/
115  /*                                                                 */  /*                                                                 */
# Line 104  err_t Paso_Solver_PCG( Line 131  err_t Paso_Solver_PCG(
131    } else {    } else {
132      maxit = *iter;      maxit = *iter;
133      tol = *resid;      tol = *resid;
134      #pragma omp parallel firstprivate(maxit,tol,convergeFlag,maxIterFlag,breakFlag) \      Performance_startMonitor(pp,PERFORMANCE_SOLVER);
135                                             private(tau_old,tau,beta,delta,gamma_1,gamma_2,alpha,norm_of_residual,num_iter)      /* initialize data */
136        #pragma omp parallel private(i0, istart, iend, ipp)
137      {      {
138         Performance_startMonitor(pp,PERFORMANCE_SOLVER);         #ifdef USE_DYNAMIC_SCHEDULING
139         /* initialize data */             #pragma omp for schedule(dynamic, 1)
140         #pragma omp for private(i0) schedule(static)             for (ipp=0; ipp < n_chunks; ++ipp) {
141         for (i0=0;i0<n;i0++) {                istart=chunk_size*ipp;
142            rs[i0]=r[i0];                iend=MIN(istart+chunk_size,n);
143            x2[i0]=x[i0];         #else
144         }             #pragma omp for schedule(static)
145         #pragma omp for private(i0) schedule(static)             for (ipp=0; ipp <np; ++ipp) {
146         for (i0=0;i0<n;i0++) {                 istart=len*ipp+MIN(ipp,rest);
147            p[i0]=0;                 iend=len*(ipp+1)+MIN(ipp+1,rest);
148            v[i0]=0;         #endif
149         }                 #pragma ivdep
150         num_iter=0;                 for (i0=istart;i0<iend;i0++) {
151         /* start of iteration */                   rs[i0]=r[i0];
152         while (!(convergeFlag || maxIterFlag || breakFlag)) {                   x2[i0]=x[i0];
153                     p[i0]=0;
154                     v[i0]=0;
155                   }
156           #ifdef USE_DYNAMIC_SCHEDULING
157             }
158           #else
159             }
160           #endif
161        }
162        num_iter=0;
163    
164        /* PGH */
165        /* without this we get a use of an unititialised var below */
166        tau = 0;
167    
168        /* start of iteration */
169        while (!(convergeFlag || maxIterFlag || breakFlag)) {
170             ++(num_iter);             ++(num_iter);
171             #pragma omp barrier  
172             #pragma omp master             /* PGH */
173             {             /* The next lines were commented out before I got here */
            sum_1 = 0;  
            sum_2 = 0;  
            sum_3 = 0;  
            sum_4 = 0;  
            sum_5 = 0;  
            }  
174             /* v=prec(r)  */             /* v=prec(r)  */
175               /* tau=v*r; */
176               /* leading to the use of an unititialised var below */
177    
178             Performance_stopMonitor(pp,PERFORMANCE_SOLVER);             Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
179             Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);             Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);
180             Paso_Solver_solvePreconditioner(A,v,r);             Paso_Solver_solvePreconditioner(A,v,r);
181             Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);             Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);
182             Performance_startMonitor(pp,PERFORMANCE_SOLVER);             Performance_startMonitor(pp,PERFORMANCE_SOLVER);
183             /* tau=v*r    */  
184             #pragma omp for private(i0) reduction(+:sum_1) schedule(static)         sum_1 = 0;
185             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)
186               {
187                      ss=0;
188                      #ifdef USE_DYNAMIC_SCHEDULING
189                          #pragma omp for schedule(dynamic, 1)
190                          for (ipp=0; ipp < n_chunks; ++ipp) {
191                             istart=chunk_size*ipp;
192                             iend=MIN(istart+chunk_size,n);
193                      #else
194                          #pragma omp for schedule(static)
195                          for (ipp=0; ipp <np; ++ipp) {
196                              istart=len*ipp+MIN(ipp,rest);
197                              iend=len*(ipp+1)+MIN(ipp+1,rest);
198                      #endif
199                  #pragma ivdep
200                              for (i0=istart;i0<iend;i0++) ss+=v[i0]*r[i0];
201                      #ifdef USE_DYNAMIC_SCHEDULING
202                        }
203                      #else
204                        }
205                      #endif
206                      #pragma omp critical
207                      {
208                         sum_1+=ss;
209                      }
210               }
211             #ifdef PASO_MPI             #ifdef PASO_MPI
212              /* 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:
213                 OMP master participates in an MPI reduction to get global sum_1 */                 OMP master participates in an MPI reduction to get global sum_1 */
214                  #pragma omp master              loc_sum[0] = sum_1;
215              {              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);  
             }  
216             #endif             #endif
217             tau_old=tau;             tau_old=tau;
218             tau=sum_1;             tau=sum_1;
219             /* p=v+beta*p */             /* p=v+beta*p */
220             if (num_iter==1) {             #pragma omp parallel private(i0, istart, iend, ipp,beta)
221                 #pragma omp for private(i0)  schedule(static)             {
222                 for (i0=0;i0<n;i0++) p[i0]=v[i0];                    #ifdef USE_DYNAMIC_SCHEDULING
223             } else {                        #pragma omp for schedule(dynamic, 1)
224                 beta=tau/tau_old;                        for (ipp=0; ipp < n_chunks; ++ipp) {
225                 #pragma omp for private(i0)  schedule(static)                           istart=chunk_size*ipp;
226                 for (i0=0;i0<n;i0++) p[i0]=v[i0]+beta*p[i0];                           iend=MIN(istart+chunk_size,n);
227                      #else
228                          #pragma omp for schedule(static)
229                          for (ipp=0; ipp <np; ++ipp) {
230                              istart=len*ipp+MIN(ipp,rest);
231                              iend=len*(ipp+1)+MIN(ipp+1,rest);
232                      #endif
233                              if (num_iter==1) {
234                      #pragma ivdep
235                                  for (i0=istart;i0<iend;i0++) p[i0]=v[i0];
236                              } else {
237                                  beta=tau/tau_old;
238                      #pragma ivdep
239                                  for (i0=istart;i0<iend;i0++) p[i0]=v[i0]+beta*p[i0];
240                              }
241                      #ifdef USE_DYNAMIC_SCHEDULING
242                        }
243                      #else
244                        }
245                      #endif
246             }             }
247             /* v=A*p */             /* v=A*p */
248             Performance_stopMonitor(pp,PERFORMANCE_SOLVER);             Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
249             Performance_startMonitor(pp,PERFORMANCE_MVM);             Performance_startMonitor(pp,PERFORMANCE_MVM);
250         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);  
251             Performance_stopMonitor(pp,PERFORMANCE_MVM);             Performance_stopMonitor(pp,PERFORMANCE_MVM);
252             Performance_startMonitor(pp,PERFORMANCE_SOLVER);             Performance_startMonitor(pp,PERFORMANCE_SOLVER);
253    
254             /* delta=p*v */             /* delta=p*v */
255             #pragma omp for private(i0) reduction(+:sum_2) schedule(static)         sum_2 = 0;
256             for (i0=0;i0<n;i0++) sum_2+=v[i0]*p[i0];             #pragma omp parallel private(i0, istart, iend, ipp,ss)
257               {
258                      ss=0;
259                      #ifdef USE_DYNAMIC_SCHEDULING
260                          #pragma omp for schedule(dynamic, 1)
261                          for (ipp=0; ipp < n_chunks; ++ipp) {
262                             istart=chunk_size*ipp;
263                             iend=MIN(istart+chunk_size,n);
264                      #else
265                          #pragma omp for schedule(static)
266                          for (ipp=0; ipp <np; ++ipp) {
267                              istart=len*ipp+MIN(ipp,rest);
268                              iend=len*(ipp+1)+MIN(ipp+1,rest);
269                      #endif
270                  #pragma ivdep
271                              for (i0=istart;i0<iend;i0++) ss+=v[i0]*p[i0];
272                      #ifdef USE_DYNAMIC_SCHEDULING
273                        }
274                      #else
275                        }
276                      #endif
277                      #pragma omp critical
278                      {
279                                 sum_2+=ss;
280                      }
281               }
282             #ifdef PASO_MPI             #ifdef PASO_MPI
283                 #pragma omp master            loc_sum[0] = sum_2;
284             {            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);  
            }  
285             #endif             #endif
286             delta=sum_2;             delta=sum_2;
287               alpha=tau/delta;
288        
289             if (! (breakFlag = (ABS(delta) <= TOLERANCE_FOR_SCALARS))) {             if (! (breakFlag = (ABS(delta) <= TOLERANCE_FOR_SCALARS))) {
                alpha=tau/delta;  
290                 /* smoother */                 /* smoother */
291                 #pragma omp for private(i0) schedule(static)             sum_3 = 0;
292                 for (i0=0;i0<n;i0++) r[i0]-=alpha*v[i0];             sum_4 = 0;
293                 #pragma omp for private(i0,d) reduction(+:sum_3,sum_4) schedule(static)                 #pragma omp parallel private(i0, istart, iend, ipp,d, ss, ss1)
294                 for (i0=0;i0<n;i0++) {                 {
295                       d=r[i0]-rs[i0];                    ss=0;
296                       sum_3+=d*d;                    ss1=0;
297                       sum_4+=d*rs[i0];                    #ifdef USE_DYNAMIC_SCHEDULING
298                          #pragma omp for schedule(dynamic, 1)
299                          for (ipp=0; ipp < n_chunks; ++ipp) {
300                             istart=chunk_size*ipp;
301                             iend=MIN(istart+chunk_size,n);
302                      #else
303                          #pragma omp for schedule(static)
304                          for (ipp=0; ipp <np; ++ipp) {
305                              istart=len*ipp+MIN(ipp,rest);
306                              iend=len*(ipp+1)+MIN(ipp+1,rest);
307                      #endif
308                  #pragma ivdep
309                              for (i0=istart;i0<iend;i0++) {
310                                    r[i0]-=alpha*v[i0];
311                                    d=r[i0]-rs[i0];
312                                    ss+=d*d;
313                                    ss1+=d*rs[i0];
314                              }
315                      #ifdef USE_DYNAMIC_SCHEDULING
316                        }
317                      #else
318                        }
319                      #endif
320                      #pragma omp critical
321                      {
322                         sum_3+=ss;
323                         sum_4+=ss1;
324                      }
325                 }                 }
326                 #ifdef PASO_MPI                 #ifdef PASO_MPI
327                     #pragma omp master                 loc_sum[0] = sum_3;
328                 {                 loc_sum[1] = sum_4;
329                   loc_sum[0] = sum_3;                 MPI_Allreduce(loc_sum, sum, 2, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
330                   loc_sum[1] = sum_4;                 sum_3=sum[0];
331                   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];  
                }  
332                  #endif                  #endif
333                  gamma_1= ( (ABS(sum_3)<= ZERO) ? 0 : -sum_4/sum_3) ;              sum_5 = 0;
334                  gamma_2= ONE-gamma_1;                  #pragma omp parallel private(i0, istart, iend, ipp, ss, gamma_1,gamma_2)
335                  #pragma omp for private(i0,x2_tmp,x_tmp,rs_tmp) schedule(static)                  {
336                  for (i0=0;i0<n;++i0) {                    gamma_1= ( (ABS(sum_3)<= ZERO) ? 0 : -sum_4/sum_3) ;
337                    rs[i0]=gamma_2*rs[i0]+gamma_1*r[i0];                    gamma_2= ONE-gamma_1;
338                    x2[i0]+=alpha*p[i0];                    ss=0;
339                    x[i0]=gamma_2*x[i0]+gamma_1*x2[i0];                    #ifdef USE_DYNAMIC_SCHEDULING
340                          #pragma omp for schedule(dynamic, 1)
341                          for (ipp=0; ipp < n_chunks; ++ipp) {
342                             istart=chunk_size*ipp;
343                             iend=MIN(istart+chunk_size,n);
344                      #else
345                          #pragma omp for schedule(static)
346                          for (ipp=0; ipp <np; ++ipp) {
347                              istart=len*ipp+MIN(ipp,rest);
348                              iend=len*(ipp+1)+MIN(ipp+1,rest);
349                      #endif
350                  #pragma ivdep
351                              for (i0=istart;i0<iend;i0++) {
352                                  rs[i0]=gamma_2*rs[i0]+gamma_1*r[i0];
353                                  x2[i0]+=alpha*p[i0];
354                                  x[i0]=gamma_2*x[i0]+gamma_1*x2[i0];
355                                  ss+=rs[i0]*rs[i0];
356                              }
357                      #ifdef USE_DYNAMIC_SCHEDULING
358                        }
359                      #else
360                        }
361                      #endif
362                      #pragma omp critical
363                      {
364                          sum_5+=ss;
365                      }
366                  }                  }
                 #pragma omp for private(i0) reduction(+:sum_5) schedule(static)  
                 for (i0=0;i0<n;++i0) sum_5+=rs[i0]*rs[i0];  
367                  #ifdef PASO_MPI                  #ifdef PASO_MPI
368                     #pragma omp master                 loc_sum[0] = sum_5;
369                 {                 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);  
                }  
370                  #endif                  #endif
371                  norm_of_residual=sqrt(sum_5);                  norm_of_residual=sqrt(sum_5);
372                  convergeFlag = norm_of_residual <= tol;                  convergeFlag = norm_of_residual <= tol;
373                  maxIterFlag = num_iter == maxit;                  maxIterFlag = num_iter == maxit;
374                  breakFlag = (ABS(tau) <= TOLERANCE_FOR_SCALARS);                  breakFlag = (ABS(tau) <= TOLERANCE_FOR_SCALARS);
375             }             }
376         }      }
377         /* end of iteration */      /* end of iteration */
378         #pragma omp master      num_iter_global=num_iter;
379         {      norm_of_residual_global=norm_of_residual;
380             num_iter_global=num_iter;      if (maxIterFlag) {
381             norm_of_residual_global=norm_of_residual;           status = SOLVER_MAXITER_REACHED;
382             if (maxIterFlag) {      } else if (breakFlag) {
383                 status = SOLVER_MAXITER_REACHED;           status = SOLVER_BREAKDOWN;
384             } else if (breakFlag) {      }
385                 status = SOLVER_BREAKDOWN;      Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
            }  
        }  
        Performance_stopMonitor(pp,PERFORMANCE_SOLVER);  
     }  /* end of parallel region */  
386      TMPMEMFREE(rs);      TMPMEMFREE(rs);
387      TMPMEMFREE(x2);      TMPMEMFREE(x2);
388      TMPMEMFREE(v);      TMPMEMFREE(v);

Legend:
Removed from v.1384  
changed lines
  Added in v.1811

  ViewVC Help
Powered by ViewVC 1.1.26