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

Revision 633 - (show annotations)
Thu Mar 23 05:37:00 2006 UTC (13 years, 11 months ago) by dhawcroft
Original Path: trunk/paso/src/Solvers/BiCGStab.c
File MIME type: text/plain
File size: 8148 byte(s)

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

## Properties

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