function [R,H,rowperm,C] = sqrB(A,nsteps,nelim,nstk,nle,Pr,B)
% SQRB  Factorization routine in SQR.
%       @(#)sqrB.m Version 1.3 11/24/97
%       Pontus Matstoms, Linkoping University.
%       Mikael Lundquist, Linköping University.
%       e-mail: pomat@mai.liu.se, milun@mai.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.

% PSTK           Stack of row indexes. Used when computing Q
% rowidx         Row index

[m,n]=size(A);
% Check the input
computeB=(nargin==7);
computeH=(nargout>=3);

if computeB, 
   B=B(Pr,:); 
   [Bm,Bn]=size(B); 
   C=zeros(n,Bn);
%else, 
%  B=0;  
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;  

if computeH,

  rowperm=zeros(m,1);
  last=m;
end

% --- 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;
    if Am>0,
      [ignore,col]=find(A((lmr+1):(lmr+Am),:)); 
      colflag(col)=ones(size(col));
    end

% - 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 ...

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

          if computeB,                     % Right-hand side front.
             BFNT=zeros(FNTm,Bn);
             BFNT(1:Am,:)=B((lmr+1):(lmr+Am),:);
	   end
	   if computeH,
	     rowidx=(lmr+1):(lmr+Am); % Rowindex of the rows from A
	   end
	     
	else
	  if computeH,
	    Fm=0;           %ML 970910
	    rowidx=[];
	  end
	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 computeB,
	      BFNT((Fm+1):(Fm+udm),:)=BSTK(:,(Bp-udm+1):Bp)';
	    end
	    if computeH,
	      rowidx=[rowidx PSTK(Bp-udm+1:Bp)];
	      %	    rowidx=[rowidx PSTK(Rp-PSTKm(spnm)+1:Rp)];
	      %	    Rp=Rp-PSTKm(spnm);
	    end
	    Bp=Bp-udm;
	    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 computeH,
      [Q,RR]=qr(FNT);
      RBF=triu(RR);
      H(iass).frontH=sparse(Q);
      H(iass).p=rowidx';
      %	 RR=qr(FNT);
      %	 RBF=triu(RR);	
      %  H(iass).frontH=sparse(tril(RR,0));
    
      if computeB
	Ctemp=Q'*BFNT
	C((rp+1):(rp+numorg),:)=Ctemp(1:numorg);
      end

    else
      if computeB,     % 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
      end
         
% - Move the first rows of frontals R to the final R.

    R(FNTcols,(rp+1):(rp+numorg))=RBF(1:numorg,1:RFn)';
    if computeH,
      rowperm((rp+1):(rp+numorg))=rowidx(1:numorg);
    end
    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 computeB,
	 BSTK(:,(Bp+1):(Bp+rdim-numorg))=RBF((numorg+1):rdim,(RFn+1):(RFn+Bn))';
       end
       if computeH,
	 PSTK((Bp+1):(Bp+rdim-numorg))=rowidx((numorg+1):rdim);
       end
       Bp=Bp+rdim-numorg; 
     else
       STKm(spnm)=0;
       STKn(spnm)=0;
    end

    if computeH,
      if FNTm>rdim,
	rowperm(last-(FNTm-rdim)+1:last)=rowidx(rdim+1:FNTm);
	last=last-(FNTm-rdim);
      end
    end
    
    lmr=lmr+Am;

  end

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

R=R';

if ~computeH & computeB
  H=C;
  clear C
end