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