/[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 3368 by gross, Mon Oct 25 04:33:31 2010 UTC revision 3369 by gross, Fri Nov 19 06:26:11 2010 UTC
# Line 43  void Paso_Solver(Paso_SystemMatrix* A,do Line 43  void Paso_Solver(Paso_SystemMatrix* A,do
43     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;
44     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;
45     double norm_max_of_residual_local,norm_max_of_residual;     double norm_max_of_residual_local,norm_max_of_residual;
46     double last_norm_max_of_residual,*scaling;     double last_norm_max_of_residual;
47  #ifdef ESYS_MPI  #ifdef ESYS_MPI
48     double loc_norm;     double loc_norm;
49  #endif  #endif
50     dim_t i,totIter=0,cntIter,method;     dim_t i,totIter=0,cntIter,method;
51     bool_t finalizeIteration;     bool_t finalizeIteration;
52     err_t errorCode=SOLVER_NO_ERROR;     err_t errorCode=SOLVER_NO_ERROR;
53     dim_t numSol = Paso_SystemMatrix_getTotalNumCols(A);     const dim_t numSol = Paso_SystemMatrix_getTotalNumCols(A);
54     dim_t numEqua = Paso_SystemMatrix_getTotalNumRows(A);     const dim_t numEqua = Paso_SystemMatrix_getTotalNumRows(A);
55     double blocktimer_precond, blocktimer_start = blocktimer_time();     double blocktimer_precond, blocktimer_start = blocktimer_time();
56     double *x0=NULL;     double *x0=NULL;
57    
# Line 94  void Paso_Solver(Paso_SystemMatrix* A,do Line 94  void Paso_Solver(Paso_SystemMatrix* A,do
94          return;          return;
95       }       }
96            
97         r=TMPMEMALLOC(numEqua,double);
98         x0=TMPMEMALLOC(numEqua,double);
99         Esys_checkPtr(r);
100         Esys_checkPtr(x0);
101         Paso_SystemMatrix_balance(A);
102         options->num_level=0;
103         options->num_inner_iter=0;
104        
105       /* ========================= */       /* ========================= */
      Performance_startMonitor(pp,PERFORMANCE_ALL);  
106       if (Esys_noError()) {       if (Esys_noError()) {
107          /* get normalization */      Performance_startMonitor(pp,PERFORMANCE_ALL);
108          scaling=Paso_SystemMatrix_borrowNormalization(A);      
109          if (Esys_noError()) {      Paso_SystemMatrix_applyBalance(A, r, b, TRUE);
110             /* get the norm of the right hand side */           /* get the norm of the right hand side */
111             norm2_of_b=0.;          norm2_of_b=0.;
112             norm_max_of_b=0.;          norm_max_of_b=0.;
113             #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)
114             {          {
115                  norm2_of_b_local=0.;                  norm2_of_b_local=0.;
116                  norm_max_of_b_local=0.;                  norm_max_of_b_local=0.;
117                  #pragma omp for private(i) schedule(static)                  #pragma omp for private(i) schedule(static)
118                  for (i = 0; i < numEqua ; ++i) {                  for (i = 0; i < numEqua ; ++i) {
119                        norm2_of_b_local += b[i] * b[i];                        norm2_of_b_local += r[i] * r[i];
120                        norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local);                        norm_max_of_b_local = MAX(ABS(r[i]),norm_max_of_b_local);
121                  }                  }
122                  #pragma omp critical                  #pragma omp critical
123                  {                  {
124                     norm2_of_b += norm2_of_b_local;                     norm2_of_b += norm2_of_b_local;
125                     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);
126                  }                  }
127             }          }
128             /* TODO: use one call */          /* TODO: use one call */
129             #ifdef ESYS_MPI          #ifdef ESYS_MPI
130             {          {
131                 loc_norm = norm2_of_b;                 loc_norm = norm2_of_b;
132                 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);
133                 loc_norm = norm_max_of_b;                 loc_norm = norm_max_of_b;
134                 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);
135             }          }
136             #endif          #endif
137             norm2_of_b=sqrt(norm2_of_b);          norm2_of_b=sqrt(norm2_of_b);
138           /* if norm2_of_b==0 we are ready: x=0 */          /* if norm2_of_b==0 we are ready: x=0 */
139           if ( IS_NAN(norm2_of_b) || IS_NAN(norm_max_of_b) ) {          if ( IS_NAN(norm2_of_b) || IS_NAN(norm_max_of_b) ) {
140              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.");
141           } else if (norm2_of_b <=0.) {          } else if (norm2_of_b <=0.) {
142    
143              #pragma omp parallel for private(i) schedule(static)              #pragma omp parallel for private(i) schedule(static)
144              for (i = 0; i < numSol; i++) x[i]=0.;              for (i = 0; i < numSol; i++) x[i]=0.;
145              if (options->verbose) printf("right hand side is identical zero.\n");              if (options->verbose) printf("right hand side is identical zero.\n");
146    
147           } else {          } else {
148    
149              if (options->verbose) {              if (options->verbose) {
150                 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);
# Line 176  void Paso_Solver(Paso_SystemMatrix* A,do Line 183  void Paso_Solver(Paso_SystemMatrix* A,do
183              Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);              Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
184              blocktimer_increment("Paso_Solver_setPreconditioner()", blocktimer_precond);              blocktimer_increment("Paso_Solver_setPreconditioner()", blocktimer_precond);
185              options->set_up_time=Esys_timer()-time_iter;              options->set_up_time=Esys_timer()-time_iter;
186              if (! Esys_noError()) return;              if (Esys_noError()) {
187    
188    
189                /* get an initial guess by evaluating the preconditioner */                /* get an initial guess by evaluating the preconditioner */
190            Paso_SystemMatrix_solvePreconditioner(A,x,b);            Paso_SystemMatrix_solvePreconditioner(A,x,r);
191    
192                /* start the iteration process :*/                totIter = 1;
193                r=TMPMEMALLOC(numEqua,double);                finalizeIteration = FALSE;
194                x0=TMPMEMALLOC(numEqua,double);                last_norm2_of_residual=norm2_of_b;
195                Esys_checkPtr(r);                last_norm_max_of_residual=norm_max_of_b;
196            Esys_checkPtr(x0);            net_time_start=Esys_timer();
               if (Esys_noError()) {  
   
                  totIter = 1;  
                  finalizeIteration = FALSE;  
                  last_norm2_of_residual=norm2_of_b;  
                  last_norm_max_of_residual=norm_max_of_b;  
          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;              if (totIter>1) Paso_SystemMatrix_applyBalance(A, r, b, TRUE); /* in the first iteration r = balance * b already */
205                      norm_max_of_residual = 0;              
                     #pragma omp parallel for private(i) schedule(static)  
                     for (i = 0; i < numEqua; i++) r[i]=b[i];  
206                      Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);                      Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);
207                      #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)                      
208                      {                      norm2_of_residual = 0;
209                norm_max_of_residual = 0;
210                #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
211                {
212                         norm2_of_residual_local = 0;                         norm2_of_residual_local = 0;
213                         norm_max_of_residual_local = 0;                         norm_max_of_residual_local = 0;
214                         #pragma omp for private(i) schedule(static)                         #pragma omp for private(i) schedule(static)
215                         for (i = 0; i < numEqua; i++) {                         for (i = 0; i < numEqua; i++) {
216                                norm2_of_residual_local+= r[i] * r[i];                                norm2_of_residual_local+= r[i] * r[i];
217                                norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);                                norm_max_of_residual_local=MAX(ABS(r[i]),norm_max_of_residual_local);
218                         }                         }
219                         #pragma omp critical                         #pragma omp critical
220                         {                         {
# Line 231  void Paso_Solver(Paso_SystemMatrix* A,do Line 232  void Paso_Solver(Paso_SystemMatrix* A,do
232                      norm2_of_residual =sqrt(norm2_of_residual);                      norm2_of_residual =sqrt(norm2_of_residual);
233                      options->residual_norm=norm2_of_residual;                      options->residual_norm=norm2_of_residual;
234    
235                    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);
236    
237                    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) {
238    
239                       if (options->verbose) printf(" divergence!\n");              if (options->verbose) printf(" divergence!\n");
240                       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.");
241    
242                    } else {               } else {
243                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 ) {  
244    
245                          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);
246                          if (options->verbose) printf(" (new tolerance = %e).\n",tol);              if (options->verbose) printf(" (new tolerance = %e).\n",tol);
247    
248                          last_norm2_of_residual=norm2_of_residual;                          last_norm2_of_residual=norm2_of_residual;
249                          last_norm_max_of_residual=norm_max_of_residual;                          last_norm_max_of_residual=norm_max_of_residual;
# Line 311  void Paso_Solver(Paso_SystemMatrix* A,do Line 310  void Paso_Solver(Paso_SystemMatrix* A,do
310                   } /* while */                   } /* while */
311                   options->net_time=Esys_timer()-net_time_start;                   options->net_time=Esys_timer()-net_time_start;
312                }                }
               MEMFREE(r);  
               MEMFREE(x0);  
               options->num_level=0;  
               options->num_inner_iter=0;  
313                options->num_iter=totIter;                options->num_iter=totIter;
314              Paso_SystemMatrix_applyBalanceInPlace(A, x, FALSE);
315             }             }
316          }      
317        }       }
318     options->time=Esys_timer()-time_iter;       MEMFREE(r);
319     Performance_stopMonitor(pp,PERFORMANCE_ALL);       MEMFREE(x0);
320     blocktimer_increment("Paso_Solver()", blocktimer_start);       options->time=Esys_timer()-time_iter;
321         Performance_stopMonitor(pp,PERFORMANCE_ALL);
322         blocktimer_increment("Paso_Solver()", blocktimer_start);
323  }  }

Legend:
Removed from v.3368  
changed lines
  Added in v.3369

  ViewVC Help
Powered by ViewVC 1.1.26