Rheolef  7.2
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  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
point u(const point &x)
#define dis_warning_macro(message)
Definition: dis_macros.h:66
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
This file is part of Rheolef.
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
int minres(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())
Definition: minres.h:100
Float beta[][pmax+1]
Definition: eta.h:25
DATA::size_type size_type
Definition: Vector.h:188
Definition: leveque.h:25
Expr1::memory_type M
Definition: vec_expr_v2.h:416