/[escript]/trunk/paso/src/Solver.c
ViewVC logotype

Diff of /trunk/paso/src/Solver.c

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

revision 700 by gross, Thu Apr 6 00:13:40 2006 UTC revision 1312 by ksteube, Mon Sep 24 06:18:44 2007 UTC
# Line 1  Line 1 
 /* $Id$ */  
1    
2    /* $Id$ */
3    
4  /*  /*******************************************************
5  ********************************************************************************   *
6  *               Copyright   2006 by ACcESS MNRF                                *   *           Copyright 2003-2007 by ACceSS MNRF
7  *                                                                              *   *       Copyright 2007 by University of Queensland
8  *                 http://www.access.edu.au                                     *   *
9  *           Primary Business: Queensland, Australia                            *   *                http://esscc.uq.edu.au
10  *     Licensed under the Open Software License version 3.0             *   *        Primary Business: Queensland, Australia
11  *        http://www.opensource.org/licenses/osl-3.0.php                        *   *  Licensed under the Open Software License version 3.0
12  ********************************************************************************   *     http://www.opensource.org/licenses/osl-3.0.php
13  */   *
14     *******************************************************/
15    
16  /**************************************************************/  /**************************************************************/
17    
# Line 26  Line 27 
27  #include "Paso.h"  #include "Paso.h"
28  #include "SystemMatrix.h"  #include "SystemMatrix.h"
29  #include "Solver.h"  #include "Solver.h"
30    #include "escript/blocktimer.h"
31    
32  /***********************************************************************************/  /***********************************************************************************/
33    
# Line 40  void Paso_Solver_free(Paso_SystemMatrix* Line 42  void Paso_Solver_free(Paso_SystemMatrix*
42  void Paso_Solver(Paso_SystemMatrix* A,double* x,double* b,Paso_Options* options,Paso_Performance* pp) {  void Paso_Solver(Paso_SystemMatrix* A,double* x,double* b,Paso_Options* options,Paso_Performance* pp) {
43      double norm2_of_b,tol,tolerance,time_iter,*r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b,      double norm2_of_b,tol,tolerance,time_iter,*r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b,
44             norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local,norm_max_of_residual_local,norm_max_of_residual,             norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local,norm_max_of_residual_local,norm_max_of_residual,
45             last_norm_max_of_residual,*scaling;             last_norm_max_of_residual,*scaling, loc_norm;
46       dim_t i,totIter,cntIter,method;       dim_t i,totIter,cntIter,method;
47       bool_t finalizeIteration;       bool_t finalizeIteration;
48       err_t errorCode;       err_t errorCode;
49       dim_t n_col = A->num_cols * A-> col_block_size;       dim_t numSol = Paso_SystemMatrix_getTotalNumCols(A);
50       dim_t n_row = A->num_rows * A-> row_block_size;       dim_t numEqua = Paso_SystemMatrix_getTotalNumRows(A);
51         double blocktimer_start = blocktimer_time();
52    
53       tolerance=MAX(options->tolerance,EPSILON);       tolerance=MAX(options->tolerance,EPSILON);
54       Paso_resetError();       Paso_resetError();
# Line 57  void Paso_Solver(Paso_SystemMatrix* A,do Line 60  void Paso_Solver(Paso_SystemMatrix* A,do
60       if (A->col_block_size != A->row_block_size) {       if (A->col_block_size != A->row_block_size) {
61          Paso_setError(TYPE_ERROR,"Iterative solver requires row and column block sizes to be equal.");          Paso_setError(TYPE_ERROR,"Iterative solver requires row and column block sizes to be equal.");
62       }       }
63       if (A->num_cols != A->num_rows) {       if (Paso_SystemMatrix_getGlobalNumCols(A) != Paso_SystemMatrix_getGlobalNumRows(A)) {
64          Paso_setError(TYPE_ERROR,"Iterative solver requires requires a square matrix.");          Paso_setError(TYPE_ERROR,"Iterative solver requires requires a square matrix.");
65          return;          return;
66       }       }
# Line 74  void Paso_Solver(Paso_SystemMatrix* A,do Line 77  void Paso_Solver(Paso_SystemMatrix* A,do
77                  norm2_of_b_local=0.;                  norm2_of_b_local=0.;
78                  norm_max_of_b_local=0.;                  norm_max_of_b_local=0.;
79                  #pragma omp for private(i) schedule(static)                  #pragma omp for private(i) schedule(static)
80                  for (i = 0; i < n_row ; ++i) {                  for (i = 0; i < numEqua ; ++i) {
81                        norm2_of_b_local += b[i] * b[i];                        norm2_of_b_local += b[i] * b[i];
82                        norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local);                        norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local);
83                  }                  }
# Line 84  void Paso_Solver(Paso_SystemMatrix* A,do Line 87  void Paso_Solver(Paso_SystemMatrix* A,do
87                     norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b);                     norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b);
88                  }                  }
89             }             }
90               /* TODO: use one call */
91               #ifdef PASO_MPI
92               {
93                   loc_norm = norm2_of_b;
94                   MPI_Allreduce(&loc_norm,&norm2_of_b, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
95                   loc_norm = norm_max_of_b;
96                   MPI_Allreduce(&loc_norm,&norm_max_of_b, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
97               }
98               #endif
99             norm2_of_b=sqrt(norm2_of_b);             norm2_of_b=sqrt(norm2_of_b);
100            
101             /* if norm2_of_b==0 we are ready: x=0 */             /* if norm2_of_b==0 we are ready: x=0 */
102             if (norm2_of_b <=0.) {             if (norm2_of_b <=0.) {
103                 #pragma omp parallel for private(i) schedule(static)                 #pragma omp parallel for private(i) schedule(static)
104                 for (i = 0; i < n_col; i++) x[i]=0.;                 for (i = 0; i < numSol; i++) x[i]=0.;
105                 if (options->verbose) printf("right hand side is identical zero.\n");                 if (options->verbose) printf("right hand side is identical zero.\n");
106             } else {             } else {
107                if (options->verbose) {                if (options->verbose) {
# Line 121  void Paso_Solver(Paso_SystemMatrix* A,do Line 133  void Paso_Solver(Paso_SystemMatrix* A,do
133                Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);                Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
134                if (! Paso_noError()) return;                if (! Paso_noError()) return;
135                        
               /* get an initial guess by evaluating the preconditioner */  
136                time_iter=Paso_timer();                time_iter=Paso_timer();
137                  Paso_SystemMatrix_allocBuffer(A);
138                  /* get an initial guess by evaluating the preconditioner */
139                #pragma omp parallel                #pragma omp parallel
140                Paso_Solver_solvePreconditioner(A,x,b);                {
141                if (! Paso_noError()) return;                   Paso_Solver_solvePreconditioner(A,x,b);
142                  }
143                /* start the iteration process :*/                /* start the iteration process :*/
144                r=TMPMEMALLOC(n_row,double);                r=TMPMEMALLOC(numEqua,double);
145                if (Paso_checkPtr(r)) return;                Paso_checkPtr(r);
146                totIter = 0;                if (Paso_noError()) {
147                finalizeIteration = FALSE;                   totIter = 0;
148                last_norm2_of_residual=norm2_of_b;                   finalizeIteration = FALSE;
149                last_norm_max_of_residual=norm_max_of_b;                   last_norm2_of_residual=norm2_of_b;
150                while (! finalizeIteration) {                   last_norm_max_of_residual=norm_max_of_b;
151                   cntIter = options->iter_max - totIter;                   while (! finalizeIteration) {
152                   finalizeIteration = TRUE;                      cntIter = options->iter_max - totIter;
153                   /*     Set initial residual. */                      finalizeIteration = TRUE;
154                   norm2_of_residual = 0;                      /*     Set initial residual. */
155                   norm_max_of_residual = 0;                      norm2_of_residual = 0;
156                   #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)                      norm_max_of_residual = 0;
157                   {                      #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
158                      #pragma omp for private(i) schedule(static)                      {
159                      for (i = 0; i < n_row; i++) r[i]=b[i];                         #pragma omp for private(i) schedule(static)
160                                     for (i = 0; i < numEqua; i++) {
161                      Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);                             r[i]=b[i];
162                                     }
163                      norm2_of_residual_local = 0;              
164                      norm_max_of_residual_local = 0;                         Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);
165                      #pragma omp for private(i) schedule(static)              
166                      for (i = 0; i < n_row; i++) {                         norm2_of_residual_local = 0;
167                             norm2_of_residual_local+= r[i] * r[i];                         norm_max_of_residual_local = 0;
168                             norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);                         #pragma omp for private(i) schedule(static)
169                           for (i = 0; i < numEqua; i++) {
170                                  norm2_of_residual_local+= r[i] * r[i];
171                                  norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);
172                           }
173                           #pragma omp critical
174                           {
175                              norm2_of_residual += norm2_of_residual_local;
176                              norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);
177                           }
178                      }                      }
179                      #pragma omp critical                      /* TODO: use one call */
180                        #ifdef PASO_MPI
181                      {                      {
182                         norm2_of_residual += norm2_of_residual_local;                          loc_norm = norm2_of_residual;
183                         norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);                          MPI_Allreduce(&loc_norm,&norm2_of_residual, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
184                            loc_norm = norm_max_of_residual;
185                            MPI_Allreduce(&loc_norm,&norm_max_of_residual, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
186                      }                      }
187                   }                      #endif
188                   norm2_of_residual =sqrt(norm2_of_residual);                      norm2_of_residual =sqrt(norm2_of_residual);
189                          
190                   if (options->verbose) printf("Step %5d: l2/lmax-norm of residual is  %e/%e",totIter,norm2_of_residual,norm_max_of_residual);                      if (options->verbose) printf("Step %5d: l2/lmax-norm of residual is  %e/%e",totIter,norm2_of_residual,norm_max_of_residual);
191                   if (totIter>0 && norm2_of_residual>=last_norm2_of_residual &&  norm_max_of_residual>=last_norm_max_of_residual) {                      if (totIter>0 && norm2_of_residual>=last_norm2_of_residual &&  norm_max_of_residual>=last_norm_max_of_residual) {
192                        if (options->verbose) printf(" divergence!\n");                           if (options->verbose) printf(" divergence!\n");
193                        Paso_setError(WARNING, "No improvement during iteration. Iterative solver gives up.");                           Paso_setError(WARNING, "No improvement during iteration. Iterative solver gives up.");
194                   } else {                      } else {
195                      /* */                         /* */
196                      if (norm2_of_residual>tolerance*norm2_of_b || norm_max_of_residual>tolerance*norm_max_of_b ) {                         if (norm2_of_residual>tolerance*norm2_of_b || norm_max_of_residual>tolerance*norm_max_of_b ) {
197                          tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b);                             tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b);
198                         if (options->verbose) printf(" (new tolerance = %e).\n",tol);                            if (options->verbose) printf(" (new tolerance = %e).\n",tol);
199                         last_norm2_of_residual=norm2_of_residual;                            last_norm2_of_residual=norm2_of_residual;
200                         last_norm_max_of_residual=norm_max_of_residual;                            last_norm_max_of_residual=norm_max_of_residual;
201                         /* call the solver */                            /* call the solver */
202                         switch (method) {                            switch (method) {
203                                case PASO_BICGSTAB:                                   case PASO_BICGSTAB:
204                                    errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol,pp);                                       errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol, pp);
205                                    break;                                       break;
206                                case PASO_PCG:                                   case PASO_PCG:
207                                    errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol,pp);                                       errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol, pp);
208                                    break;                                       break;
209                                case PASO_PRES20:                                   case PASO_PRES20:
210                                    errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20,pp);                                       errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20, pp);
211                                    break;                                       break;
212                                case PASO_GMRES:                                   case PASO_GMRES:
213                                    errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart,pp);                                       errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart, pp);
214                                    break;                                       break;
215                          }                             }
216                          totIter += cntIter;                             totIter += cntIter;
217                                            
218                          /* error handling  */                             /* error handling  */
219                          if (errorCode==NO_ERROR) {                             if (errorCode==NO_ERROR) {
220                                finalizeIteration = FALSE;                                   finalizeIteration = FALSE;
221                           } else if (errorCode==SOLVER_MAXITER_REACHED) {                              } else if (errorCode==SOLVER_MAXITER_REACHED) {
222                                 Paso_setError(WARNING,"maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");                                    Paso_setError(WARNING,"maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");
223                                 if (options->verbose) printf("Maximum number of iterations reached.!\n");                                    if (options->verbose) printf("Maximum number of iterations reached.!\n");
224                          } else if (errorCode == SOLVER_INPUT_ERROR ) {                             } else if (errorCode == SOLVER_INPUT_ERROR ) {
225                                 Paso_setError(SYSTEM_ERROR,"illegal dimension in iterative solver.");                                    Paso_setError(SYSTEM_ERROR,"illegal dimension in iterative solver.");
226                                 if (options->verbose) printf("Internal error!\n");                                    if (options->verbose) printf("Internal error!\n");
227                          } else if ( errorCode == SOLVER_BREAKDOWN ) {                             } else if ( errorCode == SOLVER_BREAKDOWN ) {
228                           if (cntIter <= 1) {                              if (cntIter <= 1) {
229                                   Paso_setError(ZERO_DIVISION_ERROR, "fatal break down in iterative solver.");                                      Paso_setError(ZERO_DIVISION_ERROR, "fatal break down in iterative solver.");
230                                   if (options->verbose) printf("Uncurable break down!\n");                                      if (options->verbose) printf("Uncurable break down!\n");
231                              } else {                                 } else {
232                               if (options->verbose) printf("Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);                                  if (options->verbose) printf("Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);
233                                  finalizeIteration = FALSE;                                     finalizeIteration = FALSE;
234                              }                                 }
235                          } else if (errorCode == SOLVER_MEMORY_ERROR) {                             } else if (errorCode == SOLVER_MEMORY_ERROR) {
236                            Paso_setError(MEMORY_ERROR,"memory allocation failed.");                               Paso_setError(MEMORY_ERROR,"memory allocation failed.");
237                               if (options->verbose) printf("Memory allocation failed!\n");                                  if (options->verbose) printf("Memory allocation failed!\n");
238                          } else if (errorCode !=SOLVER_NO_ERROR ) {                             } else if (errorCode !=SOLVER_NO_ERROR ) {
239                            Paso_setError(SYSTEM_ERROR,"unidentified error in iterative solver.");                               Paso_setError(SYSTEM_ERROR,"unidentified error in iterative solver.");
240                               if (options->verbose) printf("Unidentified error!\n");                                  if (options->verbose) printf("Unidentified error!\n");
241                          }                             }
242                     } else {                        } else {
243                        if (options->verbose) printf(". convergence! \n");                           if (options->verbose) printf(". convergence! \n");
244                          }
245                     }                     }
246                  }                   } /* while */
247                } /* while */                }
248                MEMFREE(r);                MEMFREE(r);
249                  Paso_SystemMatrix_freeBuffer(A);
250                time_iter=Paso_timer()-time_iter;                time_iter=Paso_timer()-time_iter;
251                if (options->verbose)  {                if (options->verbose)  {
252                   printf("timing: solver: %.4e sec\n",time_iter);                   printf("timing: solver: %.4e sec\n",time_iter);
# Line 228  void Paso_Solver(Paso_SystemMatrix* A,do Line 256  void Paso_Solver(Paso_SystemMatrix* A,do
256          }          }
257        }        }
258        Performance_stopMonitor(pp,PERFORMANCE_ALL);        Performance_stopMonitor(pp,PERFORMANCE_ALL);
259          blocktimer_increment("Paso_Solver()", blocktimer_start);
260  }  }

Legend:
Removed from v.700  
changed lines
  Added in v.1312

  ViewVC Help
Powered by ViewVC 1.1.26