function [] = pcgdemo(matrix_flag,precond,n)
%%%%%%%%%%%%%%%%%%     PCGDEMO       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%  For MTH 451-551  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This function demostrates the difference between the conjugate gradient
%% method with and without a Preconditioner applied to the system.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% Doug Fettig,  02/2005
%% supported by DeLoach Work Scholarship Program 
%% of the Oregon State University Honors College
%% http://oregonstate.edu/dept/honors/students/opportunities/DeLoach.htm
%% M. Peszynska, 03/2005
%% © Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% INSTRUCTIONS: 
%% calling sequence
%%      pcgdemo(matrix_flag,precond,size)
%% matrix_flag determines the type of a test matrix in the demo:
%%          1 - 1D finite difference matrix: discrete Laplacian
%%          2 - 1D finite difference matrix: Poisson problem with random coefficients
%%          3 - two dimensional finite difference matrix.
%%          4 - a random matrix (actually, its spd 'square B'*B)
%%          5 - a 'wathen' matrix (see help gallery)  
%%          6 - the 'lehmer' matrix (see help gallery)
%% precond determines which preconditioner will be applied:
%%          1 - The diagonal or Jacobi preconditioner
%%          2 - The diagonal matrix made of entries 
%%          3 - The Gauss-Seidel preconditioner (SOR with w = 1)
%%          4 - The SSOR preconditioner with \omega = 1.7 
%%          5 - The incomplete cholesky preconditioner with eps=1.e-3
%% n determines the size of the matrix: in case 2,4 it is not equal to size(A)
%% 
%% A call to PCGDEMO with all three input parameters will display a 
%% plot of the residual vs. the iterations with and without the 
%% preconditioner applied.
%%
%% A call to PCGDEMO with only the first 2 parameters will show the spectrum 
%% for the matrix and preconditioner chosen for various sizes, and residuals
%% for the largest size.
%% A call to PCGDEMO with only the first 1 parameter will show the do above for all
%% preconditioners
%% 
%% A call to PCGDEMO with no input arguments will run a default demo. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clf;

if nargin == 0                    
  matrix_flag = 1;
end;

if nargin < 2 

  for precond = 1:5  
     clf;
     show_preconditioners (matrix_flag , precond); 
     if precond== 1
         ['matrix type = ',' preconditioner = '];
     end
     [matrix_flag,precond],pause, 
  end;

elseif nargin == 2

  show_preconditioners (matrix_flag, precond)

elseif nargin > 2 

  show_residuals (matrix_flag, precond, n);

else
   return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  matrix setup

function [A] = setup_matrix ( n, matrix_type)

    switch matrix_type                                
      case 1                       %% 1 dimensional finite difference
        A = 2 * diag(ones(n,1),0) - diag(ones(n-1,1),1) - diag(ones(n-1,1),-1);
      case 2
        rand('state',0);
        cofs = 10*rand(n+1,1);
        A = zeros (n,n);
        for i=1:n A(i,i)= cofs(i)+cofs(i+1); end;
        for i=2:n A(i,i-1)= -cofs(i); A(i-1,i)=-cofs(i);end;

      case 3                      %% 2  - dimensional finite difference
        A = full(delsq(numgrid('S',n)'));
                                %% one can also use full(gallery('poisson',n));
      case 4                      %% a 'spd 'square' of a random matrix
        rand('state',0);
        Random = rand(n,n);                
        A = Random' * Random;
      case 5
        A = full(gallery('wathen',n,n));    
      case 6                      %% Initializes the demo using
        A = gallery('lehmer',n);            %% the Lehmer matrix, which is 
      otherwise
        Error = ['Incorrect entry for matrix variable,'...    
                '   Please input either 1,2,3 or 4   '];
        A = [];
        return
    end

function [P] = setup_preconditioner ( A , precond)

    D = diag (diag (A)) ;
    L = tril (A) - D ;
        
 switch precond 
    case 1        %% Jacobi or diagonal preconditioner
        P = D;
    case 2         %% 
        P = D;
        for i=1:size(A,1) P(i,i)=norm(A(i,:),2); end;
        
    case 3        %% SSOR preconditioner
        P = (D + L) * inv(D) * (D + L)';
    case 4         %% Gauss-Seidel preconditioner.
        omega = 1.7; %% should be the optimal SOR parameter which is hard to find
        P = (D + omega*L) * inv(D) * (D + omega*L)';
        P = P* omega/(2-omega);

    case 5        %% Incomplete cholesky preconditioner
        C = full( cholinc( sparse(A), 1e-3) );
        P = C' * C;      
    otherwise
        Error = ['Incorrect entry for precond variable,' %%Error messsage for 
                 '     Please input either 1,2 or 3    ']%%invalid precond entry
        P = eye ( size (A) );
        return
 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function show_preconditioners (matrix_flag, precond)

 maxn = 100;
 Cond_results = zeros (3,maxn);      %% will store the condition number(s)

 for n = 4 : maxn 

     A = setup_matrix (n, matrix_flag);
     sizeA = size(A,1);
     P = setup_preconditioner (A, precond);
     Cond_results(1 , n ) = sizeA;
     Cond_results(2 , n ) = cond ( A );            %% no preconditioner
     Cond_results(3 , n ) = cond( (inv(P) * A) );  %% with preconditioner

     if ( sizeA >= maxn)
        lastn = n; 
        break; 
     end;
  end;

  %% find and show the spectrum for the last matrix
  eigA = eig(A);
  eigPA = eig(inv(P)*A);
  clf;
  subplot(2,2,1); [m,heigA] = hist(eigA,sizeA); hist(eigA,sizeA); 
  title('Spectrum of matrix');
  subplot(2,2,2); hist(eigPA,heigA);
  title('Spectrum of preconditioned matrix');

  %% show the condition numbers
  subplot(2,2,3); 
  plot1 = loglog (Cond_results(1,:), Cond_results(2,:), 'b'); %% plot without preconditioner.
  set(plot1, 'linewidth',4);
    
  hold on;
  plot2 = loglog(Cond_results(1,:) , Cond_results(3,:), 'r'); %% plot with preconditioner.
  set(plot2, 'linewidth',3) ;
                                        
  title('Condition number vs. the size of the matrix',...
	'color', 'k');
  xlabel(' N (size of an N-by-N matrix) ', 'fontsize', 13)
  ylabel(' Condition Number ','fontsize', 13)
  label1 = [' No preconditioner  , '];
  label2 = [' With preconditioner '];
  legend(label1,label2,2);

  show_residuals (matrix_flag, precond, lastn)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function show_residuals (matrix_flag, precond, n)

  A = setup_matrix(n, matrix_flag);
  sizeA = size(A,1);
  P = setup_preconditioner (A, precond);
  x = ones(sizeA, 1);
  b = A * x;
                                 
  maxit = size(A,1);                     
  tol = 1e-6;  
                                %% Call to pcg algorithm that stores a vector
                                %% of the residuals for each iteration.
  [h,h,h,h,ResVec1] = pcg(A,b,tol,maxit,P);   %% With preconditioner. 
  [h,h,h,h,ResVec2] = pcg(A,b,tol,maxit);     %% No preconditioner.

  Iter1 = length(ResVec1); Iter2 = length(ResVec2);                                  

                               %% plot of residual with preconditioner applied.
  subplot(2,2,4);
  plot1 = semilogy(0 : Iter1 - 1 , ResVec1 / (ResVec1(1,1)), 'r') ; 
  set(plot1, 'linewidth',4) ;

  hold on
                               %% plot of residual without preconditioner. 
  plot2 = semilogy(0 : Iter2 - 1 , ResVec2  / ResVec2(1,1), 'b');
  set(plot2, 'linewidth',3) ;    
               
                                %% sets plot and display options                                         
  xlabel(' Iteration ', 'fontsize', 13)
  ylabel(' Relative Residual Value','fontsize', 13)

  title('Plot of Residual vs. the Iteration #','color', 'k');

   label1 = ['No Preconditioner,   cond(A)        = '...
                num2str( cond( full( A ) ) )];
   label2 = ['With Preconditioner, cond(Pinv*A) = '...
                num2str( cond( inv(P) * A ) ) ];
 
   legend(label2,label1,3);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  End of PCGDEMO  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
