Mercurial > octave
diff libcruft/lapack/slamc2.f @ 7789:82be108cc558
First attempt at single precision tyeps
* * *
corrections to qrupdate single precision routines
* * *
prefer demotion to single over promotion to double
* * *
Add single precision support to log2 function
* * *
Trivial PROJECT file update
* * *
Cache optimized hermitian/transpose methods
* * *
Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
| author | David Bateman <dbateman@free.fr> |
|---|---|
| date | Sun, 27 Apr 2008 22:34:17 +0200 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/slamc2.f Sun Apr 27 22:34:17 2008 +0200 @@ -0,0 +1,255 @@ + SUBROUTINE SLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX ) +* +* -- LAPACK auxiliary routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + LOGICAL RND + INTEGER BETA, EMAX, EMIN, T + REAL EPS, RMAX, RMIN +* .. +* +* Purpose +* ======= +* +* SLAMC2 determines the machine parameters specified in its argument +* list. +* +* Arguments +* ========= +* +* BETA (output) INTEGER +* The base of the machine. +* +* T (output) INTEGER +* The number of ( BETA ) digits in the mantissa. +* +* RND (output) LOGICAL +* Specifies whether proper rounding ( RND = .TRUE. ) or +* chopping ( RND = .FALSE. ) occurs in addition. This may not +* be a reliable guide to the way in which the machine performs +* its arithmetic. +* +* EPS (output) REAL +* The smallest positive number such that +* +* fl( 1.0 - EPS ) .LT. 1.0, +* +* where fl denotes the computed value. +* +* EMIN (output) INTEGER +* The minimum exponent before (gradual) underflow occurs. +* +* RMIN (output) REAL +* The smallest normalized number for the machine, given by +* BASE**( EMIN - 1 ), where BASE is the floating point value +* of BETA. +* +* EMAX (output) INTEGER +* The maximum exponent before overflow occurs. +* +* RMAX (output) REAL +* The largest positive number for the machine, given by +* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point +* value of BETA. +* +* Further Details +* =============== +* +* The computation of EPS is based on a routine PARANOIA by +* W. Kahan of the University of California at Berkeley. +* +* ===================================================================== +* +* .. Local Scalars .. + LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND + INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT, + $ NGNMIN, NGPMIN + REAL A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE, + $ SIXTH, SMALL, THIRD, TWO, ZERO +* .. +* .. External Functions .. + REAL SLAMC3 + EXTERNAL SLAMC3 +* .. +* .. External Subroutines .. + EXTERNAL SLAMC1, SLAMC4, SLAMC5 +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MIN +* .. +* .. Save statement .. + SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX, + $ LRMIN, LT +* .. +* .. Data statements .. + DATA FIRST / .TRUE. / , IWARN / .FALSE. / +* .. +* .. Executable Statements .. +* + IF( FIRST ) THEN + ZERO = 0 + ONE = 1 + TWO = 2 +* +* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of +* BETA, T, RND, EPS, EMIN and RMIN. +* +* Throughout this routine we use the function SLAMC3 to ensure +* that relevant values are stored and not held in registers, or +* are not affected by optimizers. +* +* SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. +* + CALL SLAMC1( LBETA, LT, LRND, LIEEE1 ) +* +* Start to find EPS. +* + B = LBETA + A = B**( -LT ) + LEPS = A +* +* Try some tricks to see whether or not this is the correct EPS. +* + B = TWO / 3 + HALF = ONE / 2 + SIXTH = SLAMC3( B, -HALF ) + THIRD = SLAMC3( SIXTH, SIXTH ) + B = SLAMC3( THIRD, -HALF ) + B = SLAMC3( B, SIXTH ) + B = ABS( B ) + IF( B.LT.LEPS ) + $ B = LEPS +* + LEPS = 1 +* +*+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP + 10 CONTINUE + IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN + LEPS = B + C = SLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) ) + C = SLAMC3( HALF, -C ) + B = SLAMC3( HALF, C ) + C = SLAMC3( HALF, -B ) + B = SLAMC3( HALF, C ) + GO TO 10 + END IF +*+ END WHILE +* + IF( A.LT.LEPS ) + $ LEPS = A +* +* Computation of EPS complete. +* +* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). +* Keep dividing A by BETA until (gradual) underflow occurs. This +* is detected when we cannot recover the previous A. +* + RBASE = ONE / LBETA + SMALL = ONE + DO 20 I = 1, 3 + SMALL = SLAMC3( SMALL*RBASE, ZERO ) + 20 CONTINUE + A = SLAMC3( ONE, SMALL ) + CALL SLAMC4( NGPMIN, ONE, LBETA ) + CALL SLAMC4( NGNMIN, -ONE, LBETA ) + CALL SLAMC4( GPMIN, A, LBETA ) + CALL SLAMC4( GNMIN, -A, LBETA ) + IEEE = .FALSE. +* + IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN + IF( NGPMIN.EQ.GPMIN ) THEN + LEMIN = NGPMIN +* ( Non twos-complement machines, no gradual underflow; +* e.g., VAX ) + ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN + LEMIN = NGPMIN - 1 + LT + IEEE = .TRUE. +* ( Non twos-complement machines, with gradual underflow; +* e.g., IEEE standard followers ) + ELSE + LEMIN = MIN( NGPMIN, GPMIN ) +* ( A guess; no known machine ) + IWARN = .TRUE. + END IF +* + ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN + IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN + LEMIN = MAX( NGPMIN, NGNMIN ) +* ( Twos-complement machines, no gradual underflow; +* e.g., CYBER 205 ) + ELSE + LEMIN = MIN( NGPMIN, NGNMIN ) +* ( A guess; no known machine ) + IWARN = .TRUE. + END IF +* + ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND. + $ ( GPMIN.EQ.GNMIN ) ) THEN + IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN + LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT +* ( Twos-complement machines with gradual underflow; +* no known machine ) + ELSE + LEMIN = MIN( NGPMIN, NGNMIN ) +* ( A guess; no known machine ) + IWARN = .TRUE. + END IF +* + ELSE + LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN ) +* ( A guess; no known machine ) + IWARN = .TRUE. + END IF + FIRST = .FALSE. +*** +* Comment out this if block if EMIN is ok + IF( IWARN ) THEN + FIRST = .TRUE. + WRITE( 6, FMT = 9999 )LEMIN + END IF +*** +* +* Assume IEEE arithmetic if we found denormalised numbers above, +* or if arithmetic seems to round in the IEEE style, determined +* in routine SLAMC1. A true IEEE machine should have both things +* true; however, faulty machines may have one or the other. +* + IEEE = IEEE .OR. LIEEE1 +* +* Compute RMIN by successive division by BETA. We could compute +* RMIN as BASE**( EMIN - 1 ), but some machines underflow during +* this computation. +* + LRMIN = 1 + DO 30 I = 1, 1 - LEMIN + LRMIN = SLAMC3( LRMIN*RBASE, ZERO ) + 30 CONTINUE +* +* Finally, call SLAMC5 to compute EMAX and RMAX. +* + CALL SLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX ) + END IF +* + BETA = LBETA + T = LT + RND = LRND + EPS = LEPS + EMIN = LEMIN + RMIN = LRMIN + EMAX = LEMAX + RMAX = LRMAX +* + RETURN +* + 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-', + $ ' EMIN = ', I8, / + $ ' If, after inspection, the value EMIN looks', + $ ' acceptable please comment out ', + $ / ' the IF block as marked within the code of routine', + $ ' SLAMC2,', / ' otherwise supply EMIN explicitly.', / ) +* +* End of SLAMC2 +* + END