LAPACK 3.3.1 Linear Algebra PACKage

# slamchf77.f

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