% Test input/output alternatives for the lurp code from
% http://www.math.sjsu.edu/~foster/rankrevealingcode.html.
% all the statistic, except for the last three, should
% be order of machine epsilon

disp(' ')
disp('Test input/output alternatives for the lurp code from')
disp('http://www.math.sjsu.edu/~foster/rankrevealingcode.html.')
disp('All the error statistics should be O(eps), except for the')
disp('last three (where the element growth tolerance = 2000).')
disp(' ')
ev =[];
m = 800; n =500;
A=randn(m,n);
B=randn(n,m);
C=randn(n,n);
format compact
% no output
lurp(A);lurp(B);lurp(C);
% one output
LA=lurp(A);LB=lurp(B);LC=lurp(C);
% test if one output and four outputs consistent
[L,U,P,Q]=lurp(A);
e = norm(tril(L,-1)-tril(LA,-1),1);
ev = [ev, e];
e = norm(triu(U)-triu(LA(1:n,1:n)),1 );
disp(e)
ev = [ev, e];
% permutations returned as vectors
[L,U,P,Q]=lurp(A,-1);
e=norm(A(P,Q)-L*U,1)/norm(A,1);
disp(e)
ev = [ev, e];
% set the block size to 32
[L,U,P,Q]=lurp(A,32);
e=norm(P*A*Q-L*U,1)/norm(A,1);
disp(e)
ev = [ev, e];
% set the block size to 213
[L,U,P,Q]=lurp(A,213);
e=norm(P*A*Q-L*U,1)/norm(A,1);
disp(e)
ev = [ev, e];
% test two outputs
[LA,UA]=lurp(A);[LB,UB]=lurp(B);[LC,UC]=lurp(C);
e = norm(A-LA*UA,1)/norm(A,1); 
disp(e)
ev = [ev, e];
e = norm(B-LB*UB,1)/norm(B,1);
disp(e)
ev = [ev, e];
e = norm(C-LC*UC,1)/norm(C,1);
disp(e)
ev = [ev, e];
% test tol = -2 and trep = .5
[LA,UA]=lurp(A,0,.5,-2);[LB,UB]=lurp(B,0,.5,-2);[LC,UC]=lurp(C,0,.5,-2);
e = norm(A-LA*UA,1)/norm(A,1);
disp(e)
ev = [ev, e];
e = norm(B-LB*UB,1)/norm(B,1);
disp(e)
ev = [ev, e];
e = norm(C-LC*UC,1)/norm(C,1);
disp(e)
ev = [ev, e];
% test tol = -2000, four outputs, permutations as matrices
[LA,UA,PA,QA]=lurp(A,0,0,-2000);
[LB,UB,PB,QB]=lurp(B,0,0,-2000);
[LC,UC,PC,QC]=lurp(C,0,0,-2000);
e = norm(PA*A*QA-LA*UA,1)/norm(A,1);
disp(e)
ev = [ev, e];
e = norm(PB*B*QB-LB*UB,1)/norm(B,1);
disp(e)
ev = [ev, e];
e = norm(PC*C*QC-LC*UC,1)/norm(C,1);
disp(e)
ev = [ev, e];
% test tol = -2000, four outputs, permutations as matrices, block size 1
[LA,UA,PA,QA]=lurp(A,1,0,-2000);
[LB,UB,PB,QB]=lurp(B,1,0,-2000);
[LC,UC,PC,QC]=lurp(C,1,0,-2000);
e = norm(PA*A*QA-LA*UA,1)/norm(A,1);
disp(e)
ev = [ev, e];
e = norm(PB*B*QB-LB*UB,1)/norm(B,1);
disp(e)
ev = [ev, e];
e = norm(PC*C*QC-LC*UC,1)/norm(C,1);
disp(e)
ev = [ev, e];
% test permutations output as vectors (again)
[LA,UA,PA,QA]=lurp(A,-1);
[LB,UB,PB,QB]=lurp(B,-1);
[LC,UC,PC,QC]=lurp(C,-1);
e = norm(A(PA,QA)-LA*UA,1)/norm(A,1);
disp(e)
ev = [ev, e];
e = norm(B(PB,QB)-LB*UB,1)/norm(B,1);
disp(e)
ev = [ev, e];
e = norm(C(PC,QC)-LC*UC,1)/norm(C,1);
disp(e)
ev = [ev, e];
% test leading dimension lda = 1000
[LA,UA,PA,QA]=lurp(A,-1,0,0,1000);
[LB,UB,PB,QB]=lurp(B,-1,0,0,1000);
[LC,UC,PC,QC]=lurp(C,-1,0,0,1000);
e = norm(A(PA,QA)-LA*UA,1)/norm(A,1);
disp(e)
ev = [ev, e];
e = norm(B(PB,QB)-LB*UB,1)/norm(B,1);
disp(e)
ev = [ev, e];
e = norm(C(PC,QC)-LC*UC,1)/norm(C,1);
disp(e)
ev = [ev, e];
% test minimum workspace (lwork = 1)
[LA,UA,PA,QA]=lurp(A,-1,0,0,0,0,1);
[LB,UB,PB,QB]=lurp(B,-1,0,0,0,0,1);
[LC,UC,PC,QC]=lurp(C,-1,0,0,0,0,1);
e = norm(A(PA,QA)-LA*UA,1)/norm(A,1);
disp(e)
ev = [ev, e];
e = norm(B(PB,QB)-LB*UB,1)/norm(B,1);
disp(e)
ev = [ev, e];
e = norm(C(PC,QC)-LC*UC,1)/norm(C,1);
disp(e)
ev = [ev, e];
% test relative element growth, tol = 2000 -- potential instabilities
[LA,UA]=lurp(A,0,0,2000);
[LB,UB]=lurp(B,0,0,2000);
[LC,UC]=lurp(C,0,0,2000);
e = norm(A-LA*UA,1)/norm(A,1);
disp(e)
ev = [ev, e];
e = norm(B-LB*UB,1)/norm(B,1);
disp(e)
ev = [ev, e];
e = norm(C-LC*UC,1)/norm(C,1);
disp(e)
ev = [ev, e];
disp(' ')
if ( max(ev(1:end-3)) <= max(m,n)*eps )
    disp('All tests passed that should pass.')
else
    disp('One or more routines failed a test.')
end