function [x, error, iter, flag] = cg(A, x, b, M, max_it, tol)%  -- Iterative template routine --%     Univ. of Tennessee and Oak Ridge National Laboratory%     October 1, 1993%     Details of this algorithm are described in "Templates for the%     Solution of Linear Systems: Building Blocks for Iterative%     Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra,%     Eijkhout, Pozo, Romine, and van der Vorst, SIAM Publications,%     1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps).%%  [x, error, iter, flag] = cg(A, x, b, M, max_it, tol)%% cg.m solves the symmetric positive definite linear system Ax=b % using the Conjugate Gradient method with preconditioning.%% input   A        REAL symmetric positive definite matrix%         x        REAL initial guess vector%         b        REAL right hand side vector%         M        REAL preconditioner matrix%         max_it   INTEGER maximum number of iterations%         tol      REAL error tolerance%% output  x        REAL solution vector%         error    REAL error norm%         iter     INTEGER number of iterations performed%         flag     INTEGER: 0 = solution found to tolerance%                           1 = no convergence given max_it  flag = 0;                                 % initialization  iter = 0;  bnrm2 = norm( b );  if  ( bnrm2 == 0.0 ), bnrm2 = 1.0; end  r = b - A*x;  error = norm( r ) / bnrm2;  if ( error < tol ) return, end  for iter = 1:max_it                       % begin iteration     z  = M \ r;     rho = (r'*z);     if ( iter > 1 ),                       % direction vector        beta = rho / rho_1;        p = z + beta*p;     else        p = z;     end     q = A*p;     alpha = rho / (p'*q );     x = x + alpha * p;                    % update approximation vector     r = r - alpha*q;                      % compute residual     error = norm( r ) / bnrm2;            % check convergence     if ( error <= tol ), break, end      rho_1 = rho;  end  if ( error > tol ) flag = 1; end         % no convergence% END cg.m