function [R,Pc,C] = sqr(A,B)
%SQR	Sparse orthogonal-triangular decomposition.
%       @(#)sqr.m Version 1.10 9/17/93
%       Pontus Matstoms, Linkoping University.
%       e-mail: pomat@math.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,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 exist('sqr_nemin') == 0,
   [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,
   R = sqrB(A,nsteps,nelim,nstk,nle,Pr);
   C = 0;
else
   [R,C] = sqrB(A,nsteps,nelim,nstk,nle,Pr,B);
end

