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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 682 - (hide annotations)
Mon Mar 27 02:43:09 2006 UTC (13 years, 5 months 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 jgs 150 /* $Id$ */
2    
3     /*
4 dhawcroft 631 ********************************************************************************
5 dhawcroft 633 * Copyright 2006 by ACcESS MNRF *
6 dhawcroft 631 * *
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 jgs 150 Crude modifications and translations for Paso by Matt Davies and Lutz Gross
16     */
17    
18 robwdcock 682 #include "../Paso.h"
19     #include "../SystemMatrix.h"
20 jgs 150 #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 gross 584 double * tolerance,
81     Paso_Performance* pp) {
82 jgs 150
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 gross 415 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &phat[0],ZERO, &v[0]);
174 jgs 150
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 gross 415 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &shat[0],ZERO,&t[0]);
198 jgs 150
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