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

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 24 #endif 25 26 #ifdef ESYS_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 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