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

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

```
 1 artak 1787 2 jfenwick 3981 /***************************************************************************** 3 ksteube 1811 * 4 jfenwick 5863 * Copyright (c) 2003-2016 by The University of Queensland 5 jfenwick 3981 6 ksteube 1811 * 7 * Primary Business: Queensland, Australia 8 * Licensed under the Open Software License version 3.0 9 10 * 11 jfenwick 3981 * Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 jfenwick 4657 * Development 2012-2013 by School of Earth Sciences 13 * Development from 2014 by Centre for Geoscience Computing (GeoComp) 14 jfenwick 3981 * 15 *****************************************************************************/ 16 artak 1787 17 ksteube 1811 18 artak 1979 /* MINRES iterations */ 19 artak 1787 20 #include "Solver.h" 21 #include "PasoUtil.h" 22 caltinay 6009 #include "SystemMatrix.h" 23 artak 1787 24 caltinay 4856 namespace paso { 25 artak 1787 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 gross 3140 * R (input) DOUBLE PRECISION array, dimension N. 39 caltinay 3642 * On entry, residual of initial guess x 40 artak 1787 * 41 gross 3140 * X (input/output) DOUBLE PRECISION array, dimension N. 42 artak 1787 * 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 caltinay 5963 * = NoError: Successful exit. Iterated approximate solution returned. 51 * = MaxIterReached 52 * = InputError Illegal parameter: 53 * = MemoryError : 54 * = NegativeNormError : 55 artak 1787 * 56 * ============================================================== 57 */ 58 59 caltinay 5963 SolverResult Solver_MINRES(SystemMatrix_ptr A, double* R, double* X, 60 dim_t* iter, double* tolerance, Performance* pp) 61 gross 3140 { 62 caltinay 4856 const dim_t maxit = *iter; 63 if (maxit <= 0) { 64 caltinay 5963 return InputError; 65 caltinay 4856 } 66 caltinay 4873 67 caltinay 4856 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 caltinay 4836 double tol=1., norm_scal=1.; 70 caltinay 4856 double alpha_0,alpha_1,alpha_2,alpha_3,dp = 0.0; 71 dim_t num_iter = 0; 72 caltinay 4836 const dim_t n = A->getTotalNumRows(); 73 caltinay 4856 bool convergeFlag=false; 74 caltinay 5963 SolverResult status = NoError; 75 artak 1787 76 caltinay 4856 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 artak 1787 85 caltinay 4856 // 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 caltinay 5963 status = NegativeNormError; 92 caltinay 6009 } else if (std::abs(dp) <= 0) { 93 caltinay 4856 // 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 artak 1787 105 caltinay 5963 while (!convergeFlag && status == NoError) { 106 caltinay 4856 // z <- z / gamma 107 util::scale(n, Z, 1./gamma); 108 artak 1787 109 caltinay 4856 // Az <- A*z 110 caltinay 5933 A->MatrixVector_CSR_OFFSET0(PASO_ONE, Z, PASO_ZERO, AZ); 111 artak 1787 112 caltinay 4856 // delta <- Az'.z 113 caltinay 4873 delta = util::innerProduct(n, AZ, Z, A->mpi_info); 114 artak 1787 115 caltinay 4856 // 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 artak 1787 119 caltinay 4856 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 caltinay 4873 A->solvePreconditioner(ZNEW, R); 127 caltinay 4856 128 dp = util::innerProduct(n, R, ZNEW, A->mpi_info); 129 if (dp < 0.) { 130 caltinay 5963 status = NegativeNormError; 131 caltinay 6009 } else if (std::abs(dp) == 0.) { 132 caltinay 4856 // happy break down 133 convergeFlag = true; 134 caltinay 6009 } else if (std::abs(dp) > 0.e-13 * std::abs(dp0)) { 135 caltinay 4856 // gamma <- sqrt(r'*z) 136 gamma_old = gamma; 137 gamma = sqrt(dp); 138 // QR factorisation 139 caltinay 4873 c_ancient = c_old; c_old = c; 140 caltinay 4856 s_ancient = s_old; s_old = s; 141 caltinay 4873 142 caltinay 4856 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 caltinay 4873 147 caltinay 4856 // 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 caltinay 4873 160 caltinay 4856 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 caltinay 5963 status = Breakdown; 172 caltinay 4856 } 173 caltinay 4873 util::copy(n, Z, ZNEW); 174 caltinay 4856 ++num_iter; 175 if (!convergeFlag && num_iter >= maxit) 176 caltinay 5963 status = MaxIterReached; 177 caltinay 4856 } 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 caltinay 4873 187 caltinay 4856 *iter=num_iter; 188 *tolerance=rnorm_prec/norm_scal; 189 return status; 190 artak 1787 } 191 caltinay 3642 192 caltinay 4856 } // 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