[mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv]=dh12(mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv);

function [mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv]=dh12(mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv);
%***BEGIN PROLOGUE  DH12
%***SUBSIDIARY
%***PURPOSE  Subsidiary to DHFTI, DLSEI and DWNNLS
%***LIBRARY   SLATEC
%***TYPE      doubleprecision (H12-S, DH12-D)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%      *** doubleprecision VERSION OF H12 ******
%
%     C.L.Lawson and R.J.Hanson, Jet Propulsion Laboratory, 1973 Jun 12
%     to appear in 'Solving Least Squares Problems', Prentice-Hall, 1974
%
%     Construction and/or application of a single
%     Householder transformation..     Q = I + U*(U**T)/B
%
%     MODE    = 1 or 2   to select algorithm  H1  or  H2 .
%     LPIVOT is the index of the pivot element.
%     L1,M   If L1 .LE. M   the transformation will be constructed to
%            zero elements indexed from L1 through M.   If L1 GT. M
%            THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION.
%     U(),IUE,UP    On entry to H1 U() contains the pivot vector.
%                   IUE is the storage increment between elements.
%                                       On exit from H1 U() and UP
%                   contain quantities defining the vector U of the
%                   Householder transformation.   On entry to H2 U()
%                   and UP should contain quantities previously computed
%                   by H1.  These will not be modified by H2.
%     C()    On entry to H1 or H2 C() contains a matrix which will be
%            regarded as a set of vectors to which the Householder
%            transformation is to be applied.  On exit C() contains the
%            set of transformed vectors.
%     ICE    Storage increment between elements of vectors in C().
%     ICV    Storage increment between vectors in C().
%     NCV    Number of vectors in C() to be transformed. If NCV .LE. 0
%            no operations will be done on C().
%
%***SEE ALSO  DHFTI, DLSEI, DWNNLS
%***ROUTINES CALLED  DAXPY, DDOT, DSWAP
%***REVISION HISTORY  (YYMMDD)
%   790101  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)
%   900911  Added DDOT to doubleprecision statement.  (WRB)
%***end PROLOGUE  DH12
persistent b cl clinv i i2 i3 i4 incr j kl1 kl2 klp l1m1 mml1p2 one sm ul1m1 ; 

if isempty(i), i=0; end;
if isempty(i2), i2=0; end;
if isempty(i3), i3=0; end;
if isempty(i4), i4=0; end;
if isempty(incr), incr=0; end;
if isempty(j), j=0; end;
if isempty(kl1), kl1=0; end;
if isempty(kl2), kl2=0; end;
if isempty(klp), klp=0; end;
if isempty(l1m1), l1m1=0; end;
if isempty(mml1p2), mml1p2=0; end;
if isempty(b), b=0; end;
if isempty(cl), cl=0; end;
if isempty(clinv), clinv=0; end;
if isempty(one), one=0; end;
if isempty(ul1m1), ul1m1=0; end;
if isempty(sm), sm=0; end;
u_shape=size(u);u=reshape([u(:).',zeros(1,ceil(numel(u)./prod([iue])).*prod([iue])-numel(u))],iue,[]);
c_shape=size(c);c=reshape(c,1,[]);
%     BEGIN BLOCK PERMITTING ...EXITS TO 140
%***FIRST EXECUTABLE STATEMENT  DH12
one = 1.0d0;
%
%     ...EXIT
if( 0<lpivot && lpivot<l1 && l1<=m )
cl = abs(u(1,lpivot));
if( mode~=2 )
%           ****** CONSTRUCT THE TRANSFORMATION. ******
for j = l1 : m;
cl = max(abs(u(1,j)),cl);
end; j = fix(m+1);
%     .........EXIT
if( cl<=0.0d0 )
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
return;
end;
clinv = one./cl;
sm =(u(1,lpivot).*clinv).^2;
for j = l1 : m;
sm = sm +(u(1,j).*clinv).^2;
end; j = fix(m+1);
cl = cl.*sqrt(sm);
if( u(1,lpivot)>0.0d0 )
cl = -cl;
end;
up = u(1,lpivot) - cl;
u(1,lpivot) = cl;
%        ****** APPLY THE TRANSFORMATION  I+U*(U**T)/B  TO C. ******
%
elseif( cl<=0.0d0 ) ;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
return;
%     ......EXIT
end;
%     ...EXIT
if( ncv>0 )
b = up.*u(1,lpivot);
%        B  MUST BE NONPOSITIVE HERE.  IF B = 0., RETURN.
%
if( b<0.0d0 )
b = one./b;
mml1p2 = fix(m - l1 + 2);
if( mml1p2<=20 )
i2 = fix(1 - icv + ice.*(lpivot-1));
incr = fix(ice.*(l1-lpivot));
for j = 1 : ncv;
i2 = fix(i2 + icv);
i3 = fix(i2 + incr);
i4 = fix(i3);
sm = c(i2).*up;
for i = l1 : m;
sm = sm + c(i3).*u(1,i);
i3 = fix(i3 + ice);
end; i = fix(m+1);
if( sm~=0.0d0 )
sm = sm.*b;
c(i2) = c(i2) + sm.*up;
for i = l1 : m;
c(i4) = c(i4) + sm.*u(1,i);
i4 = fix(i4 + ice);
end; i = fix(m+1);
end;
end; j = fix(ncv+1);
else;
l1m1 = fix(l1 - 1);
kl1 = fix(1 +(l1m1-1).*ice);
kl2 = fix(kl1);
klp = fix(1 +(lpivot-1).*ice);
ul1m1 = u(1,l1m1);
u(1,l1m1) = up;
if( lpivot~=l1m1 )
icv_orig=icv;    [ncv,dumvar2,icv,dumvar4,dumvar5]=dswap(ncv,c(sub2ind(size(c),max(kl1,1)):end),icv,c(sub2ind(size(c),max(klp,1)):end),icv);    icv(dumvar5~=icv_orig)=dumvar5(dumvar5~=icv_orig);   dumvar2i=find((c(sub2ind(size(c),max(kl1,1)):end))~=(dumvar2));dumvar4i=find((c(sub2ind(size(c),max(klp,1)):end))~=(dumvar4));   c(kl1-1+dumvar2i)=dumvar2(dumvar2i); c(klp-1+dumvar4i)=dumvar4(dumvar4i); 
end;
for j = 1 : ncv;
[sm ,mml1p2,u(sub2ind(size(u),1,l1m1):end),iue,c(sub2ind(size(c),max(kl1,1)):end),ice]=ddot(mml1p2,u(sub2ind(size(u),1,l1m1):end),iue,c(sub2ind(size(c),max(kl1,1)):end),ice);
sm = sm.*b;
[mml1p2,sm,u(sub2ind(size(u),1,l1m1):end),iue,c(sub2ind(size(c),max(kl1,1)):end),ice]=daxpy(mml1p2,sm,u(sub2ind(size(u),1,l1m1):end),iue,c(sub2ind(size(c),max(kl1,1)):end),ice);
kl1 = fix(kl1 + icv);
end; j = fix(ncv+1);
u(1,l1m1) = ul1m1;
%     ......EXIT
if( lpivot~=l1m1 )
kl1 = fix(kl2);
icv_orig=icv;    [ncv,dumvar2,icv,dumvar4,dumvar5]=dswap(ncv,c(sub2ind(size(c),max(kl1,1)):end),icv,c(sub2ind(size(c),max(klp,1)):end),icv);    icv(dumvar5~=icv_orig)=dumvar5(dumvar5~=icv_orig);   dumvar2i=find((c(sub2ind(size(c),max(kl1,1)):end))~=(dumvar2));dumvar4i=find((c(sub2ind(size(c),max(klp,1)):end))~=(dumvar4));   c(kl1-1+dumvar2i)=dumvar2(dumvar2i); c(klp-1+dumvar4i)=dumvar4(dumvar4i); 
end;
end;
%     ......EXIT
end;
end;
end;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
end
%DECK DHELS