function JamesWenner()

clear all;
format long
% Math 105LB Final Project
% By James Joseph Wenner
% ID#84239664

%Initialize the results file & print headers.
FID=fopen('JamesWenner.txt','w');
fprintf(1,'\nJames Joseph Wenner\nID#84239664\n');
fprintf(1,'Math 105LB Final Project\n\n\n');
fprintf(FID,'James Joseph Wenner\nID#84239664\n');
fprintf(FID,'Math 105LB Final Project\n\n\n');

%Problem 1
%Determine several eigenvalues of the given matrix.
%1: Initialize several variables and print header.
%2: Use PowerEig (for Power Method) & CheckEig to determine the largest
%   eigenvalue of A.  Print out this eigenvalue and the number of
%   iterations required.
%3: Use InvPowerEig (for Inverse Power Method) & CheckEig to determine the
%   smallest eigenvalue of A.  Print out this eigenvalue and the number of
%   iterations required.
%4: Use Deflation (for Hotelling deflation) to determine the second largest
%   eigenvalue of A. Pring out this eigenvalue and the number of iterations
%   required.
A=[4,1,1,1;1,3,-1,1;1,-1,2,0;1,1,0,2];
x0=[1;-2;0;3];
TOL=10^-7;
MAXN=1000;

fprintf(1,'Problem 1\n');
fprintf(1,'Find several eigenvalues of the given matrix A.\n\n');
fprintf(FID,'Problem 1\n');
fprintf(FID,'Find several eigenvalues of the given matrix A.\n\n');

[LARGE,LARGEN,LARGEVECT]=PowerEig(A,x0,TOL,MAXN);
LARGEVAL=CheckEig(A,LARGE,LARGEVECT,TOL);
if ~isnan(LARGEVAL)
    fprintf(1,'The largest eigenvalue of A is %.7f.\n',LARGEVAL);
    fprintf(1,'The power method required %d iterations to find',LARGEN);
    fprintf(1,' this eigenvalue.\n\n');
    fprintf(FID,'The largest eigenvalue of A is %.7f.\n',LARGEVAL);
    fprintf(FID,'The power method required %d iterations to find',LARGEN);
    fprintf(FID,' this eigenvalue.\n\n');
else
    fprintf(1,'Error in finding largest eigenvalue.\n\n');
    fprintf(FID,'Error in finding largest eigenvalue.\n\n');
end

[SMALL,SMALLN,SMALLVECT]=InvPowerEig(A,x0,TOL,MAXN);
SMALLVAL=CheckEig(A,SMALL,SMALLVECT,TOL);
if ~isnan(SMALLVAL)
    fprintf(1,'The smallest eigenvalue of A is %.7f.\n',SMALLVAL);
    fprintf(1,'The inverse power method required %d iterations ',SMALLN);
    fprintf(1,' to find this eigenvalue.\n\n');
    fprintf(FID,'The smallest eigenvalue of A is %.7f.\n',SMALLVAL);
    fprintf(FID,'The inverse power method required %d iterations ',SMALLN);
    fprintf(FID,' to find this eigenvalue.\n\n');
else
    fprintf(1,'Error in finding smallest eigenvalue.\n\n');
    fprintf(FID,'Error in finding smallest eigenvalue.\n\n');
end

[EIG2,N2]=DEFLATION(A,x0,TOL,MAXN,LARGEVAL,LARGEVECT);
if ~isnan(EIG2)
    fprintf(1,'The second largest eigenvalue of A is %.7f.\n',EIG2);
    fprintf(1,'The power method required %d iterations to find',N2);
    fprintf(1,' this eigenvalue.\n\n');
    fprintf(FID,'The second largest eigenvalue of A is %.7f.\n',EIG2);
    fprintf(FID,'The power method required %d iterations to find',N2);
    fprintf(FID,' this eigenvalue.\n\n');
else
    fprintf(1,'Error in finding second largest eigenvalue.\n\n');
    fprintf(FID,'Error in finding second largest eigenvalue.\n\n');
end

%Problem 2
%Determine all of the eigenvalues of the given matrix.
%1: Initialize several variables and print header.
%2: Use QREig (QR Algorithm) to determine all eigenvalues of A. Print out
%   these eigenvalue and the number of iterations required.
A = [0.5,0.25,0,0;0.25,0.8,0.4,0;0,0.4,0.6,0.1;0,0,0.1,1];
TOL=10^-7;
MAXN=1000;

fprintf(1,'\nProblem 2\n');
fprintf(1,'Find all eigenvalues of the given matrix.\n\n');
fprintf(FID,'\nProblem 2\n');
fprintf(FID,'Find all eigenvalues of the given matrix.\n\n');

[QREIG,QRN]=QREig(A,MAXN,TOL);
NEIG=length(QREIG);
if QRN<MAXN
    fprintf(1,'After %d iterations, the QR Algorithm gave as the ',QRN);
    fprintf(1,'the eigenvalues of A:\n   ');
    fprintf(FID,'After %d iterations, the QR Algorithm gave as the ',QRN);
    fprintf(FID,'the eigenvalues of A:\n   ');
else
    fprintf(1,'The QR Algorithm required the maximum number of ');
    fprintf(1,'iterations,\nso the accuracy can not be determined');
    fprintf(1,' for the eigenvalues:\n');
    fprintf(FID,'The QR Algorithm required the maximum number of ');
    fprintf(FID,'iterations,\nso the accuracy can not be determined');
    fprintf(FID,' for the eigenvalues:\n');
end
for i=1:NEIG
    fprintf(1,'%.7f   ',QREIG(i));
    fprintf(FID,'%.7f   ',QREIG(i));
end

fprintf(1,'\n\n\nEnd of program.\n\n\n');
fprintf(FID,'\n\n\nEnd of program.');
fclose(FID);



%[U,I,EIG]=POWEREIG(A,X,EPS,MAXCOUNT):
%Uses the Power Method to find the largest norm U of the eigenvalues of A,
%   which corresponds to the eigenvector EIG. X is the initial guess for 
%   the eigenvector, EPS is the desired tolerance, and MAXCOUNT is the
%   maximum number of iterations allowed. I is the number of iterations
%   used.
%1: Check to ensure that the initial guess is normalized.
%2: Implement the Power Method.
%   A: Check that the maximum number of iterations has not been reached.
%   B: Find y=Ax and normalize y to get x.
%   C: If both the inf-norm of the difference between two adjacent
%      eigenvectors and the error in the eigenvalue estimate is less than
%      EPS, quit and return the largest norm of an eigenvalue and its
%      eigenvector.
%   D: Otherwise, set the old eigenvector iteration as being the new guess.
function [u,i,x]=PowerEig(A,x,eps,maxcount)

if norm(x,inf)<eps
    u=NaN;
    i=inf;
    return;
end
xold=x/norm(x,inf);

for i=1:maxcount
    y=A*xold;
    u=norm(y,inf);
    if abs(u)<eps
        u=NaN;
        i=inf;
        return;
    end
    x=y/u;
    if norm(xold-x,inf)<eps && norm(A*x-u*x,2)<eps
        return;
    end
    xold=x;
end



%VAL=CheckEig(A,TEST,EVECT,TOL):
%Given a magnitude of an eigenvalue TEST and eigenvector EVECT for a matrix
%   A, determine the eigenvalue VAL of A. 
%1: If A*x-TEST*x<TOL, return TEST as the eigenvalue.
%2: If A*x+TEST*x<TOL, return -1*TEST as the eigenvalue.
%3: If neither of these cases is true, return error.
function VAL=CheckEig(A,TEST,EVECT,TOL)
TIMES=A*EVECT;
if max(max(abs(TIMES-TEST*EVECT)))<TOL
    VAL=TEST;
    return;
elseif max(max(abs(TIMES+TEST*EVECT)))<TOL
    VAL=-1*TEST;
    return;
else
    VAL=NaN;
end



%[U,I,EIG]=INVPOWEREIG(A,X,EPS,MAXCOUNT):
%Uses Inverse Power Method to find the largest norm U of the eigenvalues of
%   A, which corresponds to the eigenvector EIG. X is the initial guess for
%   the eigenvector, EPS is the desired tolerance, and MAXCOUNT is the
%   maximum number of iterations allowed. I is the number of iterations
%   used.
%1: Check to ensure that the initial guess is normalized.
%2: Calculate the inverse of A.
%3: Implement the Inverse Power Method.
%   A: Check that the maximum number of iterations has not been reached.
%   B: Find y=inv(A)*x and normalize y to get x.
%   C: If both the inf-norm of the difference between two adjacent
%      eigenvectors and the error in the eigenvalue estimate is less than
%      EPS, quit and return the smallest norm of an eigenvalue and its
%      eigenvector.
%   D: Otherwise, set the old eigenvector iteration as being the new guess.
function [u,i,x]=InvPowerEig(A,x,eps,maxcount)

if norm(x,inf)<eps
    u=NaN;
    i=inf;
    return;
end
xold=x/norm(x,inf);
AINV=inv(A);

for i=1:maxcount
    y=AINV*xold;
    u=norm(y,inf);
    if abs(u)<eps
        u=NaN;
        i=inf;
        return;
    end
    x=y/u;
    if norm(xold-x,inf)<eps && norm(A*x-x/u,2)<eps
        u=1/u;
        return;
    end
    xold=x;
end



%[EIG2,I]=DEFLATION(A,X0,TOL,MAXN,EIG1,EVECT1):
%Uses Hotelling deflation and then the Power method to find the second
%   largest eigenvalue EIG2 of the matrix A.  X0 is the initial guess for
%   the eigenvector, EIG1 is the largest eigenvalue with corresponding
%   eigenvector EVECT1, EPS is the desired tolerance, MAXCOUNT is the
%   maximum number of iterations allowed, and I is how many iterations
%   were conducted.
%1: Check that EVECT1, X0 are nonzero.
%2: Find the Hotelling deflation matrix B.
%3: Call upon PowerEig to determine the largest eigenvalue of B and the
%   corresponding eigenvector.  This eigenvalue of B is the norm of the
%   second largest eigenvalue of A.
%4: Call upon CheckEig to use the norm of the largest eigenvalue of B and
%   its corresponding eigenvector to determine the largest eigenvalue of B.
%   This eigenvalue is precisely the second largest eigenvalue of A.
function [EIG2,i]=DEFLATION(A,x0,TOL,MAXN,EIG1,EVECT1)
if norm(EVECT1,inf)<eps
    EIG2=NaN;
    return;
end
if norm(x0,inf)<eps
    EIG2=NaN;
    return;
end

B=A-(EIG1/norm(EVECT1,2))*EVECT1*EVECT1';
[EVAL2,i,EVECT2]=PowerEig(B,x0,TOL,MAXN);
EIG2=CheckEig(B,EVAL2,EVECT2,TOL);



%[EIG,I]=QREIG(A,EPS,MAXCOUNT):
%Uses the QR algorithm to find all the eigenvalues EIG of the matrix A. The
%   desired tolerance is EPS, MAXCOUNT is the maximum number of iterations
%   allowed, and I is how many iterations were conducted.
%1: Check that the maximum number of iterations has not been reached.
%2: Calculate the QR decomposition of the kth iterative matrix A.
%3: Set the (k+1)th iterative matrix as being R*Q.
%4: 
%   A: Check that the maximum number of iterations has not been reached.
%   B: Find y=inv(A)*x and normalize y to get x.
%   C: If both the inf-norm of the difference between two adjacent
%      eigenvectors and the error in the eigenvalue estimate is less than
%      EPS, quit and return the smallest norm of an eigenvalue and its
%      eigenvector.
%   D: Otherwise, set the old eigenvector iteration as being the new guess.
function [EIG,i]=QREig(A,MAX,TOL)
for i=1:MAX
    [Q,R]=qr(A);
    A=R*Q;
    EIG=diag(A);
    if max(max(abs(A-diag(EIG))))<TOL
        return;
    end
end