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