% Script demonstrating Cholesky and incomplete-Cholesky decomposition of 
% matrix A resulting from finite differences applied to 2D Poisson problem 
%     u_{xx} + u_{yy} = 4, 
%     u(0,y) = y^2, u(1,y) = 1+y^2, 
%     u(x,0) = x^2, u(x,1) = 1+x^2. 
% with h = 1/2, 1/4, 1/8, ..., 1/32. Plot sparsity of decompositions for 
% h = 1/8 case and condition number of A and C^-1*A*C^-1 for all h.

  pmax = 5;
  condA= zeros(1,pmax);
  precondA = zeros(1,pmax);
  clf
  
  for p = 1:pmax
    N = 2^p;
    m = N-1;
    h = 1/N;

%    Define A using matlab functions delsq and numgrid. Note that u(i,j)
%    refers to u(x(i),y(j)).  In converting to a vector, values are ordered
%    u(1,1),u(2,1),...,u(m,1),u(1,2),...,u(m,2),...,u(1,m),...u(m,m).

    A = delsq(numgrid('S',m+2))/h^2;  
    R = (chol(A,'lower'))';% Standard Cholesky decomp. A = GG'
    C = (ichol(A))';    % Incomplete Cholesky decomp. A approx HH'.
 
%    Plot h = 1/8 (p = 3) case

    if (p==3)  
      subplot(2,2,1)
      spy(R)
      title('Sparsity of R = chol(A), h=1/8')
      subplot(2,2,2)
      spy(C)
      title('Sparsity of C = ichol(A)'', h=1/8')
    end
    Cinv = inv(C);
    condA(p) = condest(A);
    precondA(p) = condest(Cinv*A*Cinv');
  end

  subplot(2,2,3)
  prange = 1:pmax;
  hrange = 2.^-prange;
  hmax = 2^-pmax;
  loglog(hrange,condA(prange),'b+',hrange, precondA(prange),'gx',...
         hrange,(hrange/hmax).^-2*condA(pmax),'b-')  
  xlabel('h')
  ylabel('Condition number of 2D FDA del^2') 
  legend('A','IC precond(A)')
