function Z = lusolZ(S,trans,options)

%        Z = lusolZ(S,trans,options);
% forms a nullspace operator Z from the sparse matrix S such that S'*Z = 0,
% using LU factors of either S or S' (trans = 0 or 1 respectively).
% The structure "options" contains parameters controlling how LUSOL
% computes the LU factors.
%
% Z can later be used to compute products of the form w = Z*v or s = Z'*t
% for given vectors v and t, via code of the form
%   w = lusolZv(Z,v);
%   s = lusolZt(Z,t);
%
% This function scales S using gmscale.m and calls LUSOL
% via code of the form
%   [cscale,rscale] = gmscale(S,iprint,scltol);
%   Cinv   = diag(sparse(1./cscale));
%   Rinv   = diag(sparse(1./rscale));
%   A      = Rinv*S*Cinv;
%   if trans, A = A'; end
%   Z.lu      = lusol_obj(A,options);
% The "options" structure is defined in https://github.com/nwh/lusol.

% 01 Mar 2016: First version of lusolZ.

 iprint = 0;
 scltol = 0.9;
 [cscale,rscale] = gmscale(S,iprint,scltol);
 Cinv   = diag(sparse(1./cscale));
 Rinv   = diag(sparse(1./rscale));
 A      = Rinv*S*Cinv;
 C      = diag(sparse(cscale));
 R      = diag(sparse(rscale));

 if trans, A = A'; end

 Z.lu      = lusol_obj(A,options);
 Z.trans   = trans;
 Z.C       = C;
 Z.R       = R;
 Z.Cinv    = Cinv;
 Z.Rinv    = Rinv;

 if trans, [Z.Utri,~,~] = lusolU(Z.lu); end
