[q,n,nrda,diag,irank,div,td]=dohtrl(q,n,nrda,diag,irank,div,td);

function [q,n,nrda,diag,irank,div,td]=dohtrl(q,n,nrda,diag,irank,div,td);
%***BEGIN PROLOGUE  DOHTRL
%***SUBSIDIARY
%***PURPOSE  Subsidiary to DBVSUP and DSUDS
%***LIBRARY   SLATEC
%***TYPE      doubleprecision (OHTROL-S, DOHTRL-D)
%***AUTHOR  Watts, H. A., (SNLA)
%***DESCRIPTION
%
%     For a rank deficient problem, additional orthogonal
%     HOUSEHOLDER transformations are applied to the left side
%     of Q to further reduce the triangular form.
%     Thus, after application of the routines DORTHR and DOHTRL
%     to the original matrix, the result is a nonsingular
%     triangular matrix while the remainder of the matrix
%     has been zeroed out.
%
%***SEE ALSO  DBVSUP, DSUDS
%***ROUTINES CALLED  DDOT
%***REVISION HISTORY  (YYMMDD)
%   750601  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900328  Added TYPE section.  (WRB)
%   910722  Updated AUTHOR section.  (ALS)
%***end PROLOGUE  DOHTRL
persistent dd diagk irp j k kir kirm l nmir qs sig sqd tdv ; 

if isempty(irp), irp=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(kir), kir=0; end;
if isempty(kirm), kirm=0; end;
if isempty(l), l=0; end;
if isempty(nmir), nmir=0; end;
if isempty(dd), dd=0; end;
diag_shape=size(diag);diag=reshape(diag,1,[]);
if isempty(diagk), diagk=0; end;
div_shape=size(div);div=reshape(div,1,[]);
q_shape=size(q);q=reshape([q(:).',zeros(1,ceil(numel(q)./prod([nrda])).*prod([nrda])-numel(q))],nrda,[]);
if isempty(qs), qs=0; end;
if isempty(sig), sig=0; end;
if isempty(sqd), sqd=0; end;
td_shape=size(td);td=reshape(td,1,[]);
if isempty(tdv), tdv=0; end;
%***FIRST EXECUTABLE STATEMENT  DOHTRL
nmir = fix(n - irank);
irp = fix(irank + 1);
for k = 1 : irank;
kir = fix(irp - k);
diagk = diag(kir);
sig =(diagk.*diagk) + ddot(nmir,q(sub2ind(size(q),irp,kir):end),1,q(sub2ind(size(q),irp,kir):end),1);
dd = (abs(sqrt(sig)).*sign(-diagk));
div(kir) = dd;
tdv = diagk - dd;
td(kir) = tdv;
if( k~=irank )
kirm = fix(kir - 1);
sqd = dd.*diagk - sig;
for j = 1 : kirm;
qs =((tdv.*q(kir,j))+ddot(nmir,q(sub2ind(size(q),irp,j):end),1,q(sub2ind(size(q),irp,kir):end),1))./sqd;
q(kir,j) = q(kir,j) + qs.*tdv;
for l = irp : n;
q(l,j) = q(l,j) + qs.*q(l,kir);
end; l = fix(n+1);
end; j = fix(kirm+1);
end;
end; k = fix(irank+1);
diag_shape=zeros(diag_shape);diag_shape(:)=diag(1:numel(diag_shape));diag=diag_shape;
div_shape=zeros(div_shape);div_shape(:)=div(1:numel(div_shape));div=div_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
td_shape=zeros(td_shape);td_shape(:)=td(1:numel(td_shape));td=td_shape;
end
%DECK DOMN