Home > K25104 > LinearSystems > conjugate_gradient.m

conjugate_gradient

PURPOSE ^

Implements method of conjugate gradients

SYNOPSIS ^

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

DESCRIPTION ^

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

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