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

Revision 584 - (hide annotations)
Thu Mar 9 23:03:38 2006 UTC (14 years, 1 month ago) by gross
Original Path: trunk/paso/src/Solvers/BiCGStab.c
File MIME type: text/plain
File size: 7503 byte(s)
```eigenvalues: compiles and passes tests on altix now
```
 1 jgs 150 /* \$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 gross 584 double * tolerance, 70 Paso_Performance* pp) { 71 jgs 150 72 73 /* Local variables */ 74 double *rtld=NULL,*p=NULL,*v=NULL,*t=NULL,*phat=NULL,*shat=NULL,*s=NULL; 75 double beta,norm_of_residual,sum_1,sum_2,sum_3,sum_4,norm_of_residual_global; 76 double alpha, omega, omegaNumtr, omegaDenumtr, rho, tol, rho1; 77 dim_t num_iter=0,maxit,num_iter_global; 78 dim_t i0; 79 bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE; 80 dim_t status = SOLVER_NO_ERROR; 81 82 /* adapt original routine parameters */ 83 dim_t n = A->num_cols * A-> col_block_size;; 84 double * resid = tolerance; 85 86 /* Executable Statements */ 87 88 /* allocate memory: */ 89 rtld=TMPMEMALLOC(n,double); 90 p=TMPMEMALLOC(n,double); 91 v=TMPMEMALLOC(n,double); 92 t=TMPMEMALLOC(n,double); 93 phat=TMPMEMALLOC(n,double); 94 shat=TMPMEMALLOC(n,double); 95 s=TMPMEMALLOC(n,double); 96 97 /* Test the input parameters. */ 98 99 if (n < 0) { 100 status = SOLVER_INPUT_ERROR; 101 } else if (rtld==NULL || p==NULL || v==NULL || t==NULL || phat==NULL || shat==NULL || s==NULL) { 102 status = SOLVER_MEMORY_ERROR; 103 } else { 104 105 /* now bicgstab starts : */ 106 maxit = *iter; 107 tol = *resid; 108 109 #pragma omp parallel firstprivate(maxit,tol,convergeFlag,maxIterFlag,breakFlag) \ 110 private(rho,omega,num_iter,norm_of_residual,beta,alpha,rho1) 111 { 112 num_iter =0; 113 114 /* initialize arrays */ 115 116 #pragma omp for private(i0) schedule(static) 117 for (i0 = 0; i0 < n; i0++) { 118 rtld[i0]=0; 119 p[i0]=0; 120 v[i0]=0; 121 t[i0]=0; 122 phat[i0]=0; 123 shat[i0]=0; 124 } 125 #pragma omp for private(i0) schedule(static) 126 for (i0 = 0; i0 < n; i0++) rtld[i0] = r[i0]; 127 128 /* Perform BiConjugate Gradient Stabilized iteration. */ 129 130 L10: 131 ++(num_iter); 132 #pragma omp barrier 133 #pragma omp master 134 { 135 sum_1 = 0; 136 sum_2 = 0; 137 sum_3 = 0; 138 sum_4 = 0; 139 omegaNumtr = 0.0; 140 omegaDenumtr = 0.0; 141 } 142 #pragma omp barrier 143 #pragma omp for private(i0) reduction(+:sum_1) schedule(static) 144 for (i0 = 0; i0 < n; i0++) sum_1 += rtld[i0] * r[i0]; 145 rho = sum_1; 146 147 if (! (breakFlag = (ABS(rho) <= TOLERANCE_FOR_SCALARS))) { 148 /* Compute vector P. */ 149 150 if (num_iter > 1) { 151 beta = rho / rho1 * (alpha / omega); 152 #pragma omp for private(i0) schedule(static) 153 for (i0 = 0; i0 < n; i0++) p[i0] = r[i0] + beta * (p[i0] - omega * v[i0]); 154 } else { 155 #pragma omp for private(i0) schedule(static) 156 for (i0 = 0; i0 < n; i0++) p[i0] = r[i0]; 157 } 158 159 /* Compute direction adjusting vector PHAT and scalar ALPHA. */ 160 161 Paso_Solver_solvePreconditioner(A,&phat[0], &p[0]); 162 gross 415 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &phat[0],ZERO, &v[0]); 163 jgs 150 164 #pragma omp for private(i0) reduction(+:sum_2) schedule(static) 165 for (i0 = 0; i0 < n; i0++) sum_2 += rtld[i0] * v[i0]; 166 if (! (breakFlag = (ABS(sum_2) <= TOLERANCE_FOR_SCALARS))) { 167 alpha = rho / sum_2; 168 169 #pragma omp for private(i0) reduction(+:sum_3) schedule(static) 170 for (i0 = 0; i0 < n; i0++) { 171 r[i0] -= alpha * v[i0]; 172 s[i0] = r[i0]; 173 sum_3 += s[i0] * s[i0]; 174 } 175 norm_of_residual = sqrt(sum_3); 176 177 /* Early check for tolerance. */ 178 if ( (convergeFlag = (norm_of_residual <= tol)) ) { 179 #pragma omp for private(i0) schedule(static) 180 for (i0 = 0; i0 < n; i0++) x[i0] += alpha * phat[i0]; 181 maxIterFlag = FALSE; 182 breakFlag = FALSE; 183 } else { 184 /* Compute stabilizer vector SHAT and scalar OMEGA. */ 185 Paso_Solver_solvePreconditioner(A,&shat[0], &s[0]); 186 gross 415 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &shat[0],ZERO,&t[0]); 187 jgs 150 188 #pragma omp for private(i0) reduction(+:omegaNumtr,omegaDenumtr) schedule(static) 189 for (i0 = 0; i0 < n; i0++) { 190 omegaNumtr +=t[i0] * s[i0]; 191 omegaDenumtr += t[i0] * t[i0]; 192 } 193 if (! (breakFlag = (ABS(omegaDenumtr) <= TOLERANCE_FOR_SCALARS))) { 194 omega = omegaNumtr / omegaDenumtr; 195 196 #pragma omp for private(i0) reduction(+:sum_4) schedule(static) 197 for (i0 = 0; i0 < n; i0++) { 198 x[i0] += alpha * phat[i0] + omega * shat[i0]; 199 r[i0] -= omega * t[i0]; 200 sum_4 += r[i0] * r[i0]; 201 } 202 norm_of_residual = sqrt(sum_4); 203 convergeFlag = norm_of_residual <= tol; 204 maxIterFlag = num_iter == maxit; 205 breakFlag = (ABS(omega) <= TOLERANCE_FOR_SCALARS); 206 } 207 } 208 } 209 if (!(convergeFlag || maxIterFlag || breakFlag)) { 210 rho1 = rho; 211 goto L10; 212 } 213 } 214 /* end of iteration */ 215 #pragma omp master 216 { 217 num_iter_global=num_iter; 218 norm_of_residual_global=norm_of_residual; 219 if (maxIterFlag) { 220 status = SOLVER_MAXITER_REACHED; 221 } else if (breakFlag) { 222 status = SOLVER_BREAKDOWN; 223 } 224 } 225 } /* end of parallel region */ 226 } 227 TMPMEMFREE(rtld); 228 TMPMEMFREE(p); 229 TMPMEMFREE(v); 230 TMPMEMFREE(t); 231 TMPMEMFREE(phat); 232 TMPMEMFREE(shat); 233 TMPMEMFREE(s); 234 *iter=num_iter_global; 235 *resid=norm_of_residual_global; 236 237 /* End of BICGSTAB */ 238 return status; 239 } 240 /* 241 * \$Log\$ 242 * Revision 1.2 2005/09/15 03:44:40 jgs 243 * Merge of development branch dev-02 back to main trunk on 2005-09-15 244 * 245 * Revision 1.1.2.1 2005/09/05 06:29:49 gross 246 * These files have been extracted from finley to define a stand alone libray for iterative 247 * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but 248 * has not been tested yet. 249 * 250 * 251 */

Properties

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