Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
minres.h
Go to the documentation of this file.
1# ifndef _RHEOLEF_MINRES_H
2# define _RHEOLEF_MINRES_H
3//
4// This file is part of Rheolef.
5//
6// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
7//
8// Rheolef is free software; you can redistribute it and/or modify
9// it under the terms of the GNU General Public License as published by
10// the Free Software Foundation; either version 2 of the License, or
11// (at your option) any later version.
12//
13// Rheolef is distributed in the hope that it will be useful,
14// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16// GNU General Public License for more details.
17//
18// You should have received a copy of the GNU General Public License
19// along with Rheolef; if not, write to the Free Software
20// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21//
22// =========================================================================
23// AUTHOR: Pierre.Saramito@imag.fr
24// DATE: 22 april 2009
25
26namespace rheolef {
92} // namespace rheolef
93
94#include "rheolef/solver_option.h"
95
96namespace rheolef {
97
98// [verbatim_minres_synopsis]
99template <class Matrix, class Vector, class Preconditioner>
100int minres (const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M,
101 const solver_option& sopt = solver_option())
102// [verbatim_minres_synopsis]
103// [verbatim_minres]
104{
105 // Size &max_iter, Real &tol, odiststream *p_derr = 0
106 typedef typename Vector::size_type Size;
107 typedef typename Vector::float_type Real;
108 std::string label = (sopt.label != "" ? sopt.label : "minres");
109 Vector b = M.solve(Mb);
110 Real norm_b = sqrt(fabs(dot(Mb,b)));
111 if (sopt.absolute_stopping || norm_b == Real(0.)) norm_b = 1;
112 Vector Mr = Mb - A*x;
113 Vector z = M.solve(Mr);
114 Real beta2 = dot(Mr, z);
115 Real norm_r = sqrt(fabs(beta2));
116 sopt.residue = norm_r/norm_b;
117 if (sopt.p_err) (*sopt.p_err) << "[" << label << "] #iteration residue" << std::endl
118 << "[" << label << "] 0 " << sopt.residue << std::endl;
119 if (beta2 < 0 || sopt.residue <= sopt.tol) {
120 return 0;
121 }
122 Real beta = sqrt(beta2);
123 Real eta = beta;
124 Vector Mv = Mr/beta;
125 Vector u = z/beta;
126 Real c_old = 1.;
127 Real s_old = 0.;
128 Real c = 1.;
129 Real s = 0.;
130 Vector u_old (x.ownership(), 0.);
131 Vector Mv_old (x.ownership(), 0.);
132 Vector w (x.ownership(), 0.);
133 Vector w_old (x.ownership(), 0.);
134 Vector w_old2 (x.ownership(), 0.);
135 for (sopt.n_iter = 1; sopt.n_iter <= sopt.max_iter; sopt.n_iter++) {
136 // Lanczos
137 Mr = A*u;
138 z = M.solve(Mr);
139 Real alpha = dot(Mr, u);
140 Mr = Mr - alpha*Mv - beta*Mv_old;
141 z = z - alpha*u - beta*u_old;
142 beta2 = dot(Mr, z);
143 if (beta2 < 0) {
144 dis_warning_macro ("minres: machine precision problem");
145 sopt.residue = norm_r/norm_b;
146 return 2;
147 }
148 Real beta_old = beta;
149 beta = sqrt(beta2);
150 // QR factorisation
151 Real c_old2 = c_old;
152 Real s_old2 = s_old;
153 c_old = c;
154 s_old = s;
155 Real rho0 = c_old*alpha - c_old2*s_old*beta_old;
156 Real rho2 = s_old*alpha + c_old2*c_old*beta_old;
157 Real rho1 = sqrt(sqr(rho0) + sqr(beta));
158 Real rho3 = s_old2 * beta_old;
159 // Givens rotation
160 c = rho0 / rho1;
161 s = beta / rho1;
162 // update
163 w_old2 = w_old;
164 w_old = w;
165 w = (u - rho2*w_old - rho3*w_old2)/rho1;
166 x += c*eta*w;
167 eta = -s*eta;
168 Mv_old = Mv;
169 u_old = u;
170 Mv = Mr/beta;
171 u = z/beta;
172 // check residue
173 norm_r *= s;
174 sopt.residue = norm_r/norm_b;
175 if (sopt.p_err) (*sopt.p_err) << "[" << label << "] " << sopt.n_iter << " " << sopt.residue << std::endl;
176 if (sopt.residue <= sopt.tol) return 0;
177 }
178 return 1;
179}
180// [verbatim_minres]
181
182}// namespace rheolef
183# endif // _RHEOLEF_MINRES_H
see the solver_option page for the full documentation
#define dis_warning_macro(message)
Definition dis_macros.h:66
This file is part of Rheolef.
int minres(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())
Definition minres.h:100
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
Definition eta.h:25
DATA::size_type size_type
Definition Vector.h:188
Definition leveque.h:25
Expr1::memory_type M