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

function [q,n,nrda,diag,irank,div,td]=ohtror(q,n,nrda,diag,irank,div,td);
persistent dd diagk irp j k kir kirm l nmir qs sig sqd tdv ; 

if isempty(dd), dd=0; end;
if isempty(diagk), diagk=0; end;
if isempty(qs), qs=0; end;
if isempty(sig), sig=0; end;
if isempty(sqd), sqd=0; end;
if isempty(tdv), tdv=0; end;
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;
%***BEGIN PROLOGUE  OHTROR
%***SUBSIDIARY
%***PURPOSE  Subsidiary to BVSUP
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (OHTROR-S)
%***AUTHOR  Watts, H. A., (SNLA)
%***DESCRIPTION
%
%     For a rank deficient problem, additional orthogonal
%     HOUSEHOLDER transformations are applied to the right side
%     of Q to further reduce the triangular form.
%     Thus, after application of the routines ORTHOL and OHTROR
%     to the original matrix, the result is a nonsingular
%     triangular matrix while the remainder of the matrix
%     has been zeroed out.
%
%***SEE ALSO  BVSUP
%***ROUTINES CALLED  SDOT
%***REVISION HISTORY  (YYMMDD)
%   750601  DATE WRITTEN
%   890831  Modified array declarations.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900402  Added TYPE section.  (WRB)
%   910722  Updated AUTHOR section.  (ALS)
%***end PROLOGUE  OHTROR
q_shape=size(q);q=reshape([q(:).',zeros(1,ceil(numel(q)./prod([nrda])).*prod([nrda])-numel(q))],nrda,[]);
diag_shape=size(diag);diag=reshape(diag,1,[]);
div_shape=size(div);div=reshape(div,1,[]);
td_shape=size(td);td=reshape(td,1,[]);
%***FIRST EXECUTABLE STATEMENT  OHTROR
nmir = fix(n - irank);
irp = fix(irank + 1);
for k = 1 : irank;
kir = fix(irp - k);
diagk = diag(kir);
sig =(diagk.*diagk) + sdot(nmir,q(sub2ind(size(q),kir,irp):end),nrda,q(sub2ind(size(q),kir,irp):end),nrda);
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(j,kir))+sdot(nmir,q(sub2ind(size(q),j,irp):end),nrda,q(sub2ind(size(q),kir,irp):end),nrda))./sqd;
q(j,kir) = q(j,kir) + qs.*tdv;
for l = irp : n;
q(j,l) = q(j,l) + qs.*q(kir,l);
end; l = fix(n+1);
end; j = fix(kirm+1);
end;
end; k = fix(irank+1);
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
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;
td_shape=zeros(td_shape);td_shape(:)=td(1:numel(td_shape));td=td_shape;
end
%DECK ORTBAK