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

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

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

revision 1504 by trankine, Mon Apr 14 06:43:51 2008 UTC revision 1505 by trankine, Mon Apr 14 12:57:37 2008 UTC
# Line 34  Line 34 
34  /*  free space */  /*  free space */
35    
36  void Paso_Solver_free(Paso_SystemMatrix* A) {  void Paso_Solver_free(Paso_SystemMatrix* A) {
37      Paso_Preconditioner_free(A->solver);     Paso_Preconditioner_free(A->solver);
38      A->solver=NULL;     A->solver=NULL;
39  }  }
40  /*  call the iterative solver: */  /*  call the iterative solver: */
41    
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,
43      double norm2_of_b,tol,tolerance,time_iter,*r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b,                   Paso_Options* options,Paso_Performance* pp) {
44             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, loc_norm;     double norm2_of_b,tol,tolerance,time_iter;
46       dim_t i,totIter,cntIter,method;     double *r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b;
47       bool_t finalizeIteration;     double norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local;
48       err_t errorCode;     double norm_max_of_residual_local,norm_max_of_residual;
49       dim_t numSol = Paso_SystemMatrix_getTotalNumCols(A);     double last_norm_max_of_residual,*scaling;
50       dim_t numEqua = Paso_SystemMatrix_getTotalNumRows(A);  #ifdef PASO_MPI
51       double blocktimer_start = blocktimer_time();     double loc_norm;
52    #endif
53       dim_t i,totIter,cntIter,method;
54       bool_t finalizeIteration;
55       err_t errorCode;
56       dim_t numSol = Paso_SystemMatrix_getTotalNumCols(A);
57       dim_t numEqua = Paso_SystemMatrix_getTotalNumRows(A);
58       double blocktimer_start = blocktimer_time();
59    
60       tolerance=MAX(options->tolerance,EPSILON);     tolerance=MAX(options->tolerance,EPSILON);
61       Paso_resetError();     Paso_resetError();
62       method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);     method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);
63       /* check matrix type */     /* check matrix type */
64       if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) || (A->type & MATRIX_FORMAT_SYM) ) {     if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) || (A->type & MATRIX_FORMAT_SYM) ) {
65         Paso_setError(TYPE_ERROR,"Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0.");        Paso_setError(TYPE_ERROR,"Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0.");
66       }     }
67       if (A->col_block_size != A->row_block_size) {     if (A->col_block_size != A->row_block_size) {
68          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.");
69       }     }
70       if (Paso_SystemMatrix_getGlobalNumCols(A) != Paso_SystemMatrix_getGlobalNumRows(A)) {     if (Paso_SystemMatrix_getGlobalNumCols(A) != Paso_SystemMatrix_getGlobalNumRows(A)) {
71          Paso_setError(TYPE_ERROR,"Iterative solver requires requires a square matrix.");        Paso_setError(TYPE_ERROR,"Iterative solver requires requires a square matrix.");
72          return;        return;
73       }     }
74       Performance_startMonitor(pp,PERFORMANCE_ALL);     Performance_startMonitor(pp,PERFORMANCE_ALL);
75       if (Paso_noError()) {     if (Paso_noError()) {
76          /* get normalization */        /* get normalization */
77          scaling=Paso_SystemMatrix_borrowNormalization(A);        scaling=Paso_SystemMatrix_borrowNormalization(A);
78          if (Paso_noError()) {        if (Paso_noError()) {
79             /* get the norm of the right hand side */           /* get the norm of the right hand side */
80             norm2_of_b=0.;           norm2_of_b=0.;
81             norm_max_of_b=0.;           norm_max_of_b=0.;
82             #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)
83             {           {
84                  norm2_of_b_local=0.;              norm2_of_b_local=0.;
85                  norm_max_of_b_local=0.;              norm_max_of_b_local=0.;
86                  #pragma omp for private(i) schedule(static)  #pragma omp for private(i) schedule(static)
87                  for (i = 0; i < numEqua ; ++i) {              for (i = 0; i < numEqua ; ++i) {
88                        norm2_of_b_local += b[i] * b[i];                 norm2_of_b_local += b[i] * b[i];
89                        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);
90                  }              }
91                  #pragma omp critical  #pragma omp critical
92                  {              {
93                     norm2_of_b += norm2_of_b_local;                 norm2_of_b += norm2_of_b_local;
94                     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);
95                  }              }
96             }           }
97             /* TODO: use one call */           /* TODO: use one call */
98             #ifdef PASO_MPI  #ifdef PASO_MPI
99             {           {
100                 loc_norm = norm2_of_b;              loc_norm = norm2_of_b;
101                 MPI_Allreduce(&loc_norm,&norm2_of_b, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);              MPI_Allreduce(&loc_norm,&norm2_of_b, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
102                 loc_norm = norm_max_of_b;              loc_norm = norm_max_of_b;
103                 MPI_Allreduce(&loc_norm,&norm_max_of_b, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);              MPI_Allreduce(&loc_norm,&norm_max_of_b, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
104             }           }
105             #endif  #endif
106             norm2_of_b=sqrt(norm2_of_b);           norm2_of_b=sqrt(norm2_of_b);
107            
108             /* if norm2_of_b==0 we are ready: x=0 */           /* if norm2_of_b==0 we are ready: x=0 */
109             if (norm2_of_b <=0.) {           if (norm2_of_b <=0.) {
110                 #pragma omp parallel for private(i) schedule(static)  #pragma omp parallel for private(i) schedule(static)
111                 for (i = 0; i < numSol; i++) x[i]=0.;              for (i = 0; i < numSol; i++) x[i]=0.;
112                 if (options->verbose) printf("right hand side is identical zero.\n");              if (options->verbose) printf("right hand side is identical zero.\n");
113             } else {           } else {
114                if (options->verbose) {              if (options->verbose) {
115                    printf("l2/lmax-norm of right hand side is  %e/%e.\n",norm2_of_b,norm_max_of_b);                 printf("l2/lmax-norm of right hand side is  %e/%e.\n",norm2_of_b,norm_max_of_b);
116                    printf("l2/lmax-stopping criterion is  %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance);                 printf("l2/lmax-stopping criterion is  %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance);
117                    switch (method) {                 switch (method) {
118                       case PASO_BICGSTAB:                 case PASO_BICGSTAB:
119                           printf("Iterative method is BiCGStab.\n");                    printf("Iterative method is BiCGStab.\n");
120                           break;                    break;
121                       case PASO_PCG:                 case PASO_PCG:
122                           printf("Iterative method is PCG.\n");                    printf("Iterative method is PCG.\n");
123                           break;                    break;
124                       case PASO_PRES20:                 case PASO_PRES20:
125                           printf("Iterative method is PRES20.\n");                    printf("Iterative method is PRES20.\n");
126                           break;                    break;
127                       case PASO_GMRES:                 case PASO_GMRES:
128                           if (options->restart>0) {                    if (options->restart>0) {
129                              printf("Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);                       printf("Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
130                           } else {                    } else {
131                              printf("Iterative method is GMRES(%d)\n",options->truncation);                       printf("Iterative method is GMRES(%d)\n",options->truncation);
                          }  
                          break;  
132                    }                    }
133                }                    break;
134                /* construct the preconditioner */                 }
135                }
136                /* construct the preconditioner */
137                                
138                Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);              Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
139                Paso_Solver_setPreconditioner(A,options);              Paso_Solver_setPreconditioner(A,options);
140                Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);              Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
141                if (! Paso_noError()) return;              if (! Paso_noError()) return;
142                        
143                time_iter=Paso_timer();              time_iter=Paso_timer();
144                Paso_SystemMatrix_allocBuffer(A);              Paso_SystemMatrix_allocBuffer(A);
145                /* get an initial guess by evaluating the preconditioner */              /* get an initial guess by evaluating the preconditioner */
146                #pragma omp parallel  #pragma omp parallel
147                {              {
148                   Paso_Solver_solvePreconditioner(A,x,b);                 Paso_Solver_solvePreconditioner(A,x,b);
149                }              }
150                /* start the iteration process :*/              /* start the iteration process :*/
151                r=TMPMEMALLOC(numEqua,double);              r=TMPMEMALLOC(numEqua,double);
152                Paso_checkPtr(r);              Paso_checkPtr(r);
153                if (Paso_noError()) {              if (Paso_noError()) {
154                   totIter = 0;                 totIter = 0;
155                   finalizeIteration = FALSE;                 finalizeIteration = FALSE;
156                   last_norm2_of_residual=norm2_of_b;                 last_norm2_of_residual=norm2_of_b;
157                   last_norm_max_of_residual=norm_max_of_b;                 last_norm_max_of_residual=norm_max_of_b;
158                   while (! finalizeIteration) {                 while (! finalizeIteration) {
159                      cntIter = options->iter_max - totIter;                    cntIter = options->iter_max - totIter;
160                      finalizeIteration = TRUE;                    finalizeIteration = TRUE;
161                      /*     Set initial residual. */                    /*     Set initial residual. */
162                      norm2_of_residual = 0;                    norm2_of_residual = 0;
163                      norm_max_of_residual = 0;                    norm_max_of_residual = 0;
164                      #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)
165                      {                    {
166                         #pragma omp for private(i) schedule(static)  #pragma omp for private(i) schedule(static)
167                         for (i = 0; i < numEqua; i++) {                       for (i = 0; i < numEqua; i++) {
168                             r[i]=b[i];                          r[i]=b[i];
169                         }                       }
170                            
171                         Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);                       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);
172                            
173                         norm2_of_residual_local = 0;                       norm2_of_residual_local = 0;
174                         norm_max_of_residual_local = 0;                       norm_max_of_residual_local = 0;
175                         #pragma omp for private(i) schedule(static)  #pragma omp for private(i) schedule(static)
176                         for (i = 0; i < numEqua; i++) {                       for (i = 0; i < numEqua; i++) {
177                                norm2_of_residual_local+= r[i] * r[i];                          norm2_of_residual_local+= r[i] * r[i];
178                                norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);                          norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);
179                         }                       }
180                         #pragma omp critical  #pragma omp critical
181                         {                       {
182                            norm2_of_residual += norm2_of_residual_local;                          norm2_of_residual += norm2_of_residual_local;
183                            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);
184                         }                       }
185                      }                    }
186                      /* TODO: use one call */                    /* TODO: use one call */
187                      #ifdef PASO_MPI  #ifdef PASO_MPI
188                      {                    {
189                          loc_norm = norm2_of_residual;                       loc_norm = norm2_of_residual;
190                          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);
191                          loc_norm = norm_max_of_residual;                       loc_norm = norm_max_of_residual;
192                          MPI_Allreduce(&loc_norm,&norm_max_of_residual, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);                       MPI_Allreduce(&loc_norm,&norm_max_of_residual, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
193                      }                    }
194                      #endif  #endif
195                      norm2_of_residual =sqrt(norm2_of_residual);                    norm2_of_residual =sqrt(norm2_of_residual);
196                            
197                      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);
198                      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) {
199                           if (options->verbose) printf(" divergence!\n");                       if (options->verbose) printf(" divergence!\n");
200                           Paso_setError(WARNING, "No improvement during iteration. Iterative solver gives up.");                       Paso_setError(WARNING, "No improvement during iteration. Iterative solver gives up.");
201                      } else {                    } else {
202                         /* */                       /* */
203                         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 ) {
204                             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);
205                            if (options->verbose) printf(" (new tolerance = %e).\n",tol);                          if (options->verbose) printf(" (new tolerance = %e).\n",tol);
206                            last_norm2_of_residual=norm2_of_residual;                          last_norm2_of_residual=norm2_of_residual;
207                            last_norm_max_of_residual=norm_max_of_residual;                          last_norm_max_of_residual=norm_max_of_residual;
208                            /* call the solver */                          /* call the solver */
209                            switch (method) {                          switch (method) {
210                                   case PASO_BICGSTAB:                          case PASO_BICGSTAB:
211                                       errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol, pp);                             errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol, pp);
212                                       break;                             break;
213                                   case PASO_PCG:                          case PASO_PCG:
214                                       errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol, pp);                             errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol, pp);
215                                       break;                             break;
216                                   case PASO_PRES20:                          case PASO_PRES20:
217                                       errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20, pp);                             errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20, pp);
218                                       break;                             break;
219                                   case PASO_GMRES:                          case PASO_GMRES:
220                                       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);
221                                       break;                             break;
222                             }                          }
223                             totIter += cntIter;                          totIter += cntIter;
224                                                
225                             /* error handling  */                          /* error handling  */
226                             if (errorCode==NO_ERROR) {                          if (errorCode==NO_ERROR) {
227                                   finalizeIteration = FALSE;                             finalizeIteration = FALSE;
228                              } else if (errorCode==SOLVER_MAXITER_REACHED) {                          } else if (errorCode==SOLVER_MAXITER_REACHED) {
229                                    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.");
230                                    if (options->verbose) printf("Maximum number of iterations reached.!\n");                             if (options->verbose) printf("Maximum number of iterations reached.!\n");
231                             } else if (errorCode == SOLVER_INPUT_ERROR ) {                          } else if (errorCode == SOLVER_INPUT_ERROR ) {
232                                    Paso_setError(SYSTEM_ERROR,"illegal dimension in iterative solver.");                             Paso_setError(SYSTEM_ERROR,"illegal dimension in iterative solver.");
233                                    if (options->verbose) printf("Internal error!\n");                             if (options->verbose) printf("Internal error!\n");
234                             } else if ( errorCode == SOLVER_BREAKDOWN ) {                          } else if ( errorCode == SOLVER_BREAKDOWN ) {
235                              if (cntIter <= 1) {                             if (cntIter <= 1) {
236                                      Paso_setError(ZERO_DIVISION_ERROR, "fatal break down in iterative solver.");                                Paso_setError(ZERO_DIVISION_ERROR, "fatal break down in iterative solver.");
237                                      if (options->verbose) printf("Uncurable break down!\n");                                if (options->verbose) printf("Uncurable break down!\n");
238                                 } else {                             } else {
239                                  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);
240                                     finalizeIteration = FALSE;                                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");  
241                             }                             }
242                        } else {                          } else if (errorCode == SOLVER_MEMORY_ERROR) {
243                           if (options->verbose) printf(". convergence! \n");                             Paso_setError(MEMORY_ERROR,"memory allocation failed.");
244                        }                             if (options->verbose) printf("Memory allocation failed!\n");
245                     }                          } else if (errorCode !=SOLVER_NO_ERROR ) {
246                   } /* while */                             Paso_setError(SYSTEM_ERROR,"unidentified error in iterative solver.");
247                }                             if (options->verbose) printf("Unidentified error!\n");
248                MEMFREE(r);                          }
249                Paso_SystemMatrix_freeBuffer(A);                       } else {
250                time_iter=Paso_timer()-time_iter;                          if (options->verbose) printf(". convergence! \n");
251                if (options->verbose)  {                       }
252                   printf("timing: solver: %.4e sec\n",time_iter);                    }
253                   if (totIter>0) printf("timing: per iteration step: %.4e sec\n",time_iter/totIter);                 } /* while */
254                }              }
255             }              MEMFREE(r);
256          }              Paso_SystemMatrix_freeBuffer(A);
257                time_iter=Paso_timer()-time_iter;
258                if (options->verbose)  {
259                   printf("timing: solver: %.4e sec\n",time_iter);
260                   if (totIter>0) printf("timing: per iteration step: %.4e sec\n",time_iter/totIter);
261                }
262             }
263        }        }
264        Performance_stopMonitor(pp,PERFORMANCE_ALL);     }
265        blocktimer_increment("Paso_Solver()", blocktimer_start);     Performance_stopMonitor(pp,PERFORMANCE_ALL);
266       blocktimer_increment("Paso_Solver()", blocktimer_start);
267  }  }

Legend:
Removed from v.1504  
changed lines
  Added in v.1505

  ViewVC Help
Powered by ViewVC 1.1.26