share/minpack/fortran/enorm.f

   1       double precision function enorm(n,x)
   2       integer n
   3       double precision x(n)
   4 c     **********
   5 c
   6 c     function enorm
   7 c
   8 c     given an n-vector x, this function calculates the
   9 c     euclidean norm of x.
  10 c
  11 c     the euclidean norm is computed by accumulating the sum of
  12 c     squares in three different sums. the sums of squares for the
  13 c     small and large components are scaled so that no overflows
  14 c     occur. non-destructive underflows are permitted. underflows
  15 c     and overflows do not occur in the computation of the unscaled
  16 c     sum of squares for the intermediate components.
  17 c     the definitions of small, intermediate and large components
  18 c     depend on two constants, rdwarf and rgiant. the main
  19 c     restrictions on these constants are that rdwarf**2 not
  20 c     underflow and rgiant**2 not overflow. the constants
  21 c     given here are suitable for every known computer.
  22 c
  23 c     the function statement is
  24 c
  25 c       double precision function enorm(n,x)
  26 c
  27 c     where
  28 c
  29 c       n is a positive integer input variable.
  30 c
  31 c       x is an input array of length n.
  32 c
  33 c     subprograms called
  34 c
  35 c       fortran-supplied ... dabs,dsqrt
  36 c
  37 c     argonne national laboratory. minpack project. march 1980.
  38 c     burton s. garbow, kenneth e. hillstrom, jorge j. more
  39 c
  40 c     **********
  41       integer i
  42       double precision agiant,floatn,one,rdwarf,rgiant,s1,s2,s3,xabs,
  43      *                 x1max,x3max,zero
  44       data one,zero,rdwarf,rgiant /1.0d0,0.0d0,3.834d-20,1.304d19/
  45       s1 = zero
  46       s2 = zero
  47       s3 = zero
  48       x1max = zero
  49       x3max = zero
  50       floatn = n
  51       agiant = rgiant/floatn
  52       do 90 i = 1, n
  53          xabs = dabs(x(i))
  54          if (xabs .gt. rdwarf .and. xabs .lt. agiant) go to 70
  55             if (xabs .le. rdwarf) go to 30
  56 c
  57 c              sum for large components.
  58 c
  59                if (xabs .le. x1max) go to 10
  60                   s1 = one + s1*(x1max/xabs)**2
  61                   x1max = xabs
  62                   go to 20
  63    10          continue
  64                   s1 = s1 + (xabs/x1max)**2
  65    20          continue
  66                go to 60
  67    30       continue
  68 c
  69 c              sum for small components.
  70 c
  71                if (xabs .le. x3max) go to 40
  72                   s3 = one + s3*(x3max/xabs)**2
  73                   x3max = xabs
  74                   go to 50
  75    40          continue
  76                   if (xabs .ne. zero) s3 = s3 + (xabs/x3max)**2
  77    50          continue
  78    60       continue
  79             go to 80
  80    70    continue
  81 c
  82 c           sum for intermediate components.
  83 c
  84             s2 = s2 + xabs**2
  85    80    continue
  86    90    continue
  87 c
  88 c     calculation of norm.
  89 c
  90       if (s1 .eq. zero) go to 100
  91          enorm = x1max*dsqrt(s1+(s2/x1max)/x1max)
  92          go to 130
  93   100 continue
  94          if (s2 .eq. zero) go to 110
  95             if (s2 .ge. x3max)
  96      *         enorm = dsqrt(s2*(one+(x3max/s2)*(x3max*s3)))
  97             if (s2 .lt. x3max)
  98      *         enorm = dsqrt(x3max*((s2/x3max)+(x3max*s3)))
  99             go to 120
 100   110    continue
 101             enorm = x3max*dsqrt(s3)
 102   120    continue
 103   130 continue
 104       return
 105 c
 106 c     last card of function enorm.
 107 c
 108       end