.\"
.de Id
..
.de Sp
.if n .sp
.if t .sp 0.4
..
.TH puzawa 4rheolef "rheolef-6.7" "rheolef-6.7" "rheolef-6.7"
.\" label: /*D:puzawa
.SH NAME
\fBpuzawa\fP -- Uzawa algorithm.
.\" skip: @findex puzawa
.\" 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 puzawa (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
      int &max_iter, Real &tol, const Real& rho, odiststream *p_derr=0);
.Sp
.fi
.\" end_example
.PP
.SH EXAMPLE

The simplest call to 'puzawa' has the folling form:
.\" begin_example
.Sp
.nf
    size_t max_iter = 100;
    double tol = 1e-7;
    int status = puzawa(A, x, b, EYE, max_iter, tol, 1.0, &derr);
.Sp
.fi
.\" end_example
.PP
.SH DESCRIPTION

\fBpuzawa\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 'pcg' that
computes automatically the optimal descdnt step length at
each iteration.
.PP
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 max_iter
the number of iterations performed before the tolerance was reached

.\" start item
.TP
.B tol
the residual after the final iteration
.\" end table
.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 Real, class Size>
int puzawa(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M,
    Size &max_iter, Real &tol, const Real& rho,
    odiststream *p_derr, std::string label)
{
    Vector b = M.solve(Mb);
    Real norm2_b = dot(Mb,b);
    Real norm2_r = norm2_b;
    if (norm2_b == Real(0)) norm2_b = 1;
    if (p_derr) (*p_derr) << "[" << label << "] #iteration residue" << std::endl;
    for (Size n = 0; n <= max_iter; n++) {
        Vector Mr = A*x - Mb;
        Vector r = M.solve(Mr);
        norm2_r = dot(Mr, r);
        if (p_derr) (*p_derr) << "[" << label << "] " << n << " " << sqrt(norm2_r/norm2_b) << std::endl;
        if (norm2_r <= sqr(tol)*norm2_b) {
          tol = sqrt(norm2_r/norm2_b);
          max_iter = n;
          return 0;
        }
        x  -= rho*r;
    }
    tol = sqrt(norm2_r/norm2_b);
    return 1;
}
.Sp
.fi
.\" end_example