/[escript]/trunk/paso/src/BiCGStab.c
ViewVC logotype

Diff of /trunk/paso/src/BiCGStab.c

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

revision 1552 by artak, Wed Feb 27 06:12:04 2008 UTC revision 1553 by gross, Thu May 8 09:38:07 2008 UTC
# Line 86  err_t Paso_Solver_BiCGStab( Line 86  err_t Paso_Solver_BiCGStab(
86    /* Local variables */    /* Local variables */
87    double *rtld=NULL,*p=NULL,*v=NULL,*t=NULL,*phat=NULL,*shat=NULL,*s=NULL, *buf1=NULL, *buf0=NULL;    double *rtld=NULL,*p=NULL,*v=NULL,*t=NULL,*phat=NULL,*shat=NULL,*s=NULL, *buf1=NULL, *buf0=NULL;
88    double beta,norm_of_residual,sum_1,sum_2,sum_3,sum_4,norm_of_residual_global;    double beta,norm_of_residual,sum_1,sum_2,sum_3,sum_4,norm_of_residual_global;
89    double alpha, omega, omegaNumtr, omegaDenumtr, rho, tol, rho1;    double alpha, omega, omegaNumtr, omegaDenumtr, rho, tol, rho1, loc_sum[2], sum[2];
90    dim_t num_iter=0,maxit,num_iter_global;    dim_t num_iter=0,maxit,num_iter_global;
91    dim_t i0;    dim_t i0;
92    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;
# Line 155  err_t Paso_Solver_BiCGStab( Line 155  err_t Paso_Solver_BiCGStab(
155        #pragma omp barrier        #pragma omp barrier
156        #pragma omp for private(i0) reduction(+:sum_1) schedule(static)        #pragma omp for private(i0) reduction(+:sum_1) schedule(static)
157        for (i0 = 0; i0 < n; i0++) sum_1 += rtld[i0] * r[i0];        for (i0 = 0; i0 < n; i0++) sum_1 += rtld[i0] * r[i0];
158           #ifdef PASO_MPI
159             #pragma omp master
160             {
161                loc_sum[0] = sum_1;
162                MPI_Allreduce(loc_sum, &sum_1, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
163             }
164          #endif
165        rho = sum_1;        rho = sum_1;
166                
167        if (! (breakFlag = (ABS(rho) <= TOLERANCE_FOR_SCALARS))) {        if (! (breakFlag = (ABS(rho) <= TOLERANCE_FOR_SCALARS))) {
# Line 176  err_t Paso_Solver_BiCGStab( Line 183  err_t Paso_Solver_BiCGStab(
183        
184          #pragma omp for private(i0) reduction(+:sum_2) schedule(static)          #pragma omp for private(i0) reduction(+:sum_2) schedule(static)
185      for (i0 = 0; i0 < n; i0++) sum_2 += rtld[i0] * v[i0];      for (i0 = 0; i0 < n; i0++) sum_2 += rtld[i0] * v[i0];
186            #ifdef PASO_MPI
187                #pragma omp master
188                {
189                   loc_sum[0] = sum_2;
190                   MPI_Allreduce(loc_sum, &sum_2, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
191                }
192            #endif
193          if (! (breakFlag = (ABS(sum_2) <= TOLERANCE_FOR_SCALARS))) {          if (! (breakFlag = (ABS(sum_2) <= TOLERANCE_FOR_SCALARS))) {
194         alpha = rho / sum_2;         alpha = rho / sum_2;
195    
# Line 185  err_t Paso_Solver_BiCGStab( Line 199  err_t Paso_Solver_BiCGStab(
199           s[i0] = r[i0];           s[i0] = r[i0];
200           sum_3 += s[i0] * s[i0];           sum_3 += s[i0] * s[i0];
201         }         }
202               #ifdef PASO_MPI
203                   #pragma omp master
204                   {
205                      loc_sum[0] = sum_3;
206                      MPI_Allreduce(loc_sum, &sum_3, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
207                   }
208               #endif
209         norm_of_residual = sqrt(sum_3);         norm_of_residual = sqrt(sum_3);
210                    
211         /*        Early check for tolerance. */         /*        Early check for tolerance. */
# Line 203  err_t Paso_Solver_BiCGStab( Line 224  err_t Paso_Solver_BiCGStab(
224             omegaNumtr +=t[i0] * s[i0];             omegaNumtr +=t[i0] * s[i0];
225             omegaDenumtr += t[i0] * t[i0];             omegaDenumtr += t[i0] * t[i0];
226           }           }
227                 #ifdef PASO_MPI
228                     #pragma omp master
229                     {
230                        loc_sum[0] = omegaNumtr;
231                        loc_sum[1] = omegaDenumtr;
232                        MPI_Allreduce(loc_sum, sum, 2, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
233                        omegaNumtr=sum[0];
234                        omegaDenumtr=sum[1];
235                     }
236                 #endif
237               if (! (breakFlag = (ABS(omegaDenumtr) <= TOLERANCE_FOR_SCALARS))) {               if (! (breakFlag = (ABS(omegaDenumtr) <= TOLERANCE_FOR_SCALARS))) {
238              omega = omegaNumtr / omegaDenumtr;              omega = omegaNumtr / omegaDenumtr;
239        
# Line 212  err_t Paso_Solver_BiCGStab( Line 243  err_t Paso_Solver_BiCGStab(
243                r[i0] = s[i0]-omega * t[i0];                r[i0] = s[i0]-omega * t[i0];
244                sum_4 += r[i0] * r[i0];                sum_4 += r[i0] * r[i0];
245              }              }
246                    #ifdef PASO_MPI
247                        #pragma omp master
248                        {
249                           loc_sum[0] = sum_4;
250                           MPI_Allreduce(loc_sum, &sum_4, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
251                        }
252                    #endif
253              norm_of_residual = sqrt(sum_4);              norm_of_residual = sqrt(sum_4);
254              convergeFlag = norm_of_residual <= tol;              convergeFlag = norm_of_residual <= tol;
255              maxIterFlag = num_iter == maxit;              maxIterFlag = num_iter == maxit;
Legend:
Removed from v.1552  
changed lines
  Added in v.1553
  ViewVC Help
Powered by ViewVC 1.1.26