Rheolef  7.1
an efficient C++ finite element environment
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 
26 namespace rheolef {
92 } // namespace rheolef
93 
94 #include "rheolef/solver_option.h"
95 
96 namespace rheolef {
97 
98 // [verbatim_minres_synopsis]
99 template <class Matrix, class Vector, class Preconditioner>
100 int 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  dis_warning_macro ("beta2 = " << beta2 << " < 0: stop");
121  return 0;
122  }
123  Real beta = sqrt(beta2);
124  Real eta = beta;
125  Vector Mv = Mr/beta;
126  Vector u = z/beta;
127  Real c_old = 1.;
128  Real s_old = 0.;
129  Real c = 1.;
130  Real s = 0.;
131  Vector u_old (x.ownership(), 0.);
132  Vector Mv_old (x.ownership(), 0.);
133  Vector w (x.ownership(), 0.);
134  Vector w_old (x.ownership(), 0.);
135  Vector w_old2 (x.ownership(), 0.);
136  for (sopt.n_iter = 1; sopt.n_iter <= sopt.max_iter; sopt.n_iter++) {
137  // Lanczos
138  Mr = A*u;
139  z = M.solve(Mr);
140  Real alpha = dot(Mr, u);
141  Mr = Mr - alpha*Mv - beta*Mv_old;
142  z = z - alpha*u - beta*u_old;
143  beta2 = dot(Mr, z);
144  if (beta2 < 0) {
145  dis_warning_macro ("minres: machine precision problem");
146  sopt.residue = norm_r/norm_b;
147  return 2;
148  }
149  Real beta_old = beta;
150  beta = sqrt(beta2);
151  // QR factorisation
152  Real c_old2 = c_old;
153  Real s_old2 = s_old;
154  c_old = c;
155  s_old = s;
156  Real rho0 = c_old*alpha - c_old2*s_old*beta_old;
157  Real rho2 = s_old*alpha + c_old2*c_old*beta_old;
158  Real rho1 = sqrt(sqr(rho0) + sqr(beta));
159  Real rho3 = s_old2 * beta_old;
160  // Givens rotation
161  c = rho0 / rho1;
162  s = beta / rho1;
163  // update
164  w_old2 = w_old;
165  w_old = w;
166  w = (u - rho2*w_old - rho3*w_old2)/rho1;
167  x += c*eta*w;
168  eta = -s*eta;
169  Mv_old = Mv;
170  u_old = u;
171  Mv = Mr/beta;
172  u = z/beta;
173  // check residue
174  norm_r *= s;
175  sopt.residue = norm_r/norm_b;
176  if (sopt.p_err) (*sopt.p_err) << "[" << label << "] " << sopt.n_iter << " " << sopt.residue << std::endl;
177  if (sopt.residue <= sopt.tol) return 0;
178  }
179  return 1;
180 }
181 // [verbatim_minres]
182 
183 }// namespace rheolef
184 # endif // _RHEOLEF_MINRES_H
rheolef::Vector
Definition: Vector.h:175
bdf::alpha
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
rheolef::dot
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
Definition: vec_expr_v2.h:415
rheolef::minres
int minres(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())
Definition: minres.h:100
eta
Definition: eta.h:25
mkgeo_ball.c
c
Definition: mkgeo_ball.sh:153
dis_warning_macro
#define dis_warning_macro(message)
Definition: dis_macros.h:66
rheolef::Vector::size_type
DATA::size_type size_type
Definition: Vector.h:188
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
u
Definition: leveque.h:25
u
Float u(const point &x)
Definition: transmission_error.cc:26
mkgeo_ball.b
b
Definition: mkgeo_ball.sh:152
M
Expr1::memory_type M
Definition: vec_expr_v2.h:385
rk::beta
Float beta[][pmax+1]
Definition: runge_kutta_semiimplicit.icc:60
rheolef::solver_option
see the solver_option page for the full documentation
Definition: solver_option.h:155