LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slamchf77.f
Go to the documentation of this file.
00001 *> \brief \b SLAMCHF77 deprecated
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *      REAL FUNCTION SLAMCH( CMACH )
00012 *
00013 *     .. Scalar Arguments ..
00014 *      CHARACTER          CMACH
00015 *     ..
00016 *  
00017 *
00018 *> \par Purpose:
00019 *  =============
00020 *>
00021 *> \verbatim
00022 *>
00023 *> SLAMCH determines single precision machine parameters.
00024 *> \endverbatim
00025 *
00026 *  Arguments:
00027 *  ==========
00028 *
00029 *> \param[in] CMACH
00030 *> \verbatim
00031 *>          Specifies the value to be returned by SLAMCH:
00032 *>          = 'E' or 'e',   SLAMCH := eps
00033 *>          = 'S' or 's ,   SLAMCH := sfmin
00034 *>          = 'B' or 'b',   SLAMCH := base
00035 *>          = 'P' or 'p',   SLAMCH := eps*base
00036 *>          = 'N' or 'n',   SLAMCH := t
00037 *>          = 'R' or 'r',   SLAMCH := rnd
00038 *>          = 'M' or 'm',   SLAMCH := emin
00039 *>          = 'U' or 'u',   SLAMCH := rmin
00040 *>          = 'L' or 'l',   SLAMCH := emax
00041 *>          = 'O' or 'o',   SLAMCH := rmax
00042 *>          where
00043 *>          eps   = relative machine precision
00044 *>          sfmin = safe minimum, such that 1/sfmin does not overflow
00045 *>          base  = base of the machine
00046 *>          prec  = eps*base
00047 *>          t     = number of (base) digits in the mantissa
00048 *>          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise
00049 *>          emin  = minimum exponent before (gradual) underflow
00050 *>          rmin  = underflow threshold - base**(emin-1)
00051 *>          emax  = largest exponent before overflow
00052 *>          rmax  = overflow threshold  - (base**emax)*(1-eps)
00053 *> \endverbatim
00054 *
00055 *  Authors:
00056 *  ========
00057 *
00058 *> \author Univ. of Tennessee 
00059 *> \author Univ. of California Berkeley 
00060 *> \author Univ. of Colorado Denver 
00061 *> \author NAG Ltd. 
00062 *
00063 *> \date April 2012
00064 *
00065 *> \ingroup auxOTHERauxiliary
00066 *
00067 *  =====================================================================
00068       REAL FUNCTION SLAMCH( CMACH )
00069 *
00070 *  -- LAPACK auxiliary routine (version 3.4.1) --
00071 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00072 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00073 *     April 2012
00074 *
00075 *     .. Scalar Arguments ..
00076       CHARACTER          CMACH
00077 *     ..
00078 *     .. Parameters ..
00079       REAL               ONE, ZERO
00080       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00081 *     ..
00082 *     .. Local Scalars ..
00083       LOGICAL            FIRST, LRND
00084       INTEGER            BETA, IMAX, IMIN, IT
00085       REAL               BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
00086      $                   RND, SFMIN, SMALL, T
00087 *     ..
00088 *     .. External Functions ..
00089       LOGICAL            LSAME
00090       EXTERNAL           LSAME
00091 *     ..
00092 *     .. External Subroutines ..
00093       EXTERNAL           SLAMC2
00094 *     ..
00095 *     .. Save statement ..
00096       SAVE               FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN,
00097      $                   EMAX, RMAX, PREC
00098 *     ..
00099 *     .. Data statements ..
00100       DATA               FIRST / .TRUE. /
00101 *     ..
00102 *     .. Executable Statements ..
00103 *
00104       IF( FIRST ) THEN
00105          CALL SLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX )
00106          BASE = BETA
00107          T = IT
00108          IF( LRND ) THEN
00109             RND = ONE
00110             EPS = ( BASE**( 1-IT ) ) / 2
00111          ELSE
00112             RND = ZERO
00113             EPS = BASE**( 1-IT )
00114          END IF
00115          PREC = EPS*BASE
00116          EMIN = IMIN
00117          EMAX = IMAX
00118          SFMIN = RMIN
00119          SMALL = ONE / RMAX
00120          IF( SMALL.GE.SFMIN ) THEN
00121 *
00122 *           Use SMALL plus a bit, to avoid the possibility of rounding
00123 *           causing overflow when computing  1/sfmin.
00124 *
00125             SFMIN = SMALL*( ONE+EPS )
00126          END IF
00127       END IF
00128 *
00129       IF( LSAME( CMACH, 'E' ) ) THEN
00130          RMACH = EPS
00131       ELSE IF( LSAME( CMACH, 'S' ) ) THEN
00132          RMACH = SFMIN
00133       ELSE IF( LSAME( CMACH, 'B' ) ) THEN
00134          RMACH = BASE
00135       ELSE IF( LSAME( CMACH, 'P' ) ) THEN
00136          RMACH = PREC
00137       ELSE IF( LSAME( CMACH, 'N' ) ) THEN
00138          RMACH = T
00139       ELSE IF( LSAME( CMACH, 'R' ) ) THEN
00140          RMACH = RND
00141       ELSE IF( LSAME( CMACH, 'M' ) ) THEN
00142          RMACH = EMIN
00143       ELSE IF( LSAME( CMACH, 'U' ) ) THEN
00144          RMACH = RMIN
00145       ELSE IF( LSAME( CMACH, 'L' ) ) THEN
00146          RMACH = EMAX
00147       ELSE IF( LSAME( CMACH, 'O' ) ) THEN
00148          RMACH = RMAX
00149       END IF
00150 *
00151       SLAMCH = RMACH
00152       FIRST  = .FALSE.
00153       RETURN
00154 *
00155 *     End of SLAMCH
00156 *
00157       END
00158 *
00159 ************************************************************************
00160 *
00161 *> \brief \b SLAMC1
00162 *> \details
00163 *> \b Purpose:
00164 *> \verbatim
00165 *> SLAMC1 determines the machine parameters given by BETA, T, RND, and
00166 *> IEEE1.
00167 *> \endverbatim
00168 *>
00169 *> \param[out] BETA
00170 *> \verbatim
00171 *>          The base of the machine.
00172 *> \endverbatim
00173 *>
00174 *> \param[out] T
00175 *> \verbatim
00176 *>          The number of ( BETA ) digits in the mantissa.
00177 *> \endverbatim
00178 *>
00179 *> \param[out] RND
00180 *> \verbatim
00181 *>          Specifies whether proper rounding  ( RND = .TRUE. )  or
00182 *>          chopping  ( RND = .FALSE. )  occurs in addition. This may not
00183 *>          be a reliable guide to the way in which the machine performs
00184 *>          its arithmetic.
00185 *> \endverbatim
00186 *>
00187 *> \param[out] IEEE1
00188 *> \verbatim
00189 *>          Specifies whether rounding appears to be done in the IEEE
00190 *>          'round to nearest' style.
00191 *> \endverbatim
00192 *> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
00193 *> \date April 2012
00194 *> \ingroup auxOTHERauxiliary
00195 *>
00196 *> \details \b Further \b Details
00197 *> \verbatim
00198 *>
00199 *>  The routine is based on the routine  ENVRON  by Malcolm and
00200 *>  incorporates suggestions by Gentleman and Marovich. See
00201 *>
00202 *>     Malcolm M. A. (1972) Algorithms to reveal properties of
00203 *>        floating-point arithmetic. Comms. of the ACM, 15, 949-951.
00204 *>
00205 *>     Gentleman W. M. and Marovich S. B. (1974) More on algorithms
00206 *>        that reveal properties of floating point arithmetic units.
00207 *>        Comms. of the ACM, 17, 276-277.
00208 *> \endverbatim
00209 *>
00210       SUBROUTINE SLAMC1( BETA, T, RND, IEEE1 )
00211 *
00212 *  -- LAPACK auxiliary routine (version 3.4.1) --
00213 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00214 *     November 2010
00215 *
00216 *     .. Scalar Arguments ..
00217       LOGICAL            IEEE1, RND
00218       INTEGER            BETA, T
00219 *     ..
00220 * =====================================================================
00221 *
00222 *     .. Local Scalars ..
00223       LOGICAL            FIRST, LIEEE1, LRND
00224       INTEGER            LBETA, LT
00225       REAL               A, B, C, F, ONE, QTR, SAVEC, T1, T2
00226 *     ..
00227 *     .. External Functions ..
00228       REAL               SLAMC3
00229       EXTERNAL           SLAMC3
00230 *     ..
00231 *     .. Save statement ..
00232       SAVE               FIRST, LIEEE1, LBETA, LRND, LT
00233 *     ..
00234 *     .. Data statements ..
00235       DATA               FIRST / .TRUE. /
00236 *     ..
00237 *     .. Executable Statements ..
00238 *
00239       IF( FIRST ) THEN
00240          ONE = 1
00241 *
00242 *        LBETA,  LIEEE1,  LT and  LRND  are the  local values  of  BETA,
00243 *        IEEE1, T and RND.
00244 *
00245 *        Throughout this routine  we use the function  SLAMC3  to ensure
00246 *        that relevant values are  stored and not held in registers,  or
00247 *        are not affected by optimizers.
00248 *
00249 *        Compute  a = 2.0**m  with the  smallest positive integer m such
00250 *        that
00251 *
00252 *           fl( a + 1.0 ) = a.
00253 *
00254          A = 1
00255          C = 1
00256 *
00257 *+       WHILE( C.EQ.ONE )LOOP
00258    10    CONTINUE
00259          IF( C.EQ.ONE ) THEN
00260             A = 2*A
00261             C = SLAMC3( A, ONE )
00262             C = SLAMC3( C, -A )
00263             GO TO 10
00264          END IF
00265 *+       END WHILE
00266 *
00267 *        Now compute  b = 2.0**m  with the smallest positive integer m
00268 *        such that
00269 *
00270 *           fl( a + b ) .gt. a.
00271 *
00272          B = 1
00273          C = SLAMC3( A, B )
00274 *
00275 *+       WHILE( C.EQ.A )LOOP
00276    20    CONTINUE
00277          IF( C.EQ.A ) THEN
00278             B = 2*B
00279             C = SLAMC3( A, B )
00280             GO TO 20
00281          END IF
00282 *+       END WHILE
00283 *
00284 *        Now compute the base.  a and c  are neighbouring floating point
00285 *        numbers  in the  interval  ( beta**t, beta**( t + 1 ) )  and so
00286 *        their difference is beta. Adding 0.25 to c is to ensure that it
00287 *        is truncated to beta and not ( beta - 1 ).
00288 *
00289          QTR = ONE / 4
00290          SAVEC = C
00291          C = SLAMC3( C, -A )
00292          LBETA = C + QTR
00293 *
00294 *        Now determine whether rounding or chopping occurs,  by adding a
00295 *        bit  less  than  beta/2  and a  bit  more  than  beta/2  to  a.
00296 *
00297          B = LBETA
00298          F = SLAMC3( B / 2, -B / 100 )
00299          C = SLAMC3( F, A )
00300          IF( C.EQ.A ) THEN
00301             LRND = .TRUE.
00302          ELSE
00303             LRND = .FALSE.
00304          END IF
00305          F = SLAMC3( B / 2, B / 100 )
00306          C = SLAMC3( F, A )
00307          IF( ( LRND ) .AND. ( C.EQ.A ) )
00308      $      LRND = .FALSE.
00309 *
00310 *        Try and decide whether rounding is done in the  IEEE  'round to
00311 *        nearest' style. B/2 is half a unit in the last place of the two
00312 *        numbers A and SAVEC. Furthermore, A is even, i.e. has last  bit
00313 *        zero, and SAVEC is odd. Thus adding B/2 to A should not  change
00314 *        A, but adding B/2 to SAVEC should change SAVEC.
00315 *
00316          T1 = SLAMC3( B / 2, A )
00317          T2 = SLAMC3( B / 2, SAVEC )
00318          LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
00319 *
00320 *        Now find  the  mantissa, t.  It should  be the  integer part of
00321 *        log to the base beta of a,  however it is safer to determine  t
00322 *        by powering.  So we find t as the smallest positive integer for
00323 *        which
00324 *
00325 *           fl( beta**t + 1.0 ) = 1.0.
00326 *
00327          LT = 0
00328          A = 1
00329          C = 1
00330 *
00331 *+       WHILE( C.EQ.ONE )LOOP
00332    30    CONTINUE
00333          IF( C.EQ.ONE ) THEN
00334             LT = LT + 1
00335             A = A*LBETA
00336             C = SLAMC3( A, ONE )
00337             C = SLAMC3( C, -A )
00338             GO TO 30
00339          END IF
00340 *+       END WHILE
00341 *
00342       END IF
00343 *
00344       BETA = LBETA
00345       T = LT
00346       RND = LRND
00347       IEEE1 = LIEEE1
00348       FIRST = .FALSE.
00349       RETURN
00350 *
00351 *     End of SLAMC1
00352 *
00353       END
00354 *
00355 ************************************************************************
00356 *
00357 *> \brief \b SLAMC2
00358 *> \details
00359 *> \b Purpose:
00360 *> \verbatim
00361 *> SLAMC2 determines the machine parameters specified in its argument
00362 *> list.
00363 *> \endverbatim
00364 *> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
00365 *> \date April 2012
00366 *> \ingroup auxOTHERauxiliary
00367 *>
00368 *> \param[out] BETA
00369 *> \verbatim
00370 *>          The base of the machine.
00371 *> \endverbatim
00372 *>
00373 *> \param[out] T
00374 *> \verbatim
00375 *>          The number of ( BETA ) digits in the mantissa.
00376 *> \endverbatim
00377 *>
00378 *> \param[out] RND
00379 *> \verbatim
00380 *>          Specifies whether proper rounding  ( RND = .TRUE. )  or
00381 *>          chopping  ( RND = .FALSE. )  occurs in addition. This may not
00382 *>          be a reliable guide to the way in which the machine performs
00383 *>          its arithmetic.
00384 *> \endverbatim
00385 *>
00386 *> \param[out] EPS
00387 *> \verbatim
00388 *>          The smallest positive number such that
00389 *>             fl( 1.0 - EPS ) .LT. 1.0,
00390 *>          where fl denotes the computed value.
00391 *> \endverbatim
00392 *>
00393 *> \param[out] EMIN
00394 *> \verbatim
00395 *>          The minimum exponent before (gradual) underflow occurs.
00396 *> \endverbatim
00397 *>
00398 *> \param[out] RMIN
00399 *> \verbatim
00400 *>          The smallest normalized number for the machine, given by
00401 *>          BASE**( EMIN - 1 ), where  BASE  is the floating point value
00402 *>          of BETA.
00403 *> \endverbatim
00404 *>
00405 *> \param[out] EMAX
00406 *> \verbatim
00407 *>          The maximum exponent before overflow occurs.
00408 *> \endverbatim
00409 *>
00410 *> \param[out] RMAX
00411 *> \verbatim
00412 *>          The largest positive number for the machine, given by
00413 *>          BASE**EMAX * ( 1 - EPS ), where  BASE  is the floating point
00414 *>          value of BETA.
00415 *> \endverbatim
00416 *>
00417 *> \details \b Further \b Details
00418 *> \verbatim
00419 *>
00420 *>  The computation of  EPS  is based on a routine PARANOIA by
00421 *>  W. Kahan of the University of California at Berkeley.
00422 *> \endverbatim
00423       SUBROUTINE SLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
00424 *
00425 *  -- LAPACK auxiliary routine (version 3.4.1) --
00426 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00427 *     November 2010
00428 *
00429 *     .. Scalar Arguments ..
00430       LOGICAL            RND
00431       INTEGER            BETA, EMAX, EMIN, T
00432       REAL               EPS, RMAX, RMIN
00433 *     ..
00434 * =====================================================================
00435 *
00436 *     .. Local Scalars ..
00437       LOGICAL            FIRST, IEEE, IWARN, LIEEE1, LRND
00438       INTEGER            GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
00439      $                   NGNMIN, NGPMIN
00440       REAL               A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
00441      $                   SIXTH, SMALL, THIRD, TWO, ZERO
00442 *     ..
00443 *     .. External Functions ..
00444       REAL               SLAMC3
00445       EXTERNAL           SLAMC3
00446 *     ..
00447 *     .. External Subroutines ..
00448       EXTERNAL           SLAMC1, SLAMC4, SLAMC5
00449 *     ..
00450 *     .. Intrinsic Functions ..
00451       INTRINSIC          ABS, MAX, MIN
00452 *     ..
00453 *     .. Save statement ..
00454       SAVE               FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,
00455      $                   LRMIN, LT
00456 *     ..
00457 *     .. Data statements ..
00458       DATA               FIRST / .TRUE. / , IWARN / .FALSE. /
00459 *     ..
00460 *     .. Executable Statements ..
00461 *
00462       IF( FIRST ) THEN
00463          ZERO = 0
00464          ONE = 1
00465          TWO = 2
00466 *
00467 *        LBETA, LT, LRND, LEPS, LEMIN and LRMIN  are the local values of
00468 *        BETA, T, RND, EPS, EMIN and RMIN.
00469 *
00470 *        Throughout this routine  we use the function  SLAMC3  to ensure
00471 *        that relevant values are stored  and not held in registers,  or
00472 *        are not affected by optimizers.
00473 *
00474 *        SLAMC1 returns the parameters  LBETA, LT, LRND and LIEEE1.
00475 *
00476          CALL SLAMC1( LBETA, LT, LRND, LIEEE1 )
00477 *
00478 *        Start to find EPS.
00479 *
00480          B = LBETA
00481          A = B**( -LT )
00482          LEPS = A
00483 *
00484 *        Try some tricks to see whether or not this is the correct  EPS.
00485 *
00486          B = TWO / 3
00487          HALF = ONE / 2
00488          SIXTH = SLAMC3( B, -HALF )
00489          THIRD = SLAMC3( SIXTH, SIXTH )
00490          B = SLAMC3( THIRD, -HALF )
00491          B = SLAMC3( B, SIXTH )
00492          B = ABS( B )
00493          IF( B.LT.LEPS )
00494      $      B = LEPS
00495 *
00496          LEPS = 1
00497 *
00498 *+       WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
00499    10    CONTINUE
00500          IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN
00501             LEPS = B
00502             C = SLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
00503             C = SLAMC3( HALF, -C )
00504             B = SLAMC3( HALF, C )
00505             C = SLAMC3( HALF, -B )
00506             B = SLAMC3( HALF, C )
00507             GO TO 10
00508          END IF
00509 *+       END WHILE
00510 *
00511          IF( A.LT.LEPS )
00512      $      LEPS = A
00513 *
00514 *        Computation of EPS complete.
00515 *
00516 *        Now find  EMIN.  Let A = + or - 1, and + or - (1 + BASE**(-3)).
00517 *        Keep dividing  A by BETA until (gradual) underflow occurs. This
00518 *        is detected when we cannot recover the previous A.
00519 *
00520          RBASE = ONE / LBETA
00521          SMALL = ONE
00522          DO 20 I = 1, 3
00523             SMALL = SLAMC3( SMALL*RBASE, ZERO )
00524    20    CONTINUE
00525          A = SLAMC3( ONE, SMALL )
00526          CALL SLAMC4( NGPMIN, ONE, LBETA )
00527          CALL SLAMC4( NGNMIN, -ONE, LBETA )
00528          CALL SLAMC4( GPMIN, A, LBETA )
00529          CALL SLAMC4( GNMIN, -A, LBETA )
00530          IEEE = .FALSE.
00531 *
00532          IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN
00533             IF( NGPMIN.EQ.GPMIN ) THEN
00534                LEMIN = NGPMIN
00535 *            ( Non twos-complement machines, no gradual underflow;
00536 *              e.g.,  VAX )
00537             ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
00538                LEMIN = NGPMIN - 1 + LT
00539                IEEE = .TRUE.
00540 *            ( Non twos-complement machines, with gradual underflow;
00541 *              e.g., IEEE standard followers )
00542             ELSE
00543                LEMIN = MIN( NGPMIN, GPMIN )
00544 *            ( A guess; no known machine )
00545                IWARN = .TRUE.
00546             END IF
00547 *
00548          ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
00549             IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
00550                LEMIN = MAX( NGPMIN, NGNMIN )
00551 *            ( Twos-complement machines, no gradual underflow;
00552 *              e.g., CYBER 205 )
00553             ELSE
00554                LEMIN = MIN( NGPMIN, NGNMIN )
00555 *            ( A guess; no known machine )
00556                IWARN = .TRUE.
00557             END IF
00558 *
00559          ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND.
00560      $            ( GPMIN.EQ.GNMIN ) ) THEN
00561             IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN
00562                LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
00563 *            ( Twos-complement machines with gradual underflow;
00564 *              no known machine )
00565             ELSE
00566                LEMIN = MIN( NGPMIN, NGNMIN )
00567 *            ( A guess; no known machine )
00568                IWARN = .TRUE.
00569             END IF
00570 *
00571          ELSE
00572             LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
00573 *         ( A guess; no known machine )
00574             IWARN = .TRUE.
00575          END IF
00576          FIRST = .FALSE.
00577 ***
00578 * Comment out this if block if EMIN is ok
00579          IF( IWARN ) THEN
00580             FIRST = .TRUE.
00581             WRITE( 6, FMT = 9999 )LEMIN
00582          END IF
00583 ***
00584 *
00585 *        Assume IEEE arithmetic if we found denormalised  numbers above,
00586 *        or if arithmetic seems to round in the  IEEE style,  determined
00587 *        in routine SLAMC1. A true IEEE machine should have both  things
00588 *        true; however, faulty machines may have one or the other.
00589 *
00590          IEEE = IEEE .OR. LIEEE1
00591 *
00592 *        Compute  RMIN by successive division by  BETA. We could compute
00593 *        RMIN as BASE**( EMIN - 1 ),  but some machines underflow during
00594 *        this computation.
00595 *
00596          LRMIN = 1
00597          DO 30 I = 1, 1 - LEMIN
00598             LRMIN = SLAMC3( LRMIN*RBASE, ZERO )
00599    30    CONTINUE
00600 *
00601 *        Finally, call SLAMC5 to compute EMAX and RMAX.
00602 *
00603          CALL SLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )
00604       END IF
00605 *
00606       BETA = LBETA
00607       T = LT
00608       RND = LRND
00609       EPS = LEPS
00610       EMIN = LEMIN
00611       RMIN = LRMIN
00612       EMAX = LEMAX
00613       RMAX = LRMAX
00614 *
00615       RETURN
00616 *
00617  9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
00618      $      '  EMIN = ', I8, /
00619      $      ' If, after inspection, the value EMIN looks',
00620      $      ' acceptable please comment out ',
00621      $      / ' the IF block as marked within the code of routine',
00622      $      ' SLAMC2,', / ' otherwise supply EMIN explicitly.', / )
00623 *
00624 *     End of SLAMC2
00625 *
00626       END
00627 *
00628 ************************************************************************
00629 *
00630 *> \brief \b SLAMC3
00631 *> \details
00632 *> \b Purpose:
00633 *> \verbatim
00634 *> SLAMC3  is intended to force  A  and  B  to be stored prior to doing
00635 *> the addition of  A  and  B ,  for use in situations where optimizers
00636 *> might hold one of these in a register.
00637 *> \endverbatim
00638 *>
00639 *> \param[in] A
00640 *>
00641 *> \param[in] B
00642 *> \verbatim
00643 *>          The values A and B.
00644 *> \endverbatim
00645 
00646       REAL FUNCTION SLAMC3( A, B )
00647 *
00648 *  -- LAPACK auxiliary routine (version 3.4.1) --
00649 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00650 *     November 2010
00651 *
00652 *     .. Scalar Arguments ..
00653       REAL               A, B
00654 *     ..
00655 * =====================================================================
00656 *
00657 *     .. Executable Statements ..
00658 *
00659       SLAMC3 = A + B
00660 *
00661       RETURN
00662 *
00663 *     End of SLAMC3
00664 *
00665       END
00666 *
00667 ************************************************************************
00668 *
00669 *> \brief \b SLAMC4
00670 *> \details
00671 *> \b Purpose:
00672 *> \verbatim
00673 *> SLAMC4 is a service routine for SLAMC2.
00674 *> \endverbatim
00675 *>
00676 *> \param[out] EMIN
00677 *> \verbatim
00678 *>          The minimum exponent before (gradual) underflow, computed by
00679 *>          setting A = START and dividing by BASE until the previous A
00680 *>          can not be recovered.
00681 *> \endverbatim
00682 *>
00683 *> \param[in] START
00684 *> \verbatim
00685 *>          The starting point for determining EMIN.
00686 *> \endverbatim
00687 *>
00688 *> \param[in] BASE
00689 *> \verbatim
00690 *>          The base of the machine.
00691 *> \endverbatim
00692 *>
00693       SUBROUTINE SLAMC4( EMIN, START, BASE )
00694 *
00695 *  -- LAPACK auxiliary routine (version 3.4.1) --
00696 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00697 *     November 2010
00698 *
00699 *     .. Scalar Arguments ..
00700       INTEGER            BASE
00701       INTEGER            EMIN
00702       REAL               START
00703 *     ..
00704 * =====================================================================
00705 *
00706 *     .. Local Scalars ..
00707       INTEGER            I
00708       REAL               A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
00709 *     ..
00710 *     .. External Functions ..
00711       REAL               SLAMC3
00712       EXTERNAL           SLAMC3
00713 *     ..
00714 *     .. Executable Statements ..
00715 *
00716       A = START
00717       ONE = 1
00718       RBASE = ONE / BASE
00719       ZERO = 0
00720       EMIN = 1
00721       B1 = SLAMC3( A*RBASE, ZERO )
00722       C1 = A
00723       C2 = A
00724       D1 = A
00725       D2 = A
00726 *+    WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
00727 *    $       ( D1.EQ.A ).AND.( D2.EQ.A )      )LOOP
00728    10 CONTINUE
00729       IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND.
00730      $    ( D2.EQ.A ) ) THEN
00731          EMIN = EMIN - 1
00732          A = B1
00733          B1 = SLAMC3( A / BASE, ZERO )
00734          C1 = SLAMC3( B1*BASE, ZERO )
00735          D1 = ZERO
00736          DO 20 I = 1, BASE
00737             D1 = D1 + B1
00738    20    CONTINUE
00739          B2 = SLAMC3( A*RBASE, ZERO )
00740          C2 = SLAMC3( B2 / RBASE, ZERO )
00741          D2 = ZERO
00742          DO 30 I = 1, BASE
00743             D2 = D2 + B2
00744    30    CONTINUE
00745          GO TO 10
00746       END IF
00747 *+    END WHILE
00748 *
00749       RETURN
00750 *
00751 *     End of SLAMC4
00752 *
00753       END
00754 *
00755 ************************************************************************
00756 *
00757 *> \brief \b SLAMC5
00758 *> \details
00759 *> \b Purpose:
00760 *> \verbatim
00761 *> SLAMC5 attempts to compute RMAX, the largest machine floating-point
00762 *> number, without overflow.  It assumes that EMAX + abs(EMIN) sum
00763 *> approximately to a power of 2.  It will fail on machines where this
00764 *> assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
00765 *> EMAX = 28718).  It will also fail if the value supplied for EMIN is
00766 *> too large (i.e. too close to zero), probably with overflow.
00767 *> \endverbatim
00768 *>
00769 *> \param[in] BETA
00770 *> \verbatim
00771 *>          The base of floating-point arithmetic.
00772 *> \endverbatim
00773 *>
00774 *> \param[in] P
00775 *> \verbatim
00776 *>          The number of base BETA digits in the mantissa of a
00777 *>          floating-point value.
00778 *> \endverbatim
00779 *>
00780 *> \param[in] EMIN
00781 *> \verbatim
00782 *>          The minimum exponent before (gradual) underflow.
00783 *> \endverbatim
00784 *>
00785 *> \param[in] IEEE
00786 *> \verbatim
00787 *>          A logical flag specifying whether or not the arithmetic
00788 *>          system is thought to comply with the IEEE standard.
00789 *> \endverbatim
00790 *>
00791 *> \param[out] EMAX
00792 *> \verbatim
00793 *>          The largest exponent before overflow
00794 *> \endverbatim
00795 *>
00796 *> \param[out] RMAX
00797 *> \verbatim
00798 *>          The largest machine floating-point number.
00799 *> \endverbatim
00800 *>
00801       SUBROUTINE SLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
00802 *
00803 *  -- LAPACK auxiliary routine (version 3.4.1) --
00804 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00805 *     November 2010
00806 *
00807 *     .. Scalar Arguments ..
00808       LOGICAL            IEEE
00809       INTEGER            BETA, EMAX, EMIN, P
00810       REAL               RMAX
00811 *     ..
00812 * =====================================================================
00813 *
00814 *     .. Parameters ..
00815       REAL               ZERO, ONE
00816       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00817 *     ..
00818 *     .. Local Scalars ..
00819       INTEGER            EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
00820       REAL               OLDY, RECBAS, Y, Z
00821 *     ..
00822 *     .. External Functions ..
00823       REAL               SLAMC3
00824       EXTERNAL           SLAMC3
00825 *     ..
00826 *     .. Intrinsic Functions ..
00827       INTRINSIC          MOD
00828 *     ..
00829 *     .. Executable Statements ..
00830 *
00831 *     First compute LEXP and UEXP, two powers of 2 that bound
00832 *     abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
00833 *     approximately to the bound that is closest to abs(EMIN).
00834 *     (EMAX is the exponent of the required number RMAX).
00835 *
00836       LEXP = 1
00837       EXBITS = 1
00838    10 CONTINUE
00839       TRY = LEXP*2
00840       IF( TRY.LE.( -EMIN ) ) THEN
00841          LEXP = TRY
00842          EXBITS = EXBITS + 1
00843          GO TO 10
00844       END IF
00845       IF( LEXP.EQ.-EMIN ) THEN
00846          UEXP = LEXP
00847       ELSE
00848          UEXP = TRY
00849          EXBITS = EXBITS + 1
00850       END IF
00851 *
00852 *     Now -LEXP is less than or equal to EMIN, and -UEXP is greater
00853 *     than or equal to EMIN. EXBITS is the number of bits needed to
00854 *     store the exponent.
00855 *
00856       IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
00857          EXPSUM = 2*LEXP
00858       ELSE
00859          EXPSUM = 2*UEXP
00860       END IF
00861 *
00862 *     EXPSUM is the exponent range, approximately equal to
00863 *     EMAX - EMIN + 1 .
00864 *
00865       EMAX = EXPSUM + EMIN - 1
00866       NBITS = 1 + EXBITS + P
00867 *
00868 *     NBITS is the total number of bits needed to store a
00869 *     floating-point number.
00870 *
00871       IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
00872 *
00873 *        Either there are an odd number of bits used to store a
00874 *        floating-point number, which is unlikely, or some bits are
00875 *        not used in the representation of numbers, which is possible,
00876 *        (e.g. Cray machines) or the mantissa has an implicit bit,
00877 *        (e.g. IEEE machines, Dec Vax machines), which is perhaps the
00878 *        most likely. We have to assume the last alternative.
00879 *        If this is true, then we need to reduce EMAX by one because
00880 *        there must be some way of representing zero in an implicit-bit
00881 *        system. On machines like Cray, we are reducing EMAX by one
00882 *        unnecessarily.
00883 *
00884          EMAX = EMAX - 1
00885       END IF
00886 *
00887       IF( IEEE ) THEN
00888 *
00889 *        Assume we are on an IEEE machine which reserves one exponent
00890 *        for infinity and NaN.
00891 *
00892          EMAX = EMAX - 1
00893       END IF
00894 *
00895 *     Now create RMAX, the largest machine number, which should
00896 *     be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
00897 *
00898 *     First compute 1.0 - BETA**(-P), being careful that the
00899 *     result is less than 1.0 .
00900 *
00901       RECBAS = ONE / BETA
00902       Z = BETA - ONE
00903       Y = ZERO
00904       DO 20 I = 1, P
00905          Z = Z*RECBAS
00906          IF( Y.LT.ONE )
00907      $      OLDY = Y
00908          Y = SLAMC3( Y, Z )
00909    20 CONTINUE
00910       IF( Y.GE.ONE )
00911      $   Y = OLDY
00912 *
00913 *     Now multiply by BETA**EMAX to get RMAX.
00914 *
00915       DO 30 I = 1, EMAX
00916          Y = SLAMC3( Y*BETA, ZERO )
00917    30 CONTINUE
00918 *
00919       RMAX = Y
00920       RETURN
00921 *
00922 *     End of SLAMC5
00923 *
00924       END
 All Files Functions