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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1860 - (show annotations)
Wed Oct 8 04:07:10 2008 UTC (11 years, 1 month ago) by artak
File MIME type: text/plain
File size: 6990 byte(s)
minor bug is fixed
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
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 /* TFQMR iterations */
16
17 #include "SystemMatrix.h"
18 #include "Paso.h"
19 #include "Solver.h"
20 #include "PasoUtil.h"
21
22 #ifdef _OPENMP
23 #include <omp.h>
24 #endif
25
26 #ifdef PASO_MPI
27 #include <mpi.h>
28 #endif
29
30 /*
31 *
32 * Purpose
33 * =======
34 *
35 * MINRES solves the linear system A*x = b
36 *
37 * Convergence test: norm( b - A*x )< TOL.
38 * For other measures, see the above reference.
39 *
40 * Arguments
41 * =========
42 *
43 * r (input) DOUBLE PRECISION array, dimension N.
44 * On entry, residual of inital guess x
45 *
46 * x (input/output) DOUBLE PRECISION array, dimension N.
47 * On input, the initial guess.
48 *
49 * ITER (input/output) INT
50 * On input, the maximum iterations to be performed.
51 * On output, actual number of iterations performed.
52 *
53 * INFO (output) INT
54 *
55 * = SOLVER_NO_ERROR: Successful exit. Iterated approximate solution returned.
56 * = SOLVER_MAXITER_REACHED
57 * = SOLVER_INPUT_ERROR Illegal parameter:
58 * = SOLVER_MEMORY_ERROR :
59 * = SOLVER_NEGATIVE_NORM_ERROR :
60 *
61 * ==============================================================
62 */
63
64 /* #define PASO_DYNAMIC_SCHEDULING_MVM */
65
66 #if defined PASO_DYNAMIC_SCHEDULING_MVM && defined __OPENMP
67 #define USE_DYNAMIC_SCHEDULING
68 #endif
69
70 err_t Paso_Solver_MINRES(
71 Paso_SystemMatrix * A,
72 double * r,
73 double * x,
74 dim_t *iter,
75 double *tol,
76 double * tolerance,
77 Paso_Performance* pp) {
78
79 /* Local variables */
80
81 dim_t num_iter=0,maxit;
82 bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;
83 err_t status = SOLVER_NO_ERROR;
84 dim_t n = Paso_SystemMatrix_getTotalNumRows(A);
85 double *w=NULL, *w1=NULL, *w2=NULL, *r1=NULL, *r2=NULL, *y=NULL, *v=NULL;
86
87 double Anorm,ynorm,oldb,dbar,epsln,phibar,rhs1,rhs2,rnorm,tnorm2,ynorm2,cs,sn,eps,s,alfa,denom,z,beta1,beta;
88 double gmax,gmin,oldeps,delta,gbar,gamma,phi;
89
90 double norm_of_residual;
91
92 /* */
93 /*-----------------------------------------------------------------*/
94 /* */
95 /* Start of Calculation : */
96 /* --------------------- */
97 /* */
98 /* */
99 w=TMPMEMALLOC(n,double);
100 w1=TMPMEMALLOC(n,double);
101 w2=TMPMEMALLOC(n,double);
102 r1=TMPMEMALLOC(n,double);
103 r2=TMPMEMALLOC(n,double);
104 y=TMPMEMALLOC(n,double);
105 v=TMPMEMALLOC(n,double);
106
107
108 if (w ==NULL || w1== NULL || w2== NULL || r1 == NULL || r2== NULL || y==NULL || v==NULL ) {
109 status=SOLVER_MEMORY_ERROR;
110 }
111
112 maxit = *iter;
113
114 /* Test the input parameters. */
115 if (n < 0 || maxit<=0 ) {
116 status=SOLVER_INPUT_ERROR;
117 }
118
119 Paso_Copy(n,r1,r);
120
121 Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);
122 Paso_Solver_solvePreconditioner(A,y,r1);
123 Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);
124
125 beta1=Paso_InnerProduct(n,r1,y,A->mpi_info);
126 if (beta1<0) {
127 status=SOLVER_NEGATIVE_NORM_ERROR;
128 }
129
130 if (beta1>0) {
131 beta1=sqrt(beta1);
132 }
133
134 Performance_startMonitor(pp,PERFORMANCE_SOLVER);
135
136 Paso_zeroes(n,w);
137 Paso_zeroes(n,w2);
138
139 Paso_Copy(n,r2,r1);
140
141 Anorm = 0;
142 ynorm = 0;
143 oldb = 0;
144 beta = beta1;
145 dbar = 0;
146 epsln = 0;
147 phibar = beta1;
148 rhs1 = beta1;
149 rhs2 = 0;
150 rnorm = phibar;
151 tnorm2 = 0;
152 ynorm2 = 0;
153 cs = -1;
154 sn = 0;
155 eps = 0.0001;
156
157 while (!(convergeFlag || maxIterFlag || breakFlag || (status !=SOLVER_NO_ERROR) ))
158 {
159
160 s=1/beta;
161 Paso_Update(n, 0., v, s, y);
162
163 Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
164 Performance_startMonitor(pp,PERFORMANCE_MVM);
165 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, v,ZERO,y);
166 Performance_stopMonitor(pp,PERFORMANCE_MVM);
167 Performance_startMonitor(pp,PERFORMANCE_SOLVER);
168
169 if (num_iter >= 1) {
170 Paso_Update(n, 1., y, -(beta/oldb), r1);
171 }
172
173 alfa = Paso_InnerProduct(n,v,y,A->mpi_info);
174 Paso_Update(n, 1., y, (-alfa/beta), r2);
175 Paso_Copy(n,r1,r2);
176 Paso_Copy(n,r2,y);
177
178 Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
179 Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);
180 Paso_Solver_solvePreconditioner(A,y,r2);
181 Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);
182 Performance_startMonitor(pp,PERFORMANCE_SOLVER);
183
184 oldb = beta;
185 beta = Paso_InnerProduct(n,y,r2,A->mpi_info);
186 if (beta<0) {
187 status=SOLVER_NEGATIVE_NORM_ERROR;
188 }
189
190 beta = sqrt( beta );
191 tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta;
192
193 if (num_iter==0) {
194 gmax = ABS(alfa);
195 gmin = gmax;
196 }
197
198 /* Apply previous rotation Qk-1 to get
199 [deltak epslnk+1] = [cs sn][dbark 0 ]
200 [gbar k dbar k+1] [sn -cs][alfak betak+1]. */
201
202 oldeps = epsln;
203 delta = cs * dbar + sn * alfa ;
204 gbar = sn * dbar - cs * alfa ;
205 epsln = sn * beta ;
206 dbar = - cs * beta;
207
208 gamma = sqrt(gbar*gbar+beta*beta) ;
209 gamma = MAX(gamma,eps) ;
210 cs = gbar / gamma ;
211 sn = beta / gamma ;
212 phi = cs * phibar ;
213 phibar = sn * phibar ;
214
215 /* Update x. */
216
217 denom = 1/gamma ;
218 Paso_Copy(n,w1,w2);
219 Paso_Copy(n,w2,w);
220
221 Paso_LinearCombination(n,w,denom,v,-(denom*oldeps),w1);
222 Paso_Update(n, 1., w, -(delta*denom), w2) ;
223 Paso_Update(n, 1., x, phi, w) ;
224
225 /* Go round again. */
226
227 gmax = MAX(gmax,gamma);
228 gmin = MIN(gmin,gamma);
229 z = rhs1 / gamma;
230 ynorm2 = z*z + ynorm2;
231 rhs1 = rhs2 - delta*z;
232 rhs2 = - epsln*z;
233
234 Anorm = sqrt( tnorm2 ) ;
235 ynorm = sqrt( ynorm2 ) ;
236
237 rnorm = phibar;
238
239 maxIterFlag = (num_iter > maxit);
240 norm_of_residual=rnorm/Anorm*ynorm;
241 convergeFlag=(norm_of_residual<(*tolerance));
242
243
244 if (maxIterFlag) {
245 status = SOLVER_MAXITER_REACHED;
246 } else if (breakFlag) {
247 status = SOLVER_BREAKDOWN;
248 }
249 ++(num_iter);
250 }
251 /* end of iteration */
252
253 Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
254 TMPMEMFREE(w);
255 TMPMEMFREE(w1);
256 TMPMEMFREE(w2);
257 TMPMEMFREE(r1);
258 TMPMEMFREE(r2);
259 TMPMEMFREE(y);
260 TMPMEMFREE(v);
261
262 *iter=num_iter;
263 *tol=norm_of_residual;
264
265 /* End of MINRES */
266 return status;
267 }

  ViewVC Help
Powered by ViewVC 1.1.26