OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
slamchf77.f File Reference

Go to the source code of this file.

Functions/Subroutines

real function slamch (cmach)
 SLAMCHF77 deprecated
subroutine slamc1 (beta, t, rnd, ieee1)
 SLAMC1
subroutine slamc2 (beta, t, rnd, eps, emin, rmin, emax, rmax)
 SLAMC2
real function slamc3 (a, b)
 SLAMC3
subroutine slamc4 (emin, start, base)
 SLAMC4
subroutine slamc5 (beta, p, emin, ieee, emax, rmax)
 SLAMC5

Function/Subroutine Documentation

◆ slamc1()

subroutine slamc1 ( integer beta,
integer t,
logical rnd,
logical ieee1 )

SLAMC1

Purpose:

!> SLAMC1 determines the machine parameters given by BETA, T, RND, and
!> IEEE1.
!> 
Parameters
[out]BETA
!>          The base of the machine.
!> 
[out]T
!>          The number of ( BETA ) digits in the mantissa.
!> 
[out]RND
!>          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.
!> 
[out]IEEE1
!>          Specifies whether rounding appears to be done in the IEEE
!>          'round to nearest' style.
!> 
Author
LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..

Further Details

!>
!>  The routine is based on the routine  ENVRON  by Malcolm and
!>  incorporates suggestions by Gentleman and Marovich. See
!>
!>     Malcolm M. A. (1972) Algorithms to reveal properties of
!>        floating-point arithmetic. Comms. of the ACM, 15, 949-951.
!>
!>     Gentleman W. M. and Marovich S. B. (1974) More on algorithms
!>        that reveal properties of floating point arithmetic units.
!>        Comms. of the ACM, 17, 276-277.
!> 

Definition at line 206 of file slamchf77.f.

207*
208* -- LAPACK auxiliary routine --
209* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
210*
211* .. Scalar Arguments ..
212 LOGICAL IEEE1, RND
213 INTEGER BETA, T
214* ..
215* =====================================================================
216*
217* .. Local Scalars ..
218 LOGICAL FIRST, LIEEE1, LRND
219 INTEGER LBETA, LT
220 REAL A, B, C, F, ONE, QTR, SAVEC, T1, T2
221* ..
222* .. External Functions ..
223 REAL SLAMC3
224 EXTERNAL slamc3
225* ..
226* .. Save statement ..
227 SAVE first, lieee1, lbeta, lrnd, lt
228* ..
229* .. Data statements ..
230 DATA first / .true. /
231* ..
232* .. Executable Statements ..
233*
234 IF( first ) THEN
235 one = 1
236*
237* LBETA, LIEEE1, LT and LRND are the local values of BETA,
238* IEEE1, T and RND.
239*
240* Throughout this routine we use the function SLAMC3 to ensure
241* that relevant values are stored and not held in registers, or
242* are not affected by optimizers.
243*
244* Compute a = 2.0**m with the smallest positive integer m such
245* that
246*
247* fl( a + 1.0 ) = a.
248*
249 a = 1
250 c = 1
251*
252*+ WHILE( C.EQ.ONE )LOOP
253 10 CONTINUE
254 IF( c.EQ.one ) THEN
255 a = 2*a
256 c = slamc3( a, one )
257 c = slamc3( c, -a )
258 GO TO 10
259 END IF
260*+ END WHILE
261*
262* Now compute b = 2.0**m with the smallest positive integer m
263* such that
264*
265* fl( a + b ) .gt. a.
266*
267 b = 1
268 c = slamc3( a, b )
269*
270*+ WHILE( C.EQ.A )LOOP
271 20 CONTINUE
272 IF( c.EQ.a ) THEN
273 b = 2*b
274 c = slamc3( a, b )
275 GO TO 20
276 END IF
277*+ END WHILE
278*
279* Now compute the base. a and c are neighbouring floating point
280* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
281* their difference is beta. Adding 0.25 to c is to ensure that it
282* is truncated to beta and not ( beta - 1 ).
283*
284 qtr = one / 4
285 savec = c
286 c = slamc3( c, -a )
287 lbeta = c + qtr
288*
289* Now determine whether rounding or chopping occurs, by adding a
290* bit less than beta/2 and a bit more than beta/2 to a.
291*
292 b = lbeta
293 f = slamc3( b / 2, -b / 100 )
294 c = slamc3( f, a )
295 IF( c.EQ.a ) THEN
296 lrnd = .true.
297 ELSE
298 lrnd = .false.
299 END IF
300 f = slamc3( b / 2, b / 100 )
301 c = slamc3( f, a )
302 IF( ( lrnd ) .AND. ( c.EQ.a ) )
303 $ lrnd = .false.
304*
305* Try and decide whether rounding is done in the IEEE 'round to
306* nearest' style. B/2 is half a unit in the last place of the two
307* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
308* zero, and SAVEC is odd. Thus adding B/2 to A should not change
309* A, but adding B/2 to SAVEC should change SAVEC.
310*
311 t1 = slamc3( b / 2, a )
312 t2 = slamc3( b / 2, savec )
313 lieee1 = ( t1.EQ.a ) .AND. ( t2.GT.savec ) .AND. lrnd
314*
315* Now find the mantissa, t. It should be the integer part of
316* log to the base beta of a, however it is safer to determine t
317* by powering. So we find t as the smallest positive integer for
318* which
319*
320* fl( beta**t + 1.0 ) = 1.0.
321*
322 lt = 0
323 a = 1
324 c = 1
325*
326*+ WHILE( C.EQ.ONE )LOOP
327 30 CONTINUE
328 IF( c.EQ.one ) THEN
329 lt = lt + 1
330 a = a*lbeta
331 c = slamc3( a, one )
332 c = slamc3( c, -a )
333 GO TO 30
334 END IF
335*+ END WHILE
336*
337 END IF
338*
339 beta = lbeta
340 t = lt
341 rnd = lrnd
342 ieee1 = lieee1
343 first = .false.
344 RETURN
345*
346* End of SLAMC1
347*
real function slamc3(a, b)
SLAMC3
Definition slamchf77.f:640

◆ slamc2()

subroutine slamc2 ( integer beta,
integer t,
logical rnd,
real eps,
integer emin,
real rmin,
integer emax,
real rmax )

SLAMC2

Purpose:

!> SLAMC2 determines the machine parameters specified in its argument
!> list.
!> 
Author
LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
Parameters
[out]BETA
!>          The base of the machine.
!> 
[out]T
!>          The number of ( BETA ) digits in the mantissa.
!> 
[out]RND
!>          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.
!> 
[out]EPS
!>          The smallest positive number such that
!>             fl( 1.0 - EPS ) .LT. 1.0,
!>          where fl denotes the computed value.
!> 
[out]EMIN
!>          The minimum exponent before (gradual) underflow occurs.
!> 
[out]RMIN
!>          The smallest normalized number for the machine, given by
!>          BASE**( EMIN - 1 ), where  BASE  is the floating point value
!>          of BETA.
!> 
[out]EMAX
!>          The maximum exponent before overflow occurs.
!> 
[out]RMAX
!>          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.
!> 

Definition at line 417 of file slamchf77.f.

418*
419* -- LAPACK auxiliary routine --
420* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
421*
422* .. Scalar Arguments ..
423 LOGICAL RND
424 INTEGER BETA, EMAX, EMIN, T
425 REAL EPS, RMAX, RMIN
426* ..
427* =====================================================================
428*
429* .. Local Scalars ..
430 LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND
431 INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
432 $ NGNMIN, NGPMIN
433 REAL A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
434 $ SIXTH, SMALL, THIRD, TWO, ZERO
435* ..
436* .. External Functions ..
437 REAL SLAMC3
438 EXTERNAL slamc3
439* ..
440* .. External Subroutines ..
441 EXTERNAL slamc1, slamc4, slamc5
442* ..
443* .. Intrinsic Functions ..
444 INTRINSIC abs, max, min
445* ..
446* .. Save statement ..
447 SAVE first, iwarn, lbeta, lemax, lemin, leps, lrmax,
448 $ lrmin, lt
449* ..
450* .. Data statements ..
451 DATA first / .true. / , iwarn / .false. /
452* ..
453* .. Executable Statements ..
454*
455 IF( first ) THEN
456 zero = 0
457 one = 1
458 two = 2
459*
460* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
461* BETA, T, RND, EPS, EMIN and RMIN.
462*
463* Throughout this routine we use the function SLAMC3 to ensure
464* that relevant values are stored and not held in registers, or
465* are not affected by optimizers.
466*
467* SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
468*
469 CALL slamc1( lbeta, lt, lrnd, lieee1 )
470*
471* Start to find EPS.
472*
473 b = lbeta
474 a = b**( -lt )
475 leps = a
476*
477* Try some tricks to see whether or not this is the correct EPS.
478*
479 b = two / 3
480 half = one / 2
481 sixth = slamc3( b, -half )
482 third = slamc3( sixth, sixth )
483 b = slamc3( third, -half )
484 b = slamc3( b, sixth )
485 b = abs( b )
486 IF( b.LT.leps )
487 $ b = leps
488*
489 leps = 1
490*
491*+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
492 10 CONTINUE
493 IF( ( leps.GT.b ) .AND. ( b.GT.zero ) ) THEN
494 leps = b
495 c = slamc3( half*leps, ( two**5 )*( leps**2 ) )
496 c = slamc3( half, -c )
497 b = slamc3( half, c )
498 c = slamc3( half, -b )
499 b = slamc3( half, c )
500 GO TO 10
501 END IF
502*+ END WHILE
503*
504 IF( a.LT.leps )
505 $ leps = a
506*
507* Computation of EPS complete.
508*
509* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
510* Keep dividing A by BETA until (gradual) underflow occurs. This
511* is detected when we cannot recover the previous A.
512*
513 rbase = one / lbeta
514 small = one
515 DO 20 i = 1, 3
516 small = slamc3( small*rbase, zero )
517 20 CONTINUE
518 a = slamc3( one, small )
519 CALL slamc4( ngpmin, one, lbeta )
520 CALL slamc4( ngnmin, -one, lbeta )
521 CALL slamc4( gpmin, a, lbeta )
522 CALL slamc4( gnmin, -a, lbeta )
523 ieee = .false.
524*
525 IF( ( ngpmin.EQ.ngnmin ) .AND. ( gpmin.EQ.gnmin ) ) THEN
526 IF( ngpmin.EQ.gpmin ) THEN
527 lemin = ngpmin
528* ( Non twos-complement machines, no gradual underflow;
529* e.g., VAX )
530 ELSE IF( ( gpmin-ngpmin ).EQ.3 ) THEN
531 lemin = ngpmin - 1 + lt
532 ieee = .true.
533* ( Non twos-complement machines, with gradual underflow;
534* e.g., IEEE standard followers )
535 ELSE
536 lemin = min( ngpmin, gpmin )
537* ( A guess; no known machine )
538 iwarn = .true.
539 END IF
540*
541 ELSE IF( ( ngpmin.EQ.gpmin ) .AND. ( ngnmin.EQ.gnmin ) ) THEN
542 IF( abs( ngpmin-ngnmin ).EQ.1 ) THEN
543 lemin = max( ngpmin, ngnmin )
544* ( Twos-complement machines, no gradual underflow;
545* e.g., CYBER 205 )
546 ELSE
547 lemin = min( ngpmin, ngnmin )
548* ( A guess; no known machine )
549 iwarn = .true.
550 END IF
551*
552 ELSE IF( ( abs( ngpmin-ngnmin ).EQ.1 ) .AND.
553 $ ( gpmin.EQ.gnmin ) ) THEN
554 IF( ( gpmin-min( ngpmin, ngnmin ) ).EQ.3 ) THEN
555 lemin = max( ngpmin, ngnmin ) - 1 + lt
556* ( Twos-complement machines with gradual underflow;
557* no known machine )
558 ELSE
559 lemin = min( ngpmin, ngnmin )
560* ( A guess; no known machine )
561 iwarn = .true.
562 END IF
563*
564 ELSE
565 lemin = min( ngpmin, ngnmin, gpmin, gnmin )
566* ( A guess; no known machine )
567 iwarn = .true.
568 END IF
569 first = .false.
570***
571* Comment out this if block if EMIN is ok
572 IF( iwarn ) THEN
573 first = .true.
574 WRITE( 6, fmt = 9999 )lemin
575 END IF
576***
577*
578* Assume IEEE arithmetic if we found denormalised numbers above,
579* or if arithmetic seems to round in the IEEE style, determined
580* in routine SLAMC1. A true IEEE machine should have both things
581* true; however, faulty machines may have one or the other.
582*
583 ieee = ieee .OR. lieee1
584*
585* Compute RMIN by successive division by BETA. We could compute
586* RMIN as BASE**( EMIN - 1 ), but some machines underflow during
587* this computation.
588*
589 lrmin = 1
590 DO 30 i = 1, 1 - lemin
591 lrmin = slamc3( lrmin*rbase, zero )
592 30 CONTINUE
593*
594* Finally, call SLAMC5 to compute EMAX and RMAX.
595*
596 CALL slamc5( lbeta, lt, lemin, ieee, lemax, lrmax )
597 END IF
598*
599 beta = lbeta
600 t = lt
601 rnd = lrnd
602 eps = leps
603 emin = lemin
604 rmin = lrmin
605 emax = lemax
606 rmax = lrmax
607*
608 RETURN
609*
610 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
611 $ ' EMIN = ', i8, /
612 $ ' If, after inspection, the value EMIN looks',
613 $ ' acceptable please comment out ',
614 $ / ' the IF block as marked within the code of routine',
615 $ ' SLAMC2,', / ' otherwise supply EMIN explicitly.', / )
616*
617* End of SLAMC2
618*
#define min(a, b)
Definition macros.h:20
#define max(a, b)
Definition macros.h:21
subroutine slamc5(beta, p, emin, ieee, emax, rmax)
SLAMC5
Definition slamchf77.f:793
subroutine slamc1(beta, t, rnd, ieee1)
SLAMC1
Definition slamchf77.f:207
subroutine slamc4(emin, start, base)
SLAMC4
Definition slamchf77.f:686

◆ slamc3()

real function slamc3 ( real a,
real b )

SLAMC3

Purpose:

!> SLAMC3  is intended to force  A  and  B  to be stored prior to doing
!> the addition of  A  and  B ,  for use in situations where optimizers
!> might hold one of these in a register.
!> 
Parameters
[in]A
[in]B
!>          The values A and B.
!> 

Definition at line 639 of file slamchf77.f.

640*
641* -- LAPACK auxiliary routine --
642* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
643*
644* .. Scalar Arguments ..
645 REAL A, B
646* ..
647* =====================================================================
648*
649* .. Executable Statements ..
650*
651 slamc3 = a + b
652*
653 RETURN
654*
655* End of SLAMC3
656*

◆ slamc4()

subroutine slamc4 ( integer emin,
real start,
integer base )

SLAMC4

Purpose:

!> SLAMC4 is a service routine for SLAMC2.
!> 
Parameters
[out]EMIN
!>          The minimum exponent before (gradual) underflow, computed by
!>          setting A = START and dividing by BASE until the previous A
!>          can not be recovered.
!> 
[in]START
!>          The starting point for determining EMIN.
!> 
[in]BASE
!>          The base of the machine.
!> 

Definition at line 685 of file slamchf77.f.

686*
687* -- LAPACK auxiliary routine --
688* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
689*
690* .. Scalar Arguments ..
691 INTEGER BASE
692 INTEGER EMIN
693 REAL START
694* ..
695* =====================================================================
696*
697* .. Local Scalars ..
698 INTEGER I
699 REAL A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
700* ..
701* .. External Functions ..
702 REAL SLAMC3
703 EXTERNAL slamc3
704* ..
705* .. Executable Statements ..
706*
707 a = start
708 one = 1
709 rbase = one / base
710 zero = 0
711 emin = 1
712 b1 = slamc3( a*rbase, zero )
713 c1 = a
714 c2 = a
715 d1 = a
716 d2 = a
717*+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
718* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP
719 10 CONTINUE
720 IF( ( c1.EQ.a ) .AND. ( c2.EQ.a ) .AND. ( d1.EQ.a ) .AND.
721 $ ( d2.EQ.a ) ) THEN
722 emin = emin - 1
723 a = b1
724 b1 = slamc3( a / base, zero )
725 c1 = slamc3( b1*base, zero )
726 d1 = zero
727 DO 20 i = 1, base
728 d1 = d1 + b1
729 20 CONTINUE
730 b2 = slamc3( a*rbase, zero )
731 c2 = slamc3( b2 / rbase, zero )
732 d2 = zero
733 DO 30 i = 1, base
734 d2 = d2 + b2
735 30 CONTINUE
736 GO TO 10
737 END IF
738*+ END WHILE
739*
740 RETURN
741*
742* End of SLAMC4
743*

◆ slamc5()

subroutine slamc5 ( integer beta,
integer p,
integer emin,
logical ieee,
integer emax,
real rmax )

SLAMC5

Purpose:

!> SLAMC5 attempts to compute RMAX, the largest machine floating-point
!> number, without overflow.  It assumes that EMAX + abs(EMIN) sum
!> approximately to a power of 2.  It will fail on machines where this
!> assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
!> EMAX = 28718).  It will also fail if the value supplied for EMIN is
!> too large (i.e. too close to zero), probably with overflow.
!> 
Parameters
[in]BETA
!>          The base of floating-point arithmetic.
!> 
[in]P
!>          The number of base BETA digits in the mantissa of a
!>          floating-point value.
!> 
[in]EMIN
!>          The minimum exponent before (gradual) underflow.
!> 
[in]IEEE
!>          A logical flag specifying whether or not the arithmetic
!>          system is thought to comply with the IEEE standard.
!> 
[out]EMAX
!>          The largest exponent before overflow
!> 
[out]RMAX
!>          The largest machine floating-point number.
!> 

Definition at line 792 of file slamchf77.f.

793*
794* -- LAPACK auxiliary routine --
795* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
796*
797* .. Scalar Arguments ..
798 LOGICAL IEEE
799 INTEGER BETA, EMAX, EMIN, P
800 REAL RMAX
801* ..
802* =====================================================================
803*
804* .. Parameters ..
805 REAL ZERO, ONE
806 parameter( zero = 0.0e0, one = 1.0e0 )
807* ..
808* .. Local Scalars ..
809 INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
810 REAL OLDY, RECBAS, Y, Z
811* ..
812* .. External Functions ..
813 REAL SLAMC3
814 EXTERNAL slamc3
815* ..
816* .. Intrinsic Functions ..
817 INTRINSIC mod
818* ..
819* .. Executable Statements ..
820*
821* First compute LEXP and UEXP, two powers of 2 that bound
822* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
823* approximately to the bound that is closest to abs(EMIN).
824* (EMAX is the exponent of the required number RMAX).
825*
826 lexp = 1
827 exbits = 1
828 10 CONTINUE
829 try = lexp*2
830 IF( try.LE.( -emin ) ) THEN
831 lexp = try
832 exbits = exbits + 1
833 GO TO 10
834 END IF
835 IF( lexp.EQ.-emin ) THEN
836 uexp = lexp
837 ELSE
838 uexp = try
839 exbits = exbits + 1
840 END IF
841*
842* Now -LEXP is less than or equal to EMIN, and -UEXP is greater
843* than or equal to EMIN. EXBITS is the number of bits needed to
844* store the exponent.
845*
846 IF( ( uexp+emin ).GT.( -lexp-emin ) ) THEN
847 expsum = 2*lexp
848 ELSE
849 expsum = 2*uexp
850 END IF
851*
852* EXPSUM is the exponent range, approximately equal to
853* EMAX - EMIN + 1 .
854*
855 emax = expsum + emin - 1
856 nbits = 1 + exbits + p
857*
858* NBITS is the total number of bits needed to store a
859* floating-point number.
860*
861 IF( ( mod( nbits, 2 ).EQ.1 ) .AND. ( beta.EQ.2 ) ) THEN
862*
863* Either there are an odd number of bits used to store a
864* floating-point number, which is unlikely, or some bits are
865* not used in the representation of numbers, which is possible,
866* (e.g. Cray machines) or the mantissa has an implicit bit,
867* (e.g. IEEE machines, Dec Vax machines), which is perhaps the
868* most likely. We have to assume the last alternative.
869* If this is true, then we need to reduce EMAX by one because
870* there must be some way of representing zero in an implicit-bit
871* system. On machines like Cray, we are reducing EMAX by one
872* unnecessarily.
873*
874 emax = emax - 1
875 END IF
876*
877 IF( ieee ) THEN
878*
879* Assume we are on an IEEE machine which reserves one exponent
880* for infinity and NaN.
881*
882 emax = emax - 1
883 END IF
884*
885* Now create RMAX, the largest machine number, which should
886* be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
887*
888* First compute 1.0 - BETA**(-P), being careful that the
889* result is less than 1.0 .
890*
891 recbas = one / beta
892 z = beta - one
893 y = zero
894 DO 20 i = 1, p
895 z = z*recbas
896 IF( y.LT.one )
897 $ oldy = y
898 y = slamc3( y, z )
899 20 CONTINUE
900 IF( y.GE.one )
901 $ y = oldy
902*
903* Now multiply by BETA**EMAX to get RMAX.
904*
905 DO 30 i = 1, emax
906 y = slamc3( y*beta, zero )
907 30 CONTINUE
908*
909 rmax = y
910 RETURN
911*
912* End of SLAMC5
913*

◆ slamch()

real function slamch ( character cmach)

SLAMCHF77 deprecated

Purpose:
!>
!> SLAMCH determines single precision machine parameters.
!> 
Parameters
[in]CMACH
!>          Specifies the value to be returned by SLAMCH:
!>          = 'E' or 'e',   SLAMCH := eps
!>          = 'S' or 's ,   SLAMCH := sfmin
!>          = 'B' or 'b',   SLAMCH := base
!>          = 'P' or 'p',   SLAMCH := eps*base
!>          = 'N' or 'n',   SLAMCH := t
!>          = 'R' or 'r',   SLAMCH := rnd
!>          = 'M' or 'm',   SLAMCH := emin
!>          = 'U' or 'u',   SLAMCH := rmin
!>          = 'L' or 'l',   SLAMCH := emax
!>          = 'O' or 'o',   SLAMCH := rmax
!>          where
!>          eps   = relative machine precision
!>          sfmin = safe minimum, such that 1/sfmin does not overflow
!>          base  = base of the machine
!>          prec  = eps*base
!>          t     = number of (base) digits in the mantissa
!>          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise
!>          emin  = minimum exponent before (gradual) underflow
!>          rmin  = underflow threshold - base**(emin-1)
!>          emax  = largest exponent before overflow
!>          rmax  = overflow threshold  - (base**emax)*(1-eps)
!> 
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 66 of file slamchf77.f.

67*
68* -- LAPACK auxiliary routine --
69* -- LAPACK is a software package provided by Univ. of Tennessee, --
70* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
71*
72* .. Scalar Arguments ..
73 CHARACTER CMACH
74* ..
75* .. Parameters ..
76 REAL ONE, ZERO
77 parameter( one = 1.0e+0, zero = 0.0e+0 )
78* ..
79* .. Local Scalars ..
80 LOGICAL FIRST, LRND
81 INTEGER BETA, IMAX, IMIN, IT
82 REAL BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
83 $ RND, SFMIN, SMALL, T
84* ..
85* .. External Functions ..
86 LOGICAL LSAME
87 EXTERNAL lsame
88* ..
89* .. External Subroutines ..
90 EXTERNAL slamc2
91* ..
92* .. Save statement ..
93 SAVE first, eps, sfmin, base, t, rnd, emin, rmin,
94 $ emax, rmax, prec
95* ..
96* .. Data statements ..
97 DATA first / .true. /
98* ..
99* .. Executable Statements ..
100*
101 IF( first ) THEN
102 CALL slamc2( beta, it, lrnd, eps, imin, rmin, imax, rmax )
103 base = beta
104 t = it
105 IF( lrnd ) THEN
106 rnd = one
107 eps = ( base**( 1-it ) ) / 2
108 ELSE
109 rnd = zero
110 eps = base**( 1-it )
111 END IF
112 prec = eps*base
113 emin = imin
114 emax = imax
115 sfmin = rmin
116 small = one / rmax
117 IF( small.GE.sfmin ) THEN
118*
119* Use SMALL plus a bit, to avoid the possibility of rounding
120* causing overflow when computing 1/sfmin.
121*
122 sfmin = small*( one+eps )
123 END IF
124 END IF
125*
126 IF( lsame( cmach, 'E' ) ) THEN
127 rmach = eps
128 ELSE IF( lsame( cmach, 'S' ) ) THEN
129 rmach = sfmin
130 ELSE IF( lsame( cmach, 'B' ) ) THEN
131 rmach = base
132 ELSE IF( lsame( cmach, 'P' ) ) THEN
133 rmach = prec
134 ELSE IF( lsame( cmach, 'N' ) ) THEN
135 rmach = t
136 ELSE IF( lsame( cmach, 'R' ) ) THEN
137 rmach = rnd
138 ELSE IF( lsame( cmach, 'M' ) ) THEN
139 rmach = emin
140 ELSE IF( lsame( cmach, 'U' ) ) THEN
141 rmach = rmin
142 ELSE IF( lsame( cmach, 'L' ) ) THEN
143 rmach = emax
144 ELSE IF( lsame( cmach, 'O' ) ) THEN
145 rmach = rmax
146 END IF
147*
148 slamch = rmach
149 first = .false.
150 RETURN
151*
152* End of SLAMCH
153*
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
real function slamch(cmach)
SLAMCHF77 deprecated
Definition slamchf77.f:67
subroutine slamc2(beta, t, rnd, eps, emin, rmin, emax, rmax)
SLAMC2
Definition slamchf77.f:418