Rheolef  7.1
an efficient C++ finite element environment
cg.h
Go to the documentation of this file.
1 # ifndef _RHEOLEF_CG_H
2 # define _RHEOLEF_CG_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
24 // DATE: 20 april 2009
25 
26 namespace rheolef {
86 } // namespace rheolef
87 
88 #include "rheolef/solver_option.h"
89 
90 namespace rheolef {
91 
92 // [verbatim_cg_synopsis]
93 template <class Matrix, class Vector, class Vector2, class Preconditioner>
94 int cg (const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M,
95  const solver_option& sopt = solver_option())
96 // [verbatim_cg_synopsis]
97 // [verbatim_cg]
98 {
99  typedef typename Vector::size_type Size;
100  typedef typename Vector::float_type Real;
101  std::string label = (sopt.label != "" ? sopt.label : "cg");
102  Vector b = M.solve(Mb);
103  Real norm2_b = dot(Mb,b);
104  if (sopt.absolute_stopping || norm2_b == Real(0)) norm2_b = 1;
105  Vector Mr = Mb - A*x;
106  Real norm2_r = 0;
107  if (sopt.p_err) (*sopt.p_err) << "[" << label << "] #iteration residue" << std::endl;
108  Vector p;
109  for (sopt.n_iter = 0; sopt.n_iter <= sopt.max_iter; sopt.n_iter++) {
110  Vector r = M.solve(Mr);
111  Real prev_norm2_r = norm2_r;
112  norm2_r = dot(Mr, r);
113  sopt.residue = sqrt(norm2_r/norm2_b);
114  if (sopt.p_err) (*sopt.p_err) << "[" << label << "] " << sopt.n_iter << " " << sopt.residue << std::endl;
115  if (sopt.residue <= sopt.tol) return 0;
116  if (sopt.n_iter == 0) {
117  p = r;
118  } else {
119  Real beta = norm2_r/prev_norm2_r;
120  p = r + beta*p;
121  }
122  Vector Mq = A*p;
123  Real alpha = norm2_r/dot(Mq, p);
124  x += alpha*p;
125  Mr -= alpha*Mq;
126  }
127  return 1;
128 }
129 // [verbatim_cg]
130 
131 }// namespace rheolef
132 # endif // _RHEOLEF_CG_H
rheolef::cg
int cg(const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())
Definition: cg.h:94
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
p
Definition: sphere.icc:25
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
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