/[escript]/branches/trilinos_from_5897/paso/src/MINRES.cpp
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6009 - (show 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
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
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

  ViewVC Help
Powered by ViewVC 1.1.26