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