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

Revision 1974 - (show annotations)
Thu Nov 6 02:40:10 2008 UTC (10 years, 5 months ago) by jfenwick
File MIME type: text/plain
File size: 6989 byte(s)
```Changing ZERO and ONE in Solver.h to PASO_ZERO and PASO_ONE.
Changed three static vars to #defines
```
 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 24 #endif 25 26 #ifdef PASO_MPI 27 #include 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 * tolerance, 76 Paso_Performance* pp) { 77 78 /* Local variables */ 79 80 dim_t num_iter=0,maxit; 81 bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE; 82 err_t status = SOLVER_NO_ERROR; 83 dim_t n = Paso_SystemMatrix_getTotalNumRows(A); 84 double *w=NULL, *w1=NULL, *w2=NULL, *r1=NULL, *r2=NULL, *y=NULL, *v=NULL; 85 86 double Anorm,ynorm,oldb,dbar,epsln,phibar,rhs1,rhs2,rnorm,tnorm2,ynorm2,cs,sn,eps,s,alfa,denom,z,beta1,beta; 87 double gmax,gmin,oldeps,delta,gbar,gamma,phi; 88 89 double norm_of_residual; 90 91 /* */ 92 /*-----------------------------------------------------------------*/ 93 /* */ 94 /* Start of Calculation : */ 95 /* --------------------- */ 96 /* */ 97 /* */ 98 w=TMPMEMALLOC(n,double); 99 w1=TMPMEMALLOC(n,double); 100 w2=TMPMEMALLOC(n,double); 101 r1=TMPMEMALLOC(n,double); 102 r2=TMPMEMALLOC(n,double); 103 y=TMPMEMALLOC(n,double); 104 v=TMPMEMALLOC(n,double); 105 106 107 if (w ==NULL || w1== NULL || w2== NULL || r1 == NULL || r2== NULL || y==NULL || v==NULL ) { 108 status=SOLVER_MEMORY_ERROR; 109 } 110 111 maxit = *iter; 112 113 /* Test the input parameters. */ 114 if (n < 0 || maxit<=0 ) { 115 status=SOLVER_INPUT_ERROR; 116 } 117 118 Paso_Copy(n,r1,r); 119 120 Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER); 121 Paso_Solver_solvePreconditioner(A,y,r1); 122 Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER); 123 124 beta1=Paso_InnerProduct(n,r1,y,A->mpi_info); 125 if (beta1<0) { 126 status=SOLVER_NEGATIVE_NORM_ERROR; 127 } 128 129 if (beta1>0) { 130 beta1=sqrt(beta1); 131 } 132 133 Performance_startMonitor(pp,PERFORMANCE_SOLVER); 134 135 Paso_zeroes(n,w); 136 Paso_zeroes(n,w2); 137 138 Paso_Copy(n,r2,r1); 139 140 Anorm = 0; 141 ynorm = 0; 142 oldb = 0; 143 beta = beta1; 144 dbar = 0; 145 epsln = 0; 146 phibar = beta1; 147 rhs1 = beta1; 148 rhs2 = 0; 149 rnorm = phibar; 150 tnorm2 = 0; 151 ynorm2 = 0; 152 cs = -1; 153 sn = 0; 154 eps = 0.0001; 155 156 while (!(convergeFlag || maxIterFlag || breakFlag || (status !=SOLVER_NO_ERROR) )) 157 { 158 159 s=1/beta; 160 Paso_Update(n, 0., v, s, y); 161 162 Performance_stopMonitor(pp,PERFORMANCE_SOLVER); 163 Performance_startMonitor(pp,PERFORMANCE_MVM); 164 Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(PASO_ONE, A, v,PASO_ZERO,y); 165 Performance_stopMonitor(pp,PERFORMANCE_MVM); 166 Performance_startMonitor(pp,PERFORMANCE_SOLVER); 167 168 if (num_iter >= 1) { 169 Paso_Update(n, 1., y, -(beta/oldb), r1); 170 } 171 172 alfa = Paso_InnerProduct(n,v,y,A->mpi_info); 173 Paso_Update(n, 1., y, (-alfa/beta), r2); 174 Paso_Copy(n,r1,r2); 175 Paso_Copy(n,r2,y); 176 177 Performance_stopMonitor(pp,PERFORMANCE_SOLVER); 178 Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER); 179 Paso_Solver_solvePreconditioner(A,y,r2); 180 Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER); 181 Performance_startMonitor(pp,PERFORMANCE_SOLVER); 182 183 oldb = beta; 184 beta = Paso_InnerProduct(n,y,r2,A->mpi_info); 185 if (beta<0) { 186 status=SOLVER_NEGATIVE_NORM_ERROR; 187 } 188 189 beta = sqrt( beta ); 190 tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta; 191 192 if (num_iter==0) { 193 gmax = ABS(alfa); 194 gmin = gmax; 195 } 196 197 /* Apply previous rotation Qk-1 to get 198 [deltak epslnk+1] = [cs sn][dbark 0 ] 199 [gbar k dbar k+1] [sn -cs][alfak betak+1]. */ 200 201 oldeps = epsln; 202 delta = cs * dbar + sn * alfa ; 203 gbar = sn * dbar - cs * alfa ; 204 epsln = sn * beta ; 205 dbar = - cs * beta; 206 207 gamma = sqrt(gbar*gbar+beta*beta) ; 208 gamma = MAX(gamma,eps) ; 209 cs = gbar / gamma ; 210 sn = beta / gamma ; 211 phi = cs * phibar ; 212 phibar = sn * phibar ; 213 214 /* Update x. */ 215 216 denom = 1/gamma ; 217 Paso_Copy(n,w1,w2); 218 Paso_Copy(n,w2,w); 219 220 Paso_LinearCombination(n,w,denom,v,-(denom*oldeps),w1); 221 Paso_Update(n, 1., w, -(delta*denom), w2) ; 222 Paso_Update(n, 1., x, phi, w) ; 223 224 /* Go round again. */ 225 226 gmax = MAX(gmax,gamma); 227 gmin = MIN(gmin,gamma); 228 z = rhs1 / gamma; 229 ynorm2 = z*z + ynorm2; 230 rhs1 = rhs2 - delta*z; 231 rhs2 = - epsln*z; 232 233 Anorm = sqrt( tnorm2 ) ; 234 ynorm = sqrt( ynorm2 ) ; 235 236 rnorm = phibar; 237 238 maxIterFlag = (num_iter > maxit); 239 norm_of_residual=rnorm/Anorm*ynorm; 240 convergeFlag=(norm_of_residual<(*tolerance)); 241 242 243 if (maxIterFlag) { 244 status = SOLVER_MAXITER_REACHED; 245 } else if (breakFlag) { 246 status = SOLVER_BREAKDOWN; 247 } 248 ++(num_iter); 249 } 250 /* end of iteration */ 251 252 Performance_stopMonitor(pp,PERFORMANCE_SOLVER); 253 TMPMEMFREE(w); 254 TMPMEMFREE(w1); 255 TMPMEMFREE(w2); 256 TMPMEMFREE(r1); 257 TMPMEMFREE(r2); 258 TMPMEMFREE(y); 259 TMPMEMFREE(v); 260 261 *iter=num_iter; 262 *tolerance=norm_of_residual; 263 264 /* End of MINRES */ 265 return status; 266 }