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

trunk/paso/src/Solvers/BiCGStab.c revision 682 by robwdcock, Mon Mar 27 02:43:09 2006 UTC trunk/paso/src/BiCGStab.c revision 929 by gross, Wed Jan 17 07:41:13 2007 UTC
# Line 15  Line 15 
15     Crude modifications and translations for Paso by Matt Davies and Lutz Gross     Crude modifications and translations for Paso by Matt Davies and Lutz Gross
16  */  */
17    
18  #include "../Paso.h"  #include "Paso.h"
19  #include "../SystemMatrix.h"  #include "SystemMatrix.h"
20  #include "Solver.h"  #include "Solver.h"
21  #ifdef _OPENMP  #ifdef _OPENMP
22  #include <omp.h>  #include <omp.h>
# Line 117  err_t Paso_Solver_BiCGStab( Line 117  err_t Paso_Solver_BiCGStab(
117      maxit = *iter;      maxit = *iter;
118      tol = *resid;      tol = *resid;
119        
120  #pragma omp parallel firstprivate(maxit,tol,convergeFlag,maxIterFlag,breakFlag) \      #pragma omp parallel firstprivate(maxit,tol) \
121         private(rho,omega,num_iter,norm_of_residual,beta,alpha,rho1)         private(rho,omega,num_iter,norm_of_residual,beta,alpha,rho1, convergeFlag,maxIterFlag,breakFlag)
122      {      {
123        num_iter =0;        num_iter =0;
124          convergeFlag=FALSE;
125          maxIterFlag=FALSE;
126          breakFlag=FALSE;
127    
128        /* initialize arrays */        /* initialize arrays */
129        
# Line 151  err_t Paso_Solver_BiCGStab( Line 154  err_t Paso_Solver_BiCGStab(
154      omegaDenumtr = 0.0;      omegaDenumtr = 0.0;
155        }        }
156        #pragma omp barrier        #pragma omp barrier
157          #pragma ivdep
158        #pragma omp for private(i0) reduction(+:sum_1) schedule(static)        #pragma omp for private(i0) reduction(+:sum_1) schedule(static)
159        for (i0 = 0; i0 < n; i0++) sum_1 += rtld[i0] * r[i0];        for (i0 = 0; i0 < n; i0++) sum_1 += rtld[i0] * r[i0];
160        rho = sum_1;        rho = sum_1;
# Line 160  err_t Paso_Solver_BiCGStab( Line 164  err_t Paso_Solver_BiCGStab(
164                
165      if (num_iter > 1) {      if (num_iter > 1) {
166        beta = rho / rho1 * (alpha / omega);        beta = rho / rho1 * (alpha / omega);
167              #pragma ivdep
168            #pragma omp for private(i0) schedule(static)            #pragma omp for private(i0) schedule(static)
169        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]);
170      } else {      } else {
171              #pragma ivdep
172            #pragma omp for private(i0) schedule(static)            #pragma omp for private(i0) schedule(static)
173        for (i0 = 0; i0 < n; i0++) p[i0] = r[i0];        for (i0 = 0; i0 < n; i0++) p[i0] = r[i0];
174      }      }
# Line 172  err_t Paso_Solver_BiCGStab( Line 178  err_t Paso_Solver_BiCGStab(
178          Paso_Solver_solvePreconditioner(A,&phat[0], &p[0]);          Paso_Solver_solvePreconditioner(A,&phat[0], &p[0]);
179      Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &phat[0],ZERO, &v[0]);      Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &phat[0],ZERO, &v[0]);
180        
181            // #pragma ivdep
182          #pragma omp for private(i0) reduction(+:sum_2) schedule(static)          #pragma omp for private(i0) reduction(+:sum_2) schedule(static)
183      for (i0 = 0; i0 < n; i0++) sum_2 += rtld[i0] * v[i0];      for (i0 = 0; i0 < n; i0++) sum_2 += rtld[i0] * v[i0];
184          if (! (breakFlag = (ABS(sum_2) <= TOLERANCE_FOR_SCALARS))) {          if (! (breakFlag = (ABS(sum_2) <= TOLERANCE_FOR_SCALARS))) {
185         alpha = rho / sum_2;         alpha = rho / sum_2;
186    
187               // #pragma ivdep
188             #pragma omp for private(i0) reduction(+:sum_3) schedule(static)             #pragma omp for private(i0) reduction(+:sum_3) schedule(static)
189         for (i0 = 0; i0 < n; i0++) {         for (i0 = 0; i0 < n; i0++) {
190           r[i0] -= alpha * v[i0];           r[i0] -= alpha * v[i0];
# Line 187  err_t Paso_Solver_BiCGStab( Line 195  err_t Paso_Solver_BiCGStab(
195                    
196         /*        Early check for tolerance. */         /*        Early check for tolerance. */
197         if ( (convergeFlag = (norm_of_residual <= tol)) ) {         if ( (convergeFlag = (norm_of_residual <= tol)) ) {
198                 // #pragma ivdep
199               #pragma omp for  private(i0) schedule(static)               #pragma omp for  private(i0) schedule(static)
200           for (i0 = 0; i0 < n; i0++) x[i0] += alpha * phat[i0];           for (i0 = 0; i0 < n; i0++) x[i0] += alpha * phat[i0];
201           maxIterFlag = FALSE;           maxIterFlag = FALSE;
# Line 196  err_t Paso_Solver_BiCGStab( Line 205  err_t Paso_Solver_BiCGStab(
205               Paso_Solver_solvePreconditioner(A,&shat[0], &s[0]);               Paso_Solver_solvePreconditioner(A,&shat[0], &s[0]);
206           Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &shat[0],ZERO,&t[0]);           Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &shat[0],ZERO,&t[0]);
207        
208                 // #pragma ivdep
209               #pragma omp for private(i0) reduction(+:omegaNumtr,omegaDenumtr) schedule(static)               #pragma omp for private(i0) reduction(+:omegaNumtr,omegaDenumtr) schedule(static)
210           for (i0 = 0; i0 < n; i0++) {           for (i0 = 0; i0 < n; i0++) {
211             omegaNumtr +=t[i0] * s[i0];             omegaNumtr +=t[i0] * s[i0];
# Line 204  err_t Paso_Solver_BiCGStab( Line 214  err_t Paso_Solver_BiCGStab(
214               if (! (breakFlag = (ABS(omegaDenumtr) <= TOLERANCE_FOR_SCALARS))) {               if (! (breakFlag = (ABS(omegaDenumtr) <= TOLERANCE_FOR_SCALARS))) {
215              omega = omegaNumtr / omegaDenumtr;              omega = omegaNumtr / omegaDenumtr;
216        
217                    // #pragma ivdep
218                  #pragma omp for private(i0) reduction(+:sum_4) schedule(static)                  #pragma omp for private(i0) reduction(+:sum_4) schedule(static)
219              for (i0 = 0; i0 < n; i0++) {              for (i0 = 0; i0 < n; i0++) {
220                x[i0] += alpha * phat[i0] + omega * shat[i0];                x[i0] += alpha * phat[i0] + omega * shat[i0];

Legend:
Removed from v.682  
changed lines
  Added in v.929

  ViewVC Help
Powered by ViewVC 1.1.26