function cg % CG on Ax = b % Define matrix A dim = 1000; % matrix dimension %A = hilbertMatrix(dim); A = generateMatrix(dim, 1); b = ones(dim, 1); ep = 1.0E-6; x0 = zeros(dim,1); x = cg_loop(A,b,x0,ep); end function x = cg_loop(A,b,x,ep) x1 = A\b; % "Exact" solution fprintf('#it\t||r||\t\tlog(||x-x*||_A)\n') r = A*x - b; p = -r; rr= r'*r; k = 0; while norm(r) > ep Ap = A*p; alpha = rr/(p'*Ap); x = x+alpha*p; r = r+alpha*Ap; rr1 = r'*r; beta = rr1/rr; p = -r+beta*p; rr = rr1; k = k+1; lognorm(k) = log(sqrt((x1-x)'*A*(x1-x))); fprintf('%d:\t%8f\t%8f\n', k, norm(r), lognorm(k)); end fprintf('CG converged after %d iterations\n', k); fprintf('Relative error compared to direct solve: %e\n', norm(x1-x)/norm(x1)); plot(1:k,lognorm) end function A = hilbertMatrix(dim) % A_ij = 1/(i+j-1) A = zeros(dim); for i=1:dim for j=1:dim A(i,j) = 1/(i+j-1); end end fprintf('Condition number of A: %g\n', cond(A)); end function A = generateMatrix(n, type) % Create matrix of dimension n (multiple of 5) and cond(A)=100 % Matrix types: % 1: 5 clustered eigenvalues % 2: 5 clustered and small perturbation % 3: uniformly distributed eigenvalues if (mod(n,5) ~= 0) error('Dimension should be multiple of 5'); end rng('default'); % To produce same matrix every time M = rand(n); % Random nxn matrix with elements in the range [0,10] [Q, R] = qr(M); % Q is now unitary if (type == 1) || (type == 2) rep = n/5; e = ones(rep,1); D = [1*e; 2*e; 10*e; 20*e; 100*e]; if (type == 2) for i=1:n D(i) = D(i) + 0.001*(0.5-rand())*D(i); end end elseif (type == 3) D = 1:(99/(n-1)):100; else error('Type should be 1,2 or 3') end D = diag(D); A = Q'*D*Q; end