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

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

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

revision 4855 by caltinay, Tue Apr 8 23:55:38 2014 UTC revision 4856 by caltinay, Wed Apr 9 06:46:55 2014 UTC
# Line 30  Line 30 
30  #include "Solver.h"  #include "Solver.h"
31  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
32    
33    namespace paso {
34    
35  void Paso_Solver_free(paso::SystemMatrix* A)  void Solver_free(SystemMatrix* A)
36  {  {
37      A->freePreconditioner();      A->freePreconditioner();
38  }  }
39    
40  /*  call the iterative solver: */  ///  calls the iterative solver
41  void Paso_Solver(paso::SystemMatrix_ptr A, double* x, double* b,  void Solver(SystemMatrix_ptr A, double* x, double* b, Options* options,
42                   paso::Options* options, Paso_Performance* pp)              Paso_Performance* pp)
43  {  {
44     double norm2_of_b,tol,tolerance,time_iter,net_time_start;      double norm2_of_b,tol,tolerance,time_iter,net_time_start;
45     double *r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b;      double *r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b;
46     double norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local;      double norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local;
47     double norm_max_of_residual_local,norm_max_of_residual;      double norm_max_of_residual_local,norm_max_of_residual;
48     double last_norm_max_of_residual;      double last_norm_max_of_residual;
49  #ifdef ESYS_MPI  #ifdef ESYS_MPI
50     double loc_norm;      double loc_norm;
51  #endif  #endif
52     dim_t i,totIter=0,cntIter,method;      dim_t i,totIter=0,cntIter,method;
53     bool finalizeIteration;      bool finalizeIteration;
54     err_t errorCode=SOLVER_NO_ERROR;      err_t errorCode=SOLVER_NO_ERROR;
55     const dim_t numSol = A->getTotalNumCols();      const dim_t numSol = A->getTotalNumCols();
56     const dim_t numEqua = A->getTotalNumRows();      const dim_t numEqua = A->getTotalNumRows();
57     double blocktimer_precond, blocktimer_start = blocktimer_time();      double blocktimer_precond, blocktimer_start = blocktimer_time();
58     double *x0=NULL;      double *x0=NULL;
59    
60      Esys_resetError();      Esys_resetError();
61      tolerance=options->tolerance;      tolerance=options->tolerance;
62      if (tolerance < 100.* EPSILON) {      if (tolerance < 100.* EPSILON) {
63         Esys_setError(VALUE_ERROR,"Paso_Solver: Tolerance is too small.");          Esys_setError(VALUE_ERROR,"Paso_Solver: Tolerance is too small.");
64      }      }
65      if (tolerance >1.) {      if (tolerance >1.) {
66         Esys_setError(VALUE_ERROR,"Paso_Solver: Tolerance must be less than one.");          Esys_setError(VALUE_ERROR,"Paso_Solver: Tolerance must be less than one.");
67      }      }
68      method=paso::Options::getSolver(options->method, PASO_PASO, options->symmetric, A->mpi_info);      method=Options::getSolver(options->method, PASO_PASO, options->symmetric, A->mpi_info);
69      /* check matrix type */      /* check matrix type */
70      if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) ) {      if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) ) {
71         Esys_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0.");          Esys_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0.");
72      }      }
73      if (A->col_block_size != A->row_block_size) {      if (A->col_block_size != A->row_block_size) {
74          Esys_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires row and column block sizes to be equal.");          Esys_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires row and column block sizes to be equal.");
# Line 79  void Paso_Solver(paso::SystemMatrix_ptr Line 80  void Paso_Solver(paso::SystemMatrix_ptr
80      time_iter=Esys_timer();      time_iter=Esys_timer();
81      /* this for testing only */      /* this for testing only */
82      if (method==PASO_NONLINEAR_GMRES) {      if (method==PASO_NONLINEAR_GMRES) {
83          Paso_Function_LinearSystem* F =          LinearSystem* F = new LinearSystem(A, b, options);
             Paso_Function_LinearSystem_alloc(A, b, options);  
84          A->solvePreconditioner(x, b);          A->solvePreconditioner(x, b);
85          errorCode=Paso_Solver_NewtonGMRES(F, x, options, pp);          errorCode = Solver_NewtonGMRES(F, x, options, pp);
86          if (errorCode != NO_ERROR) {          if (errorCode != NO_ERROR) {
87             Esys_setError(SYSTEM_ERROR,"Paso_Solver_NewtonGMRES: an error has occurred.");              Esys_setError(SYSTEM_ERROR,"Solver_NewtonGMRES: an error has occurred.");
88          }          }
89          Paso_Function_LinearSystem_free(F);          delete F;
90          return;          return;
91      }      }
92            
93      r=new double[numEqua];      r = new double[numEqua];
94      x0=new double[numEqua];      x0 = new double[numEqua];
95      A->balance();      A->balance();
96      options->num_level=0;      options->num_level=0;
97      options->num_inner_iter=0;      options->num_inner_iter=0;
# Line 109  void Paso_Solver(paso::SystemMatrix_ptr Line 109  void Paso_Solver(paso::SystemMatrix_ptr
109              norm_max_of_b_local=0.;              norm_max_of_b_local=0.;
110              #pragma omp for private(i) schedule(static)              #pragma omp for private(i) schedule(static)
111              for (i = 0; i < numEqua ; ++i) {              for (i = 0; i < numEqua ; ++i) {
112                        norm2_of_b_local += r[i] * r[i];                  norm2_of_b_local += r[i] * r[i];
113                        norm_max_of_b_local = MAX(ABS(r[i]),norm_max_of_b_local);                  norm_max_of_b_local = MAX(ABS(r[i]),norm_max_of_b_local);
114              }              }
115              #pragma omp critical              #pragma omp critical
116              {              {
117                     norm2_of_b += norm2_of_b_local;                  norm2_of_b += norm2_of_b_local;
118                     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);
119              }              }
120          }          }
121  #ifdef ESYS_MPI  #ifdef ESYS_MPI
# Line 136  void Paso_Solver(paso::SystemMatrix_ptr Line 136  void Paso_Solver(paso::SystemMatrix_ptr
136                  printf("right hand side is identical to zero.\n");                  printf("right hand side is identical to zero.\n");
137          } else {          } else {
138              if (options->verbose) {              if (options->verbose) {
139                 printf("Paso_Solver: l2/lmax-norm of right hand side is  %e/%e.\n",norm2_of_b,norm_max_of_b);                  printf("Paso_Solver: l2/lmax-norm of right hand side is  %e/%e.\n",norm2_of_b,norm_max_of_b);
140                 printf("Paso_Solver: l2/lmax-stopping criterion is  %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance);                  printf("Paso_Solver: l2/lmax-stopping criterion is  %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance);
141                 switch (method) {                  switch (method) {
142                 case PASO_BICGSTAB:                      case PASO_BICGSTAB:
143                    printf("Paso_Solver: Iterative method is BiCGStab.\n");                          printf("Paso_Solver: Iterative method is BiCGStab.\n");
144                    break;                      break;
145                 case PASO_PCG:                      case PASO_PCG:
146                    printf("Paso_Solver: Iterative method is PCG.\n");                          printf("Paso_Solver: Iterative method is PCG.\n");
147                    break;                      break;
148                 case PASO_TFQMR:                      case PASO_TFQMR:
149                    printf("Paso_Solver: Iterative method is TFQMR.\n");                          printf("Paso_Solver: Iterative method is TFQMR.\n");
150                    break;                      break;
151                 case PASO_MINRES:                      case PASO_MINRES:
152                    printf("Paso_Solver: Iterative method is MINRES.\n");                          printf("Paso_Solver: Iterative method is MINRES.\n");
153                    break;                      break;
154                 case PASO_PRES20:                      case PASO_PRES20:
155                    printf("Paso_Solver: Iterative method is PRES20.\n");                          printf("Paso_Solver: Iterative method is PRES20.\n");
156                    break;                      break;
157                 case PASO_GMRES:                      case PASO_GMRES:
158                    if (options->restart>0) {                          if (options->restart>0) {
159                       printf("Paso_Solver: Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);                              printf("Paso_Solver: Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
160                    } else {                          } else {
161                       printf("Paso_Solver: Iterative method is GMRES(%d)\n",options->truncation);                              printf("Paso_Solver: Iterative method is GMRES(%d)\n",options->truncation);
162                    }                          }
163                    break;                      break;
164                 }                  }
165              }              }
166    
167              /* construct the preconditioner */              // construct the preconditioner
168              blocktimer_precond = blocktimer_time();              blocktimer_precond = blocktimer_time();
169              Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);              Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
170              A->setPreconditioner(options);              A->setPreconditioner(options);
171              Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);              Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
172              blocktimer_increment("Paso_Solver_setPreconditioner()", blocktimer_precond);              blocktimer_increment("Solver_setPreconditioner()", blocktimer_precond);
173              options->set_up_time=Esys_timer()-time_iter;              options->set_up_time=Esys_timer()-time_iter;
174              if (Esys_noError()) {              if (Esys_noError()) {
175                  /* get an initial guess by evaluating the preconditioner */                  // get an initial guess by evaluating the preconditioner
176                  A->solvePreconditioner(x, r);                  A->solvePreconditioner(x, r);
177    
178                  totIter = 1;                  totIter = 1;
# Line 181  void Paso_Solver(paso::SystemMatrix_ptr Line 181  void Paso_Solver(paso::SystemMatrix_ptr
181                  last_norm_max_of_residual=norm_max_of_b;                  last_norm_max_of_residual=norm_max_of_b;
182                  net_time_start=Esys_timer();                  net_time_start=Esys_timer();
183    
184                  /* Loop */                  // main loop
185                  while (! finalizeIteration) {                  while (!finalizeIteration) {
186                      cntIter = options->iter_max - totIter;                      cntIter = options->iter_max - totIter;
187                      finalizeIteration = TRUE;                      finalizeIteration = true;
188    
189                      /* Set initial residual. */                      // Set initial residual
190                      if (totIter>1) {                      if (totIter > 1) {
191                          // in the first iteration r = balance * b already                          // in the first iteration r = balance * b already
192                          A->applyBalance(r, b, true);                          A->applyBalance(r, b, true);
193                      }                      }
194    
195                      paso::SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);                      SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);
                       
196                      norm2_of_residual = 0;                      norm2_of_residual = 0;
197                      norm_max_of_residual = 0;                      norm_max_of_residual = 0;
198                      #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)                      #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
199                      {                      {
200                         norm2_of_residual_local = 0;                          norm2_of_residual_local = 0;
201                         norm_max_of_residual_local = 0;                          norm_max_of_residual_local = 0;
202                         #pragma omp for private(i) schedule(static)                          #pragma omp for private(i) schedule(static)
203                         for (i = 0; i < numEqua; i++) {                          for (i = 0; i < numEqua; i++) {
204                                norm2_of_residual_local+= r[i] * r[i];                              norm2_of_residual_local+= r[i] * r[i];
205                                norm_max_of_residual_local=MAX(ABS(r[i]),norm_max_of_residual_local);                              norm_max_of_residual_local=MAX(ABS(r[i]),norm_max_of_residual_local);
206                         }                          }
207                         #pragma omp critical                          #pragma omp critical
208                         {                          {
209                            norm2_of_residual += norm2_of_residual_local;                              norm2_of_residual += norm2_of_residual_local;
210                            norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);                              norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);
211                         }                          }
212                      }                      }
213  #ifdef ESYS_MPI  #ifdef ESYS_MPI
214                      /* TODO: use one call */                      // TODO: use one call
215                      loc_norm = norm2_of_residual;                      loc_norm = norm2_of_residual;
216                      MPI_Allreduce(&loc_norm,&norm2_of_residual, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);                      MPI_Allreduce(&loc_norm,&norm2_of_residual, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
217                      loc_norm = norm_max_of_residual;                      loc_norm = norm_max_of_residual;
# Line 221  void Paso_Solver(paso::SystemMatrix_ptr Line 220  void Paso_Solver(paso::SystemMatrix_ptr
220                      norm2_of_residual =sqrt(norm2_of_residual);                      norm2_of_residual =sqrt(norm2_of_residual);
221                      options->residual_norm=norm2_of_residual;                      options->residual_norm=norm2_of_residual;
222    
223                       if (options->verbose) printf("Paso_Solver: Step %5d: l2/lmax-norm of residual is  %e/%e",totIter,norm2_of_residual,norm_max_of_residual);                      if (options->verbose)
224                            printf("Paso_Solver: Step %5d: l2/lmax-norm of residual is  %e/%e",totIter,norm2_of_residual,norm_max_of_residual);
225    
226                       if (totIter>1 && norm2_of_residual>=last_norm2_of_residual &&  norm_max_of_residual>=last_norm_max_of_residual) {                      if (totIter>1 && norm2_of_residual>=last_norm2_of_residual &&  norm_max_of_residual>=last_norm_max_of_residual) {
227    
228                          if (options->verbose) printf(" divergence!\n");                          if (options->verbose) printf(" divergence!\n");
229                          Esys_setError(DIVERGED, "Paso_Solver: No improvement during iteration. Iterative solver gives up.");                          Esys_setError(DIVERGED, "Paso_Solver: No improvement during iteration. Iterative solver gives up.");
230    
231                       } else {                      } else {
232                          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 ||
233                                    norm_max_of_residual>tolerance*norm_max_of_b ) {
234                          tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b);  
235                          if (options->verbose) printf(" (new tolerance = %e).\n",tol);                              tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b);
236                                if (options->verbose) printf(" (new tolerance = %e).\n",tol);
237                          last_norm2_of_residual=norm2_of_residual;  
238                          last_norm_max_of_residual=norm_max_of_residual;                              last_norm2_of_residual=norm2_of_residual;
239                                last_norm_max_of_residual=norm_max_of_residual;
240                          /* call the solver */  
241                          switch (method) {                              // call the solver
242                          case PASO_BICGSTAB:                              switch (method) {
243                             errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol, pp);                                  case PASO_BICGSTAB:
244                             break;                                      errorCode = Solver_BiCGStab(A, r, x, &cntIter, &tol, pp);
245                          case PASO_PCG:                                  break;
246                             errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol, pp);                                  case PASO_PCG:
247                             break;                                      errorCode = Solver_PCG(A, r, x, &cntIter, &tol, pp);
248                          case PASO_TFQMR:                                  break;
249                             tol=tolerance*norm2_of_residual/norm2_of_b;                                  case PASO_TFQMR:
250                             errorCode = Paso_Solver_TFQMR(A, r, x0, &cntIter, &tol, pp);                                      tol=tolerance*norm2_of_residual/norm2_of_b;
251                             #pragma omp for private(i) schedule(static)                                      errorCode = Solver_TFQMR(A, r, x0, &cntIter, &tol, pp);
252                             for (i = 0; i < numEqua; i++) {                                      #pragma omp for private(i) schedule(static)
253                              x[i]+= x0[i];                                      for (i = 0; i < numEqua; i++) {
254                             }                                          x[i]+= x0[i];
255                             break;                                      }
256                          case PASO_MINRES:                                  break;
257                             /* tol=tolerance*norm2_of_residual/norm2_of_b; */                                  case PASO_MINRES:
258                             errorCode = Paso_Solver_MINRES(A, r, x, &cntIter, &tol, pp);                                      //tol=tolerance*norm2_of_residual/norm2_of_b;
259                             break;                                      errorCode = Solver_MINRES(A, r, x, &cntIter, &tol, pp);
260                          case PASO_PRES20:                                  break;
261                             errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20, pp);                                  case PASO_PRES20:
262                             break;                                      errorCode = Solver_GMRES(A, r, x, &cntIter, &tol, 5, 20, pp);
263                          case PASO_GMRES:                                  break;
264                             errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart, pp);                                  case PASO_GMRES:
265                             break;                                      errorCode = Solver_GMRES(A, r, x, &cntIter, &tol, options->truncation, options->restart, pp);
266                          }                                  break;
267                                }
268                          totIter += cntIter;  
269                                totIter += cntIter;
270                          /* error handling  */  
271                          if (errorCode==SOLVER_NO_ERROR) {                              // error handling
272                             finalizeIteration = FALSE;                              if (errorCode == SOLVER_NO_ERROR) {
273                          } else if (errorCode==SOLVER_MAXITER_REACHED) {                                  finalizeIteration = false;
274                             Esys_setError(DIVERGED,"Paso_Solver: maximum number of iteration steps reached.\nReturned solution does not fulfil stopping criterion.");                              } else if (errorCode == SOLVER_MAXITER_REACHED) {
275                             if (options->verbose) printf("Paso_Solver: Maximum number of iterations reached.\n");                                  Esys_setError(DIVERGED, "Paso_Solver: maximum number of iteration steps reached.\nReturned solution does not fulfil stopping criterion.");
276                          } else if (errorCode == SOLVER_INPUT_ERROR ) {                                  if (options->verbose)
277                             Esys_setError(SYSTEM_ERROR,"Paso_Solver: illegal dimension in iterative solver.");                                      printf("Paso_Solver: Maximum number of iterations reached.\n");
278                             if (options->verbose) printf("Paso_Solver: Internal error!\n");                              } else if (errorCode == SOLVER_INPUT_ERROR) {
279                          } else if (errorCode == SOLVER_NEGATIVE_NORM_ERROR) {                                  Esys_setError(SYSTEM_ERROR, "Paso_Solver: illegal dimension in iterative solver.");
280                             Esys_setError(VALUE_ERROR,"Paso_Solver: negative energy norm (try other solver or preconditioner).");                                  if (options->verbose)
281                             if (options->verbose) printf("Paso_Solver: negative energy norm (try other solver or preconditioner)!\n");                                      printf("Paso_Solver: Internal error!\n");
282                          } else if ( errorCode == SOLVER_BREAKDOWN ) {                              } else if (errorCode == SOLVER_NEGATIVE_NORM_ERROR) {
283                             if (cntIter <= 1) {                                  Esys_setError(VALUE_ERROR, "Paso_Solver: negative energy norm (try other solver or preconditioner).");
284                                Esys_setError(ZERO_DIVISION_ERROR, "Paso_Solver: fatal break down in iterative solver.");                                  if (options->verbose)
285                                if (options->verbose) printf("Paso_Solver: Uncurable break down!\n");                                      printf("Paso_Solver: negative energy norm (try other solver or preconditioner)!\n");
286                             } else {                              } else if (errorCode == SOLVER_BREAKDOWN) {
287                                if (options->verbose) printf("Paso_Solver: Breakdown at iter %d (residual = %e). Restarting ...\n", totIter, tol);                                  if (cntIter <= 1) {
288                                finalizeIteration = FALSE;                                      Esys_setError(ZERO_DIVISION_ERROR, "Paso_Solver: fatal break down in iterative solver.");
289                             }                                      if (options->verbose)
290                                            printf("Paso_Solver: Uncurable break down!\n");
291                                    } else {
292                                        if (options->verbose)
293                                            printf("Paso_Solver: Breakdown at iter %d (residual = %e). Restarting ...\n", totIter, tol);
294                                        finalizeIteration = false;
295                                    }
296                                } else {
297                                    Esys_setError(SYSTEM_ERROR, "Paso_Solver: Generic error in solver.");
298                                    if (options->verbose)
299                                        printf("Paso_Solver: Generic error in solver!\n");
300                                }
301                          } else {                          } else {
302                             Esys_setError(SYSTEM_ERROR,"Paso_Solver: Generic error in solver.");                              if (options->verbose)
303                             if (options->verbose) printf("Paso_Solver: Generic error in solver!\n");                                  printf(" convergence!\n");
304                                options->converged = true;
305                          }                          }
306                        } else {                      }
307                           if (options->verbose) printf(" convergence!\n");                  } // while
308                           options->converged=TRUE;                  options->net_time = Esys_timer()-net_time_start;
                       }  
                    }  
                 } /* while */  
                 options->net_time=Esys_timer()-net_time_start;  
309              }              }
310              options->num_iter=totIter;              options->num_iter = totIter;
311              A->applyBalanceInPlace(x, false);              A->applyBalanceInPlace(x, false);
312          }          }
313      }      }
314      delete[] r;      delete[] r;
315      delete[] x0;      delete[] x0;
316      options->time=Esys_timer()-time_iter;      options->time = Esys_timer()-time_iter;
317      Performance_stopMonitor(pp,PERFORMANCE_ALL);      Performance_stopMonitor(pp, PERFORMANCE_ALL);
318      blocktimer_increment("Paso_Solver()", blocktimer_start);      blocktimer_increment("Paso_Solver()", blocktimer_start);
319  }  }
320    
321    } // namespace paso
322    

Legend:
Removed from v.4855  
changed lines
  Added in v.4856

  ViewVC Help
Powered by ViewVC 1.1.26