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

trunk/paso/src/Solvers/PCG.c revision 495 by gross, Mon Feb 6 06:32:06 2006 UTC trunk/paso/src/PCG.c revision 1579 by ksteube, Mon Jun 2 08:48:36 2008 UTC
# Line 1  Line 1 
1    
2  /* $Id$ */  /* $Id$ */
3    
4    /*******************************************************
5     *
6     *           Copyright 2003-2007 by ACceSS MNRF
7     *       Copyright 2007 by University of Queensland
8     *
9     *                http://esscc.uq.edu.au
10     *        Primary Business: Queensland, Australia
11     *  Licensed under the Open Software License version 3.0
12     *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16  /* PCG iterations */  /* PCG iterations */
17    
18  #include "SystemMatrix.h"  #include "SystemMatrix.h"
19  #include "Paso.h"  #include "Paso.h"
20  #include "Solver.h"  #include "Solver.h"
21  /* #include <math.h> */  
22  #ifdef _OPENMP  #ifdef _OPENMP
23  #include <omp.h>  #include <omp.h>
24  #endif  #endif
25    
26    #ifdef PASO_MPI
27    #include <mpi.h>
28    #endif
29    
30  /*  /*
31  *  *
32  *  Purpose  *  Purpose
# Line 47  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,
75      double * x,      double * x,
76      dim_t *iter,      dim_t *iter,
77      double * tolerance) {      double * tolerance,
78        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 = A->num_cols * A-> col_block_size;    dim_t n = Paso_SystemMatrix_getTotalNumRows(A);
87    double *resid = tolerance, *rs=NULL, *p=NULL, *v=NULL, *x2=NULL ;    double *resid = tolerance, *rs=NULL, *p=NULL, *v=NULL, *x2=NULL ;
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;    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 87  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         /* initialize data */         #ifdef USE_DYNAMIC_SCHEDULING
130         #pragma omp for private(i0) schedule(static)             #pragma omp for schedule(dynamic, 1)
131         for (i0=0;i0<n;i0++) {             for (ipp=0; ipp < n_chunks; ++ipp) {
132            rs[i0]=r[i0];                istart=chunk_size*ipp;
133            x2[i0]=x[i0];                iend=MIN(istart+chunk_size,n);
134         }         #else
135         #pragma omp for private(i0) schedule(static)             #pragma omp for schedule(static)
136         for (i0=0;i0<n;i0++) {             for (ipp=0; ipp <np; ++ipp) {
137            p[i0]=0;                 istart=len*ipp+MIN(ipp,rest);
138            v[i0]=0;                 iend=len*(ipp+1)+MIN(ipp+1,rest);
139         }         #endif
140         num_iter=0;                 #pragma ivdep
141         /* start of iteration */                 for (i0=istart;i0<iend;i0++) {
142         while (!(convergeFlag || maxIterFlag || breakFlag)) {                   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;
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);
159               Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);
160             Paso_Solver_solvePreconditioner(A,v,r);             Paso_Solver_solvePreconditioner(A,v,r);
161               Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);
162               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];             #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 omp critical
187                      {
188                         sum_1+=ss;
189                      }
190               }
191               #ifdef PASO_MPI
192                /* 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 */
194                loc_sum[0] = sum_1;
195                MPI_Allreduce(loc_sum, &sum_1, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
196               #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);
229               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);
231               Performance_stopMonitor(pp,PERFORMANCE_MVM);
232               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 omp critical
258                      {
259                                 sum_2+=ss;
260                      }
261               }
262               #ifdef PASO_MPI
263              loc_sum[0] = sum_2;
264              MPI_Allreduce(loc_sum, &sum_2, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
265               #endif
266             delta=sum_2;             delta=sum_2;
267               alpha=tau/delta;
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)
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];                    #ifdef USE_DYNAMIC_SCHEDULING
278                  }                        #pragma omp for schedule(dynamic, 1)
279                  gamma_1= ( (ABS(sum_3)<= ZERO) ? 0 : -sum_4/sum_3) ;                        for (ipp=0; ipp < n_chunks; ++ipp) {
280                  gamma_2= ONE-gamma_1;                           istart=chunk_size*ipp;
281                  #pragma omp for private(i0,x2_tmp,x_tmp,rs_tmp) schedule(static)                           iend=MIN(istart+chunk_size,n);
282                  for (i0=0;i0<n;++i0) {                    #else
283                    rs[i0]=gamma_2*rs[i0]+gamma_1*r[i0];                        #pragma omp for schedule(static)
284                    x2[i0]+=alpha*p[i0];                        for (ipp=0; ipp <np; ++ipp) {
285                    x[i0]=gamma_2*x[i0]+gamma_1*x2[i0];                            istart=len*ipp+MIN(ipp,rest);
286                              iend=len*(ipp+1)+MIN(ipp+1,rest);
287                      #endif
288                  #pragma ivdep
289                              for (i0=istart;i0<iend;i0++) {
290                                    r[i0]-=alpha*v[i0];
291                                    d=r[i0]-rs[i0];
292                                    ss+=d*d;
293                                    ss1+=d*rs[i0];
294                              }
295                      #ifdef USE_DYNAMIC_SCHEDULING
296                        }
297                      #else
298                        }
299                      #endif
300                      #pragma omp critical
301                      {
302                         sum_3+=ss;
303                         sum_4+=ss1;
304                      }
305                   }
306                   #ifdef PASO_MPI
307                   loc_sum[0] = sum_3;
308                   loc_sum[1] = sum_4;
309                   MPI_Allreduce(loc_sum, sum, 2, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
310                   sum_3=sum[0];
311                   sum_4=sum[1];
312                    #endif
313                sum_5 = 0;
314                    #pragma omp parallel private(i0, istart, iend, ipp, ss, gamma_1,gamma_2)
315                    {
316                      gamma_1= ( (ABS(sum_3)<= ZERO) ? 0 : -sum_4/sum_3) ;
317                      gamma_2= ONE-gamma_1;
318                      ss=0;
319                      #ifdef USE_DYNAMIC_SCHEDULING
320                          #pragma omp for schedule(dynamic, 1)
321                          for (ipp=0; ipp < n_chunks; ++ipp) {
322                             istart=chunk_size*ipp;
323                             iend=MIN(istart+chunk_size,n);
324                      #else
325                          #pragma omp for schedule(static)
326                          for (ipp=0; ipp <np; ++ipp) {
327                              istart=len*ipp+MIN(ipp,rest);
328                              iend=len*(ipp+1)+MIN(ipp+1,rest);
329                      #endif
330                  #pragma ivdep
331                              for (i0=istart;i0<iend;i0++) {
332                                  rs[i0]=gamma_2*rs[i0]+gamma_1*r[i0];
333                                  x2[i0]+=alpha*p[i0];
334                                  x[i0]=gamma_2*x[i0]+gamma_1*x2[i0];
335                                  ss+=rs[i0]*rs[i0];
336                              }
337                      #ifdef USE_DYNAMIC_SCHEDULING
338                        }
339                      #else
340                        }
341                      #endif
342                      #pragma omp critical
343                      {
344                          sum_5+=ss;
345                      }
346                  }                  }
347                  #pragma omp for private(i0) reduction(+:sum_5) schedule(static)                  #ifdef PASO_MPI
348                  for (i0=0;i0<n;++i0) sum_5+=rs[i0]*rs[i0];                 loc_sum[0] = sum_5;
349                   MPI_Allreduce(loc_sum, &sum_5, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
350                    #endif
351                  norm_of_residual=sqrt(sum_5);                  norm_of_residual=sqrt(sum_5);
352                  convergeFlag = norm_of_residual <= tol;                  convergeFlag = norm_of_residual <= tol;
353                  maxIterFlag = num_iter == maxit;                  maxIterFlag = num_iter == maxit;
354                  breakFlag = (ABS(tau) <= TOLERANCE_FOR_SCALARS);                  breakFlag = (ABS(tau) <= TOLERANCE_FOR_SCALARS);
355             }             }
356         }      }
357         /* end of iteration */      /* end of iteration */
358         #pragma omp master      num_iter_global=num_iter;
359         {      norm_of_residual_global=norm_of_residual;
360             num_iter_global=num_iter;      if (maxIterFlag) {
361             norm_of_residual_global=norm_of_residual;           status = SOLVER_MAXITER_REACHED;
362             if (maxIterFlag) {      } else if (breakFlag) {
363                 status = SOLVER_MAXITER_REACHED;           status = SOLVER_BREAKDOWN;
364             } else if (breakFlag) {      }
365                 status = SOLVER_BREAKDOWN;      Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
            }  
          }  
     }  /* end of parallel region */  
366      TMPMEMFREE(rs);      TMPMEMFREE(rs);
367      TMPMEMFREE(x2);      TMPMEMFREE(x2);
368      TMPMEMFREE(v);      TMPMEMFREE(v);
# Line 187  err_t Paso_Solver_PCG( Line 373  err_t Paso_Solver_PCG(
373    /*     End of PCG */    /*     End of PCG */
374    return status;    return status;
375  }  }
   
 /*  
  * $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.1579

  ViewVC Help
Powered by ViewVC 1.1.26