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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 584 - (hide annotations)
Thu Mar 9 23:03:38 2006 UTC (14 years, 1 month 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 jgs 150 /* $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 gross 584 double * tolerance,
70     Paso_Performance* pp) {
71 jgs 150
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 gross 415 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &phat[0],ZERO, &v[0]);
163 jgs 150
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 gross 415 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, &shat[0],ZERO,&t[0]);
187 jgs 150
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