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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3642 - (show annotations)
Thu Oct 27 03:41:51 2011 UTC (7 years, 5 months ago) by caltinay
File MIME type: text/plain
File size: 6864 byte(s)
Assorted spelling/comment fixes in paso.

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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 /* MINRES 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 ESYS_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 initial 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
65
66 err_t Paso_Solver_MINRES(
67 Paso_SystemMatrix * A,
68 double * R,
69 double * X,
70 dim_t *iter,
71 double * tolerance,
72 Paso_Performance* pp)
73 {
74
75 double delta,gamma=0.,gamma_old=0.,eta=0.,dp0=0., c=1.0,c_old=1.0,c_ancient=1.,s=0.0,s_old=0.0,s_ancient, norm_of_residual=0., rnorm_prec=1;
76 double tol=1., norm_scal=1.;
77 const dim_t maxit = *iter;
78 double alpha_0,alpha_1,alpha_2,alpha_3,dp = 0.0;
79 dim_t num_iter = 0;
80 double *Z=NULL, *W=NULL, *AZ=NULL, *R_old=NULL, *R_ancient=NULL, *W_old=NULL, *W_ancient=NULL, *ZNEW=NULL;
81 const dim_t n = Paso_SystemMatrix_getTotalNumRows(A);
82 bool_t convergeFlag=FALSE;
83 err_t status = SOLVER_NO_ERROR;
84 /*
85 *
86 * Start of Calculation :
87 * ---------------------
88 *
89 * */
90
91
92 /* Test the input parameters. */
93 if (n < 0 || maxit<=0 ) {
94 status=SOLVER_INPUT_ERROR;
95 }
96
97 ZNEW = TMPMEMALLOC(n,double);
98 Z = TMPMEMALLOC(n,double);
99 AZ = TMPMEMALLOC(n,double);
100 W = TMPMEMALLOC(n,double);
101 R_old = TMPMEMALLOC(n,double);
102 W_old = TMPMEMALLOC(n,double);
103 R_ancient = TMPMEMALLOC(n,double);
104 W_ancient = TMPMEMALLOC(n,double);
105
106 if (R_ancient==NULL || Z==NULL || W==NULL || AZ==NULL || R_old==NULL || W_old==NULL || W_ancient==NULL || ZNEW==NULL) {
107 status=SOLVER_MEMORY_ERROR;
108 }
109
110 if (status ==SOLVER_NO_ERROR) {
111
112 Paso_SystemMatrix_solvePreconditioner(A, Z, R); /* z <- Prec*r */
113 /* gamma <- r'*z */
114 dp=Paso_InnerProduct(n, R ,Z,A->mpi_info); /* gamma <- r'*z */
115 dp0=dp;
116 if (dp<0) {
117 status=SOLVER_NEGATIVE_NORM_ERROR;
118 } else if (! ABS(dp)>0) {
119 convergeFlag = TRUE; /* happy break down */
120 } else {
121 gamma = sqrt(dp); /* gamma <- sqrt(r'*z) */
122 eta = gamma;
123 rnorm_prec = gamma;
124 norm_of_residual=Paso_l2(n, R, A->mpi_info);
125 norm_scal=rnorm_prec/norm_of_residual;
126 tol=(*tolerance)*norm_scal;
127 }
128 }
129 while (!(convergeFlag || (status !=SOLVER_NO_ERROR) ))
130 {
131 /* z <- z / gamma */
132 Paso_Scale(n, Z,1./gamma);
133
134 /* Az <- A*z */
135 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(PASO_ONE, A, Z,PASO_ZERO,AZ);
136
137 /* delta <- Az'.z */
138 delta=Paso_InnerProduct(n,AZ,Z,A->mpi_info);
139
140 /* r_new <- Az-delta/gamma * r - gamma/gamma_old r_old */
141 if (num_iter>0) Paso_Copy(n, R_ancient, R_old); /* r__ancient <- r_old */
142 Paso_Copy(n, R_old, R); /* r_old <- r */
143 Paso_Copy(n, R, AZ); /* r <- Az */
144 Paso_AXPY(n, R, -delta/gamma, R_old); /* r <- r - delta/gamma v */
145 if (num_iter>0) Paso_AXPY(n, R, -gamma/gamma_old, R_ancient); /* r <- r - gamma/gamma_old r__ancient */
146
147 /* z <- prec*r */
148 Paso_SystemMatrix_solvePreconditioner(A, ZNEW, R);
149
150
151 dp=Paso_InnerProduct(n,R,ZNEW,A->mpi_info);
152 if (dp < 0.) {
153 status=SOLVER_NEGATIVE_NORM_ERROR;
154 } else if (ABS(dp) == 0.) {
155 convergeFlag = TRUE; /* happy break down */
156 } else if (ABS(dp) > 0.e-13 * ABS(dp0) ) {
157 /* gamma <- sqrt(r'*z) */
158 gamma_old=gamma;
159 gamma = sqrt(dp);
160 /* QR factorisation */
161
162 c_ancient = c_old; c_old = c;
163 s_ancient = s_old; s_old = s;
164
165 alpha_0 = c_old * delta - c_ancient * s_old * gamma_old;
166 alpha_1 = sqrt(alpha_0*alpha_0 + gamma*gamma);
167 alpha_2 = s_old * delta + c_ancient * c_old * gamma_old;
168 alpha_3 = s_ancient * gamma_old;
169
170 /* Givens rotation */
171
172 c = alpha_0 / alpha_1;
173 s = gamma / alpha_1;
174
175 rnorm_prec = rnorm_prec * s;
176
177 /* w_new <- (z-alpha_3 w - alpha_2 w_old)/alpha_1 */
178
179 if (num_iter>1) Paso_Copy(n, W_ancient, W_old); /* w__ancient <- w_old */
180 if (num_iter>0) Paso_Copy(n, W_old, W); /* w_old <- w */
181
182 Paso_Copy(n, W, Z);
183 if (num_iter>1) Paso_AXPY(n, W,- alpha_3,W_ancient); /* w <- w - alpha_3 w__ancient */
184 if (num_iter>0) Paso_AXPY(n, W,- alpha_2,W_old); /* w <- w - alpha_2 w_old */
185 Paso_Scale(n, W, 1.0 / alpha_1); /* w <- w / alpha_1 */
186 /* */
187 Paso_AXPY(n, X,c * eta,W); /* x <- x + c eta w */
188 eta = - s * eta;
189 convergeFlag = rnorm_prec <= tol;
190 } else {
191 status=SOLVER_BREAKDOWN;
192 }
193 Paso_Copy(n, Z, ZNEW);
194 ++(num_iter);
195 if ( !convergeFlag && (num_iter>=maxit)) status = SOLVER_MAXITER_REACHED;
196 }
197 TMPMEMFREE(Z);
198 TMPMEMFREE(ZNEW);
199 TMPMEMFREE(AZ);
200 TMPMEMFREE(R_old);
201 TMPMEMFREE(R_ancient);
202 TMPMEMFREE(W);
203 TMPMEMFREE(W_old);
204 TMPMEMFREE(W_ancient);
205
206 *iter=num_iter;
207 *tolerance=rnorm_prec/norm_scal;
208
209 /* End of MINRES */
210 return status;
211 }
212

  ViewVC Help
Powered by ViewVC 1.1.26