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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

◆ dlamc1()

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

Definition at line 131 of file dlamch.f.

132*
133* -- LAPACK auxiliary routine (version 2.1) --
134* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
135* Courant Institute, Argonne National Lab, and Rice University
136* October 31, 1992
137*
138* .. Scalar Arguments ..
139 LOGICAL IEEE1, RND
140 INTEGER BETA, T
141* ..
142*
143* Purpose
144* =======
145*
146* DLAMC1 determines the machine parameters given by BETA, T, RND, and
147* IEEE1.
148*
149* Arguments
150* =========
151*
152* BETA (output) INTEGER
153* The base of the machine.
154*
155* T (output) INTEGER
156* The number of ( BETA ) digits in the mantissa.
157*
158* RND (output) LOGICAL
159* Specifies whether proper rounding ( RND = .TRUE. ) or
160* chopping ( RND = .FALSE. ) occurs in addition. This may not
161* be a reliable guide to the way in which the machine performs
162* its arithmetic.
163*
164* IEEE1 (output) LOGICAL
165* Specifies whether rounding appears to be done in the IEEE
166* 'round to nearest' style.
167*
168* Further Details
169* ===============
170*
171* The routine is based on the routine ENVRON by Malcolm and
172* incorporates suggestions by Gentleman and Marovich. See
173*
174* Malcolm M. A. (1972) Algorithms to reveal properties of
175* floating-point arithmetic. Comms. of the ACM, 15, 949-951.
176*
177* Gentleman W. M. and Marovich S. B. (1974) More on algorithms
178* that reveal properties of floating point arithmetic units.
179* Comms. of the ACM, 17, 276-277.
180*
181* =====================================================================
182*
183* .. Local Scalars ..
184 LOGICAL FIRST, LIEEE1, LRND
185 INTEGER LBETA, LT
186 DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2
187* ..
188* .. External Functions ..
189 DOUBLE PRECISION DLAMC3
190 EXTERNAL dlamc3
191* ..
192* .. Save statement ..
193 SAVE first, lieee1, lbeta, lrnd, lt
194* ..
195* .. Data statements ..
196 DATA first / .true. /
197* ..
198* .. Executable Statements ..
199*
200 IF( first ) THEN
201 first = .false.
202 one = 1
203*
204* LBETA, LIEEE1, LT and LRND are the local values of BETA,
205* IEEE1, T and RND.
206*
207* Throughout this routine we use the function DLAMC3 to ensure
208* that relevant values are stored and not held in registers, or
209* are not affected by optimizers.
210*
211* Compute a = 2.0**m with the smallest positive integer m such
212* that
213*
214* fl( a + 1.0 ) = a.
215*
216 a = 1
217 c = 1
218*
219*+ WHILE( C.EQ.ONE )LOOP
220 10 CONTINUE
221 IF( c.EQ.one ) THEN
222 a = 2*a
223 c = dlamc3( a, one )
224 c = dlamc3( c, -a )
225 GO TO 10
226 END IF
227*+ END WHILE
228*
229* Now compute b = 2.0**m with the smallest positive integer m
230* such that
231*
232* fl( a + b ) .gt. a.
233*
234 b = 1
235 c = dlamc3( a, b )
236*
237*+ WHILE( C.EQ.A )LOOP
238 20 CONTINUE
239 IF( c.EQ.a ) THEN
240 b = 2*b
241 c = dlamc3( a, b )
242 GO TO 20
243 END IF
244*+ END WHILE
245*
246* Now compute the base. a and c are neighbouring floating point
247* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
248* their difference is beta. Adding 0.25 to c is to ensure that it
249* is truncated to beta and not ( beta - 1 ).
250*
251 qtr = one / 4
252 savec = c
253 c = dlamc3( c, -a )
254 lbeta = c + qtr
255*
256* Now determine whether rounding or chopping occurs, by adding a
257* bit less than beta/2 and a bit more than beta/2 to a.
258*
259 b = lbeta
260 f = dlamc3( b / 2, -b / 100 )
261 c = dlamc3( f, a )
262 IF( c.EQ.a ) THEN
263 lrnd = .true.
264 ELSE
265 lrnd = .false.
266 END IF
267 f = dlamc3( b / 2, b / 100 )
268 c = dlamc3( f, a )
269 IF( ( lrnd ) .AND. ( c.EQ.a ) )
270 $ lrnd = .false.
271*
272* Try and decide whether rounding is done in the IEEE 'round to
273* nearest' style. B/2 is half a unit in the last place of the two
274* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
275* zero, and SAVEC is odd. Thus adding B/2 to A should not change
276* A, but adding B/2 to SAVEC should change SAVEC.
277*
278 t1 = dlamc3( b / 2, a )
279 t2 = dlamc3( b / 2, savec )
280 lieee1 = ( t1.EQ.a ) .AND. ( t2.GT.savec ) .AND. lrnd
281*
282* Now find the mantissa, t. It should be the integer part of
283* log to the base beta of a, however it is safer to determine t
284* by powering. So we find t as the smallest positive integer for
285* which
286*
287* fl( beta**t + 1.0 ) = 1.0.
288*
289 lt = 0
290 a = 1
291 c = 1
292*
293*+ WHILE( C.EQ.ONE )LOOP
294 30 CONTINUE
295 IF( c.EQ.one ) THEN
296 lt = lt + 1
297 a = a*lbeta
298 c = dlamc3( a, one )
299 c = dlamc3( c, -a )
300 GO TO 30
301 END IF
302*+ END WHILE
303*
304 END IF
305*
306 beta = lbeta
307 t = lt
308 rnd = lrnd
309 ieee1 = lieee1
310 RETURN
311*
312* End of DLAMC1
313*
dlamc3
double precision function dlamc3(a, b)
Definition dlamch.f:578

◆ dlamc2()

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

Definition at line 318 of file dlamch.f.

319*
320* -- LAPACK auxiliary routine (version 2.1) --
321* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
322* Courant Institute, Argonne National Lab, and Rice University
323* October 31, 1992
324*
325* .. Scalar Arguments ..
326 LOGICAL RND
327 INTEGER BETA, EMAX, EMIN, T
328 DOUBLE PRECISION EPS, RMAX, RMIN
329* ..
330*
331* Purpose
332* =======
333*
334* DLAMC2 determines the machine parameters specified in its argument
335* list.
336*
337* Arguments
338* =========
339*
340* BETA (output) INTEGER
341* The base of the machine.
342*
343* T (output) INTEGER
344* The number of ( BETA ) digits in the mantissa.
345*
346* RND (output) LOGICAL
347* Specifies whether proper rounding ( RND = .TRUE. ) or
348* chopping ( RND = .FALSE. ) occurs in addition. This may not
349* be a reliable guide to the way in which the machine performs
350* its arithmetic.
351*
352* EPS (output) DOUBLE PRECISION
353* The smallest positive number such that
354*
355* fl( 1.0 - EPS ) .LT. 1.0,
356*
357* where fl denotes the computed value.
358*
359* EMIN (output) INTEGER
360* The minimum exponent before (gradual) underflow occurs.
361*
362* RMIN (output) DOUBLE PRECISION
363* The smallest normalized number for the machine, given by
364* BASE**( EMIN - 1 ), where BASE is the floating point value
365* of BETA.
366*
367* EMAX (output) INTEGER
368* The maximum exponent before overflow occurs.
369*
370* RMAX (output) DOUBLE PRECISION
371* The largest positive number for the machine, given by
372* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
373* value of BETA.
374*
375* Further Details
376* ===============
377*
378* The computation of EPS is based on a routine PARANOIA by
379* W. Kahan of the University of California at Berkeley.
380*
381* =====================================================================
382*
383* .. Local Scalars ..
384 LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND
385 INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
386 $ NGNMIN, NGPMIN
387 DOUBLE PRECISION A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
388 $ SIXTH, SMALL, THIRD, TWO, ZERO
389* ..
390* .. External Functions ..
391 DOUBLE PRECISION DLAMC3
392 EXTERNAL dlamc3
393* ..
394* .. External Subroutines ..
395 EXTERNAL dlamc1, dlamc4, dlamc5
396* ..
397* .. Intrinsic Functions ..
398 INTRINSIC abs, max, min
399* ..
400* .. Save statement ..
401 SAVE first, iwarn, lbeta, lemax, lemin, leps, lrmax,
402 $ lrmin, lt
403* ..
404* .. Data statements ..
405 DATA first / .true. / , iwarn / .false. /
406* ..
407* .. Executable Statements ..
408*
409 IF( first ) THEN
410 first = .false.
411 zero = 0
412 one = 1
413 two = 2
414*
415* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
416* BETA, T, RND, EPS, EMIN and RMIN.
417*
418* Throughout this routine we use the function DLAMC3 to ensure
419* that relevant values are stored and not held in registers, or
420* are not affected by optimizers.
421*
422* DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
423*
424 CALL dlamc1( lbeta, lt, lrnd, lieee1 )
425*
426* Start to find EPS.
427*
428 b = lbeta
429 a = b**( -lt )
430 leps = a
431*
432* Try some tricks to see whether or not this is the correct EPS.
433*
434 b = two / 3
435 half = one / 2
436 sixth = dlamc3( b, -half )
437 third = dlamc3( sixth, sixth )
438 b = dlamc3( third, -half )
439 b = dlamc3( b, sixth )
440 b = abs( b )
441 IF( b.LT.leps )
442 $ b = leps
443*
444 leps = 1
445*
446*+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
447 10 CONTINUE
448 IF( ( leps.GT.b ) .AND. ( b.GT.zero ) ) THEN
449 leps = b
450 c = dlamc3( half*leps, ( two**5 )*( leps**2 ) )
451 c = dlamc3( half, -c )
452 b = dlamc3( half, c )
453 c = dlamc3( half, -b )
454 b = dlamc3( half, c )
455 GO TO 10
456 END IF
457*+ END WHILE
458*
459 IF( a.LT.leps )
460 $ leps = a
461*
462* Computation of EPS complete.
463*
464* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
465* Keep dividing A by BETA until (gradual) underflow occurs. This
466* is detected when we cannot recover the previous A.
467*
468 rbase = one / lbeta
469 small = one
470 DO 20 i = 1, 3
471 small = dlamc3( small*rbase, zero )
472 20 CONTINUE
473 a = dlamc3( one, small )
474 CALL dlamc4( ngpmin, one, lbeta )
475 CALL dlamc4( ngnmin, -one, lbeta )
476 CALL dlamc4( gpmin, a, lbeta )
477 CALL dlamc4( gnmin, -a, lbeta )
478 ieee = .false.
479*
480 IF( ( ngpmin.EQ.ngnmin ) .AND. ( gpmin.EQ.gnmin ) ) THEN
481 IF( ngpmin.EQ.gpmin ) THEN
482 lemin = ngpmin
483* ( Non twos-complement machines, no gradual underflow;
484* e.g., VAX )
485 ELSE IF( ( gpmin-ngpmin ).EQ.3 ) THEN
486 lemin = ngpmin - 1 + lt
487 ieee = .true.
488* ( Non twos-complement machines, with gradual underflow;
489* e.g., IEEE standard followers )
490 ELSE
491 lemin = min( ngpmin, gpmin )
492* ( A guess; no known machine )
493 iwarn = .true.
494 END IF
495*
496 ELSE IF( ( ngpmin.EQ.gpmin ) .AND. ( ngnmin.EQ.gnmin ) ) THEN
497 IF( abs( ngpmin-ngnmin ).EQ.1 ) THEN
498 lemin = max( ngpmin, ngnmin )
499* ( Twos-complement machines, no gradual underflow;
500* e.g., CYBER 205 )
501 ELSE
502 lemin = min( ngpmin, ngnmin )
503* ( A guess; no known machine )
504 iwarn = .true.
505 END IF
506*
507 ELSE IF( ( abs( ngpmin-ngnmin ).EQ.1 ) .AND.
508 $ ( gpmin.EQ.gnmin ) ) THEN
509 IF( ( gpmin-min( ngpmin, ngnmin ) ).EQ.3 ) THEN
510 lemin = max( ngpmin, ngnmin ) - 1 + lt
511* ( Twos-complement machines with gradual underflow;
512* no known machine )
513 ELSE
514 lemin = min( ngpmin, ngnmin )
515* ( A guess; no known machine )
516 iwarn = .true.
517 END IF
518*
519 ELSE
520 lemin = min( ngpmin, ngnmin, gpmin, gnmin )
521* ( A guess; no known machine )
522 iwarn = .true.
523 END IF
524***
525* Comment out this if block if EMIN is ok
526 IF( iwarn ) THEN
527 first = .true.
528 WRITE( 6, fmt = 9999 )lemin
529 END IF
530***
531*
532* Assume IEEE arithmetic if we found denormalised numbers above,
533* or if arithmetic seems to round in the IEEE style, determined
534* in routine DLAMC1. A true IEEE machine should have both things
535* true; however, faulty machines may have one or the other.
536*
537 ieee = ieee .OR. lieee1
538*
539* Compute RMIN by successive division by BETA. We could compute
540* RMIN as BASE**( EMIN - 1 ), but some machines underflow during
541* this computation.
542*
543 lrmin = 1
544 DO 30 i = 1, 1 - lemin
545 lrmin = dlamc3( lrmin*rbase, zero )
546 30 CONTINUE
547*
548* Finally, call DLAMC5 to compute EMAX and RMAX.
549*
550 CALL dlamc5( lbeta, lt, lemin, ieee, lemax, lrmax )
551 END IF
552*
553 beta = lbeta
554 t = lt
555 rnd = lrnd
556 eps = leps
557 emin = lemin
558 rmin = lrmin
559 emax = lemax
560 rmax = lrmax
561*
562 RETURN
563*
564 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
565 $ ' emin = ', I8, /
566 $ ' If, after inspection, the value emin looks',
567 $ ' acceptable please comment out ',
568 $ / ' the IF block as marked within the code of routine',
569 $ ' dlamc2,', / ' otherwise supply emin explicitly.', / )
570*
571* End of DLAMC2
572*
the
end diagonal values have been computed in the(sparse) matrix id.SOL
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
dlamc4
subroutine dlamc4(emin, start, base)
Definition dlamch.f:616
dlamc1
subroutine dlamc1(beta, t, rnd, ieee1)
Definition dlamch.f:132
dlamc5
subroutine dlamc5(beta, p, emin, ieee, emax, rmax)
Definition dlamch.f:700
dlamc2
subroutine dlamc2(beta, t, rnd, eps, emin, rmin, emax, rmax)
Definition dlamch.f:319

◆ dlamc3()

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

Definition at line 577 of file dlamch.f.

578*
579* -- LAPACK auxiliary routine (version 2.1) --
580* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
581* Courant Institute, Argonne National Lab, and Rice University
582* October 31, 1992
583*
584* .. Scalar Arguments ..
585 DOUBLE PRECISION A, B
586* ..
587*
588* Purpose
589* =======
590*
591* DLAMC3 is intended to force A and B to be stored prior to doing
592* the addition of A and B , for use in situations where optimizers
593* might hold one of these in a register.
594*
595* Arguments
596* =========
597*
598* A, B (input) DOUBLE PRECISION
599* The values A and B.
600*
601* =====================================================================
602*
603* .. Executable Statements ..
604*
605 dlamc3 = a + b
606*
607 RETURN
608*
609* End of DLAMC3
610*

◆ dlamc4()

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

Definition at line 615 of file dlamch.f.

616*
617* -- LAPACK auxiliary routine (version 2.1) --
618* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
619* Courant Institute, Argonne National Lab, and Rice University
620* October 31, 1992
621*
622* .. Scalar Arguments ..
623 INTEGER BASE, EMIN
624 DOUBLE PRECISION START
625* ..
626*
627* Purpose
628* =======
629*
630* DLAMC4 is a service routine for DLAMC2.
631*
632* Arguments
633* =========
634*
635* EMIN (output) EMIN
636* The minimum exponent before (gradual) underflow, computed by
637* setting A = START and dividing by BASE until the previous A
638* can not be recovered.
639*
640* START (input) DOUBLE PRECISION
641* The starting point for determining EMIN.
642*
643* BASE (input) INTEGER
644* The base of the machine.
645*
646* =====================================================================
647*
648* .. Local Scalars ..
649 INTEGER I
650 DOUBLE PRECISION A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
651* ..
652* .. External Functions ..
653 DOUBLE PRECISION DLAMC3
654 EXTERNAL dlamc3
655* ..
656* .. Executable Statements ..
657*
658 a = start
659 one = 1
660 rbase = one / base
661 zero = 0
662 emin = 1
663 b1 = dlamc3( a*rbase, zero )
664 c1 = a
665 c2 = a
666 d1 = a
667 d2 = a
668*+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
669* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP
670 10 CONTINUE
671 IF( ( c1.EQ.a ) .AND. ( c2.EQ.a ) .AND. ( d1.EQ.a ) .AND.
672 $ ( d2.EQ.a ) ) THEN
673 emin = emin - 1
674 a = b1
675 b1 = dlamc3( a / base, zero )
676 c1 = dlamc3( b1*base, zero )
677 d1 = zero
678 DO 20 i = 1, base
679 d1 = d1 + b1
680 20 CONTINUE
681 b2 = dlamc3( a*rbase, zero )
682 c2 = dlamc3( b2 / rbase, zero )
683 d2 = zero
684 DO 30 i = 1, base
685 d2 = d2 + b2
686 30 CONTINUE
687 GO TO 10
688 END IF
689*+ END WHILE
690*
691 RETURN
692*
693* End of DLAMC4
694*

◆ dlamc5()

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

Definition at line 699 of file dlamch.f.

700*
701* -- LAPACK auxiliary routine (version 2.1) --
702* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
703* Courant Institute, Argonne National Lab, and Rice University
704* October 31, 1992
705*
706* .. Scalar Arguments ..
707 LOGICAL IEEE
708 INTEGER BETA, EMAX, EMIN, P
709 DOUBLE PRECISION RMAX
710* ..
711*
712* Purpose
713* =======
714*
715* DLAMC5 attempts to compute RMAX, the largest machine floating-point
716* number, without overflow. It assumes that EMAX + abs(EMIN) sum
717* approximately to a power of 2. It will fail on machines where this
718* assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
719* EMAX = 28718). It will also fail if the value supplied for EMIN is
720* too large (i.e. too close to zero), probably with overflow.
721*
722* Arguments
723* =========
724*
725* BETA (input) INTEGER
726* The base of floating-point arithmetic.
727*
728* P (input) INTEGER
729* The number of base BETA digits in the mantissa of a
730* floating-point value.
731*
732* EMIN (input) INTEGER
733* The minimum exponent before (gradual) underflow.
734*
735* IEEE (input) LOGICAL
736* A logical flag specifying whether or not the arithmetic
737* system is thought to comply with the IEEE standard.
738*
739* EMAX (output) INTEGER
740* The largest exponent before overflow
741*
742* RMAX (output) DOUBLE PRECISION
743* The largest machine floating-point number.
744*
745* =====================================================================
746*
747* .. Parameters ..
748 DOUBLE PRECISION ZERO, ONE
749 parameter( zero = 0.0d0, one = 1.0d0 )
750* ..
751* .. Local Scalars ..
752 INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
753 DOUBLE PRECISION OLDY, RECBAS, Y, Z
754* ..
755* .. External Functions ..
756 DOUBLE PRECISION DLAMC3
757 EXTERNAL dlamc3
758* ..
759* .. Intrinsic Functions ..
760 INTRINSIC mod
761* ..
762* .. Executable Statements ..
763*
764* First compute LEXP and UEXP, two powers of 2 that bound
765* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
766* approximately to the bound that is closest to abs(EMIN).
767* (EMAX is the exponent of the required number RMAX).
768*
769 lexp = 1
770 exbits = 1
771 10 CONTINUE
772 try = lexp*2
773 IF( try.LE.( -emin ) ) THEN
774 lexp = try
775 exbits = exbits + 1
776 GO TO 10
777 END IF
778 IF( lexp.EQ.-emin ) THEN
779 uexp = lexp
780 ELSE
781 uexp = try
782 exbits = exbits + 1
783 END IF
784*
785* Now -LEXP is less than or equal to EMIN, and -UEXP is greater
786* than or equal to EMIN. EXBITS is the number of bits needed to
787* store the exponent.
788*
789 IF( ( uexp+emin ).GT.( -lexp-emin ) ) THEN
790 expsum = 2*lexp
791 ELSE
792 expsum = 2*uexp
793 END IF
794*
795* EXPSUM is the exponent range, approximately equal to
796* EMAX - EMIN + 1 .
797*
798 emax = expsum + emin - 1
799 nbits = 1 + exbits + p
800*
801* NBITS is the total number of bits needed to store a
802* floating-point number.
803*
804 IF( ( mod( nbits, 2 ).EQ.1 ) .AND. ( beta.EQ.2 ) ) THEN
805*
806* Either there are an odd number of bits used to store a
807* floating-point number, which is unlikely, or some bits are
808* not used in the representation of numbers, which is possible,
809* (e.g. Cray machines) or the mantissa has an implicit bit,
810* (e.g. IEEE machines, Dec Vax machines), which is perhaps the
811* most likely. We have to assume the last alternative.
812* If this is true, then we need to reduce EMAX by one because
813* there must be some way of representing zero in an implicit-bit
814* system. On machines like Cray, we are reducing EMAX by one
815* unnecessarily.
816*
817 emax = emax - 1
818 END IF
819*
820 IF( ieee ) THEN
821*
822* Assume we are on an IEEE machine which reserves one exponent
823* for infinity and NaN.
824*
825 emax = emax - 1
826 END IF
827*
828* Now create RMAX, the largest machine number, which should
829* be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
830*
831* First compute 1.0 - BETA**(-P), being careful that the
832* result is less than 1.0 .
833*
834 recbas = one / beta
835 z = beta - one
836 y = zero
837 DO 20 i = 1, p
838 z = z*recbas
839 IF( y.LT.one )
840 $ oldy = y
841 y = dlamc3( y, z )
842 20 CONTINUE
843 IF( y.GE.one )
844 $ y = oldy
845*
846* Now multiply by BETA**EMAX to get RMAX.
847*
848 DO 30 i = 1, emax
849 y = dlamc3( y*beta, zero )
850 30 CONTINUE
851*
852 rmax = y
853 RETURN
854*
855* End of DLAMC5
856*

◆ dlamch()

double precision function dlamch ( character cmach)

Definition at line 1 of file dlamch.f.

2*
3* -- LAPACK auxiliary routine (version 2.1) --
4* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5* Courant Institute, Argonne National Lab, and Rice University
6* October 31, 1992
7*
8* .. Scalar Arguments ..
9 CHARACTER CMACH
10* ..
11*
12* Purpose
13* =======
14*
15* DLAMCH determines double precision machine parameters.
16*
17* Arguments
18* =========
19*
20* CMACH (input) CHARACTER*1
21* Specifies the value to be returned by DLAMCH:
22* = 'E' or 'e', DLAMCH := eps
23* = 'S' or 's , DLAMCH := sfmin
24* = 'B' or 'b', DLAMCH := base
25* = 'P' or 'p', DLAMCH := eps*base
26* = 'N' or 'n', DLAMCH := t
27* = 'R' or 'r', DLAMCH := rnd
28* = 'M' or 'm', DLAMCH := emin
29* = 'U' or 'u', DLAMCH := rmin
30* = 'L' or 'l', DLAMCH := emax
31* = 'O' or 'o', DLAMCH := rmax
32*
33* where
34*
35* eps = relative machine precision
36* sfmin = safe minimum, such that 1/sfmin does not overflow
37* base = base of the machine
38* prec = eps*base
39* t = number of (base) digits in the mantissa
40* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
41* emin = minimum exponent before (gradual) underflow
42* rmin = underflow threshold - base**(emin-1)
43* emax = largest exponent before overflow
44* rmax = overflow threshold - (base**emax)*(1-eps)
45*
46* =====================================================================
47*
48* .. Parameters ..
49 DOUBLE PRECISION ONE, ZERO
50 parameter( one = 1.0d+0, zero = 0.0d+0 )
51* ..
52* .. Local Scalars ..
53 LOGICAL FIRST, LRND
54 INTEGER BETA, IMAX, IMIN, IT
55 DOUBLE PRECISION BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
56 $ RND, SFMIN, SMALL, T
57* ..
58* .. External Functions ..
59 LOGICAL LSAME
60 EXTERNAL lsame
61* ..
62* .. External Subroutines ..
63 EXTERNAL dlamc2
64* ..
65* .. Save statement ..
66 SAVE first, eps, sfmin, base, t, rnd, emin, rmin,
67 $ emax, rmax, prec
68* ..
69* .. Data statements ..
70 DATA first / .true. /
71* ..
72* .. Executable Statements ..
73*
74 IF( first ) THEN
75 first = .false.
76 CALL dlamc2( beta, it, lrnd, eps, imin, rmin, imax, rmax )
77 base = beta
78 t = it
79 IF( lrnd ) THEN
80 rnd = one
81 eps = ( base**( 1-it ) ) / 2
82 ELSE
83 rnd = zero
84 eps = base**( 1-it )
85 END IF
86 prec = eps*base
87 emin = imin
88 emax = imax
89 sfmin = rmin
90 small = one / rmax
91 IF( small.GE.sfmin ) THEN
92*
93* Use SMALL plus a bit, to avoid the possibility of rounding
94* causing overflow when computing 1/sfmin.
95*
96 sfmin = small*( one+eps )
97 END IF
98 END IF
99*
100 IF( lsame( cmach, 'E' ) ) THEN
101 rmach = eps
102 ELSE IF( lsame( cmach, 'S' ) ) THEN
103 rmach = sfmin
104 ELSE IF( lsame( cmach, 'B' ) ) THEN
105 rmach = base
106 ELSE IF( lsame( cmach, 'P' ) ) THEN
107 rmach = prec
108 ELSE IF( lsame( cmach, 'N' ) ) THEN
109 rmach = t
110 ELSE IF( lsame( cmach, 'R' ) ) THEN
111 rmach = rnd
112 ELSE IF( lsame( cmach, 'M' ) ) THEN
113 rmach = emin
114 ELSE IF( lsame( cmach, 'U' ) ) THEN
115 rmach = rmin
116 ELSE IF( lsame( cmach, 'L' ) ) THEN
117 rmach = emax
118 ELSE IF( lsame( cmach, 'O' ) ) THEN
119 rmach = rmax
120 END IF
121*
122 dlamch = rmach
123 RETURN
124*
125* End of DLAMCH
126*
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
dlamch
double precision function dlamch(cmach)
Definition dlamch.f:2