/[escript]/trunk/paso/src/BiCGStab.c
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 584 - (show annotations)
Thu Mar 9 23:03:38 2006 UTC (13 years, 6 months 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 /* $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 <omp.h>
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 Paso_Performance* pp) {
71
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 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &phat[0],ZERO, &v[0]);
163
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 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &shat[0],ZERO,&t[0]);
187
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

  ViewVC Help
Powered by ViewVC 1.1.26