function P3()
% Checks out SparseCond

clc
randn('seed',0)
s = 1;
LU_Indices = [156 246 258 1408 ];
disp(' UF(i)    n   kappaCVL     kappaYOU')
disp('-----------------------------------------')
for i=LU_Indices
    Problem = UFget(i);
    A = Problem.A;
    [n,n] = size(A);
    y0 = randn(n,1);
    kappaCVL = SparseCondCVL(A,y0,s);
    kappaYOU = SparseCond(A,y0,s);
    fprintf('%5d  %5d  %10.3e  %10.3e\n',i,n,kappaCVL,kappaYOU)
end

function kappa = SparseCondCVL( A,y0,s )
% A is a nonsingular n-by-n matrix in sparse format
% y0 is an n-by-1 vector
% s  is a positive
% kappa is an estimate of A's 1-norm condition based on s updates
% of the form y = (A'A)\y with y = y0 at the start
    
% Compute P*inv(D)*A*Q = LU where P = I(:,p), Q = I(q,:), and D is
% diagonal in sparse format.

[L,U,p,q,D] = lu(A,[.1 .001],'vector');

% We have this: D(:,p)\A(:,q) = L*U, i.e., inv(DP')(AQ') = LU
% where P = I(:,p) and Q = I(:,q). Thus, P*inv(D)*A*Q' = LU

% Review: v = v(p) means v is overwritten by P*v.
%         v(p) = v means v is overwritten by P'*v
%         Matrix(:,p) same as Matrix*P'

% Note that
%            A  = D*P'*L*U*Q     so inv(A)  = Q'*inv(U)*inv(L)*P*inv(D)
%            A' = Q*U'*L'*P*D    so inv(A') = inv(D)*P'*inv(L')*inv(U')*Q
% This doesn't mean we compute inverses, it just clarifies the sequence
% of solve operations.

% Initializations...
normA = norm(A,1);
y = y0;
n = length(y0);
for k=1:s
    % Solve A'*d = y
    d    = y(q);
    d    = U'\d;
    d    = L'\d;
    d(p) = d;
    d    = D\d;
    % Normalize to prevent huge growth for larger values of s...
    d    = d/norm(d,1);
    % Solve Ay = d
    y = D\d;
    y = y(p);
    y = L\y;
    y = U\y;
    y(q) = y;
    %
    kappa = normA*norm(y,1);
end
 

% -3 if you compute LU twice
% -4 if LU inside loop
% -1 if you messed up the the permutation manipulations