|
|
Vatata Rm for VLC 1.0beta (Default branch)
|
|
puzawa(4rheolef) rheolef-6.1 puzawa(4rheolef) NAME
puzawa -- Uzawa algorithm. SYNOPSIS
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); EXAMPLE
The simplest call to 'puzawa' has the folling form: size_t max_iter = 100; double tol = 1e-7; int status = puzawa(A, x, b, EYE, max_iter, tol, 1.0, &derr); DESCRIPTION
puzawa solves the linear system A*x=b using the Uzawa method. The Uzawa method is a descent method in the direction opposite to the gradi- ent, 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. The return value indicates convergence within max_iter (input) iterations(0), or no convergence within max_iter iterations(1). Upon suc- cessful return, output arguments have the following values: x approximate solution to Ax = b max_iter the number of iterations performed before the tolerance was reached tol the residual after the final iteration IMPLEMENTATION
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; } rheolef-6.1 rheolef-6.1 puzawa(4rheolef)
