function Z = spqrZ(S,trans)

%        Z = spqrZ(S,trans);
% forms a nullspace operator Z from the sparse matrix S such that S'*Z = 0,
% using QR factors of either S or S' (trans = 0 or 1 respectively).
% The QR factors are computed by SPQR from SuiteSparseQR (see below).
%
% 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 = spqrZv(Z,v);
%   s = spqrZt(Z,t);
%
% This function scales S using gmscale.m and calls SPQR
% 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
%   [Q,R,P,info] = spqr(A,struct('Q','Householder'));
% The spqr code is documented in
% http://faculty.cse.tamu.edu/davis/suitesparse.html.

% 01 Mar 2016: First version of spqrZ.

 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;
 Cs     = diag(sparse(cscale));
 Rs     = diag(sparse(rscale));

 if trans, A = A'; end

 [Q,R,P,info] = spqr(A,struct('Q','Householder'));
 Z.A       = A;
 Z.Q       = Q;
 Z.R       = R;
 Z.P       = P;
 Z.info    = info;
 Z.trans   = trans;
 Z.Cs      = Cs;
 Z.Rs      = Rs;
 Z.Cinv    = Cinv;
 Z.Rinv    = Rinv;