function [R,C] = sqrB(A,nsteps,nelim,nstk,nle,Pr,B)
% SQRB  Factorization routine in SQR.
%       @(#)sqrB.m Version 1.11 3/1/93
%       Pontus Matstoms, University of Linkoping.
%       e-mail: pomat@math.liu.se
%
%       [R,C] = sqrB(A,nsteps,nelim,nstk,nle) computes the upper triangular 
%       factor R in the QR factorization of the sparse m-by-n matrix A. If a 
%       matrix B is passed to the routine, [R,C] = sqrB(A,nsteps,nelim,nstk,nle,B),
%       then the first n components of Q'B are also computed.
%

% nsteps         The number of nodes/supernodes in the elimination tree.
% nelim()        The number of columns to eliminate in the current step.
% nstk()         The number of stacked elements to merge in the current step.
% R(,)           The upper triangular matrix R.
% STK()          The stack of update matrices.
% STKcols        Column indices of stacked update matrices.
% STKm           Row dimension of stacked update matrices.
% STKn           Column dimension of stacked update matrices.
% BSTK           Right-hand side stack.

[m,n]=size(A);

% Check the input

if nargin == 6, 
   B=0; 
else, 
   B=B(Pr,:); 
   [Bm,Bn]=size(B); 
   C=zeros(n,Bn);
end

% From the vector nle (nle(i)=number of rows with leading entry i) a new vector
% NLE (NLE(i)=number of rows from A to include in the i:th frontal matrix) is
% computed.

NLE=[0 cumsum(nle)];
a=cumsum([1 nelim(1:(nsteps-1))])';
b=cumsum(nelim)+1;
NLE=NLE(b)-NLE(a);

R=spalloc(n,n,sum(symbfact(A,'col')));

sp=1;            % Pointer for STK
spc=1;           % Pointer for STKcols
spnm=0;          % Pointer for STKm and STKn
Bp=0;            % Pointer for BSTK

rp=0;                   
lmr=0;  

% --- Main loop ...

for iass=1:nsteps,

% -- Compute the dimension and the global column indices of the front.

    numorg=nelim(iass);
    numstk=nstk(iass);

    colflag=zeros(1,n);               % Global column indices in the frontal matrix.
    
% - Integer data associated with the contribution from A.
    
    Am=NLE(iass);
    FNTm=Am;
    [ignore,col]=find(A((lmr+1):(lmr+Am),:)); 
    colflag(col)=ones(size(col));

% - Integer data associated with the stack contribution.

    if ( numstk > 0 ) & ( spc>1 ),
       stkm=sum(STKm((spnm-numstk+1):spnm));
       FNTm=FNTm+stkm;
       ncol=sum(STKn((spnm-numstk+1):(spnm)));
       col=STKcols((spc-ncol):(spc-1));
       colflag(col)=ones(size(col));
    else
       stkm=0;
    end

    FNTcols=find(colflag);
    FNTn=size(FNTcols,2);
    colflag(FNTcols)=1:FNTn;          % colflag maps global to local column
                                      % indices for the frontal matrix.
                                      % local=colflag(global).

% -- Move reals to the frontal matrix.

    if FNTn > 0,
       FNT=zeros(FNTm,FNTn);

% - From A ...

       FNT(1:Am,:)=A((lmr+1):(lmr+Am),FNTcols); 
       Fm=Am;

       if nargin == 7,                % We also build the right-hand side front.
          BFNT=zeros(FNTm,Bn);
          BFNT(1:Am,:)=B((lmr+1):(lmr+Am),:);
       end

% - ... and from the stack ...

       for s=1:numstk,
         if STKm(spnm)>0,
            udm=STKm(spnm);
            udn=STKn(spnm);
            col=colflag(STKcols((spc-udn):(spc-1)));
            FNT((Fm+1):(Fm+udm),col)=reshape(STK((sp-udm*udn):(sp-1)),udm,udn);
            if nargin == 7,
               BFNT((Fm+1):(Fm+udm),:)=BSTK(:,(Bp-udm+1):Bp)';
               Bp=Bp-udm;               
            end
            Fm=Fm+udm;
            spc=spc-udn;
            sp=sp-udm*udn;
         end
         spnm=spnm-1;
       end
    else
       spnm=spnm-numstk;
    end

    numorg=min(numorg,FNTm); % Rank deficient problem may have FNTm < numorg.

% -- Factorization

    RFm=min(FNTm,FNTn);
    RFn=FNTn;
    rdim=min(RFm,RFn);

    if nargin == 7,     % A right-hand side matrix is passed to the routine.

% -    Factorize the frontal matrix and dispose computed elements of Q'B.

       RBF=triu(qr([FNT BFNT]));
       C((rp+1):(rp+numorg),:)=RBF(1:numorg,(RFn+1):(RFn+Bn)); 

    else
     
% -    Factorize the frontal matrix.

       RBF=triu(qr(FNT));
    end

% - Move the first rows of frontals R to the final R.

    R(FNTcols,(rp+1):(rp+numorg))=RBF(1:numorg,1:RFn)';

    rp=rp+numorg;

% --- Move the remaining rows to the stack.

    spnm=spnm+1;
    if rdim > numorg,
       STKm(spnm)=rdim-numorg;
       STKn(spnm)=RFn-numorg;
       STKcols(spc:(spc+STKn(spnm)-1))=FNTcols((numorg+1):RFn);
       spc=spc+STKn(spnm);
       us=STKm(spnm)*STKn(spnm);                % Size of update matrix
       STK(sp:(sp+us-1))=reshape(RBF((numorg+1):rdim,(numorg+1):RFn),1,us);
       sp=sp+us;
       if nargin == 7,
          BSTK(:,(Bp+1):(Bp+rdim-numorg))=RBF((numorg+1):rdim,(RFn+1):(RFn+Bn))';
          Bp=Bp+rdim-numorg;
       end
    else
       STKm(spnm)=0;
       STKn(spnm)=0;
    end

    lmr=lmr+Am;

end

% Up to now the transpose of R has been stored ...

R=R';

