Rheolef  7.2
an efficient C++ finite element environment
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 
25 namespace rheolef {
26 /*D:uzawa
27 NAME: @code{uzawa} -- Uzawa algorithm.
28 @findex uzawa
29 @cindex Uzawa algorithm
30 @cindex iterative solver
31 @cindex preconditioner
32 SYNOPSIS:
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 
39 EXAMPLE:
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 
49 DESCRIPTION:
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 
74 AUTHOR:
75  Pierre Saramito
76  | Pierre.Saramito@imag.fr
77  LJK-IMAG, 38041 Grenoble cedex 9, France
78 DATE:
79  20 april 2009
80 METHODS: @uzawa
81 End:
82 */
83 //<uzawa:
84 template<class Matrix, class Vector, class Preconditioner, class Real2>
85 int 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
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
DATA::size_type size_type
Definition: Vector.h:188
Expr1::memory_type M
Definition: vec_expr_v2.h:416