.\"
.de Id
..
.de Sp
.if n .sp
.if t .sp 0.4
..
.TH uzawa 4rheolef "rheolef-7.0" "rheolef-7.0" "rheolef-7.0"
.\" label: /*D:uzawa
.SH NAME
\fBuzawa\fP -- Uzawa algorithm.
.\" skip: @findex uzawa
.\" skip: @cindex Uzawa algorithm
.\" skip: @cindex iterative solver
.\" skip: @cindex preconditioner
.SH SYNOPSIS

.\" begin_example
.Sp
.nf
  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)
.Sp
.fi
.\" end_example
.PP
.SH EXAMPLE

The simplest call to 'uzawa' has the folling form:
.\" begin_example
.Sp
.nf
    solver_option sopt;
    sopt.max_iter = 100;
    sopt.tol = 1e-7;
    int status = uzawa(A, x, b, eye(), sopt);
.Sp
.fi
.\" end_example
.PP
.SH DESCRIPTION

\fBuzawa\fP 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:
.\" begin table
.\" start item
.TP
.B x
approximate solution to Ax = b

.\" start item
.TP
.B sopt.n_iter
the number of iterations performed before the tolerance was reached

.\" start item
.TP
.B sopt.residue
the residual after the final iteration
.\" end table
See also the solver_option(2).
.PP
.\" skip start:AUTHOR:

.\" skip start:DATE:

.\" skip start:METHODS:
.\" END
.SH IMPLEMENTATION
.\" begin_example
.Sp
.nf
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;
}
.Sp
.fi
.\" end_example

.\" LENGTH = 1
.SH SEE ALSO
solver_option(2)
.SH COPYRIGHT
Copyright (C) 2000-2018 Pierre Saramito <Pierre.Saramito@imag.fr> GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>.
This is free software: you are free to change and redistribute it.  There is NO WARRANTY, to the extent permitted by law.