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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 682 - (show annotations)
Mon Mar 27 02:43:09 2006 UTC (14 years ago) by robwdcock
Original Path: trunk/paso/src/Solvers/BiCGStab.c
File MIME type: text/plain
File size: 8154 byte(s)
+ NEW BUILD SYSTEM

This commit contains the new build system with cross-platform support.
Most things work are before though you can have more control.

ENVIRONMENT settings have changed:
+ You no longer require LD_LIBRARY_PATH or PYTHONPATH to point to the
esysroot for building and testing performed via scons
+ ACcESS altix users: It is recommended you change your modules to load
the latest intel compiler and other libraries required by boost to match
the setup in svn (you can override). The correct modules are as follows

module load intel_cc.9.0.026
export
MODULEPATH=${MODULEPATH}:/data/raid2/toolspp4/modulefiles/gcc-3.3.6
module load boost/1.33.0/python-2.4.1
module load python/2.4.1
module load numarray/1.3.3


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 <omp.h>
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

  ViewVC Help
Powered by ViewVC 1.1.26