home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / fortran / linpklib.arc / SNRM2.FOR < prev    next >
Encoding:
Text File  |  1984-01-20  |  4.4 KB  |  138 lines
  1.       REAL FUNCTION SNRM2 ( N, SX, INCX)
  2.       INTEGER          NEXT
  3.       REAL   SX(1),  CUTLO, CUTHI, HITEST, SUM, XMAX, ZERO, ONE
  4. C
  5. C   SMALL is used to prevent values below SMALL from being returned.
  6. C   See the end of the program for usage.
  7.       REAL SMALL
  8.       DATA SMALL / 8.43E-37 /
  9. C
  10.       DATA   ZERO, ONE /0.0E0, 1.0E0/
  11. C
  12. C     EUCLIDEAN NORM OF THE N-VECTOR STORED IN SX() 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
  19. C     FOUR PHASE METHOD     USING TWO BUILT-IN CONSTANTS THAT ARE
  20. C     HOPEFULLY APPLICABLE TO ALL MACHINES.
  21. C         CUTLO = MAXIMUM OF  SQRT(U/EPS)  OVER ALL KNOWN MACHINES.
  22. C         CUTHI = MINIMUM OF  SQRT(V)      OVER ALL KNOWN MACHINES.
  23. C     WHERE
  24. C         EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
  25. C         U   = SMALLEST POSITIVE NO.   (UNDERFLOW LIMIT)
  26. C         V   = LARGEST  NO.            (OVERFLOW  LIMIT)
  27. C
  28. C     BRIEF OUTLINE OF ALGORITHM..
  29. C
  30. C     PHASE 1    SCANS ZERO COMPONENTS.
  31. C     MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
  32. C     MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
  33. C     MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
  34. C     WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
  35. C
  36. C     VALUES FOR CUTLO AND CUTHI..
  37. C     FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
  38. C     DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
  39. C     CUTLO, S.P.   U/EPS = 2**(-102) FOR  HONEYWELL.  CLOSE SECONDS ARE
  40. C                   UNIVAC AND DEC AT 2**(-103)
  41. C                   THUS CUTLO = 2**(-51) = 4.44089E-16
  42. C     CUTHI, S.P.   V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
  43. C                   THUS CUTHI = 2**(63.5) = 1.30438E19
  44. C     CUTLO, D.P.   U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
  45. C                   THUS CUTLO = 2**(-33.5) = 8.23181D-11
  46. C     CUTHI, D.P.   SAME AS S.P.  CUTHI = 1.30438D19
  47. C     DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
  48. C     DATA CUTLO, CUTHI / 4.441E-16,  1.304E19 /
  49.       DATA CUTLO, CUTHI / 4.441E-16,  1.304E19 /
  50. C
  51.       IF(N .GT. 0) GO TO 10
  52.          SNRM2  = ZERO
  53.          GO TO 300
  54. C
  55.    10 ASSIGN 30 TO NEXT
  56.       SUM = ZERO
  57.       NN = N * INCX
  58. C                                                 BEGIN MAIN LOOP
  59.       I = 1
  60.    20    GO TO NEXT,(30, 50, 70, 110)
  61.    30 IF( ABS(SX(I)) .GT. CUTLO) GO TO 85
  62.       ASSIGN 50 TO NEXT
  63.       XMAX = ZERO
  64. C
  65. C                        PHASE 1.  SUM IS ZERO
  66. C
  67.    50 IF( SX(I) .EQ. ZERO) GO TO 200
  68.       IF( ABS(SX(I)) .GT. CUTLO) GO TO 85
  69. C
  70. C                                PREPARE FOR PHASE 2.
  71.       ASSIGN 70 TO NEXT
  72.       GO TO 105
  73. C
  74. C                                PREPARE FOR PHASE 4.
  75. C
  76.   100 I = J
  77.       ASSIGN 110 TO NEXT
  78.       SUM = (SUM / SX(I)) / SX(I)
  79.   105 XMAX = ABS(SX(I))
  80.       GO TO 115
  81. C
  82. C                   PHASE 2.  SUM IS SMALL.
  83. C                             SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
  84. C
  85.    70 IF( ABS(SX(I)) .GT. CUTLO ) GO TO 75
  86. C
  87. C                     COMMON CODE FOR PHASES 2 AND 4.
  88. C                     IN PHASE 4 SUM IS LARGE.  SCALE TO AVOID OVERFLOW.
  89. C
  90.   110 IF( ABS(SX(I)) .LE. XMAX ) GO TO 115
  91.          SUM = ONE + SUM * (XMAX / SX(I))**2
  92.          XMAX = ABS(SX(I))
  93.          GO TO 200
  94. C
  95.   115 SUM = SUM + (SX(I)/XMAX)**2
  96.       GO TO 200
  97. C
  98. C
  99. C                  PREPARE FOR PHASE 3.
  100. C
  101.    75 SUM = (SUM * XMAX) * XMAX
  102. C
  103. C
  104. C     FOR REAL OR D.P. SET HITEST = CUTHI/N
  105. C     FOR COMPLEX      SET HITEST = CUTHI/(2*N)
  106. C
  107.    85 HITEST = CUTHI/FLOAT( N )
  108. C
  109. C                   PHASE 3.  SUM IS MID-RANGE.  NO SCALING.
  110. C
  111.       DO 95 J =I,NN,INCX
  112.       IF(ABS(SX(J)) .GE. HITEST) GO TO 100
  113.    95    SUM = SUM + SX(J)**2
  114.       SNRM2 = SQRT( SUM )
  115.       GO TO 300
  116. C
  117.   200 CONTINUE
  118.       I = I + INCX
  119.       IF ( I .LE. NN ) GO TO 20
  120. C
  121. C              END OF MAIN LOOP.
  122. C
  123. C              COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
  124. C
  125.       SNRM2 = XMAX * SQRT(SUM)
  126.   300 CONTINUE
  127. C
  128. C *** 6 Jan 84, J. Fried : 
  129. C     This next line of FORTRAN was added to prevent values below 
  130. C     MS-FORTRAN's specified minimum from being returned.  Such 
  131. C     numbers have appeared during testing of the SQRDC routine 
  132. C     (test #9 in SQ).  The value of SMALL was taken from the 
  133. C     MS-FORTRAN manual.  Remember, this only applies to PC/MS-FORTRAN.
  134. C
  135.       IF ( SNRM2 .LE. SMALL ) SNRM2 = SMALL
  136.       RETURN
  137.       END
  138.