uzawa.h
Go to the documentation of this file.
# ifndef _RHEOLEF_UZAWA_H
# define _RHEOLEF_UZAWA_H
//
// This file is part of Rheolef.
//
// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
//
// Rheolef is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// Rheolef is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Rheolef; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// =========================================================================
#include "rheolef/solver_option.h"
namespace rheolef {
/*D:uzawa
NAME: @code{uzawa} -- Uzawa algorithm.
@findex uzawa
@cindex Uzawa algorithm
@cindex iterative solver
@cindex preconditioner
SYNOPSIS:
@example
template <class Matrix, class Vector, class Preconditioner, class Real>
int uzawa (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
const solver_option& sopt)
@end example
EXAMPLE:
@noindent
The simplest call to 'uzawa' has the folling form:
@example
solver_option sopt;
sopt.max_iter = 100;
sopt.tol = 1e-7;
int status = uzawa(A, x, b, eye(), sopt);
@end example
DESCRIPTION:
@noindent
@code{uzawa} solves the linear
system A*x=b using the Uzawa method. The Uzawa method is a
descent method in the direction opposite to the gradient,
with a constant step length 'rho'. The convergence is assured
when the step length 'rho' is small enough.
If matrix A is symmetric positive definite, please uses 'cg' that
computes automatically the optimal descdnt step length at
each iteration.
The return value indicates convergence within max_iter (input)
iterations (0), or no convergence within max_iter iterations (1).
Upon successful return, output arguments have the following values:
@table @code
@item x
approximate solution to Ax = b
@item sopt.n_iter
the number of iterations performed before the tolerance was reached
@item sopt.residue
the residual after the final iteration
@end table
See also the @ref{solver_option class}.
AUTHOR:
Pierre Saramito
| Pierre.Saramito@imag.fr
LJK-IMAG, 38041 Grenoble cedex 9, France
DATE:
20 april 2009
METHODS: @uzawa
End:
*/
//<uzawa:
template<class Matrix, class Vector, class Preconditioner, class Real2>
int uzawa (const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M,
const Real2& rho, const solver_option& sopt = solver_option())
{
typedef typename Vector::size_type Size;
typedef typename Vector::float_type Real;
std::string label = (sopt.label != "" ? sopt.label : "uzawa");
Vector b = M.solve(Mb);
Real norm2_b = dot(Mb,b);
Real norm2_r = norm2_b;
if (sopt.absolute_stopping || norm2_b == Real(0)) norm2_b = 1;
if (sopt.p_err) (*sopt.p_err) << "[" << label << "] #iteration residue" << std::endl;
for (sopt.n_iter = 0; sopt.n_iter <= sopt.max_iter; sopt.n_iter++) {
Vector Mr = A*x - Mb;
Vector r = M.solve(Mr);
norm2_r = dot(Mr, r);
sopt.residue = sqrt(norm2_r/norm2_b);
if (sopt.p_err) (*sopt.p_err) << "[" << label << "] " << sopt.n_iter << " " << sopt.residue << std::endl;
if (sopt.residue <= sopt.tol) return 0;
x -= rho*r;
}
return 1;
}
//>uzawa:
}// namespace rheolef
# endif // _RHEOLEF_UZAWA_H
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::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