function [R,Pc,C,H] = sqr(A,B)
%SQR	Sparse orthogonal-triangular decomposition.
%       %Z%%M% Version %I% %G%
%       Pontus Matstoms, Linkoping University.
%       Mikael Lundquist, Linköpings University
%       e-mail: pomat@mai.liu.se, milun@mai.liu.se
%
%       [R,p]=sqr(A) computes the upper triangular matrix R in the
%       QR factorization,
%                          Q'(AP)=[R 0]'.
%       The permutation matrix P is represented by the integer vector p.
%
%       [R,p,H]=sqr(A) computes also a representation of the orthogonal
%       matrix Q. This representation is returned in H. To use H on a
%       vector use the routines appH and appHT.
%
%       [R,p,C]=sqr(A,B) does also compute the n first components of C=Q'B,
%       where n is the column dimension of A. The matrix is here ordered
%       by rows in the same way as the multifrontal scheme orders the rows
%       of A.
%
%       If the GLOBAL variable sqr_nemin exists, it is chosen as amalgamation
%       parameter in sqrA.

% Call the analysis routine.

global sqr_nemin
if isempty(sqr_nemin),
   [nsteps,nelim,nstk,nle,Pc,Pr] = sqrA(A);
else
   disp(['nemin is now ',num2str(sqr_nemin)])
   [nsteps,nelim,nstk,nle,Pc,Pr] = sqrA(A,sqr_nemin);
end

% Order A by rows (Pr) and columns (Pc).

A=A(Pr,Pc);

% Call the factorization routine.

if nargin == 1,            % Only factorization
  if nargout < 3,          % Compute only R
    R = sqrB(A,nsteps,nelim,nstk,nle,Pr);
  elseif nargout ==3,      % Compute representation of Q
    [R,HH,rowperm] = sqrB(A,nsteps,nelim,nstk,nle,Pr);
    C.Pr=Pr;
    C.rowperm=rowperm;
    C.H=HH;
  else
    error('Wrong number of output parameters, see help info');    
  end
else                       % Righthand side
  if nargout <=3,           % Only righthand side
    [R,C] = sqrB(A,nsteps,nelim,nstk,nle,Pr,B);
  else
    error('Wrong number of output parameters, see help info');    
  end
end