Home > K25104 > LinearSystems > untransformed_preconditioned_CG.m

untransformed_preconditioned_CG

PURPOSE ^

Implements method of untransformed, preconditioned conjugate gradients

SYNOPSIS ^

function [ x,k ] = untransformed_preconditioned_CG(A,b,S,tol )

DESCRIPTION ^

 Implements method of untransformed, preconditioned conjugate gradients 
   to find the solution of the system Ax=b.
 Input arguments:
   A, symmetric positive definite matrix
   b, column vector
   S, preconditioner, nonsingular matrix with the same dimensions as A.
   tol, convergence is accepted when |Ax-b|<tol
 Output arguments:
   x, solution to Ax=b
   k, number of iterations

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ x,k ] =  untransformed_preconditioned_CG(A,b,S,tol )
0002 % Implements method of untransformed, preconditioned conjugate gradients
0003 %   to find the solution of the system Ax=b.
0004 % Input arguments:
0005 %   A, symmetric positive definite matrix
0006 %   b, column vector
0007 %   S, preconditioner, nonsingular matrix with the same dimensions as A.
0008 %   tol, convergence is accepted when |Ax-b|<tol
0009 % Output arguments:
0010 %   x, solution to Ax=b
0011 %   k, number of iterations
0012 
0013 [n,m]=size(A); % finding the size of A
0014 [p,q]=size(b); % finding the size of b
0015 [i,j]=size(S); % finding the size of P
0016 if n~= m;
0017     error('first input is not a square matrix');
0018 elseif q~=1;
0019     error('second input is not a column vector');
0020 elseif p~=n || i~=n || j~=m;
0021     error('input dimensions do not agree');
0022 elseif tol<=0;
0023     error('tolerance should be positive');
0024 elseif ~issymmetric(A);
0025     error('first input is not symmetric');
0026 elseif rank(S)<n
0027     error('preconditioner is singular');
0028 else
0029     [~,a] = chol(A);
0030     if a
0031         error('first input is not positive definite');
0032     end
0033 end
0034 
0035 x=zeros(n,1);   % initialize x
0036 r0=b;           % initialize r0
0037 d=S*r0;         % initialize d
0038 k=0;            % initialize counter
0039 
0040 if sqrt(r0'*r0)<tol; % if this is true then the zero vector is the solution
0041     return;
0042 else                % implement first step
0043     v=A*d;           % calculate this matrix-vector product, which will
0044                      % be needed twice more
0045     w=r0'*S*r0/(d'*v); % calculate step size in the descent direction
0046     x=x+w*d;         % update x
0047     r1=r0-w*v;       % update the residual
0048     k=k+1;           % increment counter
0049 end
0050 
0051 while sqrt(r1'*r1)>tol; % repeat the following until convergence
0052     B=r1'*S*r1/(r0'*S*r0);  % this choice of B ensures the search directions
0053                         % are conjugate with respect to A. Note that to
0054                         % calculate B we need both the current and
0055                         % previous residual, so we must keep track of
0056                         % both r0 and r1.
0057     d=S*r1+B*d;           % update the search direction
0058     r0=r1;              % update r0
0059     v=A*d;              % calculate this matrix-vector product, which will
0060                         % be needed twice more
0061     w=r0'*S*r0/(d'*v);    % calculate step size in the descent direction
0062     x=x+w*d;            % update x
0063     r1=r0-w*v;          % update r1
0064     k=k+1;              % increment counter
0065 end
0066 
0067 end

Generated on Mon 18-Jan-2016 10:25:49 by m2html © 2005