Line data    Source code

       1           0 :       double precision function pda_enorm(n,x)
       2             :       integer n
       3             :       double precision x(n)
       4             : c     **********
       5             : c
       6             : c     function pda_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 pda_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           0 :       s1 = zero
      46           0 :       s2 = zero
      47           0 :       s3 = zero
      48           0 :       x1max = zero
      49           0 :       x3max = zero
      50           0 :       floatn = n
      51           0 :       agiant = rgiant/floatn
      52           0 :       do 90 i = 1, n
      53           0 :          xabs = dabs(x(i))
      54           0 :          if (xabs .gt. rdwarf .and. xabs .lt. agiant) go to 70
      55           0 :             if (xabs .le. rdwarf) go to 30
      56             : c
      57             : c              sum for large components.
      58             : c
      59           0 :                if (xabs .le. x1max) go to 10
      60           0 :                   s1 = one + s1*(x1max/xabs)**2
      61           0 :                   x1max = xabs
      62           0 :                   go to 20
      63             :    10          continue
      64           0 :                   s1 = s1 + (xabs/x1max)**2
      65             :    20          continue
      66           0 :                go to 60
      67             :    30       continue
      68             : c
      69             : c              sum for small components.
      70             : c
      71           0 :                if (xabs .le. x3max) go to 40
      72           0 :                   s3 = one + s3*(x3max/xabs)**2
      73           0 :                   x3max = xabs
      74           0 :                   go to 50
      75             :    40          continue
      76           0 :                   if (xabs .ne. zero) s3 = s3 + (xabs/x3max)**2
      77             :    50          continue
      78             :    60       continue
      79           0 :             go to 80
      80             :    70    continue
      81             : c
      82             : c           sum for intermediate components.
      83             : c
      84           0 :             s2 = s2 + xabs**2
      85             :    80    continue
      86           0 :    90    continue
      87             : c
      88             : c     calculation of norm.
      89             : c
      90           0 :       if (s1 .eq. zero) go to 100
      91           0 :          pda_enorm = x1max*dsqrt(s1+(s2/x1max)/x1max)
      92           0 :          go to 130
      93             :   100 continue
      94           0 :          if (s2 .eq. zero) go to 110
      95           0 :             if (s2 .ge. x3max)
      96           0 :      *         pda_enorm = dsqrt(s2*(one+(x3max/s2)*(x3max*s3)))
      97           0 :             if (s2 .lt. x3max)
      98           0 :      *         pda_enorm = dsqrt(x3max*((s2/x3max)+(x3max*s3)))
      99           0 :             go to 120
     100             :   110    continue
     101           0 :             pda_enorm = x3max*dsqrt(s3)
     102             :   120    continue
     103             :   130 continue
     104           0 :       return
     105             : c
     106             : c     last card of function pda_enorm.
     107             : c
     108             :       end