INSTALL/dlamchf77.f

   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.7.0) --
  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.7.0) --
 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.7.0) --
 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.7.0) --
 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.7.0) --
 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.7.0) --
 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