/[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/BiCGStab.c revision 929 by gross, Wed Jan 17 07:41:13 2007 UTC temp/paso/src/BiCGStab.c revision 1387 by trankine, Fri Jan 11 07:45:26 2008 UTC
# Line 1  Line 1 
1    
2  /* $Id$ */  /* $Id$ */
3    
4  /*  /*******************************************************
5  ********************************************************************************   *
6  *               Copyright   2006 by ACcESS MNRF                                *   *           Copyright 2003-2007 by ACceSS MNRF
7  *                                                                              *   *       Copyright 2007 by University of Queensland
8  *                 http://www.access.edu.au                                     *   *
9  *           Primary Business: Queensland, Australia                            *   *                http://esscc.uq.edu.au
10  *     Licensed under the Open Software License version 3.0             *   *        Primary Business: Queensland, Australia
11  *        http://www.opensource.org/licenses/osl-3.0.php                        *   *  Licensed under the Open Software License version 3.0
12  ********************************************************************************   *     http://www.opensource.org/licenses/osl-3.0.php
13  */   *
14     *******************************************************/
15    
16  /*  /*
17     Crude modifications and translations for Paso by Matt Davies and Lutz Gross     Crude modifications and translations for Paso by Matt Davies and Lutz Gross
# Line 82  err_t Paso_Solver_BiCGStab( Line 84  err_t Paso_Solver_BiCGStab(
84    
85    
86    /* Local variables */    /* Local variables */
87    double *rtld=NULL,*p=NULL,*v=NULL,*t=NULL,*phat=NULL,*shat=NULL,*s=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;
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;
93    dim_t status = SOLVER_NO_ERROR;    dim_t status = SOLVER_NO_ERROR;
94      double *resid = tolerance;
95    /* adapt original routine parameters */    dim_t n = Paso_SystemMatrix_getTotalNumRows(A);
   dim_t n = A->num_cols * A-> col_block_size;;  
   double * resid = tolerance;  
96    
97    /* Executable Statements */    /* Executable Statements */
98    
# Line 104  err_t Paso_Solver_BiCGStab( Line 104  err_t Paso_Solver_BiCGStab(
104    phat=TMPMEMALLOC(n,double);    phat=TMPMEMALLOC(n,double);
105    shat=TMPMEMALLOC(n,double);    shat=TMPMEMALLOC(n,double);
106    s=TMPMEMALLOC(n,double);    s=TMPMEMALLOC(n,double);
   
107    /*     Test the input parameters. */    /*     Test the input parameters. */
108    
109    if (n < 0) {    if (n < 0) {
# Line 154  err_t Paso_Solver_BiCGStab( Line 153  err_t Paso_Solver_BiCGStab(
153      omegaDenumtr = 0.0;      omegaDenumtr = 0.0;
154        }        }
155        #pragma omp barrier        #pragma omp barrier
       #pragma ivdep  
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        rho = sum_1;        rho = sum_1;
# Line 164  err_t Paso_Solver_BiCGStab( Line 162  err_t Paso_Solver_BiCGStab(
162                
163      if (num_iter > 1) {      if (num_iter > 1) {
164        beta = rho / rho1 * (alpha / omega);        beta = rho / rho1 * (alpha / omega);
           #pragma ivdep  
165            #pragma omp for private(i0) schedule(static)            #pragma omp for private(i0) schedule(static)
166        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]);
167      } else {      } else {
           #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];        for (i0 = 0; i0 < n; i0++) p[i0] = r[i0];
170      }      }
# Line 178  err_t Paso_Solver_BiCGStab( Line 174  err_t Paso_Solver_BiCGStab(
174          Paso_Solver_solvePreconditioner(A,&phat[0], &p[0]);          Paso_Solver_solvePreconditioner(A,&phat[0], &p[0]);
175      Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &phat[0],ZERO, &v[0]);      Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &phat[0],ZERO, &v[0]);
176        
         // #pragma ivdep  
177          #pragma omp for private(i0) reduction(+:sum_2) schedule(static)          #pragma omp for private(i0) reduction(+:sum_2) schedule(static)
178      for (i0 = 0; i0 < n; i0++) sum_2 += rtld[i0] * v[i0];      for (i0 = 0; i0 < n; i0++) sum_2 += rtld[i0] * v[i0];
179          if (! (breakFlag = (ABS(sum_2) <= TOLERANCE_FOR_SCALARS))) {          if (! (breakFlag = (ABS(sum_2) <= TOLERANCE_FOR_SCALARS))) {
180         alpha = rho / sum_2;         alpha = rho / sum_2;
181    
            // #pragma ivdep  
182             #pragma omp for private(i0) reduction(+:sum_3) schedule(static)             #pragma omp for private(i0) reduction(+:sum_3) schedule(static)
183         for (i0 = 0; i0 < n; i0++) {         for (i0 = 0; i0 < n; i0++) {
184           r[i0] -= alpha * v[i0];           r[i0] -= alpha * v[i0];
# Line 195  err_t Paso_Solver_BiCGStab( Line 189  err_t Paso_Solver_BiCGStab(
189                    
190         /*        Early check for tolerance. */         /*        Early check for tolerance. */
191         if ( (convergeFlag = (norm_of_residual <= tol)) ) {         if ( (convergeFlag = (norm_of_residual <= tol)) ) {
              // #pragma ivdep  
192               #pragma omp for  private(i0) schedule(static)               #pragma omp for  private(i0) schedule(static)
193           for (i0 = 0; i0 < n; i0++) x[i0] += alpha * phat[i0];           for (i0 = 0; i0 < n; i0++) x[i0] += alpha * phat[i0];
194           maxIterFlag = FALSE;           maxIterFlag = FALSE;
# Line 205  err_t Paso_Solver_BiCGStab( Line 198  err_t Paso_Solver_BiCGStab(
198               Paso_Solver_solvePreconditioner(A,&shat[0], &s[0]);               Paso_Solver_solvePreconditioner(A,&shat[0], &s[0]);
199           Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &shat[0],ZERO,&t[0]);           Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &shat[0],ZERO,&t[0]);
200        
              // #pragma ivdep  
201               #pragma omp for private(i0) reduction(+:omegaNumtr,omegaDenumtr) schedule(static)               #pragma omp for private(i0) reduction(+:omegaNumtr,omegaDenumtr) schedule(static)
202           for (i0 = 0; i0 < n; i0++) {           for (i0 = 0; i0 < n; i0++) {
203             omegaNumtr +=t[i0] * s[i0];             omegaNumtr +=t[i0] * s[i0];
# Line 214  err_t Paso_Solver_BiCGStab( Line 206  err_t Paso_Solver_BiCGStab(
206               if (! (breakFlag = (ABS(omegaDenumtr) <= TOLERANCE_FOR_SCALARS))) {               if (! (breakFlag = (ABS(omegaDenumtr) <= TOLERANCE_FOR_SCALARS))) {
207              omega = omegaNumtr / omegaDenumtr;              omega = omegaNumtr / omegaDenumtr;
208        
                 // #pragma ivdep  
209                  #pragma omp for private(i0) reduction(+:sum_4) schedule(static)                  #pragma omp for private(i0) reduction(+:sum_4) schedule(static)
210              for (i0 = 0; i0 < n; i0++) {              for (i0 = 0; i0 < n; i0++) {
211                x[i0] += alpha * phat[i0] + omega * shat[i0];                x[i0] += alpha * phat[i0] + omega * shat[i0];
# Line 259  err_t Paso_Solver_BiCGStab( Line 250  err_t Paso_Solver_BiCGStab(
250    /*     End of BICGSTAB */    /*     End of BICGSTAB */
251    return status;    return status;
252  }  }
 /*  
  * $Log$  
  * Revision 1.2  2005/09/15 03:44:40  jgs  
  * Merge of development branch dev-02 back to main trunk on 2005-09-15  
  *  
  * Revision 1.1.2.1  2005/09/05 06:29:49  gross  
  * These files have been extracted from finley to define a stand alone libray for iterative  
  * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but  
  * has not been tested yet.  
  *  
  *  
  */  

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

  ViewVC Help
Powered by ViewVC 1.1.26