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

Go to the source code of this file.

Functions/Subroutines

double precision function dlamch (cmach)
 DLAMCHF77 deprecated
subroutine dlamc1 (beta, t, rnd, ieee1)
 DLAMC1
subroutine dlamc2 (beta, t, rnd, eps, emin, rmin, emax, rmax)
 DLAMC2
double precision function dlamc3 (a, b)
 DLAMC3
subroutine dlamc4 (emin, start, base)
 DLAMC4
subroutine dlamc5 (beta, p, emin, ieee, emax, rmax)
 DLAMC5

Function/Subroutine Documentation

◆ dlamc1()

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

DLAMC1

Purpose: 

!> DLAMC1 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 207 of file dlamchf77.f.

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

◆ dlamc2()

subroutine dlamc2 ( integer beta,
integer t,
logical rnd,
double precision eps,
integer emin,
double precision rmin,
integer emax,
double precision rmax )

DLAMC2

Purpose: 

!> DLAMC2 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 418 of file dlamchf77.f.

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

◆ dlamc3()

double precision function dlamc3 ( double precision a,
double precision b )

DLAMC3

Purpose: 

!> DLAMC3  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 640 of file dlamchf77.f.

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

◆ dlamc4()

subroutine dlamc4 ( integer emin,
double precision start,
integer base )

DLAMC4

Purpose: 

!> DLAMC4 is a service routine for DLAMC2.
!>
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 686 of file dlamchf77.f.

687*
688* -- LAPACK auxiliary routine --
689* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
690*
691* .. Scalar Arguments ..
692 INTEGER BASE, EMIN
693 DOUBLE PRECISION START
694* ..
695* =====================================================================
696*
697* .. Local Scalars ..
698 INTEGER I
699 DOUBLE PRECISION A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
700* ..
701* .. External Functions ..
702 DOUBLE PRECISION DLAMC3
703 EXTERNAL dlamc3
704* ..
705* .. Executable Statements ..
706*
707 a = start
708 one = 1
709 rbase = one / base
710 zero = 0
711 emin = 1
712 b1 = dlamc3( 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 = dlamc3( a / base, zero )
725 c1 = dlamc3( b1*base, zero )
726 d1 = zero
727 DO 20 i = 1, base
728 d1 = d1 + b1
729 20 CONTINUE
730 b2 = dlamc3( a*rbase, zero )
731 c2 = dlamc3( 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 DLAMC4
743*

◆ dlamc5()

subroutine dlamc5 ( integer beta,
integer p,
integer emin,
logical ieee,
integer emax,
double precision rmax )

DLAMC5

Purpose: 

!> DLAMC5 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 dlamchf77.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 DOUBLE PRECISION RMAX
801* ..
802* =====================================================================
803*
804* .. Parameters ..
805 DOUBLE PRECISION ZERO, ONE
806 parameter( zero = 0.0d0, one = 1.0d0 )
807* ..
808* .. Local Scalars ..
809 INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
810 DOUBLE PRECISION OLDY, RECBAS, Y, Z
811* ..
812* .. External Functions ..
813 DOUBLE PRECISION DLAMC3
814 EXTERNAL dlamc3
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 = dlamc3( 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 = dlamc3( y*beta, zero )
907 30 CONTINUE
908*
909 rmax = y
910 RETURN
911*
912* End of DLAMC5
913*

◆ dlamch()

double precision function dlamch ( character cmach)

DLAMCHF77 deprecated

Purpose:
!>
!> DLAMCHF77 determines double precision machine parameters.
!>
Parameters
[in]CMACH
!>          Specifies the value to be returned by DLAMCH:
!>          = 'E' or 'e',   DLAMCH := eps
!>          = 'S' or 's ,   DLAMCH := sfmin
!>          = 'B' or 'b',   DLAMCH := base
!>          = 'P' or 'p',   DLAMCH := eps*base
!>          = 'N' or 'n',   DLAMCH := t
!>          = 'R' or 'r',   DLAMCH := rnd
!>          = 'M' or 'm',   DLAMCH := emin
!>          = 'U' or 'u',   DLAMCH := rmin
!>          = 'L' or 'l',   DLAMCH := emax
!>          = 'O' or 'o',   DLAMCH := 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 67 of file dlamchf77.f.

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