GNU.WIKI: The GNU/Linux Knowledge Base

  [HOME] [PHP Manual] [HowTo] [ABS] [MAN1] [MAN2] [MAN3] [MAN4] [MAN5] [MAN6] [MAN7] [MAN8] [MAN9]

  [0-9] [Aa] [Bb] [Cc] [Dd] [Ee] [Ff] [Gg] [Hh] [Ii] [Jj] [Kk] [Ll] [Mm] [Nn] [Oo] [Pp] [Qq] [Rr] [Ss] [Tt] [Uu] [Vv] [Ww] [Xx] [Yy] [Zz]


       puzawa -- Uzawa algorithm.


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


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


       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   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.

       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

              the number of iterations  performed  before  the  tolerance  was

       tol    the residual after the final iteration


       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;

  All copyrights belong to their respective owners. Other content (c) 2014-2018, GNU.WIKI. Please report site errors to
Page load time: 0.185 seconds. Last modified: November 04 2018 12:49:43.