% demonstration of preconditioned conjugate gradient method:
% using M and K from 1d FEM with hat functions

 N = 256;
 h = 1/(N+1);
 x = h*[1:N]';

 [M,K] = makeMK(N);

 f1 = randn(N,1);
 f2 = exp(x);

 [a,flag,relres,iter,resvec] = pcg(K,f2,1e-6,max(200,N));

 figure(1), clf
 semilogy([0:length(resvec)-1], resvec,'b.-','linewidth',2,'markersize',20)
 hold on

 L = ichol(sparse(K));
 [a,flag,relres,iter,resvec] = pcg(K,f2,1e-6,max(200,N),L,L');
 semilogy([0:length(resvec)-1], resvec,'r.-','linewidth',2,'markersize',20)
