home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / blas / dznrm2.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  139 lines
  1.       double precision function dznrm2( n, zx, incx)
  2.       logical imag, scale
  3.       integer i, incx, ix, n, next
  4.       double precision cutlo, cuthi, hitest, sum, xmax, absx, zero, one
  5.       double complex      zx(1)
  6.       double precision dreal,dimag
  7.       double complex zdumr,zdumi
  8.       dreal(zdumr) = zdumr
  9.       dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
  10.       data         zero, one /0.0d0, 1.0d0/
  11. c
  12. c     unitary norm of the complex n-vector stored in zx() with storage
  13. c     increment incx .
  14. c     if    n .le. 0 return with result = 0.
  15. c     if n .ge. 1 then incx must be .ge. 1
  16. c
  17. c           c.l.lawson , 1978 jan 08
  18. c     modified to correct problem with negative increment, 8/21/90.
  19. c
  20. c     four phase method     using two built-in constants that are
  21. c     hopefully applicable to all machines.
  22. c         cutlo = maximum of  sqrt(u/eps)  over all known machines.
  23. c         cuthi = minimum of  sqrt(v)      over all known machines.
  24. c     where
  25. c         eps = smallest no. such that eps + 1. .gt. 1.
  26. c         u   = smallest positive no.   (underflow limit)
  27. c         v   = largest  no.            (overflow  limit)
  28. c
  29. c     brief outline of algorithm..
  30. c
  31. c     phase 1    scans zero components.
  32. c     move to phase 2 when a component is nonzero and .le. cutlo
  33. c     move to phase 3 when a component is .gt. cutlo
  34. c     move to phase 4 when a component is .ge. cuthi/m
  35. c     where m = n for x() real and m = 2*n for complex.
  36. c
  37. c     values for cutlo and cuthi..
  38. c     from the environmental parameters listed in the imsl converter
  39. c     document the limiting values are as follows..
  40. c     cutlo, s.p.   u/eps = 2**(-102) for  honeywell.  close seconds are
  41. c                   univac and dec at 2**(-103)
  42. c                   thus cutlo = 2**(-51) = 4.44089e-16
  43. c     cuthi, s.p.   v = 2**127 for univac, honeywell, and dec.
  44. c                   thus cuthi = 2**(63.5) = 1.30438e19
  45. c     cutlo, d.p.   u/eps = 2**(-67) for honeywell and dec.
  46. c                   thus cutlo = 2**(-33.5) = 8.23181d-11
  47. c     cuthi, d.p.   same as s.p.  cuthi = 1.30438d19
  48. c     data cutlo, cuthi / 8.232d-11,  1.304d19 /
  49. c     data cutlo, cuthi / 4.441e-16,  1.304e19 /
  50.       data cutlo, cuthi / 8.232d-11,  1.304d19 /
  51. c
  52.       if(n .gt. 0) go to 10
  53.          dznrm2  = zero
  54.          go to 300
  55. c
  56.    10 assign 30 to next
  57.       sum = zero
  58.       i = 1
  59.       if( incx .lt. 0 )i = (-n+1)*incx + 1
  60. c                                                 begin main loop
  61.       do 220 ix = 1,n
  62.          absx = dabs(dreal(zx(i)))
  63.          imag = .false.
  64.          go to next,(30, 50, 70, 90, 110)
  65.    30 if( absx .gt. cutlo) go to 85
  66.       assign 50 to next
  67.       scale = .false.
  68. c
  69. c                        phase 1.  sum is zero
  70. c
  71.    50 if( absx .eq. zero) go to 200
  72.       if( absx .gt. cutlo) go to 85
  73. c
  74. c                                prepare for phase 2.
  75.       assign 70 to next
  76.       go to 105
  77. c
  78. c                                prepare for phase 4.
  79. c
  80.   100 assign 110 to next
  81.       sum = (sum / absx) / absx
  82.   105 scale = .true.
  83.       xmax = absx
  84.       go to 115
  85. c
  86. c                   phase 2.  sum is small.
  87. c                             scale to avoid destructive underflow.
  88. c
  89.    70 if( absx .gt. cutlo ) go to 75
  90. c
  91. c                     common code for phases 2 and 4.
  92. c                     in phase 4 sum is large.  scale to avoid overflow.
  93. c
  94.   110 if( absx .le. xmax ) go to 115
  95.          sum = one + sum * (xmax / absx)**2
  96.          xmax = absx
  97.          go to 200
  98. c
  99.   115 sum = sum + (absx/xmax)**2
  100.       go to 200
  101. c
  102. c
  103. c                  prepare for phase 3.
  104. c
  105.    75 sum = (sum * xmax) * xmax
  106. c
  107.    85 assign 90 to next
  108.       scale = .false.
  109. c
  110. c     for real or d.p. set hitest = cuthi/n
  111. c     for complex      set hitest = cuthi/(2*n)
  112. c
  113.       hitest = cuthi/dble( 2*n )
  114. c
  115. c                   phase 3.  sum is mid-range.  no scaling.
  116. c
  117.    90 if(absx .ge. hitest) go to 100
  118.          sum = sum + absx**2
  119.   200 continue
  120. c                  control selection of real and imaginary parts.
  121. c
  122.       if(imag) go to 210
  123.          absx = dabs(dimag(zx(i)))
  124.          imag = .true.
  125.       go to next,(  50, 70, 90, 110 )
  126. c
  127.   210 continue
  128.       i = i + incx
  129.   220 continue
  130. c
  131. c              end of main loop.
  132. c              compute square root and adjust for scaling.
  133. c
  134.       dznrm2 = dsqrt(sum)
  135.       if(scale) dznrm2 = dznrm2 * xmax
  136.   300 continue
  137.       return
  138.       end
  139.