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

Revision 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 6 months ago) by ksteube
File MIME type: text/plain
File size: 7110 byte(s)
```Copyright updated in all files

```
 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 *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; 241 convergeFlag=(norm_of_residual