function[] = cgdemo(v,A,b,eps,kmax)
%%%%%%%%%%%%%%%%%%     CGDEMO       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% For MTH 451-551
%% This function demostrates the difference between the conjugate gradient
%% and the steepest descent methods graphically. Uses a call to the function
%% MYCGS.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% Doug Fettig  01/2005
%% supported by DeLoach Work Scholarship Program 
%% of the Oregon State University Honors College
%% http://oregonstate.edu/dept/honors/students/opportunities/DeLoach.htm
%% Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% INSTRUCTIONS: 
%% calling sequence
%%      cgdemo (v,A,b,eps,kmax) 
%% v is the initial guess (a vector)
%% A is a positive definite, symmetric matrix
%% b is the right hand side
%% eps, kmax determine the stopping criterium when the mycgs algorithm is called: 
%% 
%% A call to CGDEMO with no input will run a default demo,
%% finding the solution to [1,0;0,5] [x,y]'=[0,0]
%% which is equivalent to minimization of a quadratic
%% f(x,y) = 0.5x^2 + 2.5y^2
%% initial. guess of v = [1;1] is used.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

hold off

if nargin < 5                                  %%% checks if default settings are needed
    v = [1 ; 1];
    A = [0.5 0; 0 2.5];
    eps = .00001;
    kmax = 10;   
    b = [0;0];
end

if size(A,1) ~= 2                               %%% check to make sure A is 2x2.
    Error = ['A must be a 2x2 matrix']
    return
end

                                                %% prepare graphical display settings   
sol = A\b;                             
                                                
if (norm(sol,inf) < norm(v,inf) )               %% finds an appropriate scaling factor m
    m = norm(v,inf);                            %% for display and plotting purposes.
elseif (norm(sol,inf) > norm(v,inf))
    m = norm(sol,inf); 
else
    m = 1;   
end

z1 = (-m + sol(1) : m/10 : m + sol(1));          %% region of plot display: z1 is x-axis
z2 = (-m + sol(2) : m/10 : m + sol(2));          %% z2 is the y-axis
[X,Y] = meshgrid(z1,z2);                         %% initialize grid matrices
                                                 %% R sets up the 3D surface to be plotted
R = A(1,1)*X.^2 + A(2,2)*Y.^2 + 2*A(1,2)*(X.*Y) - X.*b(1) - Y.*b(2); 


for j = 1:-1:0
    [s,k,a,I] = mycgs(v,A,b,eps,kmax,j);         %% get all the iterates 
                                                 
    if k == 0                                    %% stop if mycgs terminated early
        return                                   
    end
    
    if j == 1                                    %% steepest descent first
        mesh(z1,z2,R);                           %% display 3D suface   
        view(60,60)                              %% viewpoint settings
        hold on                                             
        hidden off                                  
        grid off
        colormap(summer)                         %% color scheme
        title('Solution using steepest descent method', 'color', 'b')    
        pause
    end
    
    if j == 0                                     %% using CG 
        title('Now using Conjugate Gradient method', 'color', 'r')
        pause
    end
 
    for i = 1:k-1                                %% display vectors
        x1 = I(1,i);
        y1 = I(2,i);
        x2 = I(1,i + 1);
        y2 = I(2,i + 1);
        f1 = A(1,1)*x1.^2 + A(2,2)*y1.^2 + 2*A(1,2)*(x1.*y1) - x1.*b(1) - y1.*b(2);
        f2 = A(1,1)*x2.^2 + A(2,2)*y2.^2 + 2*A(1,2)*(x2.*y2) - x2.*b(1) - y2.*b(2);
        
        if j == 1                               %% displays vectors differently, depending on method.
           h = quiver3(x1, y1, f1, x2 - x1 , y2 - y1, f2-f1, 0, 'b*');
        else
           h = quiver3(x1, y1, f1, x2 - x1 , y2 - y1, f2-f1, 0, 'r+');
        end
        set(h,'LineWidth',3)
        pause
    end
end
