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