[dd1,dd2,dx1,dy1,dparam]=drotmg(dd1,dd2,dx1,dy1,dparam);

function [dd1,dd2,dx1,dy1,dparam]=drotmg(dd1,dd2,dx1,dy1,dparam);
persistent dflag dh11 dh12 dh21 dh22 dp1 dp2 dq1 dq2 dtemp du firstCall gam gamsq gt igo one rgamsq two zero ; if isempty(firstCall),firstCall=1;end; 

if isempty(igo), igo=0; end;
if isempty(gt), gt=zeros(1,5); end;
%***BEGIN PROLOGUE  DROTMG
%***PURPOSE  Construct a modified Givens transformation.
%***LIBRARY   SLATEC (BLAS)
%***CATEGORY  D1B10
%***TYPE      doubleprecision (SROTMG-S, DROTMG-D)
%***KEYWORDS  BLAS, LINEAR ALGEBRA, MODIFIED GIVENS ROTATION, VECTOR
%***AUTHOR  Lawson, C. L., (JPL)
%           Hanson, R. J., (SNLA)
%           Kincaid, D. R., (U. of Texas)
%           Krogh, F. T., (JPL)
%***DESCRIPTION
%
%                B L A S  Subprogram
%    Description of Parameters
%
%     --Input--
%      DD1  doubleprecision scalar
%      DD2  doubleprecision scalar
%      DX1  doubleprecision scalar
%      DX2  doubleprecision scalar
%   DPARAM  D.P. 5-vector. DPARAM(1)=DFLAG defined below.
%           Locations 2-5 contain the rotation matrix.
%
%     --Output--
%      DD1  changed to represent the effect of the transformation
%      DD2  changed to represent the effect of the transformation
%      DX1  changed to represent the effect of the transformation
%      DX2  unchanged
%
%     Construct the modified Givens transformation matrix H which zeros
%     the second component of the 2-vector  (SQRT(DD1)*DX1,SQRT(DD2)*
%     DY2)**T.
%     With DPARAM(1)=DFLAG, H has one of the following forms:
%
%     DFLAG=-1.0D0     DFLAG=0.0D0        DFLAG=1.0D0     DFLAG=-2.0D0
%
%       (DH11  DH12)    (1.0D0  DH12)    (DH11  1.0D0)    (1.0D0  0.0D0)
%     H=(          )    (          )    (          )    (          )
%       (DH21  DH22),   (DH21  1.0D0),   (-1.0D0 DH22),   (0.0D0  1.0D0).
%
%     Locations 2-5 of DPARAM contain DH11, DH21, DH12, and DH22,
%     respectively.  (Values of 1.0D0, -1.0D0, or 0.0D0 implied by the
%     value of DPARAM(1) are not stored in DPARAM.)
%
%***REFERENCES  C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
%                 Krogh, Basic linear algebra subprograms for Fortran
%                 usage, Algorithm No. 539, Transactions on Mathematical
%                 Software 5, 3 (September 1979), pp. 308-323.
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   780301  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   920316  Prologue corrected.  (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  DROTMG
if isempty(gam), gam=0; end;
if isempty(one), one=0; end;
if isempty(rgamsq), rgamsq=0; end;
if isempty(dh11), dh11=0; end;
if isempty(dh12), dh12=0; end;
if isempty(dh21), dh21=0; end;
if isempty(dh22), dh22=0; end;
if isempty(dp1), dp1=0; end;
if isempty(dp2), dp2=0; end;
if isempty(dq1), dq1=0; end;
if isempty(dq2), dq2=0; end;
if isempty(du), du=0; end;
if isempty(zero), zero=0; end;
if isempty(gamsq), gamsq=0; end;
if isempty(dflag), dflag=0; end;
if isempty(dtemp), dtemp=0; end;
if isempty(two), two=0; end;
if firstCall,   zero =[0.0d0];  end;
if firstCall,  one =[1.0d0];  end;
if firstCall,  two=[2.0d0];  end;
if firstCall,   gam =[4096.0d0];  end;
if firstCall,  gamsq =[16777216.0d0];  end;
if firstCall,  rgamsq=[5.9604645d-8];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  DROTMG
gt(:)=0;
while (1);
if( dd1>=zero )
%     CASE-DD1-NONNEGATIVE
dp2 = dd2.*dy1;
if( dp2==zero )
%       GO ZERO-H-D-AND-DX1..
dflag = -two;
dparam(1) = dflag;
return;
else;
%     REGULAR-CASE..
dp1 = dd1.*dx1;
dq2 = dp2.*dy1;
dq1 = dp1.*dx1;
%
if( ~(~(abs(dq1)>abs(dq2))) )
dh21 = -dy1./dx1;
dh12 = dp2./dp1;
%
du = one - dh12.*dh21;
%
if( du<=zero )
break;
end;
dflag = zero;
dd1 = dd1./du;
dd2 = dd2./du;
dx1 = dx1.*du;
%         GO SCALE-CHECK..
elseif( dq2<zero ) ;
break;
else;
dflag = one;
dh11 = dp1./dp2;
dh22 = dx1./dy1;
du = one + dh11.*dh22;
dtemp = dd2./du;
dd2 = dd1./du;
dd1 = dtemp;
dx1 = dy1.*du;
%         GO SCALE-CHECK
end;
end;
%     PROCEDURE..SCALE-CHECK
while( true );
if(gt(4)==0)
if(gt(3)==0)
if(gt(2)==0)
if(gt(1)==0)
if( ~(dd1<=rgamsq) )
gt(2)=1;
continue;
end;
if( dd1==zero )
gt(3)=1;
continue;
end;
igo = 300;
%              FIX-H..
%     PROCEDURE..FIX-H..
end;
gt(1)=0;
if( ~(~(dflag>=zero)) )
%
if( dflag==zero )
dh11 = one;
dh22 = one;
dflag = -one;
else;
dh21 = -one;
dh12 = one;
dflag = -one;
end;
end;
if( igo==300 )
dd1 = dd1.*gam.^2;
dx1 = dx1./gam;
dh11 = dh11./gam;
dh12 = dh12./gam;
continue;
elseif( igo==500 ) ;
dd1 = dd1./gam.^2;
dx1 = dx1.*gam;
dh11 = dh11.*gam;
dh12 = dh12.*gam;
elseif( igo==700 ) ;
dd2 = dd2.*gam.^2;
dh21 = dh21./gam;
dh22 = dh22./gam;
gt(3)=1;
continue;
elseif( igo==900 ) ;
dd2 = dd2./gam.^2;
dh21 = dh21.*gam;
dh22 = dh22.*gam;
gt(4)=1;
continue;
else;
continue;
end;
end;
gt(2)=0;
if( dd1>=gamsq )
igo = 500;
%              FIX-H..
gt(1)=1;
continue;
end;
end;
gt(3)=0;
if( abs(dd2)<=rgamsq )
if( dd2==zero )
break;
end;
igo = 700;
%              FIX-H..
gt(1)=1;
continue;
end;
end;
gt(4)=0;
if( ~(abs(dd2)>=gamsq) )
break;
end;
igo = 900;
%              FIX-H..
gt(1)=1;
continue;
end;
gt(5)=1;
end;
break;
end;
if(gt(5)==0)
%     PROCEDURE..ZERO-H-D-AND-DX1..
dflag = -one;
dh11 = zero;
dh12 = zero;
dh21 = zero;
dh22 = zero;
%
dd1 = zero;
dd2 = zero;
dx1 = zero;
end;
gt(5)=0;
%         RETURN..
if( dflag<0 )
dparam(2) = dh11;
dparam(3) = dh21;
dparam(4) = dh12;
dparam(5) = dh22;
elseif( dflag==0 ) ;
dparam(3) = dh21;
dparam(4) = dh12;
else;
dparam(2) = dh11;
dparam(5) = dh22;
end;
dparam(1) = dflag;
end %subroutine drotmg
%DECK DRSCO