puzawa(4rheolef) [debian man page]
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)
Check Out this Related Man Page
pminres(4rheolef) rheolef-6.1 pminres(4rheolef) NAME
pminres -- conjugate gradient algorithm. SYNOPSIS
template <class Matrix, class Vector, class Preconditioner, class Real> int pminres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, int &max_iter, Real &tol, odiststream *p_derr=0); EXAMPLE
The simplest call to 'pminres' has the folling form: size_t max_iter = 100; double tol = 1e-7; int status = pminres(a, x, b, EYE, max_iter, tol, &derr); DESCRIPTION
pminres solves the symmetric positive definite linear system Ax=b using the Conjugate Gradient method. 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 NOTE
pminres follows the algorithm described in "Solution of sparse indefinite systems of linear equations", C. C. Paige and M. A. Saunders, SIAM J. Numer. Anal., 12(4), 1975. For more, see http://www.stanford.edu/group/SOL/software.html and also the PhD "Iterative methods for singular linear equations and least-squares problems", S.-C. T. Choi, Stanford University, 2006, http://www.stanford.edu/group/SOL/disser- tations/sou-cheng-choi-thesis.pdf at page 60. The present implementation style is inspired from IML++ 1.2 iterative method library, http://math.nist.gov/iml++. IMPLEMENTATION
template <class Matrix, class Vector, class Preconditioner, class Real, class Size> int pminres(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M, Size &max_iter, Real &tol, odiststream *p_derr = 0, std::string label = "minres") { Vector b = M.solve(Mb); Real norm_b = sqrt(fabs(dot(Mb,b))); if (norm_b == Real(0.)) norm_b = 1; Vector Mr = Mb - A*x; Vector z = M.solve(Mr); Real beta2 = dot(Mr, z); Real norm_r = sqrt(fabs(beta2)); if (p_derr) (*p_derr) << "[" << label << "] #iteration residue" << std::endl; if (p_derr) (*p_derr) << "[" << label << "] 0 " << norm_r/norm_b << std::endl; if (beta2 < 0 || norm_r <= tol*norm_b) { tol = norm_r/norm_b; max_iter = 0; dis_warning_macro ("beta2 = " << beta2 << " < 0: stop"); return 0; } Real beta = sqrt(beta2); Real eta = beta; Vector Mv = Mr/beta; Vector u = z/beta; Real c_old = 1.; Real s_old = 0.; Real c = 1.; Real s = 0.; Vector u_old (x.ownership(), 0.); Vector Mv_old (x.ownership(), 0.); Vector w (x.ownership(), 0.); Vector w_old (x.ownership(), 0.); Vector w_old2 (x.ownership(), 0.); for (Size n = 1; n <= max_iter; n++) { // Lanczos Mr = A*u; z = M.solve(Mr); Real alpha = dot(Mr, u); Mr = Mr - alpha*Mv - beta*Mv_old; z = z - alpha*u - beta*u_old; beta2 = dot(Mr, z); if (beta2 < 0) { dis_warning_macro ("pminres: machine precision problem"); tol = norm_r/norm_b; max_iter = n; return 2; } Real beta_old = beta; beta = sqrt(beta2); // QR factorisation Real c_old2 = c_old; Real s_old2 = s_old; c_old = c; s_old = s; Real rho0 = c_old*alpha - c_old2*s_old*beta_old; Real rho2 = s_old*alpha + c_old2*c_old*beta_old; Real rho1 = sqrt(sqr(rho0) + sqr(beta)); Real rho3 = s_old2 * beta_old; // Givens rotation c = rho0 / rho1; s = beta / rho1; // update w_old2 = w_old; w_old = w; w = (u - rho2*w_old - rho3*w_old2)/rho1; x += c*eta*w; eta = -s*eta; Mv_old = Mv; u_old = u; Mv = Mr/beta; u = z/beta; // check residue norm_r *= s; if (p_derr) (*p_derr) << "[" << label << "] " << n << " " << norm_r/norm_b << std::endl; if (norm_r <= tol*norm_b) { tol = norm_r/norm_b; max_iter = n; return 0; } } tol = norm_r/norm_b; return 1; } rheolef-6.1 rheolef-6.1 pminres(4rheolef)