Annotation of rpl/lapack/lapack/dlamch.f, revision 1.1
1.1 ! bertrand 1: DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
! 2: *
! 3: * -- LAPACK auxiliary routine (version 3.2) --
! 4: * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
! 5: * November 2006
! 6: *
! 7: * .. Scalar Arguments ..
! 8: CHARACTER CMACH
! 9: * ..
! 10: *
! 11: * Purpose
! 12: * =======
! 13: *
! 14: * DLAMCH determines double precision machine parameters.
! 15: *
! 16: * Arguments
! 17: * =========
! 18: *
! 19: * CMACH (input) CHARACTER*1
! 20: * Specifies the value to be returned by DLAMCH:
! 21: * = 'E' or 'e', DLAMCH := eps
! 22: * = 'S' or 's , DLAMCH := sfmin
! 23: * = 'B' or 'b', DLAMCH := base
! 24: * = 'P' or 'p', DLAMCH := eps*base
! 25: * = 'N' or 'n', DLAMCH := t
! 26: * = 'R' or 'r', DLAMCH := rnd
! 27: * = 'M' or 'm', DLAMCH := emin
! 28: * = 'U' or 'u', DLAMCH := rmin
! 29: * = 'L' or 'l', DLAMCH := emax
! 30: * = 'O' or 'o', DLAMCH := rmax
! 31: *
! 32: * where
! 33: *
! 34: * eps = relative machine precision
! 35: * sfmin = safe minimum, such that 1/sfmin does not overflow
! 36: * base = base of the machine
! 37: * prec = eps*base
! 38: * t = number of (base) digits in the mantissa
! 39: * rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
! 40: * emin = minimum exponent before (gradual) underflow
! 41: * rmin = underflow threshold - base**(emin-1)
! 42: * emax = largest exponent before overflow
! 43: * rmax = overflow threshold - (base**emax)*(1-eps)
! 44: *
! 45: * =====================================================================
! 46: *
! 47: * .. Parameters ..
! 48: DOUBLE PRECISION ONE, ZERO
! 49: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
! 50: * ..
! 51: * .. Local Scalars ..
! 52: LOGICAL FIRST, LRND
! 53: INTEGER BETA, IMAX, IMIN, IT
! 54: DOUBLE PRECISION BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
! 55: $ RND, SFMIN, SMALL, T
! 56: * ..
! 57: * .. External Functions ..
! 58: LOGICAL LSAME
! 59: EXTERNAL LSAME
! 60: * ..
! 61: * .. External Subroutines ..
! 62: EXTERNAL DLAMC2
! 63: * ..
! 64: * .. Save statement ..
! 65: SAVE FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN,
! 66: $ EMAX, RMAX, PREC
! 67: * ..
! 68: * .. Data statements ..
! 69: DATA FIRST / .TRUE. /
! 70: * ..
! 71: * .. Executable Statements ..
! 72: *
! 73: IF( FIRST ) THEN
! 74: CALL DLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX )
! 75: BASE = BETA
! 76: T = IT
! 77: IF( LRND ) THEN
! 78: RND = ONE
! 79: EPS = ( BASE**( 1-IT ) ) / 2
! 80: ELSE
! 81: RND = ZERO
! 82: EPS = BASE**( 1-IT )
! 83: END IF
! 84: PREC = EPS*BASE
! 85: EMIN = IMIN
! 86: EMAX = IMAX
! 87: SFMIN = RMIN
! 88: SMALL = ONE / RMAX
! 89: IF( SMALL.GE.SFMIN ) THEN
! 90: *
! 91: * Use SMALL plus a bit, to avoid the possibility of rounding
! 92: * causing overflow when computing 1/sfmin.
! 93: *
! 94: SFMIN = SMALL*( ONE+EPS )
! 95: END IF
! 96: END IF
! 97: *
! 98: IF( LSAME( CMACH, 'E' ) ) THEN
! 99: RMACH = EPS
! 100: ELSE IF( LSAME( CMACH, 'S' ) ) THEN
! 101: RMACH = SFMIN
! 102: ELSE IF( LSAME( CMACH, 'B' ) ) THEN
! 103: RMACH = BASE
! 104: ELSE IF( LSAME( CMACH, 'P' ) ) THEN
! 105: RMACH = PREC
! 106: ELSE IF( LSAME( CMACH, 'N' ) ) THEN
! 107: RMACH = T
! 108: ELSE IF( LSAME( CMACH, 'R' ) ) THEN
! 109: RMACH = RND
! 110: ELSE IF( LSAME( CMACH, 'M' ) ) THEN
! 111: RMACH = EMIN
! 112: ELSE IF( LSAME( CMACH, 'U' ) ) THEN
! 113: RMACH = RMIN
! 114: ELSE IF( LSAME( CMACH, 'L' ) ) THEN
! 115: RMACH = EMAX
! 116: ELSE IF( LSAME( CMACH, 'O' ) ) THEN
! 117: RMACH = RMAX
! 118: END IF
! 119: *
! 120: DLAMCH = RMACH
! 121: FIRST = .FALSE.
! 122: RETURN
! 123: *
! 124: * End of DLAMCH
! 125: *
! 126: END
! 127: *
! 128: ************************************************************************
! 129: *
! 130: SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 )
! 131: *
! 132: * -- LAPACK auxiliary routine (version 3.2) --
! 133: * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
! 134: * November 2006
! 135: *
! 136: * .. Scalar Arguments ..
! 137: LOGICAL IEEE1, RND
! 138: INTEGER BETA, T
! 139: * ..
! 140: *
! 141: * Purpose
! 142: * =======
! 143: *
! 144: * DLAMC1 determines the machine parameters given by BETA, T, RND, and
! 145: * IEEE1.
! 146: *
! 147: * Arguments
! 148: * =========
! 149: *
! 150: * BETA (output) INTEGER
! 151: * The base of the machine.
! 152: *
! 153: * T (output) INTEGER
! 154: * The number of ( BETA ) digits in the mantissa.
! 155: *
! 156: * RND (output) LOGICAL
! 157: * Specifies whether proper rounding ( RND = .TRUE. ) or
! 158: * chopping ( RND = .FALSE. ) occurs in addition. This may not
! 159: * be a reliable guide to the way in which the machine performs
! 160: * its arithmetic.
! 161: *
! 162: * IEEE1 (output) LOGICAL
! 163: * Specifies whether rounding appears to be done in the IEEE
! 164: * 'round to nearest' style.
! 165: *
! 166: * Further Details
! 167: * ===============
! 168: *
! 169: * The routine is based on the routine ENVRON by Malcolm and
! 170: * incorporates suggestions by Gentleman and Marovich. See
! 171: *
! 172: * Malcolm M. A. (1972) Algorithms to reveal properties of
! 173: * floating-point arithmetic. Comms. of the ACM, 15, 949-951.
! 174: *
! 175: * Gentleman W. M. and Marovich S. B. (1974) More on algorithms
! 176: * that reveal properties of floating point arithmetic units.
! 177: * Comms. of the ACM, 17, 276-277.
! 178: *
! 179: * =====================================================================
! 180: *
! 181: * .. Local Scalars ..
! 182: LOGICAL FIRST, LIEEE1, LRND
! 183: INTEGER LBETA, LT
! 184: DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2
! 185: * ..
! 186: * .. External Functions ..
! 187: DOUBLE PRECISION DLAMC3
! 188: EXTERNAL DLAMC3
! 189: * ..
! 190: * .. Save statement ..
! 191: SAVE FIRST, LIEEE1, LBETA, LRND, LT
! 192: * ..
! 193: * .. Data statements ..
! 194: DATA FIRST / .TRUE. /
! 195: * ..
! 196: * .. Executable Statements ..
! 197: *
! 198: IF( FIRST ) THEN
! 199: ONE = 1
! 200: *
! 201: * LBETA, LIEEE1, LT and LRND are the local values of BETA,
! 202: * IEEE1, T and RND.
! 203: *
! 204: * Throughout this routine we use the function DLAMC3 to ensure
! 205: * that relevant values are stored and not held in registers, or
! 206: * are not affected by optimizers.
! 207: *
! 208: * Compute a = 2.0**m with the smallest positive integer m such
! 209: * that
! 210: *
! 211: * fl( a + 1.0 ) = a.
! 212: *
! 213: A = 1
! 214: C = 1
! 215: *
! 216: *+ WHILE( C.EQ.ONE )LOOP
! 217: 10 CONTINUE
! 218: IF( C.EQ.ONE ) THEN
! 219: A = 2*A
! 220: C = DLAMC3( A, ONE )
! 221: C = DLAMC3( C, -A )
! 222: GO TO 10
! 223: END IF
! 224: *+ END WHILE
! 225: *
! 226: * Now compute b = 2.0**m with the smallest positive integer m
! 227: * such that
! 228: *
! 229: * fl( a + b ) .gt. a.
! 230: *
! 231: B = 1
! 232: C = DLAMC3( A, B )
! 233: *
! 234: *+ WHILE( C.EQ.A )LOOP
! 235: 20 CONTINUE
! 236: IF( C.EQ.A ) THEN
! 237: B = 2*B
! 238: C = DLAMC3( A, B )
! 239: GO TO 20
! 240: END IF
! 241: *+ END WHILE
! 242: *
! 243: * Now compute the base. a and c are neighbouring floating point
! 244: * numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
! 245: * their difference is beta. Adding 0.25 to c is to ensure that it
! 246: * is truncated to beta and not ( beta - 1 ).
! 247: *
! 248: QTR = ONE / 4
! 249: SAVEC = C
! 250: C = DLAMC3( C, -A )
! 251: LBETA = C + QTR
! 252: *
! 253: * Now determine whether rounding or chopping occurs, by adding a
! 254: * bit less than beta/2 and a bit more than beta/2 to a.
! 255: *
! 256: B = LBETA
! 257: F = DLAMC3( B / 2, -B / 100 )
! 258: C = DLAMC3( F, A )
! 259: IF( C.EQ.A ) THEN
! 260: LRND = .TRUE.
! 261: ELSE
! 262: LRND = .FALSE.
! 263: END IF
! 264: F = DLAMC3( B / 2, B / 100 )
! 265: C = DLAMC3( F, A )
! 266: IF( ( LRND ) .AND. ( C.EQ.A ) )
! 267: $ LRND = .FALSE.
! 268: *
! 269: * Try and decide whether rounding is done in the IEEE 'round to
! 270: * nearest' style. B/2 is half a unit in the last place of the two
! 271: * numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
! 272: * zero, and SAVEC is odd. Thus adding B/2 to A should not change
! 273: * A, but adding B/2 to SAVEC should change SAVEC.
! 274: *
! 275: T1 = DLAMC3( B / 2, A )
! 276: T2 = DLAMC3( B / 2, SAVEC )
! 277: LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
! 278: *
! 279: * Now find the mantissa, t. It should be the integer part of
! 280: * log to the base beta of a, however it is safer to determine t
! 281: * by powering. So we find t as the smallest positive integer for
! 282: * which
! 283: *
! 284: * fl( beta**t + 1.0 ) = 1.0.
! 285: *
! 286: LT = 0
! 287: A = 1
! 288: C = 1
! 289: *
! 290: *+ WHILE( C.EQ.ONE )LOOP
! 291: 30 CONTINUE
! 292: IF( C.EQ.ONE ) THEN
! 293: LT = LT + 1
! 294: A = A*LBETA
! 295: C = DLAMC3( A, ONE )
! 296: C = DLAMC3( C, -A )
! 297: GO TO 30
! 298: END IF
! 299: *+ END WHILE
! 300: *
! 301: END IF
! 302: *
! 303: BETA = LBETA
! 304: T = LT
! 305: RND = LRND
! 306: IEEE1 = LIEEE1
! 307: FIRST = .FALSE.
! 308: RETURN
! 309: *
! 310: * End of DLAMC1
! 311: *
! 312: END
! 313: *
! 314: ************************************************************************
! 315: *
! 316: SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
! 317: *
! 318: * -- LAPACK auxiliary routine (version 3.2) --
! 319: * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
! 320: * November 2006
! 321: *
! 322: * .. Scalar Arguments ..
! 323: LOGICAL RND
! 324: INTEGER BETA, EMAX, EMIN, T
! 325: DOUBLE PRECISION EPS, RMAX, RMIN
! 326: * ..
! 327: *
! 328: * Purpose
! 329: * =======
! 330: *
! 331: * DLAMC2 determines the machine parameters specified in its argument
! 332: * list.
! 333: *
! 334: * Arguments
! 335: * =========
! 336: *
! 337: * BETA (output) INTEGER
! 338: * The base of the machine.
! 339: *
! 340: * T (output) INTEGER
! 341: * The number of ( BETA ) digits in the mantissa.
! 342: *
! 343: * RND (output) LOGICAL
! 344: * Specifies whether proper rounding ( RND = .TRUE. ) or
! 345: * chopping ( RND = .FALSE. ) occurs in addition. This may not
! 346: * be a reliable guide to the way in which the machine performs
! 347: * its arithmetic.
! 348: *
! 349: * EPS (output) DOUBLE PRECISION
! 350: * The smallest positive number such that
! 351: *
! 352: * fl( 1.0 - EPS ) .LT. 1.0,
! 353: *
! 354: * where fl denotes the computed value.
! 355: *
! 356: * EMIN (output) INTEGER
! 357: * The minimum exponent before (gradual) underflow occurs.
! 358: *
! 359: * RMIN (output) DOUBLE PRECISION
! 360: * The smallest normalized number for the machine, given by
! 361: * BASE**( EMIN - 1 ), where BASE is the floating point value
! 362: * of BETA.
! 363: *
! 364: * EMAX (output) INTEGER
! 365: * The maximum exponent before overflow occurs.
! 366: *
! 367: * RMAX (output) DOUBLE PRECISION
! 368: * The largest positive number for the machine, given by
! 369: * BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
! 370: * value of BETA.
! 371: *
! 372: * Further Details
! 373: * ===============
! 374: *
! 375: * The computation of EPS is based on a routine PARANOIA by
! 376: * W. Kahan of the University of California at Berkeley.
! 377: *
! 378: * =====================================================================
! 379: *
! 380: * .. Local Scalars ..
! 381: LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND
! 382: INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
! 383: $ NGNMIN, NGPMIN
! 384: DOUBLE PRECISION A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
! 385: $ SIXTH, SMALL, THIRD, TWO, ZERO
! 386: * ..
! 387: * .. External Functions ..
! 388: DOUBLE PRECISION DLAMC3
! 389: EXTERNAL DLAMC3
! 390: * ..
! 391: * .. External Subroutines ..
! 392: EXTERNAL DLAMC1, DLAMC4, DLAMC5
! 393: * ..
! 394: * .. Intrinsic Functions ..
! 395: INTRINSIC ABS, MAX, MIN
! 396: * ..
! 397: * .. Save statement ..
! 398: SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,
! 399: $ LRMIN, LT
! 400: * ..
! 401: * .. Data statements ..
! 402: DATA FIRST / .TRUE. / , IWARN / .FALSE. /
! 403: * ..
! 404: * .. Executable Statements ..
! 405: *
! 406: IF( FIRST ) THEN
! 407: ZERO = 0
! 408: ONE = 1
! 409: TWO = 2
! 410: *
! 411: * LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
! 412: * BETA, T, RND, EPS, EMIN and RMIN.
! 413: *
! 414: * Throughout this routine we use the function DLAMC3 to ensure
! 415: * that relevant values are stored and not held in registers, or
! 416: * are not affected by optimizers.
! 417: *
! 418: * DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
! 419: *
! 420: CALL DLAMC1( LBETA, LT, LRND, LIEEE1 )
! 421: *
! 422: * Start to find EPS.
! 423: *
! 424: B = LBETA
! 425: A = B**( -LT )
! 426: LEPS = A
! 427: *
! 428: * Try some tricks to see whether or not this is the correct EPS.
! 429: *
! 430: B = TWO / 3
! 431: HALF = ONE / 2
! 432: SIXTH = DLAMC3( B, -HALF )
! 433: THIRD = DLAMC3( SIXTH, SIXTH )
! 434: B = DLAMC3( THIRD, -HALF )
! 435: B = DLAMC3( B, SIXTH )
! 436: B = ABS( B )
! 437: IF( B.LT.LEPS )
! 438: $ B = LEPS
! 439: *
! 440: LEPS = 1
! 441: *
! 442: *+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
! 443: 10 CONTINUE
! 444: IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN
! 445: LEPS = B
! 446: C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
! 447: C = DLAMC3( HALF, -C )
! 448: B = DLAMC3( HALF, C )
! 449: C = DLAMC3( HALF, -B )
! 450: B = DLAMC3( HALF, C )
! 451: GO TO 10
! 452: END IF
! 453: *+ END WHILE
! 454: *
! 455: IF( A.LT.LEPS )
! 456: $ LEPS = A
! 457: *
! 458: * Computation of EPS complete.
! 459: *
! 460: * Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
! 461: * Keep dividing A by BETA until (gradual) underflow occurs. This
! 462: * is detected when we cannot recover the previous A.
! 463: *
! 464: RBASE = ONE / LBETA
! 465: SMALL = ONE
! 466: DO 20 I = 1, 3
! 467: SMALL = DLAMC3( SMALL*RBASE, ZERO )
! 468: 20 CONTINUE
! 469: A = DLAMC3( ONE, SMALL )
! 470: CALL DLAMC4( NGPMIN, ONE, LBETA )
! 471: CALL DLAMC4( NGNMIN, -ONE, LBETA )
! 472: CALL DLAMC4( GPMIN, A, LBETA )
! 473: CALL DLAMC4( GNMIN, -A, LBETA )
! 474: IEEE = .FALSE.
! 475: *
! 476: IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN
! 477: IF( NGPMIN.EQ.GPMIN ) THEN
! 478: LEMIN = NGPMIN
! 479: * ( Non twos-complement machines, no gradual underflow;
! 480: * e.g., VAX )
! 481: ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
! 482: LEMIN = NGPMIN - 1 + LT
! 483: IEEE = .TRUE.
! 484: * ( Non twos-complement machines, with gradual underflow;
! 485: * e.g., IEEE standard followers )
! 486: ELSE
! 487: LEMIN = MIN( NGPMIN, GPMIN )
! 488: * ( A guess; no known machine )
! 489: IWARN = .TRUE.
! 490: END IF
! 491: *
! 492: ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
! 493: IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
! 494: LEMIN = MAX( NGPMIN, NGNMIN )
! 495: * ( Twos-complement machines, no gradual underflow;
! 496: * e.g., CYBER 205 )
! 497: ELSE
! 498: LEMIN = MIN( NGPMIN, NGNMIN )
! 499: * ( A guess; no known machine )
! 500: IWARN = .TRUE.
! 501: END IF
! 502: *
! 503: ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND.
! 504: $ ( GPMIN.EQ.GNMIN ) ) THEN
! 505: IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN
! 506: LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
! 507: * ( Twos-complement machines with gradual underflow;
! 508: * no known machine )
! 509: ELSE
! 510: LEMIN = MIN( NGPMIN, NGNMIN )
! 511: * ( A guess; no known machine )
! 512: IWARN = .TRUE.
! 513: END IF
! 514: *
! 515: ELSE
! 516: LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
! 517: * ( A guess; no known machine )
! 518: IWARN = .TRUE.
! 519: END IF
! 520: FIRST = .FALSE.
! 521: ***
! 522: * Comment out this if block if EMIN is ok
! 523: IF( IWARN ) THEN
! 524: FIRST = .TRUE.
! 525: WRITE( 6, FMT = 9999 )LEMIN
! 526: END IF
! 527: ***
! 528: *
! 529: * Assume IEEE arithmetic if we found denormalised numbers above,
! 530: * or if arithmetic seems to round in the IEEE style, determined
! 531: * in routine DLAMC1. A true IEEE machine should have both things
! 532: * true; however, faulty machines may have one or the other.
! 533: *
! 534: IEEE = IEEE .OR. LIEEE1
! 535: *
! 536: * Compute RMIN by successive division by BETA. We could compute
! 537: * RMIN as BASE**( EMIN - 1 ), but some machines underflow during
! 538: * this computation.
! 539: *
! 540: LRMIN = 1
! 541: DO 30 I = 1, 1 - LEMIN
! 542: LRMIN = DLAMC3( LRMIN*RBASE, ZERO )
! 543: 30 CONTINUE
! 544: *
! 545: * Finally, call DLAMC5 to compute EMAX and RMAX.
! 546: *
! 547: CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )
! 548: END IF
! 549: *
! 550: BETA = LBETA
! 551: T = LT
! 552: RND = LRND
! 553: EPS = LEPS
! 554: EMIN = LEMIN
! 555: RMIN = LRMIN
! 556: EMAX = LEMAX
! 557: RMAX = LRMAX
! 558: *
! 559: RETURN
! 560: *
! 561: 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
! 562: $ ' EMIN = ', I8, /
! 563: $ ' If, after inspection, the value EMIN looks',
! 564: $ ' acceptable please comment out ',
! 565: $ / ' the IF block as marked within the code of routine',
! 566: $ ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / )
! 567: *
! 568: * End of DLAMC2
! 569: *
! 570: END
! 571: *
! 572: ************************************************************************
! 573: *
! 574: DOUBLE PRECISION FUNCTION DLAMC3( A, B )
! 575: *
! 576: * -- LAPACK auxiliary routine (version 3.2) --
! 577: * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
! 578: * November 2006
! 579: *
! 580: * .. Scalar Arguments ..
! 581: DOUBLE PRECISION A, B
! 582: * ..
! 583: *
! 584: * Purpose
! 585: * =======
! 586: *
! 587: * DLAMC3 is intended to force A and B to be stored prior to doing
! 588: * the addition of A and B , for use in situations where optimizers
! 589: * might hold one of these in a register.
! 590: *
! 591: * Arguments
! 592: * =========
! 593: *
! 594: * A (input) DOUBLE PRECISION
! 595: * B (input) DOUBLE PRECISION
! 596: * The values A and B.
! 597: *
! 598: * =====================================================================
! 599: *
! 600: * .. Executable Statements ..
! 601: *
! 602: DLAMC3 = A + B
! 603: *
! 604: RETURN
! 605: *
! 606: * End of DLAMC3
! 607: *
! 608: END
! 609: *
! 610: ************************************************************************
! 611: *
! 612: SUBROUTINE DLAMC4( EMIN, START, BASE )
! 613: *
! 614: * -- LAPACK auxiliary routine (version 3.2) --
! 615: * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
! 616: * November 2006
! 617: *
! 618: * .. Scalar Arguments ..
! 619: INTEGER BASE, EMIN
! 620: DOUBLE PRECISION START
! 621: * ..
! 622: *
! 623: * Purpose
! 624: * =======
! 625: *
! 626: * DLAMC4 is a service routine for DLAMC2.
! 627: *
! 628: * Arguments
! 629: * =========
! 630: *
! 631: * EMIN (output) INTEGER
! 632: * The minimum exponent before (gradual) underflow, computed by
! 633: * setting A = START and dividing by BASE until the previous A
! 634: * can not be recovered.
! 635: *
! 636: * START (input) DOUBLE PRECISION
! 637: * The starting point for determining EMIN.
! 638: *
! 639: * BASE (input) INTEGER
! 640: * The base of the machine.
! 641: *
! 642: * =====================================================================
! 643: *
! 644: * .. Local Scalars ..
! 645: INTEGER I
! 646: DOUBLE PRECISION A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
! 647: * ..
! 648: * .. External Functions ..
! 649: DOUBLE PRECISION DLAMC3
! 650: EXTERNAL DLAMC3
! 651: * ..
! 652: * .. Executable Statements ..
! 653: *
! 654: A = START
! 655: ONE = 1
! 656: RBASE = ONE / BASE
! 657: ZERO = 0
! 658: EMIN = 1
! 659: B1 = DLAMC3( A*RBASE, ZERO )
! 660: C1 = A
! 661: C2 = A
! 662: D1 = A
! 663: D2 = A
! 664: *+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
! 665: * $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP
! 666: 10 CONTINUE
! 667: IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND.
! 668: $ ( D2.EQ.A ) ) THEN
! 669: EMIN = EMIN - 1
! 670: A = B1
! 671: B1 = DLAMC3( A / BASE, ZERO )
! 672: C1 = DLAMC3( B1*BASE, ZERO )
! 673: D1 = ZERO
! 674: DO 20 I = 1, BASE
! 675: D1 = D1 + B1
! 676: 20 CONTINUE
! 677: B2 = DLAMC3( A*RBASE, ZERO )
! 678: C2 = DLAMC3( B2 / RBASE, ZERO )
! 679: D2 = ZERO
! 680: DO 30 I = 1, BASE
! 681: D2 = D2 + B2
! 682: 30 CONTINUE
! 683: GO TO 10
! 684: END IF
! 685: *+ END WHILE
! 686: *
! 687: RETURN
! 688: *
! 689: * End of DLAMC4
! 690: *
! 691: END
! 692: *
! 693: ************************************************************************
! 694: *
! 695: SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
! 696: *
! 697: * -- LAPACK auxiliary routine (version 3.2) --
! 698: * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
! 699: * November 2006
! 700: *
! 701: * .. Scalar Arguments ..
! 702: LOGICAL IEEE
! 703: INTEGER BETA, EMAX, EMIN, P
! 704: DOUBLE PRECISION RMAX
! 705: * ..
! 706: *
! 707: * Purpose
! 708: * =======
! 709: *
! 710: * DLAMC5 attempts to compute RMAX, the largest machine floating-point
! 711: * number, without overflow. It assumes that EMAX + abs(EMIN) sum
! 712: * approximately to a power of 2. It will fail on machines where this
! 713: * assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
! 714: * EMAX = 28718). It will also fail if the value supplied for EMIN is
! 715: * too large (i.e. too close to zero), probably with overflow.
! 716: *
! 717: * Arguments
! 718: * =========
! 719: *
! 720: * BETA (input) INTEGER
! 721: * The base of floating-point arithmetic.
! 722: *
! 723: * P (input) INTEGER
! 724: * The number of base BETA digits in the mantissa of a
! 725: * floating-point value.
! 726: *
! 727: * EMIN (input) INTEGER
! 728: * The minimum exponent before (gradual) underflow.
! 729: *
! 730: * IEEE (input) LOGICAL
! 731: * A logical flag specifying whether or not the arithmetic
! 732: * system is thought to comply with the IEEE standard.
! 733: *
! 734: * EMAX (output) INTEGER
! 735: * The largest exponent before overflow
! 736: *
! 737: * RMAX (output) DOUBLE PRECISION
! 738: * The largest machine floating-point number.
! 739: *
! 740: * =====================================================================
! 741: *
! 742: * .. Parameters ..
! 743: DOUBLE PRECISION ZERO, ONE
! 744: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
! 745: * ..
! 746: * .. Local Scalars ..
! 747: INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
! 748: DOUBLE PRECISION OLDY, RECBAS, Y, Z
! 749: * ..
! 750: * .. External Functions ..
! 751: DOUBLE PRECISION DLAMC3
! 752: EXTERNAL DLAMC3
! 753: * ..
! 754: * .. Intrinsic Functions ..
! 755: INTRINSIC MOD
! 756: * ..
! 757: * .. Executable Statements ..
! 758: *
! 759: * First compute LEXP and UEXP, two powers of 2 that bound
! 760: * abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
! 761: * approximately to the bound that is closest to abs(EMIN).
! 762: * (EMAX is the exponent of the required number RMAX).
! 763: *
! 764: LEXP = 1
! 765: EXBITS = 1
! 766: 10 CONTINUE
! 767: TRY = LEXP*2
! 768: IF( TRY.LE.( -EMIN ) ) THEN
! 769: LEXP = TRY
! 770: EXBITS = EXBITS + 1
! 771: GO TO 10
! 772: END IF
! 773: IF( LEXP.EQ.-EMIN ) THEN
! 774: UEXP = LEXP
! 775: ELSE
! 776: UEXP = TRY
! 777: EXBITS = EXBITS + 1
! 778: END IF
! 779: *
! 780: * Now -LEXP is less than or equal to EMIN, and -UEXP is greater
! 781: * than or equal to EMIN. EXBITS is the number of bits needed to
! 782: * store the exponent.
! 783: *
! 784: IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
! 785: EXPSUM = 2*LEXP
! 786: ELSE
! 787: EXPSUM = 2*UEXP
! 788: END IF
! 789: *
! 790: * EXPSUM is the exponent range, approximately equal to
! 791: * EMAX - EMIN + 1 .
! 792: *
! 793: EMAX = EXPSUM + EMIN - 1
! 794: NBITS = 1 + EXBITS + P
! 795: *
! 796: * NBITS is the total number of bits needed to store a
! 797: * floating-point number.
! 798: *
! 799: IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
! 800: *
! 801: * Either there are an odd number of bits used to store a
! 802: * floating-point number, which is unlikely, or some bits are
! 803: * not used in the representation of numbers, which is possible,
! 804: * (e.g. Cray machines) or the mantissa has an implicit bit,
! 805: * (e.g. IEEE machines, Dec Vax machines), which is perhaps the
! 806: * most likely. We have to assume the last alternative.
! 807: * If this is true, then we need to reduce EMAX by one because
! 808: * there must be some way of representing zero in an implicit-bit
! 809: * system. On machines like Cray, we are reducing EMAX by one
! 810: * unnecessarily.
! 811: *
! 812: EMAX = EMAX - 1
! 813: END IF
! 814: *
! 815: IF( IEEE ) THEN
! 816: *
! 817: * Assume we are on an IEEE machine which reserves one exponent
! 818: * for infinity and NaN.
! 819: *
! 820: EMAX = EMAX - 1
! 821: END IF
! 822: *
! 823: * Now create RMAX, the largest machine number, which should
! 824: * be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
! 825: *
! 826: * First compute 1.0 - BETA**(-P), being careful that the
! 827: * result is less than 1.0 .
! 828: *
! 829: RECBAS = ONE / BETA
! 830: Z = BETA - ONE
! 831: Y = ZERO
! 832: DO 20 I = 1, P
! 833: Z = Z*RECBAS
! 834: IF( Y.LT.ONE )
! 835: $ OLDY = Y
! 836: Y = DLAMC3( Y, Z )
! 837: 20 CONTINUE
! 838: IF( Y.GE.ONE )
! 839: $ Y = OLDY
! 840: *
! 841: * Now multiply by BETA**EMAX to get RMAX.
! 842: *
! 843: DO 30 I = 1, EMAX
! 844: Y = DLAMC3( Y*BETA, ZERO )
! 845: 30 CONTINUE
! 846: *
! 847: RMAX = Y
! 848: RETURN
! 849: *
! 850: * End of DLAMC5
! 851: *
! 852: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlamch.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack