INSTALL/slamchf77.f

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