  function P4B()
% Illustrates GMRES with incomplete LU preconditioner obtained via LDMinc.
close all
% An example that isn't really diagonnaly dominant but OK..
fid = fopen('LU_156.txt','r');
i=[]; j = []; vals=[];
while ~feof(fid)
    s = fgetl(fid); t = str2num(s);
    i = [i; t(1)]; j = [j; t(2)]; vals = [vals; t(3)];
end
fclose(fid);
A = sparse(i,j,vals);
[n,n] = size(A);
A = A + 0.007*norm(A,1)*eye(n,n);
spy(A)
title('The Matrix A','fontsize',14)
shg
b = randn(n,1);
fprintf('n = %1d, nnz(A) = %1d\n\n',n,nnz(A))

% Key parameters for the preconditioned GMRES method...
% Maximum number of iterations...
itMax = 1000;
% Stop if || Ax - b||/||b| <= tol
tol = .00000000001;

% Without any preconditioner...
L = speye(n);
U = speye(n);
[x,flag,relres,itLast] = gmres(A,b,[],tol,itMax,L,U);
fprintf('  Preconditioner    Storage      relres        itLast     \n')
fprintf('------------------------------------------------------------\n')
fprintf(' %15s    %5d       %10.3e     %4d      \n','None',0,relres,itLast(2))

tYOU = 0;
for p = [5]
   % Incomplete LDM...
   tic
   [L,d,M]  = LDMinc(A,p,p);
   tYOU = tYOU + toc;
   % Form LU factors because that is what gmres wants...
   U = diag(d)*M';
   figure
   subplot(1,2,1)
   spy(L)
   title(sprintf('L    (p = %1d)',p),'fontsize',14)
   subplot(1,2,2)
   spy(U) 
   title(sprintf('U    (q = %1d)',p),'fontsize',14)
   [x,flag,relres,itLast] = gmres(A,b,[],tol,itMax,L,U);
   s = sprintf('LDMinc(A,%1d,%1d)',p,p);
   fprintf(' %15s    %7d     %10.3e     %4d       \n',s,nnz(L)+nnz(M),relres,itLast(2))  
end

tic
for p = [5]
   % Incomplete LDM...
   [L,d,M]  = LDMincCVL(A,p,p);
end
tCVL = toc;

fprintf('\nn = %5d   tYOU = %7.3f',n,tYOU)
fprintf('\nn = %5d   tCVL = %7.3f',n,tCVL)


