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

trunk/esys2/paso/src/Solvers/Solver.c revision 153 by jgs, Tue Oct 25 01:51:20 2005 UTC trunk/paso/src/Solvers/Solver.c revision 633 by dhawcroft, Thu Mar 23 05:37:00 2006 UTC
# Line 1  Line 1 
1  /* $Id$ */  /* $Id$ */
2    
3    
4    /*
5    ********************************************************************************
6    *               Copyright   2006 by ACcESS MNRF                                *
7    *                                                                              *
8    *                 http://www.access.edu.au                                     *
9    *           Primary Business: Queensland, Australia                            *
10    *     Licensed under the Open Software License version 3.0             *
11    *        http://www.opensource.org/licenses/osl-3.0.php                        *
12    ********************************************************************************
13    */
14    
15  /**************************************************************/  /**************************************************************/
16    
17  /* Paso: SystemMatrix: controls iterative linear system solvers  */  /* Paso: SystemMatrix: controls iterative linear system solvers  */
# Line 15  Line 27 
27  #include "SystemMatrix.h"  #include "SystemMatrix.h"
28  #include "Solver.h"  #include "Solver.h"
29    
 #if PTR_OFFSET !=0 || INDEX_OFFSET!=0  
 #error Paso library usage requires PTR_OFFSET=0 and INDEX_OFFSET=0  
 #endif  
   
   
30  /***********************************************************************************/  /***********************************************************************************/
31    
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->iterative);      Paso_Preconditioner_free(A->solver);
36      A->iterative=NULL;      A->solver=NULL;
37  }  }
38  /*  call the iterative solver: */  /*  call the iterative solver: */
39    
40  void Paso_Solver(Paso_SystemMatrix* A,double* x,double* b,Paso_Options* options) {  void Paso_Solver(Paso_SystemMatrix* A,double* x,double* b,Paso_Options* options,Paso_Performance* pp) {
41      double norm2_of_b,tol,tolerance,time_iter,time_prec,*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,
42             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,
43             last_norm_max_of_residual,*scaling;             last_norm_max_of_residual,*scaling;
44      dim_t i,totIter,cntIter,method;       dim_t i,totIter,cntIter,method;
45      bool_t finalizeIteration;       bool_t finalizeIteration;
46      err_t errorCode;       err_t errorCode;
47      dim_t n_col = A->num_cols * A-> col_block_size;       dim_t n_col = A->num_cols * A-> col_block_size;
48      dim_t n_row = A->num_rows * A-> row_block_size;       dim_t n_row = A->num_rows * A-> row_block_size;
49    
50      tolerance=MAX(options->tolerance,EPSILON);       tolerance=MAX(options->tolerance,EPSILON);
51      Paso_resetError();       Paso_resetError();
52      method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);       method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);
53      /* check matrix type */       /* check matrix type */
54      if (A->type!=CSR) {       if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) || (A->type & MATRIX_FORMAT_SYM) ) {
55        Paso_setError(TYPE_ERROR,"Iterative solver requires CSR format.");         Paso_setError(TYPE_ERROR,"Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0.");
56      }       }
57      if (A->col_block_size != A->row_block_size) {       if (A->col_block_size != A->row_block_size) {
58         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.");
59       }       }
60       if (A->num_cols != A->num_rows) {       if (A->num_cols != A->num_rows) {
61         Paso_setError(TYPE_ERROR,"Iterative solver requires requires a square matrix.");          Paso_setError(TYPE_ERROR,"Iterative solver requires requires a square matrix.");
62         return;          return;
63       }       }
64         Performance_startMonitor(pp,PERFORMANCE_ALL);
65       if (Paso_noError()) {       if (Paso_noError()) {
66          /* get normalization */          /* get normalization */
67          scaling=Paso_SystemMatrix_borrowNormalization(A);          scaling=Paso_SystemMatrix_borrowNormalization(A);
68          if (! Paso_noError()) return;          if (Paso_noError()) {
69          /* get the norm of the right hand side */             /* get the norm of the right hand side */
70          norm2_of_b=0.;             norm2_of_b=0.;
71          norm_max_of_b=0.;             norm_max_of_b=0.;
72          #pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local)             #pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local)
73          {             {
74                  norm2_of_b_local=0.;                  norm2_of_b_local=0.;
75                  norm_max_of_b_local=0.;                  norm_max_of_b_local=0.;
76                  #pragma omp for private(i) schedule(static)                  #pragma omp for private(i) schedule(static)
# Line 75  void Paso_Solver(Paso_SystemMatrix* A,do Line 83  void Paso_Solver(Paso_SystemMatrix* A,do
83                     norm2_of_b += norm2_of_b_local;                     norm2_of_b += norm2_of_b_local;
84                     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);
85                  }                  }
86          }             }
87          norm2_of_b=sqrt(norm2_of_b);             norm2_of_b=sqrt(norm2_of_b);
88            
89          /* if norm2_of_b==0 we are ready: x=0 */             /* if norm2_of_b==0 we are ready: x=0 */
90          if (norm2_of_b <=0.) {             if (norm2_of_b <=0.) {
91              #pragma omp parallel for private(i) schedule(static)                 #pragma omp parallel for private(i) schedule(static)
92              for (i = 0; i < n_col; i++) x[i]=0.;                 for (i = 0; i < n_col; i++) x[i]=0.;
93              if (options->verbose) printf("right hand side is identical zero.\n");                 if (options->verbose) printf("right hand side is identical zero.\n");
94              return;             } else {
95          }                if (options->verbose) {
96          if (options->verbose) {                    printf("l2/lmax-norm of right hand side is  %e/%e.\n",norm2_of_b,norm_max_of_b);
97              printf("l2/lmax-norm of right hand side is  %e/%e.\n",norm2_of_b,norm_max_of_b);                    printf("l2/lmax-stopping criterion is  %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance);
98              printf("l2/lmax-stopping criterion is  %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance);                    switch (method) {
99              switch (method) {                       case PASO_BICGSTAB:
100                 case PASO_BICGSTAB:                           printf("Iterative method is BiCGStab.\n");
101                     printf("Iterative method is BiCGStab.\n");                           break;
102                     break;                       case PASO_PCG:
103                 case PASO_PCG:                           printf("Iterative method is PCG.\n");
104                     printf("Iterative method is PCG.\n");                           break;
105                     break;                       case PASO_PRES20:
106                 case PASO_PRES20:                           printf("Iterative method is PRES20.\n");
107                     printf("Iterative method is PRES20.\n");                           break;
108                     break;                       case PASO_GMRES:
109                 case PASO_GMRES:                           if (options->restart>0) {
110                     if (options->restart>0) {                              printf("Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
111                        printf("Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);                           } else {
112                                printf("Iterative method is GMRES(%d)\n",options->truncation);
113                             }
114                             break;
115                      }
116                  }
117                  /* construct the preconditioner */
118                  
119                  Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
120                  Paso_Solver_setPreconditioner(A,options);
121                  Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
122                  if (! Paso_noError()) return;
123              
124                  /* get an initial guess by evaluating the preconditioner */
125                  time_iter=Paso_timer();
126                  #pragma omp parallel
127                  Paso_Solver_solvePreconditioner(A,x,b);
128                  if (! Paso_noError()) return;
129                  /* start the iteration process :*/
130                  r=TMPMEMALLOC(n_row,double);
131                  if (Paso_checkPtr(r)) return;
132                  totIter = 0;
133                  finalizeIteration = FALSE;
134                  last_norm2_of_residual=norm2_of_b;
135                  last_norm_max_of_residual=norm_max_of_b;
136                  while (! finalizeIteration) {
137                     cntIter = options->iter_max - totIter;
138                     finalizeIteration = TRUE;
139                     /*     Set initial residual. */
140                     norm2_of_residual = 0;
141                     norm_max_of_residual = 0;
142                     #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
143                     {
144                        #pragma omp for private(i) schedule(static)
145                        for (i = 0; i < n_row; i++) r[i]=b[i];
146              
147                        Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);
148              
149                        norm2_of_residual_local = 0;
150                        norm_max_of_residual_local = 0;
151                        #pragma omp for private(i) schedule(static)
152                        for (i = 0; i < n_row; i++) {
153                               norm2_of_residual_local+= r[i] * r[i];
154                               norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);
155                        }
156                        #pragma omp critical
157                        {
158                           norm2_of_residual += norm2_of_residual_local;
159                           norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);
160                        }
161                     }
162                     norm2_of_residual =sqrt(norm2_of_residual);
163              
164                     if (options->verbose) printf("Step %5d: l2/lmax-norm of residual is  %e/%e",totIter,norm2_of_residual,norm_max_of_residual);
165                     if (totIter>0 && norm2_of_residual>=last_norm2_of_residual &&  norm_max_of_residual>=last_norm_max_of_residual) {
166                          if (options->verbose) printf(" divergence!\n");
167                          Paso_setError(WARNING, "No improvement during iteration. Iterative solver gives up.");
168                     } else {
169                        /* */
170                        if (norm2_of_residual>tolerance*norm2_of_b || norm_max_of_residual>tolerance*norm_max_of_b ) {
171                            tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b);
172                           if (options->verbose) printf(" (new tolerance = %e).\n",tol);
173                           last_norm2_of_residual=norm2_of_residual;
174                           last_norm_max_of_residual=norm_max_of_residual;
175                           /* call the solver */
176                           switch (method) {
177                                  case PASO_BICGSTAB:
178                                      errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol,pp);
179                                      break;
180                                  case PASO_PCG:
181                                      errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol,pp);
182                                      break;
183                                  case PASO_PRES20:
184                                      errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20,pp);
185                                      break;
186                                  case PASO_GMRES:
187                                      errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart,pp);
188                                      break;
189                            }
190                            totIter += cntIter;
191                      
192                            /* error handling  */
193                            if (errorCode==NO_ERROR) {
194                                  finalizeIteration = FALSE;
195                             } else if (errorCode==SOLVER_MAXITER_REACHED) {
196                                   Paso_setError(WARNING,"maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");
197                                   if (options->verbose) printf("Maximum number of iterations reached.!\n");
198                            } else if (errorCode == SOLVER_INPUT_ERROR ) {
199                                   Paso_setError(SYSTEM_ERROR,"illegal dimension in iterative solver.");
200                                   if (options->verbose) printf("Internal error!\n");
201                            } else if ( errorCode == SOLVER_BREAKDOWN ) {
202                             if (cntIter <= 1) {
203                                     Paso_setError(ZERO_DIVISION_ERROR, "fatal break down in iterative solver.");
204                                     if (options->verbose) printf("Uncurable break down!\n");
205                                } else {
206                                 if (options->verbose) printf("Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);
207                                    finalizeIteration = FALSE;
208                                }
209                            } else if (errorCode == SOLVER_MEMORY_ERROR) {
210                              Paso_setError(MEMORY_ERROR,"memory allocation failed.");
211                                 if (options->verbose) printf("Memory allocation failed!\n");
212                            } else if (errorCode !=SOLVER_NO_ERROR ) {
213                              Paso_setError(SYSTEM_ERROR,"unidentified error in iterative solver.");
214                                 if (options->verbose) printf("Unidentified error!\n");
215                            }
216                     } else {                     } else {
217                        printf("Iterative method is GMRES(%d)\n",options->truncation);                        if (options->verbose) printf(". convergence! \n");
218                     }                     }
219                     break;                  }
220              }                } /* while */
221          }                MEMFREE(r);
222          /* construct the preconditioner */                time_iter=Paso_timer()-time_iter;
223                          if (options->verbose)  {
224          time_prec=Paso_timer();                   printf("timing: solver: %.4e sec\n",time_iter);
225          Paso_Solver_setPreconditioner(A,options);                   if (totIter>0) printf("timing: per iteration step: %.4e sec\n",time_iter/totIter);
         time_prec=Paso_timer()-time_prec;  
         if (! Paso_noError()) return;  
       
         /* get an initial guess by evaluating the preconditioner */  
         time_iter=Paso_timer();  
         #pragma omp parallel  
         Paso_Solver_solvePreconditioner(A,x,b);  
         if (! Paso_noError()) return;  
         /* start the iteration process :*/  
         r=TMPMEMALLOC(n_row,double);  
         if (Paso_checkPtr(r)) return;  
         totIter = 0;  
         finalizeIteration = FALSE;  
         last_norm2_of_residual=norm2_of_b;  
         last_norm_max_of_residual=norm_max_of_b;  
         while (! finalizeIteration) {  
            cntIter = options->iter_max - totIter;  
            finalizeIteration = TRUE;  
            /*     Set initial residual. */  
            norm2_of_residual = 0;  
            norm_max_of_residual = 0;  
            #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)  
            {  
               #pragma omp for private(i) schedule(static)  
               for (i = 0; i < n_row; i++) r[i]=b[i];  
       
               Paso_SystemMatrix_MatrixVector(DBLE(-1), A, x, DBLE(1), r);  
       
               norm2_of_residual_local = 0;  
               norm_max_of_residual_local = 0;  
               #pragma omp for private(i) schedule(static)  
               for (i = 0; i < n_row; i++) {  
                      norm2_of_residual_local+= r[i] * r[i];  
                      norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);  
               }  
               #pragma omp critical  
               {  
                  norm2_of_residual += norm2_of_residual_local;  
                  norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);  
226                }                }
227             }             }
            norm2_of_residual =sqrt(norm2_of_residual);  
       
            if (options->verbose) printf("Step %5d: l2/lmax-norm of residual is  %e/%e",totIter,norm2_of_residual,norm_max_of_residual);  
            if (totIter>0 && norm2_of_residual>=last_norm2_of_residual &&  norm_max_of_residual>=last_norm_max_of_residual) {  
                 if (options->verbose) printf(" divergence!\n");  
                 Paso_setError(WARNING, "No improvement during iteration. Iterative solver gives up.");  
            } else {  
               /* */  
               if (norm2_of_residual>tolerance*norm2_of_b || norm_max_of_residual>tolerance*norm_max_of_b ) {  
                   tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b);  
                  if (options->verbose) printf(" (new tolerance = %e).\n",tol);  
                  last_norm2_of_residual=norm2_of_residual;  
                  last_norm_max_of_residual=norm_max_of_residual;  
                  /* call the solver */  
                  switch (method) {  
                         case PASO_BICGSTAB:  
                             errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol);  
                             break;  
                         case PASO_PCG:  
                             errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol);  
                             break;  
                         case PASO_PRES20:  
                             errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20);  
                             break;  
                         case PASO_GMRES:  
                             errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart);  
                             break;  
                   }  
                   totIter += cntIter;  
               
                   /* error handling  */  
                   if (errorCode==NO_ERROR) {  
                         finalizeIteration = FALSE;  
                   } else if (errorCode==SOLVER_MAXITER_REACHED) {  
                          Paso_setError(WARNING,"maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");  
                          if (options->verbose) printf("Maximum number of iterations reached.!\n");  
                   } else if (errorCode == SOLVER_INPUT_ERROR ) {  
                          Paso_setError(SYSTEM_ERROR,"illegal dimension in iterative solver.");  
                          if (options->verbose) printf("Internal error!\n");  
                   } else if ( errorCode == SOLVER_BREAKDOWN ) {  
                   if (cntIter <= 1) {  
                            Paso_setError(ZERO_DIVISION_ERROR, "fatal break down in iterative solver.");  
                            if (options->verbose) printf("Uncurable break down!\n");  
                       } else {  
                       if (options->verbose) printf("Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);  
                           finalizeIteration = FALSE;  
                       }  
                   } else if (errorCode == SOLVER_MEMORY_ERROR) {  
                    Paso_setError(MEMORY_ERROR,"memory allocation failed.");  
                        if (options->verbose) printf("Memory allocation failed!\n");  
                   } else if (errorCode !=SOLVER_NO_ERROR ) {  
                    Paso_setError(SYSTEM_ERROR,"unidentified error in iterative solver.");  
                        if (options->verbose) printf("Unidentified error!\n");  
                   }  
              } else {  
                 if (options->verbose) printf(". convergence! \n");  
              }  
           }  
         } /* while */  
         MEMFREE(r);  
         time_iter=Paso_timer()-time_iter;  
         if (options->verbose)  {  
            printf("timing: solver: %.4e sec\n",time_prec+time_iter);  
            printf("timing: preconditioner: %.4e sec\n",time_prec);  
            if (totIter>0) printf("timing: per iteration step: %.4e sec\n",time_iter/totIter);  
228          }          }
229     }        }
230          Performance_stopMonitor(pp,PERFORMANCE_ALL);
231  }  }

Legend:
Removed from v.153  
changed lines
  Added in v.633

  ViewVC Help
Powered by ViewVC 1.1.26