/[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 1553 by gross, Thu May 8 09:38:07 2008 UTC revision 1556 by gross, Mon May 12 00:54:58 2008 UTC
# Line 116  err_t Paso_Solver_BiCGStab( Line 116  err_t Paso_Solver_BiCGStab(
116      maxit = *iter;      maxit = *iter;
117      tol = *resid;      tol = *resid;
118        
119      #pragma omp parallel firstprivate(maxit,tol) \      num_iter =0;
120         private(rho,omega,num_iter,norm_of_residual,beta,alpha,rho1, convergeFlag,maxIterFlag,breakFlag)      convergeFlag=FALSE;
121      {      maxIterFlag=FALSE;
122        num_iter =0;      breakFlag=FALSE;
       convergeFlag=FALSE;  
       maxIterFlag=FALSE;  
       breakFlag=FALSE;  
123    
124        /* initialize arrays */      /* initialize arrays */
125      
126        #pragma omp for private(i0) schedule(static)      #pragma omp parallel for private(i0) schedule(static)
127        for (i0 = 0; i0 < n; i0++) {      for (i0 = 0; i0 < n; i0++) {
128      rtld[i0]=0;      rtld[i0]=0;
129      p[i0]=0;      p[i0]=0;
130      v[i0]=0;      v[i0]=0;
131      t[i0]=0;      t[i0]=0;
132      phat[i0]=0;      phat[i0]=0;
133      shat[i0]=0;      shat[i0]=0;
134        }          rtld[i0] = r[i0];
135        #pragma omp for private(i0) schedule(static)      }
       for (i0 = 0; i0 < n; i0++) rtld[i0] = r[i0];  
136        
137        /*     Perform BiConjugate Gradient Stabilized iteration. */      /*     Perform BiConjugate Gradient Stabilized iteration. */
138        
139      L10:      L10:
140        ++(num_iter);        ++(num_iter);
       #pragma omp barrier  
       #pragma omp master  
       {  
141      sum_1 = 0;      sum_1 = 0;
142      sum_2 = 0;      sum_2 = 0;
143      sum_3 = 0;      sum_3 = 0;
144      sum_4 = 0;      sum_4 = 0;
145      omegaNumtr = 0.0;      omegaNumtr = 0.0;
146      omegaDenumtr = 0.0;        omegaDenumtr = 0.0;
147        }        #pragma omp parallel for private(i0) reduction(+:sum_1) schedule(static)
       #pragma omp barrier  
       #pragma omp for private(i0) reduction(+:sum_1) schedule(static)  
148        for (i0 = 0; i0 < n; i0++) sum_1 += rtld[i0] * r[i0];        for (i0 = 0; i0 < n; i0++) sum_1 += rtld[i0] * r[i0];
149         #ifdef PASO_MPI        #ifdef PASO_MPI
150           #pragma omp master            loc_sum[0] = sum_1;
151           {            MPI_Allreduce(loc_sum, &sum_1, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
             loc_sum[0] = sum_1;  
             MPI_Allreduce(loc_sum, &sum_1, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);  
          }  
152        #endif        #endif
153        rho = sum_1;        rho = sum_1;
154                
# Line 169  err_t Paso_Solver_BiCGStab( Line 157  err_t Paso_Solver_BiCGStab(
157                
158      if (num_iter > 1) {      if (num_iter > 1) {
159        beta = rho / rho1 * (alpha / omega);        beta = rho / rho1 * (alpha / omega);
160            #pragma omp for private(i0) schedule(static)            #pragma omp parallel for private(i0) schedule(static)
161        for (i0 = 0; i0 < n; i0++) p[i0] = r[i0] + beta * (p[i0] - omega * v[i0]);        for (i0 = 0; i0 < n; i0++) p[i0] = r[i0] + beta * (p[i0] - omega * v[i0]);
162      } else {      } else {
163            #pragma omp for private(i0) schedule(static)            #pragma omp parallel for private(i0) schedule(static)
164        for (i0 = 0; i0 < n; i0++) p[i0] = r[i0];        for (i0 = 0; i0 < n; i0++) p[i0] = r[i0];
165      }      }
166        
# Line 181  err_t Paso_Solver_BiCGStab( Line 169  err_t Paso_Solver_BiCGStab(
169          Paso_Solver_solvePreconditioner(A,&phat[0], &p[0]);          Paso_Solver_solvePreconditioner(A,&phat[0], &p[0]);
170      Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &phat[0],ZERO, &v[0]);      Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &phat[0],ZERO, &v[0]);
171        
172          #pragma omp for private(i0) reduction(+:sum_2) schedule(static)          #pragma omp parallel for private(i0) reduction(+:sum_2) schedule(static)
173      for (i0 = 0; i0 < n; i0++) sum_2 += rtld[i0] * v[i0];      for (i0 = 0; i0 < n; i0++) sum_2 += rtld[i0] * v[i0];
174          #ifdef PASO_MPI          #ifdef PASO_MPI
175              #pragma omp master             loc_sum[0] = sum_2;
176              {              MPI_Allreduce(loc_sum, &sum_2, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
                loc_sum[0] = sum_2;  
                MPI_Allreduce(loc_sum, &sum_2, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);  
             }  
177          #endif          #endif
178          if (! (breakFlag = (ABS(sum_2) <= TOLERANCE_FOR_SCALARS))) {          if (! (breakFlag = (ABS(sum_2) <= TOLERANCE_FOR_SCALARS))) {
179         alpha = rho / sum_2;         alpha = rho / sum_2;
180    
181             #pragma omp for private(i0) reduction(+:sum_3) schedule(static)             #pragma omp parallel for private(i0) reduction(+:sum_3) schedule(static)
182         for (i0 = 0; i0 < n; i0++) {         for (i0 = 0; i0 < n; i0++) {
183           r[i0] -= alpha * v[i0];           r[i0] -= alpha * v[i0];
184           s[i0] = r[i0];           s[i0] = r[i0];
185           sum_3 += s[i0] * s[i0];           sum_3 += s[i0] * s[i0];
186         }         }
187             #ifdef PASO_MPI             #ifdef PASO_MPI
188                 #pragma omp master                 loc_sum[0] = sum_3;
189                 {                 MPI_Allreduce(loc_sum, &sum_3, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
                   loc_sum[0] = sum_3;  
                   MPI_Allreduce(loc_sum, &sum_3, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);  
                }  
190             #endif             #endif
191         norm_of_residual = sqrt(sum_3);         norm_of_residual = sqrt(sum_3);
192                    
193         /*        Early check for tolerance. */         /*        Early check for tolerance. */
194         if ( (convergeFlag = (norm_of_residual <= tol)) ) {         if ( (convergeFlag = (norm_of_residual <= tol)) ) {
195               #pragma omp for  private(i0) schedule(static)               #pragma omp parallel for  private(i0) schedule(static)
196           for (i0 = 0; i0 < n; i0++) x[i0] += alpha * phat[i0];           for (i0 = 0; i0 < n; i0++) x[i0] += alpha * phat[i0];
197           maxIterFlag = FALSE;           maxIterFlag = FALSE;
198           breakFlag = FALSE;           breakFlag = FALSE;
# Line 219  err_t Paso_Solver_BiCGStab( Line 201  err_t Paso_Solver_BiCGStab(
201               Paso_Solver_solvePreconditioner(A,&shat[0], &s[0]);               Paso_Solver_solvePreconditioner(A,&shat[0], &s[0]);
202           Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &shat[0],ZERO,&t[0]);           Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &shat[0],ZERO,&t[0]);
203        
204               #pragma omp for private(i0) reduction(+:omegaNumtr,omegaDenumtr) schedule(static)               #pragma omp parallel for private(i0) reduction(+:omegaNumtr,omegaDenumtr) schedule(static)
205           for (i0 = 0; i0 < n; i0++) {           for (i0 = 0; i0 < n; i0++) {
206             omegaNumtr +=t[i0] * s[i0];             omegaNumtr +=t[i0] * s[i0];
207             omegaDenumtr += t[i0] * t[i0];             omegaDenumtr += t[i0] * t[i0];
208           }           }
209               #ifdef PASO_MPI               #ifdef PASO_MPI
210                   #pragma omp master                  loc_sum[0] = omegaNumtr;
211                   {                  loc_sum[1] = omegaDenumtr;
212                      loc_sum[0] = omegaNumtr;                  MPI_Allreduce(loc_sum, sum, 2, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
213                      loc_sum[1] = omegaDenumtr;                  omegaNumtr=sum[0];
214                      MPI_Allreduce(loc_sum, sum, 2, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);                  omegaDenumtr=sum[1];
                     omegaNumtr=sum[0];  
                     omegaDenumtr=sum[1];  
                  }  
215               #endif               #endif
216               if (! (breakFlag = (ABS(omegaDenumtr) <= TOLERANCE_FOR_SCALARS))) {               if (! (breakFlag = (ABS(omegaDenumtr) <= TOLERANCE_FOR_SCALARS))) {
217              omega = omegaNumtr / omegaDenumtr;              omega = omegaNumtr / omegaDenumtr;
218        
219                  #pragma omp for private(i0) reduction(+:sum_4) schedule(static)                  #pragma omp parallel for private(i0) reduction(+:sum_4) schedule(static)
220              for (i0 = 0; i0 < n; i0++) {              for (i0 = 0; i0 < n; i0++) {
221                x[i0] += alpha * phat[i0] + omega * shat[i0];                x[i0] += alpha * phat[i0] + omega * shat[i0];
222                r[i0] = s[i0]-omega * t[i0];                r[i0] = s[i0]-omega * t[i0];
223                sum_4 += r[i0] * r[i0];                sum_4 += r[i0] * r[i0];
224              }              }
225                  #ifdef PASO_MPI                  #ifdef PASO_MPI
226                      #pragma omp master                     loc_sum[0] = sum_4;
227                      {                      MPI_Allreduce(loc_sum, &sum_4, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
                        loc_sum[0] = sum_4;  
                        MPI_Allreduce(loc_sum, &sum_4, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);  
                     }  
228                  #endif                  #endif
229              norm_of_residual = sqrt(sum_4);              norm_of_residual = sqrt(sum_4);
230              convergeFlag = norm_of_residual <= tol;              convergeFlag = norm_of_residual <= tol;
# Line 263  err_t Paso_Solver_BiCGStab( Line 239  err_t Paso_Solver_BiCGStab(
239      }      }
240        }        }
241        /* end of iteration */        /* end of iteration */
242        #pragma omp master        num_iter_global=num_iter;
243        {        norm_of_residual_global=norm_of_residual;
244      num_iter_global=num_iter;        if (maxIterFlag) {
     norm_of_residual_global=norm_of_residual;  
     if (maxIterFlag) {  
245          status = SOLVER_MAXITER_REACHED;          status = SOLVER_MAXITER_REACHED;
246      } else if (breakFlag) {        } else if (breakFlag) {
247          status = SOLVER_BREAKDOWN;          status = SOLVER_BREAKDOWN;
     }  
248        }        }
     }  /* end of parallel region */  
249    }    }
250    TMPMEMFREE(rtld);    TMPMEMFREE(rtld);
251    TMPMEMFREE(p);    TMPMEMFREE(p);

Legend:
Removed from v.1553  
changed lines
  Added in v.1556

  ViewVC Help
Powered by ViewVC 1.1.26