% This Matlab m file tests the code available from
% http://www.math.sjsu.edu/~foster/rankrevealingcode.html
% all the statistic, except for the lunp statistic, should
% be order of machine epsilon.

disp(' ')
disp('Test the code from')
disp('http://www.math.sjsu.edu/~foster/rankrevealingcode.html.')
disp('All the error statistics should be O(eps), except')
disp('for the first one (Gaussian elimination with no pivoting).')
disp(' ')
clear ev
m = 500;
n= 500;
A=randn(m,n);
[L,U]=lunp(A);
e = norm(A-L*U,1)/norm(A,1);
ev(1) = e;
disp(['lunp  : ',num2str(e)])
[L,U]=lucp(A);
e = norm(A-L*U,1)/norm(A,1);
ev(2) = e;
disp(['lucp  : ',num2str(e)])
[L,U,P,Q]=lupan3(A);
e = norm(P*A*Q-L*U,1)/norm(A,1);
ev(3) = e;
disp(['lupan3: ',num2str(e)])
[L,U]=lupp(A);
e = norm(A-L*U,1)/norm(A,1);
ev(4) = e;
disp(['lupp  : ',num2str(e)])
[L,U]=lurp(A);
e = norm(A-L*U,1)/norm(A,1);
ev(5) = e;
disp(['lurp  : ',num2str(e)])
[L,U]=lutp(A);
e = norm(A-L*U,1)/norm(A,1);
ev(6) = e;
disp(['lutp  : ',num2str(e)])
R =qrnp(A);
e = norm(A'*A-R'*R,1)/norm(A'*A,1);
ev(7) = e;
disp(['qrnp  : ',num2str(e)])
R=qrp(A);
e = norm(A'*A-R'*R,1)/norm(A'*A,1);
ev(8) = e;
disp(['qrp   : ',num2str(e)])
R=qrtp(A);
e = norm(A'*A-R'*R,1)/norm(A'*A,1);
ev(9) = e;
disp(['qrtp  : ',num2str(e)])
nb = 0;trep=0;rcond1=0;lda=0;job=3;
[Q, R, P,time, rank1] = qrxp(A,nb,trep,rcond1,lda,job);
e=norm(A*P-Q*R,1)/norm(A,1);
ev(10) = e;
disp(['qrxp  : ',num2str(e)])
[Q, R, P,time, rank1] = qryp(A,nb,trep,rcond1,lda,job);
e=norm(A*P-Q*R,1)/norm(A,1);
ev(11) = e;
disp(['qryp  : ',num2str(e)])
s = svdmex(A);
s2 = svd(A);
e = norm(s - s2) / norm(s);
ev(12) = e;
disp(['svdmex: ',num2str(e)])
disp(' ')
if ( max(ev(2:end)) <= max(m,n)*eps )
    disp('All tests passed that should pass.')
else
    disp('One or more routines failed a test.')
end
