LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlamchf77.f
Go to the documentation of this file.
1*> \brief \b DLAMCHF77 deprecated
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
12*
13* .. Scalar Arguments ..
14* CHARACTER CMACH
15* ..
16*
17*
18*> \par Purpose:
19* =============
20*>
21*> \verbatim
22*>
23*> DLAMCHF77 determines double precision machine parameters.
24*> \endverbatim
25*
26* Arguments:
27* ==========
28*
29*> \param[in] CMACH
30*> \verbatim
31*> Specifies the value to be returned by DLAMCH:
32*> = 'E' or 'e', DLAMCH := eps
33*> = 'S' or 's , DLAMCH := sfmin
34*> = 'B' or 'b', DLAMCH := base
35*> = 'P' or 'p', DLAMCH := eps*base
36*> = 'N' or 'n', DLAMCH := t
37*> = 'R' or 'r', DLAMCH := rnd
38*> = 'M' or 'm', DLAMCH := emin
39*> = 'U' or 'u', DLAMCH := rmin
40*> = 'L' or 'l', DLAMCH := emax
41*> = 'O' or 'o', DLAMCH := rmax
42*> where
43*> eps = relative machine precision
44*> sfmin = safe minimum, such that 1/sfmin does not overflow
45*> base = base of the machine
46*> prec = eps*base
47*> t = number of (base) digits in the mantissa
48*> rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
49*> emin = minimum exponent before (gradual) underflow
50*> rmin = underflow threshold - base**(emin-1)
51*> emax = largest exponent before overflow
52*> rmax = overflow threshold - (base**emax)*(1-eps)
53*> \endverbatim
54*
55* Authors:
56* ========
57*
58*> \author Univ. of Tennessee
59*> \author Univ. of California Berkeley
60*> \author Univ. of Colorado Denver
61*> \author NAG Ltd.
62*
63
64*> \ingroup auxOTHERauxiliary
65*
66* =====================================================================
67 DOUBLE PRECISION FUNCTION dlamch( CMACH )
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*
155 END
156*
157************************************************************************
158*
159*> \brief \b DLAMC1
160*> \details
161*> \b Purpose:
162*> \verbatim
163*> DLAMC1 determines the machine parameters given by BETA, T, RND, and
164*> IEEE1.
165*> \endverbatim
166*>
167*> \param[out] BETA
168*> \verbatim
169*> The base of the machine.
170*> \endverbatim
171*>
172*> \param[out] T
173*> \verbatim
174*> The number of ( BETA ) digits in the mantissa.
175*> \endverbatim
176*>
177*> \param[out] RND
178*> \verbatim
179*> Specifies whether proper rounding ( RND = .TRUE. ) or
180*> chopping ( RND = .FALSE. ) occurs in addition. This may not
181*> be a reliable guide to the way in which the machine performs
182*> its arithmetic.
183*> \endverbatim
184*>
185*> \param[out] IEEE1
186*> \verbatim
187*> Specifies whether rounding appears to be done in the IEEE
188*> 'round to nearest' style.
189*> \endverbatim
190*> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
191*> \ingroup auxOTHERauxiliary
192*>
193*> \details \b Further \b Details
194*> \verbatim
195*>
196*> The routine is based on the routine ENVRON by Malcolm and
197*> incorporates suggestions by Gentleman and Marovich. See
198*>
199*> Malcolm M. A. (1972) Algorithms to reveal properties of
200*> floating-point arithmetic. Comms. of the ACM, 15, 949-951.
201*>
202*> Gentleman W. M. and Marovich S. B. (1974) More on algorithms
203*> that reveal properties of floating point arithmetic units.
204*> Comms. of the ACM, 17, 276-277.
205*> \endverbatim
206*>
207 SUBROUTINE dlamc1( BETA, T, RND, IEEE1 )
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*
349 END
350*
351************************************************************************
352*
353*> \brief \b DLAMC2
354*> \details
355*> \b Purpose:
356*> \verbatim
357*> DLAMC2 determines the machine parameters specified in its argument
358*> list.
359*> \endverbatim
360*> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
361*> \ingroup auxOTHERauxiliary
362*>
363*> \param[out] BETA
364*> \verbatim
365*> The base of the machine.
366*> \endverbatim
367*>
368*> \param[out] T
369*> \verbatim
370*> The number of ( BETA ) digits in the mantissa.
371*> \endverbatim
372*>
373*> \param[out] RND
374*> \verbatim
375*> Specifies whether proper rounding ( RND = .TRUE. ) or
376*> chopping ( RND = .FALSE. ) occurs in addition. This may not
377*> be a reliable guide to the way in which the machine performs
378*> its arithmetic.
379*> \endverbatim
380*>
381*> \param[out] EPS
382*> \verbatim
383*> The smallest positive number such that
384*> fl( 1.0 - EPS ) .LT. 1.0,
385*> where fl denotes the computed value.
386*> \endverbatim
387*>
388*> \param[out] EMIN
389*> \verbatim
390*> The minimum exponent before (gradual) underflow occurs.
391*> \endverbatim
392*>
393*> \param[out] RMIN
394*> \verbatim
395*> The smallest normalized number for the machine, given by
396*> BASE**( EMIN - 1 ), where BASE is the floating point value
397*> of BETA.
398*> \endverbatim
399*>
400*> \param[out] EMAX
401*> \verbatim
402*> The maximum exponent before overflow occurs.
403*> \endverbatim
404*>
405*> \param[out] RMAX
406*> \verbatim
407*> The largest positive number for the machine, given by
408*> BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
409*> value of BETA.
410*> \endverbatim
411*>
412*> \details \b Further \b Details
413*> \verbatim
414*>
415*> The computation of EPS is based on a routine PARANOIA by
416*> W. Kahan of the University of California at Berkeley.
417*> \endverbatim
418 SUBROUTINE dlamc2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
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*
620 END
621*
622************************************************************************
623*
624*> \brief \b DLAMC3
625*> \details
626*> \b Purpose:
627*> \verbatim
628*> DLAMC3 is intended to force A and B to be stored prior to doing
629*> the addition of A and B , for use in situations where optimizers
630*> might hold one of these in a register.
631*> \endverbatim
632*>
633*> \param[in] A
634*>
635*> \param[in] B
636*> \verbatim
637*> The values A and B.
638*> \endverbatim
639
640 DOUBLE PRECISION FUNCTION dlamc3( A, B )
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*
658 END
659*
660************************************************************************
661*
662*> \brief \b DLAMC4
663*> \details
664*> \b Purpose:
665*> \verbatim
666*> DLAMC4 is a service routine for DLAMC2.
667*> \endverbatim
668*>
669*> \param[out] EMIN
670*> \verbatim
671*> The minimum exponent before (gradual) underflow, computed by
672*> setting A = START and dividing by BASE until the previous A
673*> can not be recovered.
674*> \endverbatim
675*>
676*> \param[in] START
677*> \verbatim
678*> The starting point for determining EMIN.
679*> \endverbatim
680*>
681*> \param[in] BASE
682*> \verbatim
683*> The base of the machine.
684*> \endverbatim
685*>
686 SUBROUTINE dlamc4( EMIN, START, BASE )
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*
744 END
745*
746************************************************************************
747*
748*> \brief \b DLAMC5
749*> \details
750*> \b Purpose:
751*> \verbatim
752*> DLAMC5 attempts to compute RMAX, the largest machine floating-point
753*> number, without overflow. It assumes that EMAX + abs(EMIN) sum
754*> approximately to a power of 2. It will fail on machines where this
755*> assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
756*> EMAX = 28718). It will also fail if the value supplied for EMIN is
757*> too large (i.e. too close to zero), probably with overflow.
758*> \endverbatim
759*>
760*> \param[in] BETA
761*> \verbatim
762*> The base of floating-point arithmetic.
763*> \endverbatim
764*>
765*> \param[in] P
766*> \verbatim
767*> The number of base BETA digits in the mantissa of a
768*> floating-point value.
769*> \endverbatim
770*>
771*> \param[in] EMIN
772*> \verbatim
773*> The minimum exponent before (gradual) underflow.
774*> \endverbatim
775*>
776*> \param[in] IEEE
777*> \verbatim
778*> A logical flag specifying whether or not the arithmetic
779*> system is thought to comply with the IEEE standard.
780*> \endverbatim
781*>
782*> \param[out] EMAX
783*> \verbatim
784*> The largest exponent before overflow
785*> \endverbatim
786*>
787*> \param[out] RMAX
788*> \verbatim
789*> The largest machine floating-point number.
790*> \endverbatim
791*>
792 SUBROUTINE dlamc5( BETA, P, EMIN, IEEE, EMAX, RMAX )
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*
914 END
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
dlamc3
double precision function dlamc3(A, B)
DLAMC3
Definition: dlamch.f:169
dlamc2
subroutine dlamc2(BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX)
DLAMC2
Definition: dlamchf77.f:419
dlamc5
subroutine dlamc5(BETA, P, EMIN, IEEE, EMAX, RMAX)
DLAMC5
Definition: dlamchf77.f:793
dlamc1
subroutine dlamc1(BETA, T, RND, IEEE1)
DLAMC1
Definition: dlamchf77.f:208
dlamc4
subroutine dlamc4(EMIN, START, BASE)
DLAMC4
Definition: dlamchf77.f:687
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53