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.