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