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