clc;clear;
for Ex=1:3
    if Ex==1
        Problem ='bcsstk08';
        load([Problem '.mat']);
    elseif Ex==2
        Problem ='nasa1824';
        load([Problem '.mat']);
    elseif Ex==3
        Problem ='pds02';
        load([Problem '.mat']);
    end
    [m n]=size(A);
    str=sprintf('\nFor problem %s: n = %d', Problem, n); disp(str);
    xstar=ones(n, 1); b=A*xstar;
    tic
    Ainv=inv(A);
    x=Ainv*b;
    AinvSolveTT = toc;
    
    rsduNorm = norm(A*x-b,2); errNorm = norm(x-xstar,2);
    str=sprintf('Ainv solve  time = %6.4f; ||Ax-b|| = %11.4e; ||x-x^*|| = %11.4e', AinvSolveTT, rsduNorm, errNorm);disp(str);
    
    tic
    x=A\b;
    AbackslashbSolveTT = toc;
    
    rsduNorm = norm(A*x-b,2); errNorm = norm(x-xstar,2);
    str=sprintf('A backslash time = %6.4f; ||Ax-b|| = %11.4e; ||x-x^*|| = %11.4e', AbackslashbSolveTT, rsduNorm, errNorm);disp(str);
end
str=sprintf('\nHit enter to continue');disp(str); pause
clc;clear;
for Part=1:2
    %
    %Th e value of ToK forces the right number of loops
    %
    if Part==1
        ToK = 7;
    else
        ToK = 4;
    end
    for k=1:ToK
        if Part==1
            n=50*(2^k)
            A=ArrowMatrix(n);
        else
            if k==1
                load('bcsstk08.mat');
            elseif k==2
                load('nasa2910.mat');
            elseif k==3
                load('c-36.mat');
            else
                load('pds02.mat');
            end
            n = size(A,1);
        end
        nnzA=nnz(A);
        str=sprintf('\nA has dimesnion %d and %d nonzeros:', n, nnzA);disp(str)
        %Display the sparsity pattern
        spy(A)
        str=sprintf('\nHit enter to continue');disp(str); pause
        %
        %Perform the experiments with inv(A)
        %
        tic;Ainv=inv(A);AinvTime=toc;
        nnzAinv=nnz(Ainv);
        densityAinv=100*nnzAinv/(n*n);
        str=sprintf('inv(A) gives A^{-1} with %d nonzeros (fill factor %g and density %2.0f%%) in %g secs', nnzAinv, nnzAinv/nnzA,densityAinv,AinvTime);disp(str)
        str=sprintf('Time =  %g secs', AinvTime);disp(str)
        %
        %Perform the experiments with [L U P]=lu(A)
        %
        tic;[L U P]=lu(A);LUtime=toc;
        nnzLU=nnz(L)+nnz(U)-n;
        densityLU=100*nnzLU/(n*n);
        str=sprintf('[L U P]=lu(A) gives L and U with %d nonzeros (fill factor %g and density %2.0f%%) in %g secs', nnzLU,nnzLU/nnzA,densityLU,LUtime);disp(str)
        str=sprintf('Time =  %g secs', LUtime);disp(str)
        
        %
        %Generate the sparse version of A (since it is full in the case of the Arrowhead experiemnts
        %
        SparseA=sparse(A);
        %
        %Perform the experiments with [L U P Q]=lu(A)
        %
        tic;[L U P Q]=lu(SparseA);LUtime=toc;
        nnzLU=nnz(L)+nnz(U)-n;
        densityLU=100*nnzLU/(n*n);
        str=sprintf('[L U P Q]=lu(SparseA) gives L and U with %d nonzeros (fill factor %g and density %2.0f%%) in %g secs', nnzLU,nnzLU/nnzA,densityLU,LUtime);disp(str)
        str=sprintf('Time =  %g secs', LUtime);disp(str)
        %
        %If the density of the inverse is low, illustrate its sparsity pattern
        %
        if densityAinv < 5
            SparseAinv = sparse(Ainv);
            spy(SparseAinv)
            str=sprintf('\nHit enter to continue');disp(str)
            pause
        end
        
    end
end
str=sprintf('\nHit enter to continue');disp(str); pause
clc;clear;
for Ex=2:6
    eps=1e-1 ;
    n = 1000;
    rng(1);
    if Ex == 1
        Problem ='A = [2 -1; -1 3]';
        A = [2 -1; -1 3];
    elseif Ex==2
        Problem ='A=eye(n)';
        A = eye(n);
    elseif Ex==3
        Problem ='662_bus';
        load([Problem '.mat']);
    elseif Ex==4
        Problem ='bcsstk08';
        load([Problem '.mat']);
    elseif Ex==5
        Problem ='I + eps*SpSymmRand';
        n=1000;
        A=speye(n) + eps*SpSymmRand(n, 10*n);
    elseif Ex==6
        Problem ='I + uu^T';
        n=1000;
        u = rand(n, 1);u=u/norm(u,2);
        A=eye(n) + u*u';
    end
    n = size(A, 1);
    fprintf('\nProblem: %s with n = %d\n', Problem, n)
    AnalyseEigs(A)
    fprintf('Press any key to continue\n');pause
    
    if Ex == 1
        xsol=ones(n, 1);
    else
        xsol=rand(n, 1);
    end  
    b=A*xsol;

    x=zeros(n, 1);

    tic;
    [x m] = MyCG(A, b);
    CGtime = toc;

    CG_r_error = norm(A*x-b,2);
    CG_x_error = norm(x-xsol,2);
    fprintf('CG: %d iterations in %g sec for ||Ax-b||=%g and ||x-x^*||=%g\n', m, CGtime, CG_r_error, CG_x_error)

    tic; x=A\b;BackslashTime=toc;
    Backslash_r_error=norm(A*x-b,2);
    Backslash_x_error=norm(x-xsol,2);
    fprintf('Backslash: %g sec for ||Ax-b||=%g and ||x-x^*||=%g\n', BackslashTime, Backslash_r_error, Backslash_x_error)

end