% P4
% Checks out the function [un,vn,sigman] = triple(L,U,P)
  clc
% Might want to reduce the following during debugging...
  n = 200;         % Problem dimension.
  nTrials = 100;   % Number of random trials for each kappa. 
  nRepeat = 1000;  % Repetition factor for benchmarking
  
  for kappa = logspace(3,12,4)
%    Check performance on examples with condition kappa...
     u_sum = 0; v_sum = 0; s_sum = 0; s_minmax = 0; uv_error = 0;
     for eg = 1:nTrials
%       Generate a random A with 2-norm condition kappa...
        U = orth(randn(n,n)); 
        V = orth(randn(n,n));
        s = rand(n,1); s = -sort(-s); s(1) = 1; s(n) = 1/kappa; 
        A = U*diag(s)*V';
%       Apply and assess the condition estimator...
        [L,U,P] = lu(A);
        [un,vn,sigman] = Triple(L,U,P);
        u_sum = u_sum + abs(un'*U(:,n));
        v_sum = v_sum + abs(vn'*V(:,n));
        s_sum = s_sum + sigman;
        s_minmax = max(s_minmax,sigman);
        uv_error = max(uv_error,abs(norm(un)-1)+abs(norm(vn)-1));
     end
     
 %   Estimate the efficiency of Triple as a multiple of the
 %   time required to solve an n-by-n triangular system...
     b = rand(n,1);
     tic
     for k=1:nRepeat
        x = U\b;
     end
     T1 = toc;
     tic
     for k=1:nRepeat
        [un,vn,sigman] = Triple(L,U,P);
     end
     T2 = toc;
     
 %   Display results for this kappa...
     fprintf('\n\nTrue sigman       =  %10.6e\n',1/kappa)
     fprintf('Average sigman    =  %10.6e\n',s_sum/nTrials)
     fprintf('Worst sigman      =  %10.6e\n',s_minmax)
     fprintf('Average u-cosine  =%10.6f\n',u_sum/nTrials)
     fprintf('Average v-cosine  =%10.6f\n',v_sum/nTrials)
     fprintf('uv_error          =  %10.6e\n',uv_error)
     fprintf('Efficiency Metric = %6.3f\n',T2/T1)
  end