LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
slamchf77.f
Go to the documentation of this file.
1*> \brief \b SLAMCHF77 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* REAL FUNCTION SLAMCH( CMACH )
12*
13* .. Scalar Arguments ..
14* CHARACTER CMACH
15* ..
16*
17*
18*> \par Purpose:
19* =============
20*>
21*> \verbatim
22*>
23*> SLAMCH determines single precision machine parameters.
24*> \endverbatim
25*
26* Arguments:
27* ==========
28*
29*> \param[in] CMACH
30*> \verbatim
31*> Specifies the value to be returned by SLAMCH:
32*> = 'E' or 'e', SLAMCH := eps
33*> = 'S' or 's , SLAMCH := sfmin
34*> = 'B' or 'b', SLAMCH := base
35*> = 'P' or 'p', SLAMCH := eps*base
36*> = 'N' or 'n', SLAMCH := t
37*> = 'R' or 'r', SLAMCH := rnd
38*> = 'M' or 'm', SLAMCH := emin
39*> = 'U' or 'u', SLAMCH := rmin
40*> = 'L' or 'l', SLAMCH := emax
41*> = 'O' or 'o', SLAMCH := 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*> \ingroup auxOTHERauxiliary
64*
65* =====================================================================
66 REAL function slamch( cmach )
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*
154 END
155*
156************************************************************************
157*
158*> \brief \b SLAMC1
159*> \details
160*> \b Purpose:
161*> \verbatim
162*> SLAMC1 determines the machine parameters given by BETA, T, RND, and
163*> IEEE1.
164*> \endverbatim
165*>
166*> \param[out] BETA
167*> \verbatim
168*> The base of the machine.
169*> \endverbatim
170*>
171*> \param[out] T
172*> \verbatim
173*> The number of ( BETA ) digits in the mantissa.
174*> \endverbatim
175*>
176*> \param[out] RND
177*> \verbatim
178*> Specifies whether proper rounding ( RND = .TRUE. ) or
179*> chopping ( RND = .FALSE. ) occurs in addition. This may not
180*> be a reliable guide to the way in which the machine performs
181*> its arithmetic.
182*> \endverbatim
183*>
184*> \param[out] IEEE1
185*> \verbatim
186*> Specifies whether rounding appears to be done in the IEEE
187*> 'round to nearest' style.
188*> \endverbatim
189*> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
190*> \ingroup auxOTHERauxiliary
191*>
192*> \details \b Further \b Details
193*> \verbatim
194*>
195*> The routine is based on the routine ENVRON by Malcolm and
196*> incorporates suggestions by Gentleman and Marovich. See
197*>
198*> Malcolm M. A. (1972) Algorithms to reveal properties of
199*> floating-point arithmetic. Comms. of the ACM, 15, 949-951.
200*>
201*> Gentleman W. M. and Marovich S. B. (1974) More on algorithms
202*> that reveal properties of floating point arithmetic units.
203*> Comms. of the ACM, 17, 276-277.
204*> \endverbatim
205*>
206 SUBROUTINE slamc1( BETA, T, RND, IEEE1 )
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*
348 END
349*
350************************************************************************
351*
352*> \brief \b SLAMC2
353*> \details
354*> \b Purpose:
355*> \verbatim
356*> SLAMC2 determines the machine parameters specified in its argument
357*> list.
358*> \endverbatim
359*> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
360*> \ingroup auxOTHERauxiliary
361*>
362*> \param[out] BETA
363*> \verbatim
364*> The base of the machine.
365*> \endverbatim
366*>
367*> \param[out] T
368*> \verbatim
369*> The number of ( BETA ) digits in the mantissa.
370*> \endverbatim
371*>
372*> \param[out] RND
373*> \verbatim
374*> Specifies whether proper rounding ( RND = .TRUE. ) or
375*> chopping ( RND = .FALSE. ) occurs in addition. This may not
376*> be a reliable guide to the way in which the machine performs
377*> its arithmetic.
378*> \endverbatim
379*>
380*> \param[out] EPS
381*> \verbatim
382*> The smallest positive number such that
383*> fl( 1.0 - EPS ) .LT. 1.0,
384*> where fl denotes the computed value.
385*> \endverbatim
386*>
387*> \param[out] EMIN
388*> \verbatim
389*> The minimum exponent before (gradual) underflow occurs.
390*> \endverbatim
391*>
392*> \param[out] RMIN
393*> \verbatim
394*> The smallest normalized number for the machine, given by
395*> BASE**( EMIN - 1 ), where BASE is the floating point value
396*> of BETA.
397*> \endverbatim
398*>
399*> \param[out] EMAX
400*> \verbatim
401*> The maximum exponent before overflow occurs.
402*> \endverbatim
403*>
404*> \param[out] RMAX
405*> \verbatim
406*> The largest positive number for the machine, given by
407*> BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
408*> value of BETA.
409*> \endverbatim
410*>
411*> \details \b Further \b Details
412*> \verbatim
413*>
414*> The computation of EPS is based on a routine PARANOIA by
415*> W. Kahan of the University of California at Berkeley.
416*> \endverbatim
417 SUBROUTINE slamc2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
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*
619 END
620*
621************************************************************************
622*
623*> \brief \b SLAMC3
624*> \details
625*> \b Purpose:
626*> \verbatim
627*> SLAMC3 is intended to force A and B to be stored prior to doing
628*> the addition of A and B , for use in situations where optimizers
629*> might hold one of these in a register.
630*> \endverbatim
631*>
632*> \param[in] A
633*>
634*> \param[in] B
635*> \verbatim
636*> The values A and B.
637*> \endverbatim
638
639 REAL function slamc3( a, b )
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*
657 END
658*
659************************************************************************
660*
661*> \brief \b SLAMC4
662*> \details
663*> \b Purpose:
664*> \verbatim
665*> SLAMC4 is a service routine for SLAMC2.
666*> \endverbatim
667*>
668*> \param[out] EMIN
669*> \verbatim
670*> The minimum exponent before (gradual) underflow, computed by
671*> setting A = START and dividing by BASE until the previous A
672*> can not be recovered.
673*> \endverbatim
674*>
675*> \param[in] START
676*> \verbatim
677*> The starting point for determining EMIN.
678*> \endverbatim
679*>
680*> \param[in] BASE
681*> \verbatim
682*> The base of the machine.
683*> \endverbatim
684*>
685 SUBROUTINE slamc4( EMIN, START, BASE )
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*
744 END
745*
746************************************************************************
747*
748*> \brief \b SLAMC5
749*> \details
750*> \b Purpose:
751*> \verbatim
752*> SLAMC5 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 slamc5( 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 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*
914 END
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
slamc3
real function slamc3(A, B)
SLAMC3
Definition: slamch.f:169
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
slamc1
subroutine slamc1(BETA, T, RND, IEEE1)
SLAMC1
Definition: slamchf77.f:207
slamc5
subroutine slamc5(BETA, P, EMIN, IEEE, EMAX, RMAX)
SLAMC5
Definition: slamchf77.f:793
slamc2
subroutine slamc2(BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX)
SLAMC2
Definition: slamchf77.f:418
slamc4
subroutine slamc4(EMIN, START, BASE)
SLAMC4
Definition: slamchf77.f:686