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