Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
uzawa.h
Go to the documentation of this file.
1# ifndef _RHEOLEF_UZAWA_H
2# define _RHEOLEF_UZAWA_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#include "rheolef/solver_option.h"
24
25namespace rheolef {
26/*D:uzawa
27NAME: @code{uzawa} -- Uzawa algorithm.
28@findex uzawa
29@cindex Uzawa algorithm
30@cindex iterative solver
31@cindex preconditioner
32SYNOPSIS:
33@example
34 template <class Matrix, class Vector, class Preconditioner, class Real>
35 int uzawa (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
36 const solver_option& sopt)
37@end example
38
39EXAMPLE:
40 @noindent
41 The simplest call to 'uzawa' has the folling form:
42 @example
43 solver_option sopt;
44 sopt.max_iter = 100;
45 sopt.tol = 1e-7;
46 int status = uzawa(A, x, b, eye(), sopt);
47 @end example
48
49DESCRIPTION:
50 @noindent
51 @code{uzawa} solves the linear
52 system A*x=b using the Uzawa method. The Uzawa method is a
53 descent method in the direction opposite to the gradient,
54 with a constant step length 'rho'. The convergence is assured
55 when the step length 'rho' is small enough.
56 If matrix A is symmetric positive definite, please uses 'cg' that
57 computes automatically the optimal descdnt step length at
58 each iteration.
59 The return value indicates convergence within max_iter (input)
60 iterations (0), or no convergence within max_iter iterations (1).
61 Upon successful return, output arguments have the following values:
62 @table @code
63 @item x
64 approximate solution to Ax = b
65
66 @item sopt.n_iter
67 the number of iterations performed before the tolerance was reached
68
69 @item sopt.residue
70 the residual after the final iteration
71 @end table
72 See also the @ref{solver_option class}.
73
74AUTHOR:
75 Pierre Saramito
76 | Pierre.Saramito@imag.fr
77 LJK-IMAG, 38041 Grenoble cedex 9, France
78DATE:
79 20 april 2009
80METHODS: @uzawa
81End:
82*/
83//<uzawa:
84template<class Matrix, class Vector, class Preconditioner, class Real2>
85int uzawa (const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M,
86 const Real2& rho, const solver_option& sopt = solver_option())
87{
88 typedef typename Vector::size_type Size;
89 typedef typename Vector::float_type Real;
90 std::string label = (sopt.label != "" ? sopt.label : "uzawa");
91 Vector b = M.solve(Mb);
92 Real norm2_b = dot(Mb,b);
93 Real norm2_r = norm2_b;
94 if (sopt.absolute_stopping || norm2_b == Real(0)) norm2_b = 1;
95 if (sopt.p_err) (*sopt.p_err) << "[" << label << "] #iteration residue" << std::endl;
96 for (sopt.n_iter = 0; sopt.n_iter <= sopt.max_iter; sopt.n_iter++) {
97 Vector Mr = A*x - Mb;
98 Vector r = M.solve(Mr);
99 norm2_r = dot(Mr, r);
100 sopt.residue = sqrt(norm2_r/norm2_b);
101 if (sopt.p_err) (*sopt.p_err) << "[" << label << "] " << sopt.n_iter << " " << sopt.residue << std::endl;
102 if (sopt.residue <= sopt.tol) return 0;
103 x -= rho*r;
104 }
105 return 1;
106}
107//>uzawa:
108}// namespace rheolef
109# endif // _RHEOLEF_UZAWA_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 uzawa(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M, const Real2 &rho, const solver_option &sopt=solver_option())
Definition uzawa.h:85
DATA::size_type size_type
Definition Vector.h:188
Expr1::memory_type M