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*
14*> \par Purpose:
15*  =============
16*>
17*> \verbatim
18*>
19*> DLAMCHF77 determines double precision machine parameters.
20*> \endverbatim
21*
22*  Arguments:
23*  ==========
24*
25*> \param[in] CMACH
26*> \verbatim
27*>          Specifies the value to be returned by DLAMCH:
28*>          = 'E' or 'e',   DLAMCH := eps
29*>          = 'S' or 's ,   DLAMCH := sfmin
30*>          = 'B' or 'b',   DLAMCH := base
31*>          = 'P' or 'p',   DLAMCH := eps*base
32*>          = 'N' or 'n',   DLAMCH := t
33*>          = 'R' or 'r',   DLAMCH := rnd
34*>          = 'M' or 'm',   DLAMCH := emin
35*>          = 'U' or 'u',   DLAMCH := rmin
36*>          = 'L' or 'l',   DLAMCH := emax
37*>          = 'O' or 'o',   DLAMCH := rmax
38*>          where
39*>          eps   = relative machine precision
40*>          sfmin = safe minimum, such that 1/sfmin does not overflow
41*>          base  = base of the machine
42*>          prec  = eps*base
43*>          t     = number of (base) digits in the mantissa
44*>          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise
45*>          emin  = minimum exponent before (gradual) underflow
46*>          rmin  = underflow threshold - base**(emin-1)
47*>          emax  = largest exponent before overflow
48*>          rmax  = overflow threshold  - (base**emax)*(1-eps)
49*> \endverbatim
50*
51*  Authors:
52*  ========
53*
54*> \author Univ. of Tennessee
55*> \author Univ. of California Berkeley
56*> \author Univ. of Colorado Denver
57*> \author NAG Ltd.
58*
59*> \date April 2012
60*
61*> \ingroup auxOTHERauxiliary
62*
63*  =====================================================================
64      DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
65*
66*  -- LAPACK auxiliary routine (version 3.4.1) --
67*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
68*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
69*     April 2012
70*
71*     .. Scalar Arguments ..
72      CHARACTER          CMACH
73*     ..
74*     .. Parameters ..
75      DOUBLE PRECISION   ONE, ZERO
76      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
77*     ..
78*     .. Local Scalars ..
79      LOGICAL            FIRST, LRND
80      INTEGER            BETA, IMAX, IMIN, IT
81      DOUBLE PRECISION   BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
82     $                   RND, SFMIN, SMALL, T
83*     ..
84*     .. External Functions ..
85      LOGICAL            LSAME
86      EXTERNAL           LSAME
87*     ..
88*     .. External Subroutines ..
89      EXTERNAL           DLAMC2
90*     ..
91*     .. Save statement ..
92      SAVE               FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN,
93     $                   EMAX, RMAX, PREC
94*     ..
95*     .. Data statements ..
96      DATA               FIRST / .TRUE. /
97*     ..
98*     .. Executable Statements ..
99*
100      IF( FIRST ) THEN
101         CALL DLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX )
102         BASE = BETA
103         T = IT
104         IF( LRND ) THEN
105            RND = ONE
106            EPS = ( BASE**( 1-IT ) ) / 2
107         ELSE
108            RND = ZERO
109            EPS = BASE**( 1-IT )
110         END IF
111         PREC = EPS*BASE
112         EMIN = IMIN
113         EMAX = IMAX
114         SFMIN = RMIN
115         SMALL = ONE / RMAX
116         IF( SMALL.GE.SFMIN ) THEN
117*
118*           Use SMALL plus a bit, to avoid the possibility of rounding
119*           causing overflow when computing  1/sfmin.
120*
121            SFMIN = SMALL*( ONE+EPS )
122         END IF
123      END IF
124*
125      IF( LSAME( CMACH, 'E' ) ) THEN
126         RMACH = EPS
127      ELSE IF( LSAME( CMACH, 'S' ) ) THEN
128         RMACH = SFMIN
129      ELSE IF( LSAME( CMACH, 'B' ) ) THEN
130         RMACH = BASE
131      ELSE IF( LSAME( CMACH, 'P' ) ) THEN
132         RMACH = PREC
133      ELSE IF( LSAME( CMACH, 'N' ) ) THEN
134         RMACH = T
135      ELSE IF( LSAME( CMACH, 'R' ) ) THEN
136         RMACH = RND
137      ELSE IF( LSAME( CMACH, 'M' ) ) THEN
138         RMACH = EMIN
139      ELSE IF( LSAME( CMACH, 'U' ) ) THEN
140         RMACH = RMIN
141      ELSE IF( LSAME( CMACH, 'L' ) ) THEN
142         RMACH = EMAX
143      ELSE IF( LSAME( CMACH, 'O' ) ) THEN
144         RMACH = RMAX
145      END IF
146*
147      DLAMCH = RMACH
148      FIRST  = .FALSE.
149      RETURN
150*
151*     End of DLAMCH
152*
153      END
154*
155************************************************************************
156*
157*> \brief \b DLAMC1
158*> \details
159*> \b Purpose:
160*> \verbatim
161*> DLAMC1 determines the machine parameters given by BETA, T, RND, and
162*> IEEE1.
163*> \endverbatim
164*>
165*> \param[out] BETA
166*> \verbatim
167*>          The base of the machine.
168*> \endverbatim
169*>
170*> \param[out] T
171*> \verbatim
172*>          The number of ( BETA ) digits in the mantissa.
173*> \endverbatim
174*>
175*> \param[out] RND
176*> \verbatim
177*>          Specifies whether proper rounding  ( RND = .TRUE. )  or
178*>          chopping  ( RND = .FALSE. )  occurs in addition. This may not
179*>          be a reliable guide to the way in which the machine performs
180*>          its arithmetic.
181*> \endverbatim
182*>
183*> \param[out] IEEE1
184*> \verbatim
185*>          Specifies whether rounding appears to be done in the IEEE
186*>          'round to nearest' style.
187*> \endverbatim
188*> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
189*> \date April 2012
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 DLAMC1( BETA, T, RND, IEEE1 )
207*
208*  -- LAPACK auxiliary routine (version 3.4.1) --
209*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
210*     November 2010
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*> \date April 2012
362*> \ingroup auxOTHERauxiliary
363*>
364*> \param[out] BETA
365*> \verbatim
366*>          The base of the machine.
367*> \endverbatim
368*>
369*> \param[out] T
370*> \verbatim
371*>          The number of ( BETA ) digits in the mantissa.
372*> \endverbatim
373*>
374*> \param[out] RND
375*> \verbatim
376*>          Specifies whether proper rounding  ( RND = .TRUE. )  or
377*>          chopping  ( RND = .FALSE. )  occurs in addition. This may not
378*>          be a reliable guide to the way in which the machine performs
379*>          its arithmetic.
380*> \endverbatim
381*>
382*> \param[out] EPS
383*> \verbatim
384*>          The smallest positive number such that
385*>             fl( 1.0 - EPS ) .LT. 1.0,
386*>          where fl denotes the computed value.
387*> \endverbatim
388*>
389*> \param[out] EMIN
390*> \verbatim
391*>          The minimum exponent before (gradual) underflow occurs.
392*> \endverbatim
393*>
394*> \param[out] RMIN
395*> \verbatim
396*>          The smallest normalized number for the machine, given by
397*>          BASE**( EMIN - 1 ), where  BASE  is the floating point value
398*>          of BETA.
399*> \endverbatim
400*>
401*> \param[out] EMAX
402*> \verbatim
403*>          The maximum exponent before overflow occurs.
404*> \endverbatim
405*>
406*> \param[out] RMAX
407*> \verbatim
408*>          The largest positive number for the machine, given by
409*>          BASE**EMAX * ( 1 - EPS ), where  BASE  is the floating point
410*>          value of BETA.
411*> \endverbatim
412*>
413*> \details \b Further \b Details
414*> \verbatim
415*>
416*>  The computation of  EPS  is based on a routine PARANOIA by
417*>  W. Kahan of the University of California at Berkeley.
418*> \endverbatim
419      SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
420*
421*  -- LAPACK auxiliary routine (version 3.4.1) --
422*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
423*     November 2010
424*
425*     .. Scalar Arguments ..
426      LOGICAL            RND
427      INTEGER            BETA, EMAX, EMIN, T
428      DOUBLE PRECISION   EPS, RMAX, RMIN
429*     ..
430* =====================================================================
431*
432*     .. Local Scalars ..
433      LOGICAL            FIRST, IEEE, IWARN, LIEEE1, LRND
434      INTEGER            GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
435     $                   NGNMIN, NGPMIN
436      DOUBLE PRECISION   A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
437     $                   SIXTH, SMALL, THIRD, TWO, ZERO
438*     ..
439*     .. External Functions ..
440      DOUBLE PRECISION   DLAMC3
441      EXTERNAL           DLAMC3
442*     ..
443*     .. External Subroutines ..
444      EXTERNAL           DLAMC1, DLAMC4, DLAMC5
445*     ..
446*     .. Intrinsic Functions ..
447      INTRINSIC          ABS, MAX, MIN
448*     ..
449*     .. Save statement ..
450      SAVE               FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,
451     $                   LRMIN, LT
452*     ..
453*     .. Data statements ..
454      DATA               FIRST / .TRUE. / , IWARN / .FALSE. /
455*     ..
456*     .. Executable Statements ..
457*
458      IF( FIRST ) THEN
459         ZERO = 0
460         ONE = 1
461         TWO = 2
462*
463*        LBETA, LT, LRND, LEPS, LEMIN and LRMIN  are the local values of
464*        BETA, T, RND, EPS, EMIN and RMIN.
465*
466*        Throughout this routine  we use the function  DLAMC3  to ensure
467*        that relevant values are stored  and not held in registers,  or
468*        are not affected by optimizers.
469*
470*        DLAMC1 returns the parameters  LBETA, LT, LRND and LIEEE1.
471*
472         CALL DLAMC1( LBETA, LT, LRND, LIEEE1 )
473*
474*        Start to find EPS.
475*
476         B = LBETA
477         A = B**( -LT )
478         LEPS = A
479*
480*        Try some tricks to see whether or not this is the correct  EPS.
481*
482         B = TWO / 3
483         HALF = ONE / 2
484         SIXTH = DLAMC3( B, -HALF )
485         THIRD = DLAMC3( SIXTH, SIXTH )
486         B = DLAMC3( THIRD, -HALF )
487         B = DLAMC3( B, SIXTH )
488         B = ABS( B )
489         IF( B.LT.LEPS )
490     $      B = LEPS
491*
492         LEPS = 1
493*
494*+       WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
495   10    CONTINUE
496         IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN
497            LEPS = B
498            C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
499            C = DLAMC3( HALF, -C )
500            B = DLAMC3( HALF, C )
501            C = DLAMC3( HALF, -B )
502            B = DLAMC3( HALF, C )
503            GO TO 10
504         END IF
505*+       END WHILE
506*
507         IF( A.LT.LEPS )
508     $      LEPS = A
509*
510*        Computation of EPS complete.
511*
512*        Now find  EMIN.  Let A = + or - 1, and + or - (1 + BASE**(-3)).
513*        Keep dividing  A by BETA until (gradual) underflow occurs. This
514*        is detected when we cannot recover the previous A.
515*
516         RBASE = ONE / LBETA
517         SMALL = ONE
518         DO 20 I = 1, 3
519            SMALL = DLAMC3( SMALL*RBASE, ZERO )
520   20    CONTINUE
521         A = DLAMC3( ONE, SMALL )
522         CALL DLAMC4( NGPMIN, ONE, LBETA )
523         CALL DLAMC4( NGNMIN, -ONE, LBETA )
524         CALL DLAMC4( GPMIN, A, LBETA )
525         CALL DLAMC4( GNMIN, -A, LBETA )
526         IEEE = .FALSE.
527*
528         IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN
529            IF( NGPMIN.EQ.GPMIN ) THEN
530               LEMIN = NGPMIN
531*            ( Non twos-complement machines, no gradual underflow;
532*              e.g.,  VAX )
533            ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
534               LEMIN = NGPMIN - 1 + LT
535               IEEE = .TRUE.
536*            ( Non twos-complement machines, with gradual underflow;
537*              e.g., IEEE standard followers )
538            ELSE
539               LEMIN = MIN( NGPMIN, GPMIN )
540*            ( A guess; no known machine )
541               IWARN = .TRUE.
542            END IF
543*
544         ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
545            IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
546               LEMIN = MAX( NGPMIN, NGNMIN )
547*            ( Twos-complement machines, no gradual underflow;
548*              e.g., CYBER 205 )
549            ELSE
550               LEMIN = MIN( NGPMIN, NGNMIN )
551*            ( A guess; no known machine )
552               IWARN = .TRUE.
553            END IF
554*
555         ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND.
556     $            ( GPMIN.EQ.GNMIN ) ) THEN
557            IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN
558               LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
559*            ( Twos-complement machines with gradual underflow;
560*              no known machine )
561            ELSE
562               LEMIN = MIN( NGPMIN, NGNMIN )
563*            ( A guess; no known machine )
564               IWARN = .TRUE.
565            END IF
566*
567         ELSE
568            LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
569*         ( A guess; no known machine )
570            IWARN = .TRUE.
571         END IF
572         FIRST = .FALSE.
573***
574* Comment out this if block if EMIN is ok
575         IF( IWARN ) THEN
576            FIRST = .TRUE.
577            WRITE( 6, FMT = 9999 )LEMIN
578         END IF
579***
580*
581*        Assume IEEE arithmetic if we found denormalised  numbers above,
582*        or if arithmetic seems to round in the  IEEE style,  determined
583*        in routine DLAMC1. A true IEEE machine should have both  things
584*        true; however, faulty machines may have one or the other.
585*
586         IEEE = IEEE .OR. LIEEE1
587*
588*        Compute  RMIN by successive division by  BETA. We could compute
589*        RMIN as BASE**( EMIN - 1 ),  but some machines underflow during
590*        this computation.
591*
592         LRMIN = 1
593         DO 30 I = 1, 1 - LEMIN
594            LRMIN = DLAMC3( LRMIN*RBASE, ZERO )
595   30    CONTINUE
596*
597*        Finally, call DLAMC5 to compute EMAX and RMAX.
598*
599         CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )
600      END IF
601*
602      BETA = LBETA
603      T = LT
604      RND = LRND
605      EPS = LEPS
606      EMIN = LEMIN
607      RMIN = LRMIN
608      EMAX = LEMAX
609      RMAX = LRMAX
610*
611      RETURN
612*
613 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
614     $      '  EMIN = ', I8, /
615     $      ' If, after inspection, the value EMIN looks',
616     $      ' acceptable please comment out ',
617     $      / ' the IF block as marked within the code of routine',
618     $      ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / )
619*
620*     End of DLAMC2
621*
622      END
623*
624************************************************************************
625*
626*> \brief \b DLAMC3
627*> \details
628*> \b Purpose:
629*> \verbatim
630*> DLAMC3  is intended to force  A  and  B  to be stored prior to doing
631*> the addition of  A  and  B ,  for use in situations where optimizers
632*> might hold one of these in a register.
633*> \endverbatim
634*>
635*> \param[in] A
636*>
637*> \param[in] B
638*> \verbatim
639*>          The values A and B.
640*> \endverbatim
641
642      DOUBLE PRECISION FUNCTION DLAMC3( A, B )
643*
644*  -- LAPACK auxiliary routine (version 3.4.1) --
645*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
646*     November 2010
647*
648*     .. Scalar Arguments ..
649      DOUBLE PRECISION   A, B
650*     ..
651* =====================================================================
652*
653*     .. Executable Statements ..
654*
655      DLAMC3 = A + B
656*
657      RETURN
658*
659*     End of DLAMC3
660*
661      END
662*
663************************************************************************
664*
665*> \brief \b DLAMC4
666*> \details
667*> \b Purpose:
668*> \verbatim
669*> DLAMC4 is a service routine for DLAMC2.
670*> \endverbatim
671*>
672*> \param[out] EMIN
673*> \verbatim
674*>          The minimum exponent before (gradual) underflow, computed by
675*>          setting A = START and dividing by BASE until the previous A
676*>          can not be recovered.
677*> \endverbatim
678*>
679*> \param[in] START
680*> \verbatim
681*>          The starting point for determining EMIN.
682*> \endverbatim
683*>
684*> \param[in] BASE
685*> \verbatim
686*>          The base of the machine.
687*> \endverbatim
688*>
689      SUBROUTINE DLAMC4( EMIN, START, BASE )
690*
691*  -- LAPACK auxiliary routine (version 3.4.1) --
692*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
693*     November 2010
694*
695*     .. Scalar Arguments ..
696      INTEGER            BASE, EMIN
697      DOUBLE PRECISION   START
698*     ..
699* =====================================================================
700*
701*     .. Local Scalars ..
702      INTEGER            I
703      DOUBLE PRECISION   A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
704*     ..
705*     .. External Functions ..
706      DOUBLE PRECISION   DLAMC3
707      EXTERNAL           DLAMC3
708*     ..
709*     .. Executable Statements ..
710*
711      A = START
712      ONE = 1
713      RBASE = ONE / BASE
714      ZERO = 0
715      EMIN = 1
716      B1 = DLAMC3( A*RBASE, ZERO )
717      C1 = A
718      C2 = A
719      D1 = A
720      D2 = A
721*+    WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
722*    $       ( D1.EQ.A ).AND.( D2.EQ.A )      )LOOP
723   10 CONTINUE
724      IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND.
725     $    ( D2.EQ.A ) ) THEN
726         EMIN = EMIN - 1
727         A = B1
728         B1 = DLAMC3( A / BASE, ZERO )
729         C1 = DLAMC3( B1*BASE, ZERO )
730         D1 = ZERO
731         DO 20 I = 1, BASE
732            D1 = D1 + B1
733   20    CONTINUE
734         B2 = DLAMC3( A*RBASE, ZERO )
735         C2 = DLAMC3( B2 / RBASE, ZERO )
736         D2 = ZERO
737         DO 30 I = 1, BASE
738            D2 = D2 + B2
739   30    CONTINUE
740         GO TO 10
741      END IF
742*+    END WHILE
743*
744      RETURN
745*
746*     End of DLAMC4
747*
748      END
749*
750************************************************************************
751*
752*> \brief \b DLAMC5
753*> \details
754*> \b Purpose:
755*> \verbatim
756*> DLAMC5 attempts to compute RMAX, the largest machine floating-point
757*> number, without overflow.  It assumes that EMAX + abs(EMIN) sum
758*> approximately to a power of 2.  It will fail on machines where this
759*> assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
760*> EMAX = 28718).  It will also fail if the value supplied for EMIN is
761*> too large (i.e. too close to zero), probably with overflow.
762*> \endverbatim
763*>
764*> \param[in] BETA
765*> \verbatim
766*>          The base of floating-point arithmetic.
767*> \endverbatim
768*>
769*> \param[in] P
770*> \verbatim
771*>          The number of base BETA digits in the mantissa of a
772*>          floating-point value.
773*> \endverbatim
774*>
775*> \param[in] EMIN
776*> \verbatim
777*>          The minimum exponent before (gradual) underflow.
778*> \endverbatim
779*>
780*> \param[in] IEEE
781*> \verbatim
782*>          A logical flag specifying whether or not the arithmetic
783*>          system is thought to comply with the IEEE standard.
784*> \endverbatim
785*>
786*> \param[out] EMAX
787*> \verbatim
788*>          The largest exponent before overflow
789*> \endverbatim
790*>
791*> \param[out] RMAX
792*> \verbatim
793*>          The largest machine floating-point number.
794*> \endverbatim
795*>
796      SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
797*
798*  -- LAPACK auxiliary routine (version 3.4.1) --
799*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
800*     November 2010
801*
802*     .. Scalar Arguments ..
803      LOGICAL            IEEE
804      INTEGER            BETA, EMAX, EMIN, P
805      DOUBLE PRECISION   RMAX
806*     ..
807* =====================================================================
808*
809*     .. Parameters ..
810      DOUBLE PRECISION   ZERO, ONE
811      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
812*     ..
813*     .. Local Scalars ..
814      INTEGER            EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
815      DOUBLE PRECISION   OLDY, RECBAS, Y, Z
816*     ..
817*     .. External Functions ..
818      DOUBLE PRECISION   DLAMC3
819      EXTERNAL           DLAMC3
820*     ..
821*     .. Intrinsic Functions ..
822      INTRINSIC          MOD
823*     ..
824*     .. Executable Statements ..
825*
826*     First compute LEXP and UEXP, two powers of 2 that bound
827*     abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
828*     approximately to the bound that is closest to abs(EMIN).
829*     (EMAX is the exponent of the required number RMAX).
830*
831      LEXP = 1
832      EXBITS = 1
833   10 CONTINUE
834      TRY = LEXP*2
835      IF( TRY.LE.( -EMIN ) ) THEN
836         LEXP = TRY
837         EXBITS = EXBITS + 1
838         GO TO 10
839      END IF
840      IF( LEXP.EQ.-EMIN ) THEN
841         UEXP = LEXP
842      ELSE
843         UEXP = TRY
844         EXBITS = EXBITS + 1
845      END IF
846*
847*     Now -LEXP is less than or equal to EMIN, and -UEXP is greater
848*     than or equal to EMIN. EXBITS is the number of bits needed to
849*     store the exponent.
850*
851      IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
852         EXPSUM = 2*LEXP
853      ELSE
854         EXPSUM = 2*UEXP
855      END IF
856*
857*     EXPSUM is the exponent range, approximately equal to
858*     EMAX - EMIN + 1 .
859*
860      EMAX = EXPSUM + EMIN - 1
861      NBITS = 1 + EXBITS + P
862*
863*     NBITS is the total number of bits needed to store a
864*     floating-point number.
865*
866      IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
867*
868*        Either there are an odd number of bits used to store a
869*        floating-point number, which is unlikely, or some bits are
870*        not used in the representation of numbers, which is possible,
871*        (e.g. Cray machines) or the mantissa has an implicit bit,
872*        (e.g. IEEE machines, Dec Vax machines), which is perhaps the
873*        most likely. We have to assume the last alternative.
874*        If this is true, then we need to reduce EMAX by one because
875*        there must be some way of representing zero in an implicit-bit
876*        system. On machines like Cray, we are reducing EMAX by one
877*        unnecessarily.
878*
879         EMAX = EMAX - 1
880      END IF
881*
882      IF( IEEE ) THEN
883*
884*        Assume we are on an IEEE machine which reserves one exponent
885*        for infinity and NaN.
886*
887         EMAX = EMAX - 1
888      END IF
889*
890*     Now create RMAX, the largest machine number, which should
891*     be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
892*
893*     First compute 1.0 - BETA**(-P), being careful that the
894*     result is less than 1.0 .
895*
896      RECBAS = ONE / BETA
897      Z = BETA - ONE
898      Y = ZERO
899      DO 20 I = 1, P
900         Z = Z*RECBAS
901         IF( Y.LT.ONE )
902     $      OLDY = Y
903         Y = DLAMC3( Y, Z )
904   20 CONTINUE
905      IF( Y.GE.ONE )
906     $   Y = OLDY
907*
908*     Now multiply by BETA**EMAX to get RMAX.
909*
910      DO 30 I = 1, EMAX
911         Y = DLAMC3( Y*BETA, ZERO )
912   30 CONTINUE
913*
914      RMAX = Y
915      RETURN
916*
917*     End of DLAMC5
918*
919      END
920