subroutine DLARTG7( f, h, g, r, w, rtmin )
   use LA_CONSTANTS, &
   only: wp=>dp, zero=>dzero, half=>dhalf, one=>done, &
         safmin=>dsafmin, safmax=>dsafmax
!
!  .. Scalar Arguments ..
   real(wp) :: f, g, h, r, w
!  ..
!  .. Local Scalars ..
   real(wp) :: fs, gs, hs, u, rtmin, fhh, fhl, gh, gl, p0, e0
!  ..
!  .. Intrinsic Functions ..
   intrinsic :: abs, sign, sqrt
!  ..
!  .. Executable Statements ..
!
   if( g == zero ) then
      r = f
      w = h
   else if( f == zero ) then
      r = g
      w = zero
   else if( f > rtmin .and. g > rtmin ) then
      CALL SQUARE(f,h,fhh,fhl)
      CALL TWO_PROD(g,g,gh,gl)
      CALL ADD(fhh,fhl,gh,gl,p0,e0)
      IF (p0 .EQ. ZERO) THEN
         r = zero
         w = zero
      ELSE
         CALL USER_SQRT(p0,e0,r,w)
      ENDIF
   else
      u = min( safmax, max( safmin, f, g ) )
      fs = f / u
      gs = g / u
      hs = h / u
      CALL SQUARE(fs,hs,fhh,fhl)
      CALL TWO_PROD(gs,gs,gh,gl)
      CALL ADD(fhh,fhl,gh,gl,p0,e0)
      IF (p0 .EQ. ZERO) THEN
         r = zero
         w = zero
      ELSE
         CALL USER_SQRT(p0,e0,r,w)
      ENDIF
      r = r*u
      w = w*u
   end if
   return
end subroutine

DOUBLE PRECISION FUNCTION USERDNRM2( N, X, INCX )
  IMPLICIT NONE
  
  INTEGER N, I, INCX
  DOUBLE PRECISION X(*)
  DOUBLE PRECISION TMP1, TMP2, TMP3, TMP4
  
  DOUBLE PRECISION   ZERO
  PARAMETER          ( ZERO = 0.0D0 )
  
  TMP1 = ZERO
  TMP2 = ZERO
  
  IF (INCX .EQ. 1) THEN
     
     DO I=1,N
        CALL TWO_PROD(X(I),X(I),TMP3,TMP4)
        CALL ADD(TMP1,TMP2,TMP3,TMP4,TMP1,TMP2)
     ENDDO
     IF (TMP1 .EQ. ZERO) THEN
        TMP3 = ZERO
     ELSE
        CALL USER_SQRT(TMP1,TMP2,TMP3,TMP4)
     ENDIF
     USERDNRM2 = TMP3
  ELSE
     
     DO I=1,N
        CALL TWO_PROD(X(1+INCX*(I-1)),X(1+INCX*(I-1)),TMP3,TMP4)
        CALL ADD(TMP1,TMP2,TMP3,TMP4,TMP1,TMP2)
     ENDDO
     IF (TMP1 .EQ. ZERO) THEN
        TMP3 = ZERO
     ELSE
        CALL USER_SQRT(TMP1,TMP2,TMP3,TMP4)
     ENDIF
     USERDNRM2 = TMP3
  ENDIF
  
  RETURN
  
END FUNCTION USERDNRM2

SUBROUTINE FAST_TWO_SUM(A0,B0,S0,E0)
  IMPLICIT NONE
  DOUBLE PRECISION A0,B0,S0,E0
  S0=A0+B0
  E0=B0-(S0-A0)
  RETURN
END SUBROUTINE FAST_TWO_SUM

SUBROUTINE TWO_SUM(A0,B0,S0,E0)
  IMPLICIT NONE
  DOUBLE PRECISION A0,B0,S0,E0,V0
  S0=A0+B0
  V0=S0-A0
  E0=(A0-(S0-V0))+(B0-V0)
  RETURN
END SUBROUTINE TWO_SUM

SUBROUTINE ADD(AH,AL,BH,BL,CH,CL)
  IMPLICIT NONE
  DOUBLE PRECISION AH,AL,BH,BL,CH,CL,SH,EH
  CALL TWO_SUM(AH,BH,SH,EH)
  EH=EH+AL+BL
  CALL FAST_TWO_SUM(SH,EH,CH,CL)
  RETURN
END SUBROUTINE ADD

SUBROUTINE ADD2(AH,AL,BH,BL,CH)
  IMPLICIT NONE
  DOUBLE PRECISION AH,AL,BH,BL,CH,SH,EH
  CALL TWO_SUM(AH,BH,SH,EH)
  EH=EH+AL+BL
  CH=SH+EH
  RETURN
END SUBROUTINE ADD2

SUBROUTINE TWO_PROD(A0,B0,P0,E0)
  IMPLICIT NONE
  DOUBLE PRECISION A0,B0,AH,AL,BH,BL,P0,E0
  DOUBLE PRECISION DFMA0
  EXTERNAL DFMA0
  P0=A0*B0
  E0=DFMA0(A0,B0,-P0)
  RETURN
END SUBROUTINE TWO_PROD

SUBROUTINE SQUARE(AH,AL,CH,CL)
  IMPLICIT NONE
  DOUBLE PRECISION AH,AL,CH,CL,P1,P2
  DOUBLE PRECISION TWO
  PARAMETER          ( TWO = 2.0D0 )
  DOUBLE PRECISION DFMA0
  EXTERNAL DFMA0
  CALL TWO_PROD(AH,AH,P1,P2)
  P2=DFMA0(TWO*AL,AH,P2)
  CALL FAST_TWO_SUM(P1,P2,CH,CL)
  RETURN
END SUBROUTINE SQUARE

SUBROUTINE USER_SQRT(AH,AL,BH,BL)
  IMPLICIT NONE
  DOUBLE PRECISION AH,AL,BH,BL,X,CH,CL,DH
  DOUBLE PRECISION ZERO, HALF, ONE
  PARAMETER          ( ZERO=0.0D0, HALF=0.50D0, ONE=1.0D0 )
  
  X=SQRT(AH)
  CALL TWO_PROD(X,X,CH,CL)
  CALL ADD2(AH,AL,-CH,-CL,DH)
  CALL FAST_TWO_SUM(X,HALF*DH/X,BH,BL)
  
  RETURN
END SUBROUTINE USER_SQRT