share/odepack/fortran/dnrm2.f

   1 *DECK DNRM2
   2       DOUBLE PRECISION FUNCTION DNRM2 (N, DX, INCX)
   3 C***BEGIN PROLOGUE  DNRM2
   4 C***PURPOSE  Compute the Euclidean length (L2 norm) of a vector.
   5 C***CATEGORY  D1A3B
   6 C***TYPE      DOUBLE PRECISION (SNRM2-S, DNRM2-D, SCNRM2-C)
   7 C***KEYWORDS  BLAS, EUCLIDEAN LENGTH, EUCLIDEAN NORM, L2,
   8 C             LINEAR ALGEBRA, UNITARY, VECTOR
   9 C***AUTHOR  Lawson, C. L., (JPL)
  10 C           Hanson, R. J., (SNLA)
  11 C           Kincaid, D. R., (U. of Texas)
  12 C           Krogh, F. T., (JPL)
  13 C***DESCRIPTION
  14 C
  15 C                B L A S  Subprogram
  16 C    Description of parameters
  17 C
  18 C     --Input--
  19 C        N  number of elements in input vector(s)
  20 C       DX  double precision vector with N elements
  21 C     INCX  storage spacing between elements of DX
  22 C
  23 C     --Output--
  24 C    DNRM2  double precision result (zero if N .LE. 0)
  25 C
  26 C     Euclidean norm of the N-vector stored in DX with storage
  27 C     increment INCX.
  28 C     If N .LE. 0, return with result = 0.
  29 C     If N .GE. 1, then INCX must be .GE. 1
  30 C
  31 C     Four phase method using two built-in constants that are
  32 C     hopefully applicable to all machines.
  33 C         CUTLO = maximum of  SQRT(U/EPS)  over all known machines.
  34 C         CUTHI = minimum of  SQRT(V)      over all known machines.
  35 C     where
  36 C         EPS = smallest no. such that EPS + 1. .GT. 1.
  37 C         U   = smallest positive no.   (underflow limit)
  38 C         V   = largest  no.            (overflow  limit)
  39 C
  40 C     Brief outline of algorithm.
  41 C
  42 C     Phase 1 scans zero components.
  43 C     move to phase 2 when a component is nonzero and .LE. CUTLO
  44 C     move to phase 3 when a component is .GT. CUTLO
  45 C     move to phase 4 when a component is .GE. CUTHI/M
  46 C     where M = N for X() real and M = 2*N for complex.
  47 C
  48 C     Values for CUTLO and CUTHI.
  49 C     From the environmental parameters listed in the IMSL converter
  50 C     document the limiting values are as follows:
  51 C     CUTLO, S.P.   U/EPS = 2**(-102) for  Honeywell.  Close seconds are
  52 C                   Univac and DEC at 2**(-103)
  53 C                   Thus CUTLO = 2**(-51) = 4.44089E-16
  54 C     CUTHI, S.P.   V = 2**127 for Univac, Honeywell, and DEC.
  55 C                   Thus CUTHI = 2**(63.5) = 1.30438E19
  56 C     CUTLO, D.P.   U/EPS = 2**(-67) for Honeywell and DEC.
  57 C                   Thus CUTLO = 2**(-33.5) = 8.23181D-11
  58 C     CUTHI, D.P.   same as S.P.  CUTHI = 1.30438D19
  59 C     DATA CUTLO, CUTHI /8.232D-11,  1.304D19/
  60 C     DATA CUTLO, CUTHI /4.441E-16,  1.304E19/
  61 C
  62 C***REFERENCES  C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
  63 C                 Krogh, Basic linear algebra subprograms for Fortran
  64 C                 usage, Algorithm No. 539, Transactions on Mathematical
  65 C                 Software 5, 3 (September 1979), pp. 308-323.
  66 C***ROUTINES CALLED  (NONE)
  67 C***REVISION HISTORY  (YYMMDD)
  68 C   791001  DATE WRITTEN
  69 C   890531  Changed all specific intrinsics to generic.  (WRB)
  70 C   890831  Modified array declarations.  (WRB)
  71 C   890831  REVISION DATE from Version 3.2
  72 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  73 C   920501  Reformatted the REFERENCES section.  (WRB)
  74 C***END PROLOGUE  DNRM2
  75       INTEGER NEXT
  76       DOUBLE PRECISION DX(*), CUTLO, CUTHI, HITEST, SUM, XMAX, ZERO,
  77      +                 ONE
  78       SAVE CUTLO, CUTHI, ZERO, ONE
  79       DATA ZERO, ONE /0.0D0, 1.0D0/
  80 C
  81       DATA CUTLO, CUTHI /8.232D-11,  1.304D19/
  82 C***FIRST EXECUTABLE STATEMENT  DNRM2
  83       IF (N .GT. 0) GO TO 10
  84          DNRM2  = ZERO
  85          GO TO 300
  86 C
  87    10 ASSIGN 30 TO NEXT
  88       SUM = ZERO
  89       NN = N * INCX
  90 C
  91 C                                                 BEGIN MAIN LOOP
  92 C
  93       I = 1
  94    20    GO TO NEXT,(30, 50, 70, 110)
  95    30 IF (ABS(DX(I)) .GT. CUTLO) GO TO 85
  96       ASSIGN 50 TO NEXT
  97       XMAX = ZERO
  98 C
  99 C                        PHASE 1.  SUM IS ZERO
 100 C
 101    50 IF (DX(I) .EQ. ZERO) GO TO 200
 102       IF (ABS(DX(I)) .GT. CUTLO) GO TO 85
 103 C
 104 C                                PREPARE FOR PHASE 2.
 105 C
 106       ASSIGN 70 TO NEXT
 107       GO TO 105
 108 C
 109 C                                PREPARE FOR PHASE 4.
 110 C
 111   100 I = J
 112       ASSIGN 110 TO NEXT
 113       SUM = (SUM / DX(I)) / DX(I)
 114   105 XMAX = ABS(DX(I))
 115       GO TO 115
 116 C
 117 C                   PHASE 2.  SUM IS SMALL.
 118 C                             SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
 119 C
 120    70 IF (ABS(DX(I)) .GT. CUTLO) GO TO 75
 121 C
 122 C                     COMMON CODE FOR PHASES 2 AND 4.
 123 C                     IN PHASE 4 SUM IS LARGE.  SCALE TO AVOID OVERFLOW.
 124 C
 125   110 IF (ABS(DX(I)) .LE. XMAX) GO TO 115
 126          SUM = ONE + SUM * (XMAX / DX(I))**2
 127          XMAX = ABS(DX(I))
 128          GO TO 200
 129 C
 130   115 SUM = SUM + (DX(I)/XMAX)**2
 131       GO TO 200
 132 C
 133 C                  PREPARE FOR PHASE 3.
 134 C
 135    75 SUM = (SUM * XMAX) * XMAX
 136 C
 137 C     FOR REAL OR D.P. SET HITEST = CUTHI/N
 138 C     FOR COMPLEX      SET HITEST = CUTHI/(2*N)
 139 C
 140    85 HITEST = CUTHI / N
 141 C
 142 C                   PHASE 3.  SUM IS MID-RANGE.  NO SCALING.
 143 C
 144       DO 95 J = I,NN,INCX
 145       IF (ABS(DX(J)) .GE. HITEST) GO TO 100
 146    95    SUM = SUM + DX(J)**2
 147       DNRM2 = SQRT(SUM)
 148       GO TO 300
 149 C
 150   200 CONTINUE
 151       I = I + INCX
 152       IF (I .LE. NN) GO TO 20
 153 C
 154 C              END OF MAIN LOOP.
 155 C
 156 C              COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
 157 C
 158       DNRM2 = XMAX * SQRT(SUM)
 159   300 CONTINUE
 160       RETURN
 161       END