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, tmp
!  ..
!  .. 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)
      CALL USER_SQRT(p0,e0,r,w)
   else
      u = min( safmax, max( safmin, f, g ) )
      CALL USER_DIV2(f,h,u,fs,hs)
      CALL USER_DIV4(g,u,gs,tmp)
      CALL SQUARE(fs,hs,fhh,fhl)
      CALL SQUARE(gs,tmp,gh,gl)
      CALL ADD(fhh,fhl,gh,gl,p0,e0)
      CALL USER_SQRT(p0,e0,r,w)
      CALL PROD2(u,r,w,r,w)
   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
     CALL USER_SQRT(TMP1,TMP2,TMP3,TMP4)
     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
     CALL USER_SQRT(TMP1,TMP2,TMP3,TMP4)
     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 ADDP(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 TWO_SUM(SH,EH,CH,CL)
  RETURN
END SUBROUTINE ADDP

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

SUBROUTINE PROD(AH,AL,BH,BL,CH,CL)
  IMPLICIT NONE
  DOUBLE PRECISION AH,AL,BH,BL,CH,CL,SH,EH
  DOUBLE PRECISION DFMA0
  EXTERNAL DFMA0
  CALL TWO_PROD(AH,BH,SH,EH)
  EH = DFMA0(AL,BH,EH)
  EH = DFMA0(AH,BL,EH)
  CALL FAST_TWO_SUM(SH,EH,CH,CL)
  RETURN
END SUBROUTINE PROD

SUBROUTINE PROD2(AH,BH,BL,CH,CL)
  IMPLICIT NONE
  DOUBLE PRECISION AH,BH,BL,CH,CL,SH,EH
  DOUBLE PRECISION DFMA0
  EXTERNAL DFMA0
  CALL TWO_PROD(AH,BH,SH,EH)
  EH = DFMA0(AH,BL,EH)
  CALL FAST_TWO_SUM(SH,EH,CH,CL)
  RETURN
END SUBROUTINE PROD2

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 )
  DOUBLE PRECISION DFMA0
  EXTERNAL DFMA0

  IF (AH .EQ. ZERO) THEN
     BH = ZERO
     BL = ZERO
  ELSE
     X = SQRT(AH)
     DH=HALF*(DFMA0(-X,X,AH)+AL)/X
     CALL FAST_TWO_SUM(X,DH,BH,BL)
  ENDIF
  RETURN
END SUBROUTINE USER_SQRT

SUBROUTINE USER_DIV(AH,AL,BH,BL,CH,CL)
  IMPLICIT NONE
  DOUBLE PRECISION AH,AL,BH,BL,X,CH,CL,EH,EL,GH
  DOUBLE PRECISION ZERO, HALF, ONE
  PARAMETER          ( ZERO=0.0D0, HALF=0.50D0, ONE=1.0D0 )
  DOUBLE PRECISION DFMA0
  EXTERNAL DFMA0

  X = AH/BH
  GH = DFMA0(-X,BL,DFMA0(-X,BH,AH)+AL)/BH
  CALL FAST_TWO_SUM(X,GH,CH,CL)
  RETURN
END SUBROUTINE USER_DIV

SUBROUTINE USER_DIV2(AH,AL,BH,CH,CL)
  IMPLICIT NONE
  DOUBLE PRECISION BH,X,GH,AH,AL,CH,CL,EH,EL
  DOUBLE PRECISION ZERO, HALF, ONE
  PARAMETER          ( ZERO=0.0D0, HALF=0.50D0, ONE=1.0D0 )
  DOUBLE PRECISION DFMA0
  EXTERNAL DFMA0
  
  X = AH/BH
  GH = (DFMA0(-X,BH,AH)+AL)/BH
  CALL FAST_TWO_SUM(X,GH,CH,CL)
  RETURN
END SUBROUTINE USER_DIV2

SUBROUTINE USER_DIV4(AH,BH,CH,CL)
  IMPLICIT NONE
  DOUBLE PRECISION AH,BH,X,GH,EH,EL,CH,CL
  DOUBLE PRECISION ZERO, HALF, ONE
  PARAMETER          ( ZERO=0.0D0, HALF=0.50D0, ONE=1.0D0 )
  DOUBLE PRECISION DFMA0
  EXTERNAL DFMA0
  
  X = AH/BH
  GH = DFMA0(-X,BH,AH)/BH
  CALL FAST_TWO_SUM(X,GH,CH,CL)
  RETURN
END SUBROUTINE USER_DIV4