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

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

Parent Directory Parent Directory | Revision Log Revision Log


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


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     #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 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