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

revision 1571 by gross, Sat May 24 22:28:33 2008 UTC 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 78  err_t Paso_Solver_PCG( Line 77  err_t Paso_Solver_PCG(
77      Paso_Performance* pp) {      Paso_Performance* pp) {
78    
79    /* Local variables */    /* Local variables */
80    dim_t num_iter=0,maxit,num_iter_global, chunk_size=-1, len,rest, n_chunks, np, ipp;    dim_t num_iter=0,maxit,num_iter_global, chunk_size=-1, len,rest, np, ipp;
81    register double ss,ss1;    register double ss,ss1;
82    dim_t i0, istart, iend;    dim_t i0, istart, iend;
83    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;
# Line 86  err_t Paso_Solver_PCG( Line 85  err_t Paso_Solver_PCG(
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    char* chksz_chr;  #endif
91    np=omp_get_max_threads();    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  #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");      chksz_chr=getenv("PASO_CHUNK_SIZE_PCG");
103      if (chksz_chr!=NULL) sscanf(chksz_chr, "%d",&chunk_size);      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);      chunk_size=MIN(MAX(1,chunk_size),n/np);
106      n_chunks=n/chunk_size;      n_chunks=n/chunk_size;
107      if (n_chunks*chunk_size<n) n_chunks+=1;      if (n_chunks*chunk_size<n) n_chunks+=1;
108  #else  #else
109        np=omp_get_max_threads();
110      len=n/np;      len=n/np;
111      rest=n-len*np;      rest=n-len*np;
112  #endif  #endif
# Line 137  err_t Paso_Solver_PCG( Line 146  err_t Paso_Solver_PCG(
146                 istart=len*ipp+MIN(ipp,rest);                 istart=len*ipp+MIN(ipp,rest);
147                 iend=len*(ipp+1)+MIN(ipp+1,rest);                 iend=len*(ipp+1)+MIN(ipp+1,rest);
148         #endif         #endif
149             #pragma ivdep                 #pragma ivdep
150             for (i0=istart;i0<iend;i0++) {                 for (i0=istart;i0<iend;i0++) {
151                   rs[i0]=r[i0];                   rs[i0]=r[i0];
152                   x2[i0]=x[i0];                   x2[i0]=x[i0];
153                   p[i0]=0;                   p[i0]=0;
154                   v[i0]=0;                   v[i0]=0;
155             }                 }
156         #ifdef USE_DYNAMIC_SCHEDULING         #ifdef USE_DYNAMIC_SCHEDULING
157           }           }
158         #else         #else
# Line 151  err_t Paso_Solver_PCG( Line 160  err_t Paso_Solver_PCG(
160         #endif         #endif
161      }      }
162      num_iter=0;      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 */      /* start of iteration */
169      while (!(convergeFlag || maxIterFlag || breakFlag)) {      while (!(convergeFlag || maxIterFlag || breakFlag)) {
170             ++(num_iter);             ++(num_iter);
171    
172               /* PGH */
173               /* The next lines were commented out before I got here */
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         sum_1 = 0;         sum_1 = 0;
185             #pragma omp parallel private(i0, istart, iend, ipp, ss)             #pragma omp parallel private(i0, istart, iend, ipp, ss)
186             {             {
# Line 171  err_t Paso_Solver_PCG( Line 191  err_t Paso_Solver_PCG(
191                           istart=chunk_size*ipp;                           istart=chunk_size*ipp;
192                           iend=MIN(istart+chunk_size,n);                           iend=MIN(istart+chunk_size,n);
193                    #else                    #else
194                        #pragma omp for schedule(static)                        #pragma omp for schedule(static)
195                        for (ipp=0; ipp <np; ++ipp) {                        for (ipp=0; ipp <np; ++ipp) {
196                            istart=len*ipp+MIN(ipp,rest);                            istart=len*ipp+MIN(ipp,rest);
197                            iend=len*(ipp+1)+MIN(ipp+1,rest);                            iend=len*(ipp+1)+MIN(ipp+1,rest);
198                    #endif                    #endif
199                    #pragma ivdep                #pragma ivdep
200                    for (i0=istart;i0<iend;i0++) ss+=v[i0]*r[i0];                            for (i0=istart;i0<iend;i0++) ss+=v[i0]*r[i0];
201                    #ifdef USE_DYNAMIC_SCHEDULING                    #ifdef USE_DYNAMIC_SCHEDULING
202                      }                      }
203                    #else                    #else
204                      }                      }
205                    #endif                    #endif
206                    #pragma critical                    #pragma omp critical
207                    {                    {
208                      sum_1+=ss;                       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              loc_saum[0] = sum_1;              loc_sum[0] = sum_1;
215              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);
216             #endif             #endif
217             tau_old=tau;             tau_old=tau;
# Line 210  err_t Paso_Solver_PCG( Line 230  err_t Paso_Solver_PCG(
230                            istart=len*ipp+MIN(ipp,rest);                            istart=len*ipp+MIN(ipp,rest);
231                            iend=len*(ipp+1)+MIN(ipp+1,rest);                            iend=len*(ipp+1)+MIN(ipp+1,rest);
232                    #endif                    #endif
233                    if (num_iter==1) {                            if (num_iter==1) {
234                        #pragma ivdep                    #pragma ivdep
235                        for (i0=istart;i0<iend;i0++) p[i0]=v[i0];                                for (i0=istart;i0<iend;i0++) p[i0]=v[i0];
236                    } else {                            } else {
237                        beta=tau/tau_old;                                beta=tau/tau_old;
238                        #pragma ivdep                    #pragma ivdep
239                        for (i0=istart;i0<iend;i0++) p[i0]=v[i0]+beta*p[i0];                                for (i0=istart;i0<iend;i0++) p[i0]=v[i0]+beta*p[i0];
240                    }                            }
241                    #ifdef USE_DYNAMIC_SCHEDULING                    #ifdef USE_DYNAMIC_SCHEDULING
242                      }                      }
243                    #else                    #else
# Line 247  err_t Paso_Solver_PCG( Line 267  err_t Paso_Solver_PCG(
267                            istart=len*ipp+MIN(ipp,rest);                            istart=len*ipp+MIN(ipp,rest);
268                            iend=len*(ipp+1)+MIN(ipp+1,rest);                            iend=len*(ipp+1)+MIN(ipp+1,rest);
269                    #endif                    #endif
270                    #pragma ivdep                #pragma ivdep
271                    for (i0=istart;i0<iend;i0++) ss+=v[i0]*p[i0];                            for (i0=istart;i0<iend;i0++) ss+=v[i0]*p[i0];
272                    #ifdef USE_DYNAMIC_SCHEDULING                    #ifdef USE_DYNAMIC_SCHEDULING
273                      }                      }
274                    #else                    #else
275                      }                      }
276                    #endif                    #endif
277                    #pragma critical                    #pragma omp critical
278                    {                    {
279                       sum_2+=ss;                               sum_2+=ss;
280                    }                    }
281             }             }
282             #ifdef PASO_MPI             #ifdef PASO_MPI
# Line 264  err_t Paso_Solver_PCG( Line 284  err_t Paso_Solver_PCG(
284            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);
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))) {
290                 /* smoother */                 /* smoother */
291             sum_3 = 0;             sum_3 = 0;
292             sum_4 = 0;             sum_4 = 0;
293                 #pragma omp parallel private(i0, istart, iend, ipp,d, ss, ss1, alpha)                 #pragma omp parallel private(i0, istart, iend, ipp,d, ss, ss1)
294                 {                 {
295                    ss=0;                    ss=0;
296                    ss1=0;                    ss1=0;
                   alpha=tau/delta;  
297                    #ifdef USE_DYNAMIC_SCHEDULING                    #ifdef USE_DYNAMIC_SCHEDULING
298                        #pragma omp for schedule(dynamic, 1)                        #pragma omp for schedule(dynamic, 1)
299                        for (ipp=0; ipp < n_chunks; ++ipp) {                        for (ipp=0; ipp < n_chunks; ++ipp) {
# Line 286  err_t Paso_Solver_PCG( Line 305  err_t Paso_Solver_PCG(
305                            istart=len*ipp+MIN(ipp,rest);                            istart=len*ipp+MIN(ipp,rest);
306                            iend=len*(ipp+1)+MIN(ipp+1,rest);                            iend=len*(ipp+1)+MIN(ipp+1,rest);
307                    #endif                    #endif
308                    #pragma ivdep                #pragma ivdep
309                    for (i0=istart;i0<iend;i0++) {                            for (i0=istart;i0<iend;i0++) {
310                          r[i0]-=alpha*v[i0];                                  r[i0]-=alpha*v[i0];
311                          d=r[i0]-rs[i0];                                  d=r[i0]-rs[i0];
312                          ss+=d*d;                                  ss+=d*d;
313                          ss1+=d*rs[i0];                                  ss1+=d*rs[i0];
314                    }                            }
315                    #ifdef USE_DYNAMIC_SCHEDULING                    #ifdef USE_DYNAMIC_SCHEDULING
316                      }                      }
317                    #else                    #else
318                      }                      }
319                    #endif                    #endif
320                    #pragma critical                    #pragma omp critical
321                    {                    {
322                        sum_3+=ss;                       sum_3+=ss;
323                        sum_4+=ss1;                       sum_4+=ss1;
324                    }                    }
325                 }                 }
326                 #ifdef PASO_MPI                 #ifdef PASO_MPI
# Line 328  err_t Paso_Solver_PCG( Line 347  err_t Paso_Solver_PCG(
347                            istart=len*ipp+MIN(ipp,rest);                            istart=len*ipp+MIN(ipp,rest);
348                            iend=len*(ipp+1)+MIN(ipp+1,rest);                            iend=len*(ipp+1)+MIN(ipp+1,rest);
349                    #endif                    #endif
350                    #pragma ivdep                #pragma ivdep
351                    for (i0=istart;i0<iend;i0++) {                            for (i0=istart;i0<iend;i0++) {
352                        rs[i0]=gamma_2*rs[i0]+gamma_1*r[i0];                                rs[i0]=gamma_2*rs[i0]+gamma_1*r[i0];
353                        x2[i0]+=alpha*p[i0];                                x2[i0]+=alpha*p[i0];
354                        x[i0]=gamma_2*x[i0]+gamma_1*x2[i0];                                x[i0]=gamma_2*x[i0]+gamma_1*x2[i0];
355                        ss+=rs[i0]*rs[i0];                                ss+=rs[i0]*rs[i0];
356                    }                            }
357                    #ifdef USE_DYNAMIC_SCHEDULING                    #ifdef USE_DYNAMIC_SCHEDULING
358                      }                      }
359                    #else                    #else
# Line 342  err_t Paso_Solver_PCG( Line 361  err_t Paso_Solver_PCG(
361                    #endif                    #endif
362                    #pragma omp critical                    #pragma omp critical
363                    {                    {
364                       sum_5+=ss;                        sum_5+=ss;
365                    }                    }
366                  }                  }
367                  #ifdef PASO_MPI                  #ifdef PASO_MPI

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

  ViewVC Help
Powered by ViewVC 1.1.26