# Contents of /branches/trilinos_from_5897/paso/src/MINRES.cpp

Revision 6009 - (show annotations)
Wed Mar 2 04:13:26 2016 UTC (3 years, 2 months ago) by caltinay
File size: 5838 byte(s)
```Much needed sync with trunk...

```
 1 2 /***************************************************************************** 3 * 4 * Copyright (c) 2003-2016 by The University of Queensland 5 6 * 7 * Primary Business: Queensland, Australia 8 * Licensed under the Open Software License version 3.0 9 10 * 11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 * Development 2012-2013 by School of Earth Sciences 13 * Development from 2014 by Centre for Geoscience Computing (GeoComp) 14 * 15 *****************************************************************************/ 16 17 18 /* MINRES iterations */ 19 20 #include "Solver.h" 21 #include "PasoUtil.h" 22 #include "SystemMatrix.h" 23 24 namespace paso { 25 26 /* 27 * 28 * Purpose 29 * ======= 30 * 31 * MINRES solves the linear system A*x = b 32 * 33 * Convergence test: norm( b - A*x )< TOL. 34 * 35 * Arguments 36 * ========= 37 * 38 * R (input) DOUBLE PRECISION array, dimension N. 39 * On entry, residual of initial guess x 40 * 41 * X (input/output) DOUBLE PRECISION array, dimension N. 42 * On input, the initial guess. 43 * 44 * ITER (input/output) INT 45 * On input, the maximum iterations to be performed. 46 * On output, actual number of iterations performed. 47 * 48 * INFO (output) INT 49 * 50 * = NoError: Successful exit. Iterated approximate solution returned. 51 * = MaxIterReached 52 * = InputError Illegal parameter: 53 * = MemoryError : 54 * = NegativeNormError : 55 * 56 * ============================================================== 57 */ 58 59 SolverResult Solver_MINRES(SystemMatrix_ptr A, double* R, double* X, 60 dim_t* iter, double* tolerance, Performance* pp) 61 { 62 const dim_t maxit = *iter; 63 if (maxit <= 0) { 64 return InputError; 65 } 66 67 double delta,gamma=0.,gamma_old=0.,eta=0.,dp0=0., c=1.0,c_old=1.0; 68 double c_ancient=1.,s=0.0,s_old=0.0,s_ancient, norm_of_residual=0., rnorm_prec=1; 69 double tol=1., norm_scal=1.; 70 double alpha_0,alpha_1,alpha_2,alpha_3,dp = 0.0; 71 dim_t num_iter = 0; 72 const dim_t n = A->getTotalNumRows(); 73 bool convergeFlag=false; 74 SolverResult status = NoError; 75 76 double* ZNEW = new double[n]; 77 double* Z = new double[n]; 78 double* AZ = new double[n]; 79 double* W = new double[n]; 80 double* R_old = new double[n]; 81 double* W_old = new double[n]; 82 double* R_ancient = new double[n]; 83 double* W_ancient = new double[n]; 84 85 // z <- Prec*r 86 A->solvePreconditioner(Z, R); 87 // gamma <- r'*z 88 dp = util::innerProduct(n, R ,Z,A->mpi_info); 89 dp0 = dp; 90 if (dp < 0) { 91 status = NegativeNormError; 92 } else if (std::abs(dp) <= 0) { 93 // happy break down 94 convergeFlag = true; 95 } else { 96 // gamma <- sqrt(r'*z) 97 gamma = sqrt(dp); 98 eta = gamma; 99 rnorm_prec = gamma; 100 norm_of_residual=util::l2(n, R, A->mpi_info); 101 norm_scal=rnorm_prec/norm_of_residual; 102 tol=(*tolerance)*norm_scal; 103 } 104 105 while (!convergeFlag && status == NoError) { 106 // z <- z / gamma 107 util::scale(n, Z, 1./gamma); 108 109 // Az <- A*z 110 A->MatrixVector_CSR_OFFSET0(PASO_ONE, Z, PASO_ZERO, AZ); 111 112 // delta <- Az'.z 113 delta = util::innerProduct(n, AZ, Z, A->mpi_info); 114 115 // r_new <- Az-delta/gamma * r - gamma/gamma_old r_old 116 if (num_iter>0) 117 util::copy(n, R_ancient, R_old); // r__ancient <- r_old 118 119 util::copy(n, R_old, R); // r_old <- r 120 util::copy(n, R, AZ); // r <- Az 121 util::AXPY(n, R, -delta/gamma, R_old); // r <- r - delta/gamma v 122 if (num_iter > 0) 123 util::AXPY(n, R, -gamma/gamma_old, R_ancient); // r <- r - gamma/gamma_old r__ancient 124 125 // z <- prec*r 126 A->solvePreconditioner(ZNEW, R); 127 128 dp = util::innerProduct(n, R, ZNEW, A->mpi_info); 129 if (dp < 0.) { 130 status = NegativeNormError; 131 } else if (std::abs(dp) == 0.) { 132 // happy break down 133 convergeFlag = true; 134 } else if (std::abs(dp) > 0.e-13 * std::abs(dp0)) { 135 // gamma <- sqrt(r'*z) 136 gamma_old = gamma; 137 gamma = sqrt(dp); 138 // QR factorisation 139 c_ancient = c_old; c_old = c; 140 s_ancient = s_old; s_old = s; 141 142 alpha_0 = c_old * delta - c_ancient * s_old * gamma_old; 143 alpha_1 = sqrt(alpha_0*alpha_0 + gamma*gamma); 144 alpha_2 = s_old * delta + c_ancient * c_old * gamma_old; 145 alpha_3 = s_ancient * gamma_old; 146 147 // Givens rotation 148 c = alpha_0 / alpha_1; 149 s = gamma / alpha_1; 150 151 rnorm_prec = rnorm_prec * s; 152 153 // w_new <- (z-alpha_3 w - alpha_2 w_old)/alpha_1 154 155 if (num_iter > 1) 156 util::copy(n, W_ancient, W_old); // w__ancient <- w_old 157 if (num_iter > 0) 158 util::copy(n, W_old, W); // w_old <- w 159 160 util::copy(n, W, Z); 161 if (num_iter > 1) 162 util::AXPY(n, W, -alpha_3, W_ancient); // w <- w - alpha_3 w__ancient 163 if (num_iter > 0) 164 util::AXPY(n, W, -alpha_2, W_old); // w <- w - alpha_2 w_old 165 util::scale(n, W, 1.0 / alpha_1); // w <- w / alpha_1 166 167 util::AXPY(n, X, c * eta, W); // x <- x + c eta w 168 eta = - s * eta; 169 convergeFlag = rnorm_prec <= tol; 170 } else { 171 status = Breakdown; 172 } 173 util::copy(n, Z, ZNEW); 174 ++num_iter; 175 if (!convergeFlag && num_iter >= maxit) 176 status = MaxIterReached; 177 } 178 delete[] Z; 179 delete[] ZNEW; 180 delete[] AZ; 181 delete[] R_old; 182 delete[] R_ancient; 183 delete[] W; 184 delete[] W_old; 185 delete[] W_ancient; 186 187 *iter=num_iter; 188 *tolerance=rnorm_prec/norm_scal; 189 return status; 190 } 191 192 } // namespace paso 193

## Properties

Name Value
svn:mergeinfo /branches/4.0fordebian/paso/src/MINRES.cpp:5567-5588 /branches/amg_from_3530/paso/src/MINRES.cpp:3531-3826 /branches/lapack2681/paso/src/MINRES.cpp:2682-2741 /branches/pasowrap/paso/src/MINRES.cpp:3661-3674 /branches/py3_attempt2/paso/src/MINRES.cpp:3871-3891 /branches/restext/paso/src/MINRES.cpp:2610-2624 /branches/ripleygmg_from_3668/paso/src/MINRES.cpp:3669-3791 /branches/stage3.0/paso/src/MINRES.cpp:2569-2590 /branches/symbolic_from_3470/paso/src/MINRES.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/paso/src/MINRES.cpp:3517-3974 /release/3.0/paso/src/MINRES.cpp:2591-2601 /release/4.0/paso/src/MINRES.cpp:5380-5406 /trunk/paso/src/MINRES.cpp:4257-4344,5898-6007 /trunk/ripley/test/python/paso/src/MINRES.cpp:3480-3515