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

Revision 1425 - (show annotations)
Wed Feb 27 06:12:04 2008 UTC (12 years, 7 months ago) by artak
File MIME type: text/plain
File size: 7602 byte(s)
```Added new test for BiCGStab in run_simplesolve.py
algorithmic change in BICGSTAB: from r=r-wt to r=s-wt

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

Properties

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