/[escript]/branches/domexper/paso/src/Solver.c
ViewVC logotype

Diff of /branches/domexper/paso/src/Solver.c

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

revision 3233 by jfenwick, Fri Oct 1 01:53:46 2010 UTC revision 3234 by jfenwick, Mon Oct 4 01:46:30 2010 UTC
# Line 32  Line 32 
32  /*  free space */  /*  free space */
33    
34  void Paso_Solver_free(Paso_SystemMatrix* A) {  void Paso_Solver_free(Paso_SystemMatrix* A) {
35     Paso_Preconditioner_free(A->solver);     Paso_SystemMatrix_freePreconditioner(A);
    A->solver=NULL;  
36  }  }
37  /*  call the iterative solver: */  /*  call the iterative solver: */
38    
# Line 86  void Paso_Solver(Paso_SystemMatrix* A,do Line 85  void Paso_Solver(Paso_SystemMatrix* A,do
85       if (method==PASO_NONLINEAR_GMRES) {       if (method==PASO_NONLINEAR_GMRES) {
86          Paso_Function* F=NULL;          Paso_Function* F=NULL;
87          F=Paso_Function_LinearSystem_alloc(A,b,options);          F=Paso_Function_LinearSystem_alloc(A,b,options);
88          Paso_Solver_solvePreconditioner(A,x,b);          Paso_SystemMatrix_solvePreconditioner(A,x,b);
89          errorCode=Paso_Solver_NewtonGMRES(F,x,options,pp);          errorCode=Paso_Solver_NewtonGMRES(F,x,options,pp);
90          if (errorCode!=NO_ERROR) {          if (errorCode!=NO_ERROR) {
91             Esys_setError(SYSTEM_ERROR,"Paso_Solver_NewtonGMRES: an error has occured.");             Esys_setError(SYSTEM_ERROR,"Paso_Solver_NewtonGMRES: an error has occured.");
# Line 133  void Paso_Solver(Paso_SystemMatrix* A,do Line 132  void Paso_Solver(Paso_SystemMatrix* A,do
132           if ( IS_NAN(norm2_of_b) || IS_NAN(norm_max_of_b) ) {           if ( IS_NAN(norm2_of_b) || IS_NAN(norm_max_of_b) ) {
133              Esys_setError(VALUE_ERROR,"Paso_Solver: Matrix or right hand side contains undefined values.");              Esys_setError(VALUE_ERROR,"Paso_Solver: Matrix or right hand side contains undefined values.");
134           } else if (norm2_of_b <=0.) {           } else if (norm2_of_b <=0.) {
135    
136              #pragma omp parallel for private(i) schedule(static)              #pragma omp parallel for private(i) schedule(static)
137              for (i = 0; i < numSol; i++) x[i]=0.;              for (i = 0; i < numSol; i++) x[i]=0.;
138              if (options->verbose) printf("right hand side is identical zero.\n");              if (options->verbose) printf("right hand side is identical zero.\n");
139    
140           } else {           } else {
141    
142              if (options->verbose) {              if (options->verbose) {
143                 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);
144                 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);
# Line 164  void Paso_Solver(Paso_SystemMatrix* A,do Line 166  void Paso_Solver(Paso_SystemMatrix* A,do
166                    }                    }
167                    break;                    break;
168                 }                 }
169    
170              }              }
171    
172              /* construct the preconditioner */              /* construct the preconditioner */
173              blocktimer_precond = blocktimer_time();              blocktimer_precond = blocktimer_time();
174              Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);              Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
175              Paso_Solver_setPreconditioner(A,options);          Paso_SystemMatrix_setPreconditioner(A,options);
176              Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);              Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
177              blocktimer_increment("Paso_Solver_setPreconditioner()", blocktimer_precond);              blocktimer_increment("Paso_Solver_setPreconditioner()", blocktimer_precond);
178              options->set_up_time=Esys_timer()-time_iter;              options->set_up_time=Esys_timer()-time_iter;
179              if (! Esys_noError()) return;              if (! Esys_noError()) return;
180    
181    
182                /* get an initial guess by evaluating the preconditioner */                /* get an initial guess by evaluating the preconditioner */
183                Paso_Solver_solvePreconditioner(A,x,b);            Paso_SystemMatrix_solvePreconditioner(A,x,b);
184    
185                /* start the iteration process :*/                /* start the iteration process :*/
186                r=TMPMEMALLOC(numEqua,double);                r=TMPMEMALLOC(numEqua,double);
187                x0=TMPMEMALLOC(numEqua,double);                x0=TMPMEMALLOC(numEqua,double);
188                Esys_checkPtr(r);                Esys_checkPtr(r);
189            Esys_checkPtr(x0);            Esys_checkPtr(x0);
190                if (Esys_noError()) {                if (Esys_noError()) {
191    
192                   totIter = 1;                   totIter = 1;
193                   finalizeIteration = FALSE;                   finalizeIteration = FALSE;
194                   last_norm2_of_residual=norm2_of_b;                   last_norm2_of_residual=norm2_of_b;
195                   last_norm_max_of_residual=norm_max_of_b;                   last_norm_max_of_residual=norm_max_of_b;
196           net_time_start=Esys_timer();           net_time_start=Esys_timer();
197    
198                   /* Loop */                   /* Loop */
199                   while (! finalizeIteration) {                   while (! finalizeIteration) {
200                      cntIter = options->iter_max - totIter;                      cntIter = options->iter_max - totIter;
201                      finalizeIteration = TRUE;                      finalizeIteration = TRUE;
202    
203                      /*     Set initial residual. */                      /*     Set initial residual. */
204                      norm2_of_residual = 0;                      norm2_of_residual = 0;
205                      norm_max_of_residual = 0;                      norm_max_of_residual = 0;
# Line 222  void Paso_Solver(Paso_SystemMatrix* A,do Line 232  void Paso_Solver(Paso_SystemMatrix* A,do
232                      options->residual_norm=norm2_of_residual;                      options->residual_norm=norm2_of_residual;
233    
234                    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) printf("Paso_Solver: Step %5d: l2/lmax-norm of residual is  %e/%e",totIter,norm2_of_residual,norm_max_of_residual);
235    
236                    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) {
237    
238                       if (options->verbose) printf(" divergence!\n");                       if (options->verbose) printf(" divergence!\n");
239                       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.");
240    
241                    } else {                    } else {
242    
243                       /* */                       /* */
244                       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 ) {
245    
246                          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);
247                          if (options->verbose) printf(" (new tolerance = %e).\n",tol);                          if (options->verbose) printf(" (new tolerance = %e).\n",tol);
248    
249                          last_norm2_of_residual=norm2_of_residual;                          last_norm2_of_residual=norm2_of_residual;
250                          last_norm_max_of_residual=norm_max_of_residual;                          last_norm_max_of_residual=norm_max_of_residual;
251    
252                          /* call the solver */                          /* call the solver */
253                          switch (method) {                          switch (method) {
254                          case PASO_BICGSTAB:                          case PASO_BICGSTAB:
# Line 249  void Paso_Solver(Paso_SystemMatrix* A,do Line 266  void Paso_Solver(Paso_SystemMatrix* A,do
266                             }                             }
267                             break;                             break;
268                          case PASO_MINRES:                          case PASO_MINRES:
269                             tol=tolerance*norm2_of_residual/norm2_of_b;                             /* tol=tolerance*norm2_of_residual/norm2_of_b; */
270                             errorCode = Paso_Solver_MINRES(A, r, x0, &cntIter, &tol, pp);                             errorCode = Paso_Solver_MINRES(A, r, x, &cntIter, &tol, pp);
                            #pragma omp for private(i) schedule(static)  
                            for (i = 0; i < numEqua; i++) {  
                             x[i]+= x0[i];  
                            }  
271                             break;                             break;
272                          case PASO_PRES20:                          case PASO_PRES20:
273                             errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20, pp);                             errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20, pp);
# Line 263  void Paso_Solver(Paso_SystemMatrix* A,do Line 276  void Paso_Solver(Paso_SystemMatrix* A,do
276                             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);
277                             break;                             break;
278                          }                          }
279    
280                          totIter += cntIter;                          totIter += cntIter;
281    
282                          /* error handling  */                          /* error handling  */
283                          if (errorCode==SOLVER_NO_ERROR) {                          if (errorCode==SOLVER_NO_ERROR) {
284                             finalizeIteration = FALSE;                             finalizeIteration = FALSE;
# Line 273  void Paso_Solver(Paso_SystemMatrix* A,do Line 288  void Paso_Solver(Paso_SystemMatrix* A,do
288                          } else if (errorCode == SOLVER_INPUT_ERROR ) {                          } else if (errorCode == SOLVER_INPUT_ERROR ) {
289                             Esys_setError(SYSTEM_ERROR,"Paso_Solver: illegal dimension in iterative solver.");                             Esys_setError(SYSTEM_ERROR,"Paso_Solver: illegal dimension in iterative solver.");
290                             if (options->verbose) printf("Paso_Solver: Internal error!\n");                             if (options->verbose) printf("Paso_Solver: Internal error!\n");
291                } else if (errorCode == SOLVER_NEGATIVE_NORM_ERROR) {
292                   Esys_setError(VALUE_ERROR,"Paso_Solver: negative energy norm (try other solver or preconditioner).");
293                   if (options->verbose) printf("Paso_Solver: negative energy norm (try other solver or preconditioner)!\n");
294                          } else if ( errorCode == SOLVER_BREAKDOWN ) {                          } else if ( errorCode == SOLVER_BREAKDOWN ) {
295                             if (cntIter <= 1) {                             if (cntIter <= 1) {
296                                Esys_setError(ZERO_DIVISION_ERROR, "Paso_Solver: fatal break down in iterative solver.");                                Esys_setError(ZERO_DIVISION_ERROR, "Paso_Solver: fatal break down in iterative solver.");
297                                if (options->verbose) printf("Paso_Solver: Uncurable break down!\n");                                if (options->verbose) printf("Paso_Solver: Uncurable break down!\n");
298                             } else {                             } else {
299                                if (options->verbose) printf("Paso_Solver: Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);                                if (options->verbose) printf("Paso_Solver: Breakdown at iter %d (residual = %e). Restarting ...\n", totIter, tol);
300                                finalizeIteration = FALSE;                                finalizeIteration = FALSE;
301                             }                             }
302                            } else {
303                   Esys_setError(SYSTEM_ERROR,"Paso_Solver:generic error in solver.");
304                   if (options->verbose) printf("Paso_Solver: generic error in solver!\n");
305                          }                          }
306                        } else {                        } else {
307                           if (options->verbose) printf(". convergence! \n");                           if (options->verbose) printf(". convergence! \n");

Legend:
Removed from v.3233  
changed lines
  Added in v.3234

  ViewVC Help
Powered by ViewVC 1.1.26