ENORM Function

  real(dp) FUNCTION ENORM(n,x)

    !  Given an N-vector X, this function calculates the
    !  Euclidean norm of X.
    !
    !  The Euclidean norm is computed by accumulating the sum of
    !  squares in three different sums. The sums of squares for the
    !  small and large components are scaled so that no overflows
    !  occur. Non-destructive underflows are permitted. Underflows
    !  and overflows do not occur in the computation of the unscaled
    !  sum of squares for the intermediate components.
    !  The definitions of small, intermediate and large components
    !  depend on two constants, RDWARF and RGIANT. The main
    !  restrictions on these constants are that RDWARF**2 not
    !  underflow and RGIANT**2 not overflow. The constants
    !  given here are suitable for every known computer.
    !
    !  The function statement is
    !
    !  real function enorm(n,x)
    !
    !  where
    !
    !  N is a positive integer input variable.
    !
    !  X is an input array of length N.
    !
    !  Argonne National Laboratory. Minpack project. March 1980.
    !  Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More

    IMPLICIT NONE

    INTEGER n,i

    real(dp) x(n)
    real(dp) agiant,floatn,one,rdwarf,rgiant,s1,s2,s3,xabs, &
         x1max,x3max,zero

    one = 1.0D0
    zero = 0.0D0
    rdwarf = 3.834d-20
    rgiant = 1.304d19

    s1 = zero
    s2 = zero
    s3 = zero
    x1max = zero
    x3max = zero
    floatn = dble(n)
    agiant = rgiant/floatn
    do i = 1, n
       xabs = abs(x(i))
       if ((xabs > rdwarf) .and. (xabs < agiant)) goto 70
       if (xabs <= rdwarf) goto 30

       !  Sum for large components.

       if (xabs <= x1max) goto 10
       s1 = one + s1*(x1max/xabs)**2
       x1max = xabs
       goto 20

10     continue
       s1 = s1 + (xabs/x1max)**2

20     continue
       goto 60

30     continue

       !  Sum for small components.

       if (xabs <= x3max) goto 40
       s3 = one + s3*(x3max/xabs)**2
       x3max = xabs
       goto 50

40     continue
       if (xabs  /=  zero) s3 = s3 + (xabs/x3max)**2

50     continue
60     continue
       goto 80

70     continue

       !  Sum for intermediate components.

       s2 = s2 + xabs**2

80     continue
    end do

    !  Calculation of norm.

    if (s1 == zero) goto 100
    enorm = x1max*sqrt(s1+(s2/x1max)/x1max)
    goto 130

100 continue
    if (s2 == zero) goto 110
    if (s2 >= x3max) &
         enorm = sqrt(s2*(one+(x3max/s2)*(x3max*s3)))
    if (s2 < x3max) &
         enorm = sqrt(x3max*((s2/x3max)+(x3max*s3)))
    goto 120
110 continue
    enorm = x3max*sqrt(s3)

120 continue
130 continue

    return
  end FUNCTION ENORM