Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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
26namespace rheolef {
86} // namespace rheolef
87
88#include "rheolef/solver_option.h"
89
90namespace rheolef {
91
92// [verbatim_cg_synopsis]
93template <class Matrix, class Vector, class Vector2, class Preconditioner>
94int 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
see the solver_option page for the full documentation
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)
int cg(const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())
Definition cg.h:94
Definition sphere.icc:25
DATA::size_type size_type
Definition Vector.h:188
Expr1::memory_type M