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