function hw4_driver
% HW4_DRIVER  Driver code for HW4 programming assignment.
%    Note: We use the lu function in MATLAB, but use our own
%          implementation for everything else. 
%          DO NOT USE the \ OPERATOR.

mindim = 2;
maxdim = 12;

% Compute condition numbers for both methods
cycles = 500;

conds = zeros (2, maxdim - mindim + 1);
times = zeros (2, maxdim - mindim + 1);

for i = mindim:maxdim
    A = gen_hilbert( i);

    tic;
    % Run cond1_estimate for many times for proper timing
    for j = 1 : cycles
        conds(1, i-1) = cond1_estimate(A);
    end
    times(1,i-1) = toc/cycles;

    % Run cond1_estimate for many times for proper timing
    tic;
    for j = 1 : cycles
        conds(2, i-1) = cond1_explicit(A);
    end
    times(2,i-1) = toc/cycles;
end

plotresults( conds, times);
end

function kappa_est = cond1_estimate( A )
%COND1_estimate   Estimate condition number of a given matrix in 1-norm.

m = size(A, 1);

% Step 1: LU factorization of A ...
[L, U, P] = lu( A );

% Step 2: Calculate U'*v = c for unknown v by changing c entry by entry.
% Implement this step by yourself.
v = zeros(m, 1);
v(1) = abs ( 1 / U(1, 1) );

for i = 2:m
    assert(false);
end

% Step 2: Calculate (L'*P)*y = v for unknown y
% Implement this step using back_substitution.
assert(false);

% Step 3: Solve A z = y using PA=LU
% Implement this step using forward_substitution and/or back_substitution.
assert(false);

% Return the estimated condition number ...
norm1_Ainv_est = norm1(z) / norm1(y);
kappa_est = norm1( A) * norm1_Ainv_est;

end

function kappa = cond1_explicit( A )
%COND1_EXPLICIT   Compute condition number of a given matrix A by
%     explicitly computing the inverse of A using LU factorization
%     with partial pivoting.

[L, U, P] = lu( A);

% Inverse of A
Ainv = zeros(size(A));

% Implement the procedure for computing inverse of A and save into Ainv
%   using L, U, forward_substitution, and back_substitution here.
assert(false);

% Return the computed condition number.
kappa = norm1( A) * norm1( Ainv);
end

function n = norm1(A)
% NORM1   Compute 1-norm of a given matrix

n = 0;
for j=1:size(A,2)
    n = max( n, sum( abs(A(:,j))));
end
end

function v = forward_substitution (L, c)
% Forward substitution for lower triangular matrix.

v = zeros (size(c));

for i = 1:size(c, 1)
    v(i) = c(i) - L(i, 1:i-1) * v(1:i-1);
    v(i) = v(i) / L(i, i);
end
end

function v = back_substitution (U, c)
% Back substitution for upper triangular matrix.

m = size(c, 1);
v = zeros(size(c));

for i = m:-1:1
    v(i) = c(i) - U(i, i+1:m) * v(i+1:m);
    v(i) = v(i) / U(i, i);
end
end

function H = gen_hilbert(n)
% GEN_HILBERT      Generate Hilbert matrix.
% H = GEN_HILBERT(N), where H is n-by-n Hilbert matrix.

H = zeros(n, n);

for i = 1:n
    for j = 1:n
        H(i,j) = 1/(i+j-1);
    end
end
end

function plotresults( conds, times)

mindim = 2; maxdim = size(conds,2)+1;

% Plot the estimate and computed results ...
figure (1);
hold off;
semilogy(mindim:maxdim, conds(1, :),'-bs', mindim:maxdim, conds(2, :),'--rx');
title('Condition numbers of Hilbert matrices','fontsize',12,'fontweight','b');
xlabel('Dimension of matrices','fontsize',10,'fontweight','b');
ylabel('Condition numbers','fontsize',10,'fontweight','b');
legend('Estimate','Computed','Location','NorthWest');
print -depsc conds.eps

% Plot the running times for estimate and computed algorithms ...
figure (2);
hold off;
plot(mindim:maxdim, times(1, :),'-bs',mindim:maxdim, times(2, :),'--rx');
title('Running times of algorithms','fontsize',12,'fontweight','b');
xlabel('Dimensions of matrices','fontsize',10,'fontweight','b');
ylabel('Running times (seconds)','fontsize',10,'fontweight','b');

legend('Estimate','Computed','Location','NorthWest');
print -depsc errors.eps
end