NAME

uzawa -- Uzawa algorithm.

SYNOPSIS

  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)

EXAMPLE

The simplest call to 'uzawa' has the folling form:

    solver_option sopt;


    sopt.max_iter = 100;


    sopt.tol = 1e-7;


    int status = uzawa(A, x, b, eye(), sopt);

DESCRIPTION

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:

x

approximate solution to Ax = b

sopt.n_iter

the number of iterations performed before the tolerance was reached

sopt.residue

the residual after the final iteration See also the solver_option(2).

IMPLEMENTATION

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;
}

SEE ALSO

solver_option(2)

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.