function testlusolZ( trans )

%        testlusolZ( trans );
% tests w = lusolZv(Z,v) and s = lusolZt(Z,t) for a range of matrices S
% stored in ../data/.
%
% For each of "ntest" random vectors v, the product w = Z*v is formed.
% The vectors z = S'*w should be zero.
% errZv returns the average norm(z,inf).
% timZv returns the average time for computing w.
%
% For each of "ntest" random vectors u, the product t = S*u is formed.
% The vectors s = Z'*t should be zero.    
% errZt returns the average norm(s,inf).
% timZt returns the average time for computing s.

% 01 Mar 2016: First version of testlusolZ.

 models =  { ...
     'iIT341 '; ...
     'iAF692 '; ...
     'iSB619 '; ...
     'iJN746 '; ...
     'iJR904 '; ...
     'iNJ661 '; ...
     'iND750 '; ...
     'iAF1260'; ...
     'Recon1 '; ...
     'ThMa   '; ...
     'GlcAer '; ...
     };
 nmodel         = length(models);
 factol         = 2.0;
 options        = lusol();
 options.pivot  = 'TPP';
 options.Ltol1  = factol;

 fprintf('\n   factol = %7.2f', factol)
 if ~trans
     fprintf('     Scaled S = LU     %s\n\n', options.pivot)
 else
     fprintf('     Scaled S'' = LU     %s\n\n', options.pivot)
 end
 fprintf('   model')
 fprintf(' |   rank(S)')
 fprintf(' |    nnz(L)    nnz(U)')
 fprintf(' |           errZv      timZv')
 fprintf(' |           errZt      timZt\n')
 fprintf('--------')
 fprintf('------------')
 fprintf('----------------------')
 fprintf('-----------------------------')
 fprintf('-----------------------------\n')
 
 for i = 1:nmodel
     name  = models{i};
     eval(['load data/' name])
     S     = eval([ 'model_' name '.S']);
     [m,n] = size(S);
   
     % Factorize S or S'
     options.nzinit = round(20*nnz(S)/factol);

     Z     = lusolZ(S,trans,options);
     r     = Z.lu.rank;
     nL    = Z.lu.stats.lenL;
     nU    = Z.lu.stats.lenU;
     nZ    = m - r;

     ntest = 50;
     errZv = 0;
     errZt = 0;
    
     tic;
     for j = 1:ntest
         v     = rand(nZ,1) - 0.5;
         w     = lusolZv(Z,v);
         z     = S'*w;
         errZv = errZv + norm(z,inf);
     end
     errZv = errZv / ntest;
     timZv = toc / ntest;
     
     tic;
     for j = 1:ntest
         u     = rand(n,1) - 0.5;
         t     = S*u;
         s     = lusolZt(Z,t);
         errZt = errZt + norm(s,inf);
     end
     errZt = errZt / ntest;
     timZt = toc / ntest;
     fprintf(' %7s', name)
     fprintf(' |   %7d', r)
     fprintf(' |   %7d   %7d', nL, nU)
     fprintf(' |    %8.6e   %8.4f'  , errZv, timZv)
     fprintf(' |    %8.6e   %8.4f\n', errZt, timZt)
 end