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

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

Parent Directory Parent Directory | Revision Log Revision Log


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 * http://www.uq.edu.au
6 ksteube 1811 *
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 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

  ViewVC Help
Powered by ViewVC 1.1.26