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 853