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