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