# Contents of /trunk/paso/src/Solvers/BiCGStab.c

Revision 150 - (show annotations)
Thu Sep 15 03:44:45 2005 UTC (14 years, 6 months ago) by jgs
Original Path: trunk/esys2/paso/src/Solvers/BiCGStab.c
File MIME type: text/plain
File size: 7453 byte(s)
```Merge of development branch dev-02 back to main trunk on 2005-09-15

```
 1 /* \$Id\$ */ 2 3 /* 4 Crude modifications and translations for Paso by Matt Davies and Lutz Gross 5 */ 6 7 #include "Paso.h" 8 #include "SystemMatrix.h" 9 #include "Solver.h" 10 #ifdef _OPENMP 11 #include 12 #endif 13 14 /* -- Iterative template routine -- 15 * Univ. of Tennessee and Oak Ridge National Laboratory 16 * October 1, 1993 17 * Details of this algorithm are described in "Templates for the 18 * Solution of Linear Systems: Building Blocks for Iterative 19 * Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra, 20 * Eijkhout, Pozo, Romine, and van der Vorst, SIAM Publications, 21 * 1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps). 22 * 23 * Purpose 24 * ======= 25 * 26 * BICGSTAB solves the linear system A*x = b using the 27 * BiConjugate Gradient Stabilized iterative method with 28 * preconditioning. 29 * 30 * Convergence test: norm( b - A*x )< TOL. 31 * For other measures, see the above reference. 32 * 33 * Arguments 34 * ========= 35 * 36 * A (input) 37 * 38 * R (input) DOUBLE PRECISION array, dimension N. 39 * On entry, residual of inital guess X 40 * 41 * X (input/output) DOUBLE PRECISION array, dimension N. 42 * On input, the initial guess. 43 * 44 * ITER (input/output) INT 45 * On input, the maximum iterations to be performed. 46 * On output, actual number of iterations performed. 47 * 48 * RESID (input/output) DOUBLE PRECISION 49 * On input, the allowable convergence measure for 50 * norm( b - A*x ) 51 * On output, the final value of this measure. 52 * 53 * return value 54 * 55 * = SOLVER_NO_ERROR: Successful exit. Iterated approximate solution returned. 56 * = SOLVER_MAXITER_REACHED 57 * = SOLVER_INPUT_ERROR Illegal parameter: 58 * = SOLVER_BREAKDOWN: If parameters RHO or OMEGA become smaller 59 * = SOLVER_MEMORY_ERROR : If parameters RHO or OMEGA become smaller 60 * 61 * ============================================================== 62 */ 63 64 err_t Paso_Solver_BiCGStab( 65 Paso_SystemMatrix * A, 66 double * r, 67 double * x, 68 dim_t *iter, 69 double * tolerance) { 70 71 72 /* Local variables */ 73 double *rtld=NULL,*p=NULL,*v=NULL,*t=NULL,*phat=NULL,*shat=NULL,*s=NULL; 74 double beta,norm_of_residual,sum_1,sum_2,sum_3,sum_4,norm_of_residual_global; 75 double alpha, omega, omegaNumtr, omegaDenumtr, rho, tol, rho1; 76 dim_t num_iter=0,maxit,num_iter_global; 77 dim_t i0; 78 bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE; 79 dim_t status = SOLVER_NO_ERROR; 80 81 /* adapt original routine parameters */ 82 dim_t n = A->num_cols * A-> col_block_size;; 83 double * resid = tolerance; 84 85 /* Executable Statements */ 86 87 /* allocate memory: */ 88 rtld=TMPMEMALLOC(n,double); 89 p=TMPMEMALLOC(n,double); 90 v=TMPMEMALLOC(n,double); 91 t=TMPMEMALLOC(n,double); 92 phat=TMPMEMALLOC(n,double); 93 shat=TMPMEMALLOC(n,double); 94 s=TMPMEMALLOC(n,double); 95 96 /* Test the input parameters. */ 97 98 if (n < 0) { 99 status = SOLVER_INPUT_ERROR; 100 } else if (rtld==NULL || p==NULL || v==NULL || t==NULL || phat==NULL || shat==NULL || s==NULL) { 101 status = SOLVER_MEMORY_ERROR; 102 } else { 103 104 /* now bicgstab starts : */ 105 maxit = *iter; 106 tol = *resid; 107 108 #pragma omp parallel firstprivate(maxit,tol,convergeFlag,maxIterFlag,breakFlag) \ 109 private(rho,omega,num_iter,norm_of_residual,beta,alpha,rho1) 110 { 111 num_iter =0; 112 113 /* initialize arrays */ 114 115 #pragma omp for private(i0) schedule(static) 116 for (i0 = 0; i0 < n; i0++) { 117 rtld[i0]=0; 118 p[i0]=0; 119 v[i0]=0; 120 t[i0]=0; 121 phat[i0]=0; 122 shat[i0]=0; 123 } 124 #pragma omp for private(i0) schedule(static) 125 for (i0 = 0; i0 < n; i0++) rtld[i0] = r[i0]; 126 127 /* Perform BiConjugate Gradient Stabilized iteration. */ 128 129 L10: 130 ++(num_iter); 131 #pragma omp barrier 132 #pragma omp master 133 { 134 sum_1 = 0; 135 sum_2 = 0; 136 sum_3 = 0; 137 sum_4 = 0; 138 omegaNumtr = 0.0; 139 omegaDenumtr = 0.0; 140 } 141 #pragma omp barrier 142 #pragma omp for private(i0) reduction(+:sum_1) schedule(static) 143 for (i0 = 0; i0 < n; i0++) sum_1 += rtld[i0] * r[i0]; 144 rho = sum_1; 145 146 if (! (breakFlag = (ABS(rho) <= TOLERANCE_FOR_SCALARS))) { 147 /* Compute vector P. */ 148 149 if (num_iter > 1) { 150 beta = rho / rho1 * (alpha / omega); 151 #pragma omp for private(i0) schedule(static) 152 for (i0 = 0; i0 < n; i0++) p[i0] = r[i0] + beta * (p[i0] - omega * v[i0]); 153 } else { 154 #pragma omp for private(i0) schedule(static) 155 for (i0 = 0; i0 < n; i0++) p[i0] = r[i0]; 156 } 157 158 /* Compute direction adjusting vector PHAT and scalar ALPHA. */ 159 160 Paso_Solver_solvePreconditioner(A,&phat[0], &p[0]); 161 Paso_SystemMatrix_MatrixVector(ONE, A, &phat[0],ZERO, &v[0]); 162 163 #pragma omp for private(i0) reduction(+:sum_2) schedule(static) 164 for (i0 = 0; i0 < n; i0++) sum_2 += rtld[i0] * v[i0]; 165 if (! (breakFlag = (ABS(sum_2) <= TOLERANCE_FOR_SCALARS))) { 166 alpha = rho / sum_2; 167 168 #pragma omp for private(i0) reduction(+:sum_3) schedule(static) 169 for (i0 = 0; i0 < n; i0++) { 170 r[i0] -= alpha * v[i0]; 171 s[i0] = r[i0]; 172 sum_3 += s[i0] * s[i0]; 173 } 174 norm_of_residual = sqrt(sum_3); 175 176 /* Early check for tolerance. */ 177 if ( (convergeFlag = (norm_of_residual <= tol)) ) { 178 #pragma omp for private(i0) schedule(static) 179 for (i0 = 0; i0 < n; i0++) x[i0] += alpha * phat[i0]; 180 maxIterFlag = FALSE; 181 breakFlag = FALSE; 182 } else { 183 /* Compute stabilizer vector SHAT and scalar OMEGA. */ 184 Paso_Solver_solvePreconditioner(A,&shat[0], &s[0]); 185 Paso_SystemMatrix_MatrixVector(ONE, A, &shat[0],ZERO,&t[0]); 186 187 #pragma omp for private(i0) reduction(+:omegaNumtr,omegaDenumtr) schedule(static) 188 for (i0 = 0; i0 < n; i0++) { 189 omegaNumtr +=t[i0] * s[i0]; 190 omegaDenumtr += t[i0] * t[i0]; 191 } 192 if (! (breakFlag = (ABS(omegaDenumtr) <= TOLERANCE_FOR_SCALARS))) { 193 omega = omegaNumtr / omegaDenumtr; 194 195 #pragma omp for private(i0) reduction(+:sum_4) schedule(static) 196 for (i0 = 0; i0 < n; i0++) { 197 x[i0] += alpha * phat[i0] + omega * shat[i0]; 198 r[i0] -= omega * t[i0]; 199 sum_4 += r[i0] * r[i0]; 200 } 201 norm_of_residual = sqrt(sum_4); 202 convergeFlag = norm_of_residual <= tol; 203 maxIterFlag = num_iter == maxit; 204 breakFlag = (ABS(omega) <= TOLERANCE_FOR_SCALARS); 205 } 206 } 207 } 208 if (!(convergeFlag || maxIterFlag || breakFlag)) { 209 rho1 = rho; 210 goto L10; 211 } 212 } 213 /* end of iteration */ 214 #pragma omp master 215 { 216 num_iter_global=num_iter; 217 norm_of_residual_global=norm_of_residual; 218 if (maxIterFlag) { 219 status = SOLVER_MAXITER_REACHED; 220 } else if (breakFlag) { 221 status = SOLVER_BREAKDOWN; 222 } 223 } 224 } /* end of parallel region */ 225 } 226 TMPMEMFREE(rtld); 227 TMPMEMFREE(p); 228 TMPMEMFREE(v); 229 TMPMEMFREE(t); 230 TMPMEMFREE(phat); 231 TMPMEMFREE(shat); 232 TMPMEMFREE(s); 233 *iter=num_iter_global; 234 *resid=norm_of_residual_global; 235 236 /* End of BICGSTAB */ 237 return status; 238 } 239 /* 240 * \$Log\$ 241 * Revision 1.2 2005/09/15 03:44:40 jgs 242 * Merge of development branch dev-02 back to main trunk on 2005-09-15 243 * 244 * Revision 1.1.2.1 2005/09/05 06:29:49 gross 245 * These files have been extracted from finley to define a stand alone libray for iterative 246 * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but 247 * has not been tested yet. 248 * 249 * 250 */

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision