amos/d1mach.f

   1       DOUBLE PRECISION FUNCTION D1MACH(I)
   2       INTEGER I
   3 C
   4 C  DOUBLE-PRECISION MACHINE CONSTANTS
   5 C  D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
   6 C  D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
   7 C  D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
   8 C  D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
   9 C  D1MACH( 5) = LOG10(B)
  10 C
  11       INTEGER SMALL(2)
  12       INTEGER LARGE(2)
  13       INTEGER RIGHT(2)
  14       INTEGER DIVER(2)
  15       INTEGER LOG10(2)
  16       INTEGER SC, CRAY1(38), J
  17       COMMON /D9MACH/ CRAY1
  18       SAVE SMALL, LARGE, RIGHT, DIVER, LOG10, SC
  19       DOUBLE PRECISION DMACH(5)
  20       EQUIVALENCE (DMACH(1),SMALL(1))
  21       EQUIVALENCE (DMACH(2),LARGE(1))
  22       EQUIVALENCE (DMACH(3),RIGHT(1))
  23       EQUIVALENCE (DMACH(4),DIVER(1))
  24       EQUIVALENCE (DMACH(5),LOG10(1))
  25 C  THIS VERSION ADAPTS AUTOMATICALLY TO MOST CURRENT MACHINES.
  26 C  R1MACH CAN HANDLE AUTO-DOUBLE COMPILING, BUT THIS VERSION OF
  27 C  D1MACH DOES NOT, BECAUSE WE DO NOT HAVE QUAD CONSTANTS FOR
  28 C  MANY MACHINES YET.
  29 C  TO COMPILE ON OLDER MACHINES, ADD A C IN COLUMN 1
  30 C  ON THE NEXT LINE
  31       DATA SC/0/
  32 C  AND REMOVE THE C FROM COLUMN 1 IN ONE OF THE SECTIONS BELOW.
  33 C  CONSTANTS FOR EVEN OLDER MACHINES CAN BE OBTAINED BY
  34 C          mail netlib@research.bell-labs.com
  35 C          send old1mach from blas
  36 C  PLEASE SEND CORRECTIONS TO dmg OR ehg@bell-labs.com.
  37 C
  38 C     MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES.
  39 C      DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 /
  40 C      DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 /
  41 C      DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 /
  42 C      DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 /
  43 C      DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 /, SC/987/
  44 C
  45 C     MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING
  46 C     32-BIT INTEGERS.
  47 C      DATA SMALL(1),SMALL(2) /    8388608,           0 /
  48 C      DATA LARGE(1),LARGE(2) / 2147483647,          -1 /
  49 C      DATA RIGHT(1),RIGHT(2) /  612368384,           0 /
  50 C      DATA DIVER(1),DIVER(2) /  620756992,           0 /
  51 C      DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 /, SC/987/
  52 C
  53 C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES.
  54 C      DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 /
  55 C      DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 /
  56 C      DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 /
  57 C      DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 /
  58 C      DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 /, SC/987/
  59 C
  60 C     ON FIRST CALL, IF NO DATA UNCOMMENTED, TEST MACHINE TYPES.
  61       IF (SC .NE. 987) THEN
  62          DMACH(1) = 1.D13
  63          IF (      SMALL(1) .EQ. 1117925532
  64      *       .AND. SMALL(2) .EQ. -448790528) THEN
  65 *           *** IEEE BIG ENDIAN ***
  66             SMALL(1) = 1048576
  67             SMALL(2) = 0
  68             LARGE(1) = 2146435071
  69             LARGE(2) = -1
  70             RIGHT(1) = 1017118720
  71             RIGHT(2) = 0
  72             DIVER(1) = 1018167296
  73             DIVER(2) = 0
  74             LOG10(1) = 1070810131
  75             LOG10(2) = 1352628735
  76          ELSE IF ( SMALL(2) .EQ. 1117925532
  77      *       .AND. SMALL(1) .EQ. -448790528) THEN
  78 *           *** IEEE LITTLE ENDIAN ***
  79             SMALL(2) = 1048576
  80             SMALL(1) = 0
  81             LARGE(2) = 2146435071
  82             LARGE(1) = -1
  83             RIGHT(2) = 1017118720
  84             RIGHT(1) = 0
  85             DIVER(2) = 1018167296
  86             DIVER(1) = 0
  87             LOG10(2) = 1070810131
  88             LOG10(1) = 1352628735
  89          ELSE IF ( SMALL(1) .EQ. -2065213935
  90      *       .AND. SMALL(2) .EQ. 10752) THEN
  91 *               *** VAX WITH D_FLOATING ***
  92             SMALL(1) = 128
  93             SMALL(2) = 0
  94             LARGE(1) = -32769
  95             LARGE(2) = -1
  96             RIGHT(1) = 9344
  97             RIGHT(2) = 0
  98             DIVER(1) = 9472
  99             DIVER(2) = 0
 100             LOG10(1) = 546979738
 101             LOG10(2) = -805796613
 102          ELSE IF ( SMALL(1) .EQ. 1267827943
 103      *       .AND. SMALL(2) .EQ. 704643072) THEN
 104 *               *** IBM MAINFRAME ***
 105             SMALL(1) = 1048576
 106             SMALL(2) = 0
 107             LARGE(1) = 2147483647
 108             LARGE(2) = -1
 109             RIGHT(1) = 856686592
 110             RIGHT(2) = 0
 111             DIVER(1) = 873463808
 112             DIVER(2) = 0
 113             LOG10(1) = 1091781651
 114             LOG10(2) = 1352628735
 115          ELSE IF ( SMALL(1) .EQ. 1120022684
 116      *       .AND. SMALL(2) .EQ. -448790528) THEN
 117 *           *** CONVEX C-1 ***
 118             SMALL(1) = 1048576
 119             SMALL(2) = 0
 120             LARGE(1) = 2147483647
 121             LARGE(2) = -1
 122             RIGHT(1) = 1019215872
 123             RIGHT(2) = 0
 124             DIVER(1) = 1020264448
 125             DIVER(2) = 0
 126             LOG10(1) = 1072907283
 127             LOG10(2) = 1352628735
 128          ELSE IF ( SMALL(1) .EQ. 815547074
 129      *       .AND. SMALL(2) .EQ. 58688) THEN
 130 *           *** VAX G-FLOATING ***
 131             SMALL(1) = 16
 132             SMALL(2) = 0
 133             LARGE(1) = -32769
 134             LARGE(2) = -1
 135             RIGHT(1) = 15552
 136             RIGHT(2) = 0
 137             DIVER(1) = 15568
 138             DIVER(2) = 0
 139             LOG10(1) = 1142112243
 140             LOG10(2) = 2046775455
 141          ELSE
 142             DMACH(2) = 1.D27 + 1
 143             DMACH(3) = 1.D27
 144             LARGE(2) = LARGE(2) - RIGHT(2)
 145             IF (LARGE(2) .EQ. 64 .AND. SMALL(2) .EQ. 0) THEN
 146                CRAY1(1) = 67291416
 147                DO 10 J = 1, 20
 148                   CRAY1(J+1) = CRAY1(J) + CRAY1(J)
 149  10               CONTINUE
 150                CRAY1(22) = CRAY1(21) + 321322
 151                DO 20 J = 22, 37
 152                   CRAY1(J+1) = CRAY1(J) + CRAY1(J)
 153  20               CONTINUE
 154                IF (CRAY1(38) .EQ. SMALL(1)) THEN
 155 *                  *** CRAY ***
 156                   CALL I1MCRY(SMALL(1), J, 8285, 8388608, 0)
 157                   SMALL(2) = 0
 158                   CALL I1MCRY(LARGE(1), J, 24574, 16777215, 16777215)
 159                   CALL I1MCRY(LARGE(2), J, 0, 16777215, 16777214)
 160                   CALL I1MCRY(RIGHT(1), J, 16291, 8388608, 0)
 161                   RIGHT(2) = 0
 162                   CALL I1MCRY(DIVER(1), J, 16292, 8388608, 0)
 163                   DIVER(2) = 0
 164                   CALL I1MCRY(LOG10(1), J, 16383, 10100890, 8715215)
 165                   CALL I1MCRY(LOG10(2), J, 0, 16226447, 9001388)
 166                ELSE
 167                   WRITE(*,9000)
 168                   STOP 779
 169                   END IF
 170             ELSE
 171                WRITE(*,9000)
 172                STOP 779
 173                END IF
 174             END IF
 175          SC = 987
 176          END IF
 177 *    SANITY CHECK
 178       IF (DMACH(4) .GE. 1.0D0) STOP 778
 179       IF (I .LT. 1 .OR. I .GT. 5) THEN
 180          WRITE(*,*) 'D1MACH(I): I =',I,' is out of bounds.'
 181          STOP
 182          END IF
 183       D1MACH = DMACH(I)
 184       RETURN
 185  9000 FORMAT(/' Adjust D1MACH by uncommenting data statements'/
 186      *' appropriate for your machine.')
 187 * /* Standard C source for D1MACH -- remove the * in column 1 */
 188 *#include <stdio.h>
 189 *#include <float.h>
 190 *#include <math.h>
 191 *double d1mach_(long *i)
 192 *{
 193 *       switch(*i){
 194 *         case 1: return DBL_MIN;
 195 *         case 2: return DBL_MAX;
 196 *         case 3: return DBL_EPSILON/FLT_RADIX;
 197 *         case 4: return DBL_EPSILON;
 198 *         case 5: return log10(FLT_RADIX);
 199 *         }
 200 *       fprintf(stderr, "invalid argument: d1mach(%ld)\n", *i);
 201 *       exit(1); return 0; /* some compilers demand return values */
 202 *}
 203       END
 204       SUBROUTINE I1MCRY(A, A1, B, C, D)
 205 **** SPECIAL COMPUTATION FOR OLD CRAY MACHINES ****
 206       INTEGER A, A1, B, C, D
 207       A1 = 16777216*B + C
 208       A = 16777216*A1 + D
 209       END