Home > K25104 > LinearSystems > transformed_preconditioned_CG.m

transformed_preconditioned_CG

PURPOSE ^

Implements method of transformed, preconditioned conjugate gradients

SYNOPSIS ^

function [ x,k ] = transformed_preconditioned_CG(A,b,P,tol )

DESCRIPTION ^

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

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