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

trunk/paso/src/Solvers/PCG.c revision 495 by gross, Mon Feb 6 06:32:06 2006 UTC branches/doubleplusgood/paso/src/PCG.cpp revision 4261 by jfenwick, Wed Feb 27 06:09:33 2013 UTC
# Line 1  Line 1
1  /* \$Id\$ */
2    /*****************************************************************************
3    *
4    * Copyright (c) 2003-2013 by University of Queensland
5    * http://www.uq.edu.au
6    *
7    * Primary Business: Queensland, Australia
10    *
11    * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15
16
17  /* PCG iterations */  /* PCG iterations */
# Line 6  Line 19
19  #include "SystemMatrix.h"  #include "SystemMatrix.h"
20  #include "Paso.h"  #include "Paso.h"
21  #include "Solver.h"  #include "Solver.h"
22  /* #include <math.h> */
23  #ifdef _OPENMP  #ifdef _OPENMP
24  #include <omp.h>  #include <omp.h>
25  #endif  #endif
26
27    #ifdef ESYS_MPI
28    #include <mpi.h>
29    #endif
30
31  /*  /*
32  *  *
33  *  Purpose  *  Purpose
34  *  =======  *  =======
35  *  *
36  *  PCG solves the linear system A*x = b using the  *  PCG solves the linear system A*x = b using the
37  *  preconditioned conjugate gradient method plus a smoother  *  preconditioned conjugate gradient method plus a smoother.
38  *  A has to be symmetric.  *  A has to be symmetric.
39  *  *
40  *  Convergence test: norm( b - A*x )< TOL.  *  Convergence test: norm( b - A*x )< TOL.
# Line 27  Line 44
44  *  =========  *  =========
45  *  *
46  *  r       (input) DOUBLE PRECISION array, dimension N.  *  r       (input) DOUBLE PRECISION array, dimension N.
47  *          On entry, residual of inital guess x  *          On entry, residual of initial guess x.
48  *  *
49  *  x       (input/output) DOUBLE PRECISION array, dimension N.  *  x       (input/output) DOUBLE PRECISION array, dimension N.
50  *          On input, the initial guess.  *          On input, the initial guess.
# Line 39  Line 56
56  *  INFO    (output) INT  *  INFO    (output) INT
57  *  *
58  *          = SOLVER_NO_ERROR: Successful exit. Iterated approximate solution returned.  *          = SOLVER_NO_ERROR: Successful exit. Iterated approximate solution returned.
59  *          = SOLVEr_MAXITER_REACHED  *          = SOLVER_MAXITER_REACHED
60  *          = SOLVER_INPUT_ERROR Illegal parameter:  *          = SOLVER_INPUT_ERROR Illegal parameter:
61  *          = SOLVEr_BREAKDOWN: If parameters rHO or OMEGA become smaller  *          = SOLVER_BREAKDOWN: If parameters RHO or OMEGA become smaller
62  *          = SOLVER_MEMORY_ERROR : If parameters rHO or OMEGA become smaller  *          = SOLVER_MEMORY_ERROR : If parameters RHO or OMEGA become smaller
63  *  *
64  *  ==============================================================  *  ==============================================================
65  */  */
66
67    /* #define PASO_DYNAMIC_SCHEDULING_MVM */
68
69    #if defined PASO_DYNAMIC_SCHEDULING_MVM && defined __OPENMP
70    #define USE_DYNAMIC_SCHEDULING
71    #endif
72
73  err_t Paso_Solver_PCG(  err_t Paso_Solver_PCG(
74      Paso_SystemMatrix * A,      Paso_SystemMatrix * A,
75      double * r,      double * r,
76      double * x,      double * x,
77      dim_t *iter,      dim_t *iter,
78      double * tolerance) {      double * tolerance,
79        Paso_Performance* pp) {
80
81    /* Local variables */    /* Local variables */
82    dim_t num_iter=0,maxit,num_iter_global;    dim_t num_iter=0,maxit,num_iter_global, len,rest, np, ipp;
83    dim_t i0;  #ifdef USE_DYNAMIC_SCHEDULING
84      dim_t chunk_size=-1;
85    #endif
86      register double ss,ss1;
87      dim_t i0, istart, iend;
88    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;
89    err_t status = SOLVER_NO_ERROR;    err_t status = SOLVER_NO_ERROR;
90    dim_t n = A->num_cols * A-> col_block_size;    dim_t n = Paso_SystemMatrix_getTotalNumRows(A);
91    double *resid = tolerance, *rs=NULL, *p=NULL, *v=NULL, *x2=NULL ;    double *resid = tolerance, *rs=NULL, *p=NULL, *v=NULL, *x2=NULL ;
92    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;
93    double norm_of_residual,norm_of_residual_global;  #ifdef ESYS_MPI
94    register double r_tmp,d,rs_tmp,x2_tmp,x_tmp;    double loc_sum[2], sum[2];
95    #endif
96      double norm_of_residual=0,norm_of_residual_global;
97      register double d;
98
99      /* There should not be any executable code before this ifdef */
100
101    #ifdef USE_DYNAMIC_SCHEDULING
102
103        /* Watch out for these declarations (see above) */
104        char* chksz_chr;
105        dim_t n_chunks;
106
107        chksz_chr=getenv("PASO_CHUNK_SIZE_PCG");
108        if (chksz_chr!=NULL) sscanf(chksz_chr, "%d",&chunk_size);
110        chunk_size=MIN(MAX(1,chunk_size),n/np);
111        n_chunks=n/chunk_size;
112        if (n_chunks*chunk_size<n) n_chunks+=1;
113    #else
115        len=n/np;
116        rest=n-len*np;
117    #endif
118  /*                                                                 */  /*                                                                 */
119  /*-----------------------------------------------------------------*/  /*-----------------------------------------------------------------*/
120  /*                                                                 */  /*                                                                 */
# Line 87  err_t Paso_Solver_PCG( Line 136  err_t Paso_Solver_PCG(
136    } else {    } else {
137      maxit = *iter;      maxit = *iter;
138      tol = *resid;      tol = *resid;
139      #pragma omp parallel firstprivate(maxit,tol,convergeFlag,maxIterFlag,breakFlag) \      Performance_startMonitor(pp,PERFORMANCE_SOLVER);
140                                             private(tau_old,tau,beta,delta,gamma_1,gamma_2,alpha,norm_of_residual,num_iter)      /* initialize data */
141        #pragma omp parallel private(i0, istart, iend, ipp)
142      {      {
143         /* initialize data */         #ifdef USE_DYNAMIC_SCHEDULING
144         #pragma omp for private(i0) schedule(static)             #pragma omp for schedule(dynamic, 1)
145         for (i0=0;i0<n;i0++) {             for (ipp=0; ipp < n_chunks; ++ipp) {
146            rs[i0]=r[i0];                istart=chunk_size*ipp;
147            x2[i0]=x[i0];                iend=MIN(istart+chunk_size,n);
148         }         #else
149         #pragma omp for private(i0) schedule(static)             #pragma omp for schedule(static)
150         for (i0=0;i0<n;i0++) {             for (ipp=0; ipp <np; ++ipp) {
151            p[i0]=0;                 istart=len*ipp+MIN(ipp,rest);
152            v[i0]=0;                 iend=len*(ipp+1)+MIN(ipp+1,rest);
153         }         #endif
154         num_iter=0;                 #pragma ivdep
155         /* start of iteration */                 for (i0=istart;i0<iend;i0++) {
156         while (!(convergeFlag || maxIterFlag || breakFlag)) {                   rs[i0]=r[i0];
157                     x2[i0]=x[i0];
158                     p[i0]=0;
159                     v[i0]=0;
160                   }
161           #ifdef USE_DYNAMIC_SCHEDULING
162             }
163           #else
164             }
165           #endif
166        }
167        num_iter=0;
168
169        /* PGH */
170        /* without this we get a use of an uninitialised var below */
171        tau = 0;
172
173        /* start of iteration */
174        while (!(convergeFlag || maxIterFlag || breakFlag)) {
175             ++(num_iter);             ++(num_iter);
176             #pragma omp barrier
177             #pragma omp master             /* PGH */
178               /* The next lines were commented out before I got here */
179               /* v=prec(r)  */
180               /* tau=v*r; */
181               /* leading to the use of an uninitialised var below */
182
183               Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
184               Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);
185           Paso_SystemMatrix_solvePreconditioner(A,v,r);
186               Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);
187               Performance_startMonitor(pp,PERFORMANCE_SOLVER);
188
189           sum_1 = 0;
190               #pragma omp parallel private(i0, istart, iend, ipp, ss)
191             {             {
192             sum_1 = 0;                    ss=0;
193             sum_2 = 0;                    #ifdef USE_DYNAMIC_SCHEDULING
194             sum_3 = 0;                        #pragma omp for schedule(dynamic, 1)
195             sum_4 = 0;                        for (ipp=0; ipp < n_chunks; ++ipp) {
196             sum_5 = 0;                           istart=chunk_size*ipp;
197                             iend=MIN(istart+chunk_size,n);
198                      #else
199                          #pragma omp for schedule(static)
200                          for (ipp=0; ipp <np; ++ipp) {
201                              istart=len*ipp+MIN(ipp,rest);
202                              iend=len*(ipp+1)+MIN(ipp+1,rest);
203                      #endif
204                  #pragma ivdep
205                              for (i0=istart;i0<iend;i0++) ss+=v[i0]*r[i0];
206                      #ifdef USE_DYNAMIC_SCHEDULING
207                        }
208                      #else
209                        }
210                      #endif
211                      #pragma omp critical
212                      {
213                         sum_1+=ss;
214                      }
215             }             }
216             /* v=prec(r)  */             #ifdef ESYS_MPI
217             Paso_Solver_solvePreconditioner(A,v,r);              /* In case we have many MPI processes, each of which may have several OMP threads:
218             /* tau=v*r    */                 OMP master participates in an MPI reduction to get global sum_1 */
219             #pragma omp for private(i0) reduction(+:sum_1) schedule(static)              loc_sum[0] = sum_1;
220             for (i0=0;i0<n;i0++) sum_1+=v[i0]*r[i0];              MPI_Allreduce(loc_sum, &sum_1, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
221               #endif
222             tau_old=tau;             tau_old=tau;
223             tau=sum_1;             tau=sum_1;
224             /* p=v+beta*p */             /* p=v+beta*p */
225             if (num_iter==1) {             #pragma omp parallel private(i0, istart, iend, ipp,beta)
226                 #pragma omp for private(i0)  schedule(static)             {
227                 for (i0=0;i0<n;i0++) p[i0]=v[i0];                    #ifdef USE_DYNAMIC_SCHEDULING
228             } else {                        #pragma omp for schedule(dynamic, 1)
229                 beta=tau/tau_old;                        for (ipp=0; ipp < n_chunks; ++ipp) {
230                 #pragma omp for private(i0)  schedule(static)                           istart=chunk_size*ipp;
231                 for (i0=0;i0<n;i0++) p[i0]=v[i0]+beta*p[i0];                           iend=MIN(istart+chunk_size,n);
232                      #else
233                          #pragma omp for schedule(static)
234                          for (ipp=0; ipp <np; ++ipp) {
235                              istart=len*ipp+MIN(ipp,rest);
236                              iend=len*(ipp+1)+MIN(ipp+1,rest);
237                      #endif
238                              if (num_iter==1) {
239                      #pragma ivdep
240                                  for (i0=istart;i0<iend;i0++) p[i0]=v[i0];
241                              } else {
242                                  beta=tau/tau_old;
243                      #pragma ivdep
244                                  for (i0=istart;i0<iend;i0++) p[i0]=v[i0]+beta*p[i0];
245                              }
246                      #ifdef USE_DYNAMIC_SCHEDULING
247                        }
248                      #else
249                        }
250                      #endif
251             }             }
252             /* v=A*p */             /* v=A*p */
253         Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, p,ZERO,v);             Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
254               Performance_startMonitor(pp,PERFORMANCE_MVM);
255           Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(PASO_ONE, A, p,PASO_ZERO,v);
256               Performance_stopMonitor(pp,PERFORMANCE_MVM);
257               Performance_startMonitor(pp,PERFORMANCE_SOLVER);
258
259             /* delta=p*v */             /* delta=p*v */
260             #pragma omp for private(i0) reduction(+:sum_2) schedule(static)         sum_2 = 0;
261             for (i0=0;i0<n;i0++) sum_2+=v[i0]*p[i0];             #pragma omp parallel private(i0, istart, iend, ipp,ss)
262               {
263                      ss=0;
264                      #ifdef USE_DYNAMIC_SCHEDULING
265                          #pragma omp for schedule(dynamic, 1)
266                          for (ipp=0; ipp < n_chunks; ++ipp) {
267                             istart=chunk_size*ipp;
268                             iend=MIN(istart+chunk_size,n);
269                      #else
270                          #pragma omp for schedule(static)
271                          for (ipp=0; ipp <np; ++ipp) {
272                              istart=len*ipp+MIN(ipp,rest);
273                              iend=len*(ipp+1)+MIN(ipp+1,rest);
274                      #endif
275                  #pragma ivdep
276                              for (i0=istart;i0<iend;i0++) ss+=v[i0]*p[i0];
277                      #ifdef USE_DYNAMIC_SCHEDULING
278                        }
279                      #else
280                        }
281                      #endif
282                      #pragma omp critical
283                      {
284                                 sum_2+=ss;
285                      }
286               }
287               #ifdef ESYS_MPI
288              loc_sum[0] = sum_2;
289              MPI_Allreduce(loc_sum, &sum_2, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
290               #endif
291             delta=sum_2;             delta=sum_2;
292               alpha=tau/delta;
293
294             if (! (breakFlag = (ABS(delta) <= TOLERANCE_FOR_SCALARS))) {             if (! (breakFlag = (ABS(delta) <= TOLERANCE_FOR_SCALARS))) {
alpha=tau/delta;
295                 /* smoother */                 /* smoother */
296                 #pragma omp for private(i0) schedule(static)             sum_3 = 0;
297                 for (i0=0;i0<n;i0++) r[i0]-=alpha*v[i0];             sum_4 = 0;
298                 #pragma omp for private(i0,d) reduction(+:sum_3,sum_4) schedule(static)                 #pragma omp parallel private(i0, istart, iend, ipp,d, ss, ss1)
299                 for (i0=0;i0<n;i0++) {                 {
300                       d=r[i0]-rs[i0];                    ss=0;
301                       sum_3+=d*d;                    ss1=0;
302                       sum_4+=d*rs[i0];                    #ifdef USE_DYNAMIC_SCHEDULING
303                  }                        #pragma omp for schedule(dynamic, 1)
304                  gamma_1= ( (ABS(sum_3)<= ZERO) ? 0 : -sum_4/sum_3) ;                        for (ipp=0; ipp < n_chunks; ++ipp) {
305                  gamma_2= ONE-gamma_1;                           istart=chunk_size*ipp;
306                  #pragma omp for private(i0,x2_tmp,x_tmp,rs_tmp) schedule(static)                           iend=MIN(istart+chunk_size,n);
307                  for (i0=0;i0<n;++i0) {                    #else
308                    rs[i0]=gamma_2*rs[i0]+gamma_1*r[i0];                        #pragma omp for schedule(static)
309                    x2[i0]+=alpha*p[i0];                        for (ipp=0; ipp <np; ++ipp) {
310                    x[i0]=gamma_2*x[i0]+gamma_1*x2[i0];                            istart=len*ipp+MIN(ipp,rest);
311                              iend=len*(ipp+1)+MIN(ipp+1,rest);
312                      #endif
313                  #pragma ivdep
314                              for (i0=istart;i0<iend;i0++) {
315                                    r[i0]-=alpha*v[i0];
316                                    d=r[i0]-rs[i0];
317                                    ss+=d*d;
318                                    ss1+=d*rs[i0];
319                              }
320                      #ifdef USE_DYNAMIC_SCHEDULING
321                        }
322                      #else
323                        }
324                      #endif
325                      #pragma omp critical
326                      {
327                         sum_3+=ss;
328                         sum_4+=ss1;
329                      }
330                   }
331                   #ifdef ESYS_MPI
332                   loc_sum[0] = sum_3;
333                   loc_sum[1] = sum_4;
334                   MPI_Allreduce(loc_sum, sum, 2, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
335                   sum_3=sum[0];
336                   sum_4=sum[1];
337                    #endif
338                sum_5 = 0;
339                    #pragma omp parallel private(i0, istart, iend, ipp, ss, gamma_1,gamma_2)
340                    {
341                      gamma_1= ( (ABS(sum_3)<= PASO_ZERO) ? 0 : -sum_4/sum_3) ;
342                      gamma_2= PASO_ONE-gamma_1;
343                      ss=0;
344                      #ifdef USE_DYNAMIC_SCHEDULING
345                          #pragma omp for schedule(dynamic, 1)
346                          for (ipp=0; ipp < n_chunks; ++ipp) {
347                             istart=chunk_size*ipp;
348                             iend=MIN(istart+chunk_size,n);
349                      #else
350                          #pragma omp for schedule(static)
351                          for (ipp=0; ipp <np; ++ipp) {
352                              istart=len*ipp+MIN(ipp,rest);
353                              iend=len*(ipp+1)+MIN(ipp+1,rest);
354                      #endif
355                  #pragma ivdep
356                              for (i0=istart;i0<iend;i0++) {
357                                  rs[i0]=gamma_2*rs[i0]+gamma_1*r[i0];
358                                  x2[i0]+=alpha*p[i0];
359                                  x[i0]=gamma_2*x[i0]+gamma_1*x2[i0];
360                                  ss+=rs[i0]*rs[i0];
361                              }
362                      #ifdef USE_DYNAMIC_SCHEDULING
363                        }
364                      #else
365                        }
366                      #endif
367                      #pragma omp critical
368                      {
369                          sum_5+=ss;
370                      }
371                  }                  }
372                  #pragma omp for private(i0) reduction(+:sum_5) schedule(static)                  #ifdef ESYS_MPI
373                  for (i0=0;i0<n;++i0) sum_5+=rs[i0]*rs[i0];                 loc_sum[0] = sum_5;
374                   MPI_Allreduce(loc_sum, &sum_5, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
375                    #endif
376                  norm_of_residual=sqrt(sum_5);                  norm_of_residual=sqrt(sum_5);
377                  convergeFlag = norm_of_residual <= tol;                  convergeFlag = norm_of_residual <= tol;
378                  maxIterFlag = num_iter == maxit;                  maxIterFlag = num_iter > maxit;
379                  breakFlag = (ABS(tau) <= TOLERANCE_FOR_SCALARS);                  breakFlag = (ABS(tau) <= TOLERANCE_FOR_SCALARS);
380             }             }
381         }      }
382         /* end of iteration */      /* end of iteration */
383         #pragma omp master      num_iter_global=num_iter;
384         {      norm_of_residual_global=norm_of_residual;
385             num_iter_global=num_iter;      if (maxIterFlag) {
386             norm_of_residual_global=norm_of_residual;           status = SOLVER_MAXITER_REACHED;
387             if (maxIterFlag) {      } else if (breakFlag) {
388                 status = SOLVER_MAXITER_REACHED;           status = SOLVER_BREAKDOWN;
389             } else if (breakFlag) {      }
390                 status = SOLVER_BREAKDOWN;      Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
}
}
}  /* end of parallel region */
391      TMPMEMFREE(rs);      TMPMEMFREE(rs);
392      TMPMEMFREE(x2);      TMPMEMFREE(x2);
393      TMPMEMFREE(v);      TMPMEMFREE(v);
# Line 188  err_t Paso_Solver_PCG( Line 399  err_t Paso_Solver_PCG(
399    return status;    return status;
400  }  }
401
/*
* \$Log\$
* Revision 1.2  2005/09/15 03:44:40  jgs
* Merge of development branch dev-02 back to main trunk on 2005-09-15
*
* Revision 1.1.2.1  2005/09/05 06:29:49  gross
* These files have been extracted from finley to define a stand alone libray for iterative
* linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
* has not been tested yet.
*
*
*/

Legend:
 Removed from v.495 changed lines Added in v.4261