1*          where
2*
3*          eps   = relative machine precision
4*          sfmin = safe minimum, such that 1/sfmin does not overflow
5*          base  = base of the machine
6*          prec  = eps*base
7*          t     = number of (base) digits in the mantissa
8*          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise
9*          emin  = minimum exponent before (gradual) underflow
10*          rmin  = underflow threshold - base**(emin-1)
11*          emax  = largest exponent before overflow
12*          rmax  = overflow threshold  - (base**emax)*(1-eps)
13*
14* =====================================================================
15*
16*     .. Parameters ..
17      DOUBLE PRECISION   ONE, ZERO
18      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
19*     ..
20*     .. Local Scalars ..
21      LOGICAL            FIRST, LRND
22      INTEGER            BETA, IMAX, IMIN, IT
23      DOUBLE PRECISION   BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
24     $                   RND, SFMIN, SMALL, T
25*     ..
26*     .. External Functions ..
27      LOGICAL            LSAME
28      EXTERNAL           LSAME
29*     ..
30*     .. External Subroutines ..
31      EXTERNAL           DLAMC2
32*     ..
33*     .. Save statement ..
34      SAVE               FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN,
35     $                   EMAX, RMAX, PREC
36*     ..
37*     .. Data statements ..
38      DATA               FIRST / .TRUE. /
39*     ..
40*     .. Executable Statements ..
41*
42      IF( FIRST ) THEN
43         CALL DLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX )
44         BASE = BETA
45         T = IT
46         IF( LRND ) THEN
47            RND = ONE
48            EPS = ( BASE**( 1-IT ) ) / 2
49         ELSE
50            RND = ZERO
51            EPS = BASE**( 1-IT )
52         END IF
53         PREC = EPS*BASE
54         EMIN = IMIN
55         EMAX = IMAX
56         SFMIN = RMIN
57         SMALL = ONE / RMAX
58         IF( SMALL.GE.SFMIN ) THEN
59*
60*           Use SMALL plus a bit, to avoid the possibility of rounding
61*           causing overflow when computing  1/sfmin.
62*
63            SFMIN = SMALL*( ONE+EPS )
64         END IF
65      END IF
66*
67      IF( LSAME( CMACH, 'E' ) ) THEN
68         RMACH = EPS
69      ELSE IF( LSAME( CMACH, 'S' ) ) THEN
70         RMACH = SFMIN
71      ELSE IF( LSAME( CMACH, 'B' ) ) THEN
72         RMACH = BASE
73      ELSE IF( LSAME( CMACH, 'P' ) ) THEN
74         RMACH = PREC
75      ELSE IF( LSAME( CMACH, 'N' ) ) THEN
76         RMACH = T
77      ELSE IF( LSAME( CMACH, 'R' ) ) THEN
78         RMACH = RND
79      ELSE IF( LSAME( CMACH, 'M' ) ) THEN
80         RMACH = EMIN
81      ELSE IF( LSAME( CMACH, 'U' ) ) THEN
82         RMACH = RMIN
83      ELSE IF( LSAME( CMACH, 'L' ) ) THEN
84         RMACH = EMAX
85      ELSE IF( LSAME( CMACH, 'O' ) ) THEN
86         RMACH = RMAX
87      END IF
88*
89      DLAMCH = RMACH
90      FIRST  = .FALSE.
91      RETURN
92*
93*     End of DLAMCH
94*
95      END
96*
97************************************************************************
98*
99      SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 )
100*
101*  -- LAPACK auxiliary routine (version 3.2) --
102*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
103*     November 2006
104*
105*     .. Scalar Arguments ..
106      LOGICAL            IEEE1, RND
107      INTEGER            BETA, T
108*     ..
109*
110*  Purpose
111*  =======
112*
113*  DLAMC1 determines the machine parameters given by BETA, T, RND, and
114*  IEEE1.
115*
116*  Arguments
117*  =========
118*
119*  BETA    (output) INTEGER
120*          The base of the machine.
121*
122*  T       (output) INTEGER
123*          The number of ( BETA ) digits in the mantissa.
124*
125*  RND     (output) LOGICAL
126*          Specifies whether proper rounding  ( RND = .TRUE. )  or
127*          chopping  ( RND = .FALSE. )  occurs in addition. This may not
128*          be a reliable guide to the way in which the machine performs
129*          its arithmetic.
130*
131*  IEEE1   (output) LOGICAL
132*          Specifies whether rounding appears to be done in the IEEE
133*          'round to nearest' style.
134*
135*  Further Details
136*  ===============
137*
138*  The routine is based on the routine  ENVRON  by Malcolm and
139*  incorporates suggestions by Gentleman and Marovich. See
140*
141*     Malcolm M. A. (1972) Algorithms to reveal properties of
142*        floating-point arithmetic. Comms. of the ACM, 15, 949-951.
143*
144*     Gentleman W. M. and Marovich S. B. (1974) More on algorithms
145*        that reveal properties of floating point arithmetic units.
146*        Comms. of the ACM, 17, 276-277.
147*
148* =====================================================================
149*
150*     .. Local Scalars ..
151      LOGICAL            FIRST, LIEEE1, LRND
152      INTEGER            LBETA, LT
153      DOUBLE PRECISION   A, B, C, F, ONE, QTR, SAVEC, T1, T2
154*     ..
155*     .. External Functions ..
156      DOUBLE PRECISION   DLAMC3
157      EXTERNAL           DLAMC3
158*     ..
159*     .. Save statement ..
160      SAVE               FIRST, LIEEE1, LBETA, LRND, LT
161*     ..
162*     .. Data statements ..
163      DATA               FIRST / .TRUE. /
164*     ..
165*     .. Executable Statements ..
166*
167      IF( FIRST ) THEN
168         ONE = 1
169*
170*        LBETA,  LIEEE1,  LT and  LRND  are the  local values  of  BETA,
171*        IEEE1, T and RND.
172*
173*        Throughout this routine  we use the function  DLAMC3  to ensure
174*        that relevant values are  stored and not held in registers,  or
175*        are not affected by optimizers.
176*
177*        Compute  a = 2.0**m  with the  smallest positive integer m such
178*        that
179*
180*           fl( a + 1.0 ) = a.
181*
182         A = 1
183         C = 1
184*
185*+       WHILE( C.EQ.ONE )LOOP
186   10    CONTINUE
187         IF( C.EQ.ONE ) THEN
188            A = 2*A
189            C = DLAMC3( A, ONE )
190            C = DLAMC3( C, -A )
191            GO TO 10
192         END IF
193*+       END WHILE
194*
195*        Now compute  b = 2.0**m  with the smallest positive integer m
196*        such that
197*
198*           fl( a + b ) .gt. a.
199*
200         B = 1
201         C = DLAMC3( A, B )
202*
203*+       WHILE( C.EQ.A )LOOP
204   20    CONTINUE
205         IF( C.EQ.A ) THEN
206            B = 2*B
207            C = DLAMC3( A, B )
208            GO TO 20
209         END IF
210*+       END WHILE
211*
212*        Now compute the base.  a and c  are neighbouring floating point
213*        numbers  in the  interval  ( beta**t, beta**( t + 1 ) )  and so
214*        their difference is beta. Adding 0.25 to c is to ensure that it
215*        is truncated to beta and not ( beta - 1 ).
216*
217         QTR = ONE / 4
218         SAVEC = C
219         C = DLAMC3( C, -A )
220         LBETA = C + QTR
221*
222*        Now determine whether rounding or chopping occurs,  by adding a
223*        bit  less  than  beta/2  and a  bit  more  than  beta/2  to  a.
224*
225         B = LBETA
226         F = DLAMC3( B / 2, -B / 100 )
227         C = DLAMC3( F, A )
228         IF( C.EQ.A ) THEN
229            LRND = .TRUE.
230         ELSE
231            LRND = .FALSE.
232         END IF
233         F = DLAMC3( B / 2, B / 100 )
234         C = DLAMC3( F, A )
235         IF( ( LRND ) .AND. ( C.EQ.A ) )
236     $      LRND = .FALSE.
237*
238*        Try and decide whether rounding is done in the  IEEE  'round to
239*        nearest' style. B/2 is half a unit in the last place of the two
240*        numbers A and SAVEC. Furthermore, A is even, i.e. has last  bit
241*        zero, and SAVEC is odd. Thus adding B/2 to A should not  change
242*        A, but adding B/2 to SAVEC should change SAVEC.
243*
244         T1 = DLAMC3( B / 2, A )
245         T2 = DLAMC3( B / 2, SAVEC )
246         LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
247*
248*        Now find  the  mantissa, t.  It should  be the  integer part of
249*        log to the base beta of a,  however it is safer to determine  t
250*        by powering.  So we find t as the smallest positive integer for
251*        which
252*
253*           fl( beta**t + 1.0 ) = 1.0.
254*
255         LT = 0
256         A = 1
257         C = 1
258*
259*+       WHILE( C.EQ.ONE )LOOP
260   30    CONTINUE
261         IF( C.EQ.ONE ) THEN
262            LT = LT + 1
263            A = A*LBETA
264            C = DLAMC3( A, ONE )
265            C = DLAMC3( C, -A )
266            GO TO 30
267         END IF
268*+       END WHILE
269*
270      END IF
271*
272      BETA = LBETA
273      T = LT
274      RND = LRND
275      IEEE1 = LIEEE1
276      FIRST = .FALSE.
277      RETURN
278*
279*     End of DLAMC1
280*
281      END
282*
283************************************************************************
284*
285      SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
286*
287*  -- LAPACK auxiliary routine (version 3.2) --
288*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
289*     November 2006
290*
291*     .. Scalar Arguments ..
292      LOGICAL            RND
293      INTEGER            BETA, EMAX, EMIN, T
294      DOUBLE PRECISION   EPS, RMAX, RMIN
295*     ..
296*
297*  Purpose
298*  =======
299*
300*  DLAMC2 determines the machine parameters specified in its argument
301*  list.
302*
303*  Arguments
304*  =========
305*
306*  BETA    (output) INTEGER
307*          The base of the machine.
308*
309*  T       (output) INTEGER
310*          The number of ( BETA ) digits in the mantissa.
311*
312*  RND     (output) LOGICAL
313*          Specifies whether proper rounding  ( RND = .TRUE. )  or
314*          chopping  ( RND = .FALSE. )  occurs in addition. This may not
315*          be a reliable guide to the way in which the machine performs
316*          its arithmetic.
317*
318*  EPS     (output) DOUBLE PRECISION
319*          The smallest positive number such that
320*
321*             fl( 1.0 - EPS ) .LT. 1.0,
322*
323*          where fl denotes the computed value.
324*
325*  EMIN    (output) INTEGER
326*          The minimum exponent before (gradual) underflow occurs.
327*
328*  RMIN    (output) DOUBLE PRECISION
329*          The smallest normalized number for the machine, given by
330*          BASE**( EMIN - 1 ), where  BASE  is the floating point value
331*          of BETA.
332*
333*  EMAX    (output) INTEGER
334*          The maximum exponent before overflow occurs.
335*
336*  RMAX    (output) DOUBLE PRECISION
337*          The largest positive number for the machine, given by
338*          BASE**EMAX * ( 1 - EPS ), where  BASE  is the floating point
339*          value of BETA.
340*
341*  Further Details
342*  ===============
343*
344*  The computation of  EPS  is based on a routine PARANOIA by
345*  W. Kahan of the University of California at Berkeley.
346*
347* =====================================================================
348*
349*     .. Local Scalars ..
350      LOGICAL            FIRST, IEEE, IWARN, LIEEE1, LRND
351      INTEGER            GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
352     $                   NGNMIN, NGPMIN
353      DOUBLE PRECISION   A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
354     $                   SIXTH, SMALL, THIRD, TWO, ZERO
355*     ..
356*     .. External Functions ..
357      DOUBLE PRECISION   DLAMC3
358      EXTERNAL           DLAMC3
359*     ..
360*     .. External Subroutines ..
361      EXTERNAL           DLAMC1, DLAMC4, DLAMC5
362*     ..
363*     .. Intrinsic Functions ..
364      INTRINSIC          ABS, MAX, MIN
365*     ..
366*     .. Save statement ..
367      SAVE               FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,
368     $                   LRMIN, LT
369*     ..
370*     .. Data statements ..
371      DATA               FIRST / .TRUE. / , IWARN / .FALSE. /
372*     ..
373*     .. Executable Statements ..
374*
375      IF( FIRST ) THEN
376         ZERO = 0
377         ONE = 1
378         TWO = 2
379*
380*        LBETA, LT, LRND, LEPS, LEMIN and LRMIN  are the local values of
381*        BETA, T, RND, EPS, EMIN and RMIN.
382*
383*        Throughout this routine  we use the function  DLAMC3  to ensure
384*        that relevant values are stored  and not held in registers,  or
385*        are not affected by optimizers.
386*
387*        DLAMC1 returns the parameters  LBETA, LT, LRND and LIEEE1.
388*
389         CALL DLAMC1( LBETA, LT, LRND, LIEEE1 )
390*
391*        Start to find EPS.
392*
393         B = LBETA
394         A = B**( -LT )
395         LEPS = A
396*
397*        Try some tricks to see whether or not this is the correct  EPS.
398*
399         B = TWO / 3
400         HALF = ONE / 2
401         SIXTH = DLAMC3( B, -HALF )
402         THIRD = DLAMC3( SIXTH, SIXTH )
403         B = DLAMC3( THIRD, -HALF )
404         B = DLAMC3( B, SIXTH )
405         B = ABS( B )
406         IF( B.LT.LEPS )
407     $      B = LEPS
408*
409         LEPS = 1
410*
411*+       WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
412   10    CONTINUE
413         IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN
414            LEPS = B
415            C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
416            C = DLAMC3( HALF, -C )
417            B = DLAMC3( HALF, C )
418            C = DLAMC3( HALF, -B )
419            B = DLAMC3( HALF, C )
420            GO TO 10
421         END IF
422*+       END WHILE
423*
424         IF( A.LT.LEPS )
425     $      LEPS = A
426*
427*        Computation of EPS complete.
428*
429*        Now find  EMIN.  Let A = + or - 1, and + or - (1 + BASE**(-3)).
430*        Keep dividing  A by BETA until (gradual) underflow occurs. This
431*        is detected when we cannot recover the previous A.
432*
433         RBASE = ONE / LBETA
434         SMALL = ONE
435         DO 20 I = 1, 3
436            SMALL = DLAMC3( SMALL*RBASE, ZERO )
437   20    CONTINUE
438         A = DLAMC3( ONE, SMALL )
439         CALL DLAMC4( NGPMIN, ONE, LBETA )
440         CALL DLAMC4( NGNMIN, -ONE, LBETA )
441         CALL DLAMC4( GPMIN, A, LBETA )
442         CALL DLAMC4( GNMIN, -A, LBETA )
443         IEEE = .FALSE.
444*
445         IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN
446            IF( NGPMIN.EQ.GPMIN ) THEN
447               LEMIN = NGPMIN
448*            ( Non twos-complement machines, no gradual underflow;
449*              e.g.,  VAX )
450            ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
451               LEMIN = NGPMIN - 1 + LT
452               IEEE = .TRUE.
453*            ( Non twos-complement machines, with gradual underflow;
454*              e.g., IEEE standard followers )
455            ELSE
456               LEMIN = MIN( NGPMIN, GPMIN )
457*            ( A guess; no known machine )
458               IWARN = .TRUE.
459            END IF
460*
461         ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
462            IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
463               LEMIN = MAX( NGPMIN, NGNMIN )
464*            ( Twos-complement machines, no gradual underflow;
465*              e.g., CYBER 205 )
466            ELSE
467               LEMIN = MIN( NGPMIN, NGNMIN )
468*            ( A guess; no known machine )
469               IWARN = .TRUE.
470            END IF
471*
472         ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND.
473     $            ( GPMIN.EQ.GNMIN ) ) THEN
474            IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN
475               LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
476*            ( Twos-complement machines with gradual underflow;
477*              no known machine )
478            ELSE
479               LEMIN = MIN( NGPMIN, NGNMIN )
480*            ( A guess; no known machine )
481               IWARN = .TRUE.
482            END IF
483*
484         ELSE
485            LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
486*         ( A guess; no known machine )
487            IWARN = .TRUE.
488         END IF
489         FIRST = .FALSE.
490***
491* Comment out this if block if EMIN is ok
492         IF( IWARN ) THEN
493            FIRST = .TRUE.
494            WRITE( 6, FMT = 9999 )LEMIN
495         END IF
496***
497*
498*        Assume IEEE arithmetic if we found denormalised  numbers above,
499*        or if arithmetic seems to round in the  IEEE style,  determined
500*        in routine DLAMC1. A true IEEE machine should have both  things
501*        true; however, faulty machines may have one or the other.
502*
503         IEEE = IEEE .OR. LIEEE1
504*
505*        Compute  RMIN by successive division by  BETA. We could compute
506*        RMIN as BASE**( EMIN - 1 ),  but some machines underflow during
507*        this computation.
508*
509         LRMIN = 1
510         DO 30 I = 1, 1 - LEMIN
511            LRMIN = DLAMC3( LRMIN*RBASE, ZERO )
512   30    CONTINUE
513*
514*        Finally, call DLAMC5 to compute EMAX and RMAX.
515*
516         CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )
517      END IF
518*
519      BETA = LBETA
520      T = LT
521      RND = LRND
522      EPS = LEPS
523      EMIN = LEMIN
524      RMIN = LRMIN
525      EMAX = LEMAX
526      RMAX = LRMAX
527*
528      RETURN
529*
530 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
531     $      '  EMIN = ', I8, /
532     $      ' If, after inspection, the value EMIN looks',
533     $      ' acceptable please comment out ',
534     $      / ' the IF block as marked within the code of routine',
535     $      ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / )
536*
537*     End of DLAMC2
538*
539      END
540*
541************************************************************************
542*
543      DOUBLE PRECISION FUNCTION DLAMC3( A, B )
544*
545*  -- LAPACK auxiliary routine (version 3.2) --
546*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
547*     November 2006
548*
549*     .. Scalar Arguments ..
550      DOUBLE PRECISION   A, B
551*     ..
552*
553*  Purpose
554*  =======
555*
556*  DLAMC3  is intended to force  A  and  B  to be stored prior to doing
557*  the addition of  A  and  B ,  for use in situations where optimizers
558*  might hold one of these in a register.
559*
560*  Arguments
561*  =========
562*
563*  A       (input) DOUBLE PRECISION
564*  B       (input) DOUBLE PRECISION
565*          The values A and B.
566*
567* =====================================================================
568*
569*     .. Executable Statements ..
570*
571      DLAMC3 = A + B
572*
573      RETURN
574*
575*     End of DLAMC3
576*
577      END
578*
579************************************************************************
580*
581      SUBROUTINE DLAMC4( EMIN, START, BASE )
582*
583*  -- LAPACK auxiliary routine (version 3.2) --
584*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
585*     November 2006
586*
587*     .. Scalar Arguments ..
588      INTEGER            BASE, EMIN
589      DOUBLE PRECISION   START
590*     ..
591*
592*  Purpose
593*  =======
594*
595*  DLAMC4 is a service routine for DLAMC2.
596*
597*  Arguments
598*  =========
599*
600*  EMIN    (output) INTEGER
601*          The minimum exponent before (gradual) underflow, computed by
602*          setting A = START and dividing by BASE until the previous A
603*          can not be recovered.
604*
605*  START   (input) DOUBLE PRECISION
606*          The starting point for determining EMIN.
607*
608*  BASE    (input) INTEGER
609*          The base of the machine.
610*
611* =====================================================================
612*
613*     .. Local Scalars ..
614      INTEGER            I
615      DOUBLE PRECISION   A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
616*     ..
617*     .. External Functions ..
618      DOUBLE PRECISION   DLAMC3
619      EXTERNAL           DLAMC3
620*     ..
621*     .. Executable Statements ..
622*
623      A = START
624      ONE = 1
625      RBASE = ONE / BASE
626      ZERO = 0
627      EMIN = 1
628      B1 = DLAMC3( A*RBASE, ZERO )
629      C1 = A
630      C2 = A
631      D1 = A
632      D2 = A
633*+    WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
634*    $       ( D1.EQ.A ).AND.( D2.EQ.A )      )LOOP
635   10 CONTINUE
636      IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND.
637     $    ( D2.EQ.A ) ) THEN
638         EMIN = EMIN - 1
639         A = B1
640         B1 = DLAMC3( A / BASE, ZERO )
641         C1 = DLAMC3( B1*BASE, ZERO )
642         D1 = ZERO
643         DO 20 I = 1, BASE
644            D1 = D1 + B1
645   20    CONTINUE
646         B2 = DLAMC3( A*RBASE, ZERO )
647         C2 = DLAMC3( B2 / RBASE, ZERO )
648         D2 = ZERO
649         DO 30 I = 1, BASE
650            D2 = D2 + B2
651   30    CONTINUE
652         GO TO 10
653      END IF
654*+    END WHILE
655*
656      RETURN
657*
658*     End of DLAMC4
659*
660      END
661*
662************************************************************************
663*
664      SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
665*
666*  -- LAPACK auxiliary routine (version 3.2) --
667*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
668*     November 2006
669*
670*     .. Scalar Arguments ..
671      LOGICAL            IEEE
672      INTEGER            BETA, EMAX, EMIN, P
673      DOUBLE PRECISION   RMAX
674*     ..
675*
676*  Purpose
677*  =======
678*
679*  DLAMC5 attempts to compute RMAX, the largest machine floating-point
680*  number, without overflow.  It assumes that EMAX + abs(EMIN) sum
681*  approximately to a power of 2.  It will fail on machines where this
682*  assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
683*  EMAX = 28718).  It will also fail if the value supplied for EMIN is
684*  too large (i.e. too close to zero), probably with overflow.
685*
686*  Arguments
687*  =========
688*
689*  BETA    (input) INTEGER
690*          The base of floating-point arithmetic.
691*
692*  P       (input) INTEGER
693*          The number of base BETA digits in the mantissa of a
694*          floating-point value.
695*
696*  EMIN    (input) INTEGER
697*          The minimum exponent before (gradual) underflow.
698*
699*  IEEE    (input) LOGICAL
700*          A logical flag specifying whether or not the arithmetic
701*          system is thought to comply with the IEEE standard.
702*
703*  EMAX    (output) INTEGER
704*          The largest exponent before overflow
705*
706*  RMAX    (output) DOUBLE PRECISION
707*          The largest machine floating-point number.
708*
709* =====================================================================
710*
711*     .. Parameters ..
712      DOUBLE PRECISION   ZERO, ONE
713      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
714*     ..
715*     .. Local Scalars ..
716      INTEGER            EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
717      DOUBLE PRECISION   OLDY, RECBAS, Y, Z
718*     ..
719*     .. External Functions ..
720      DOUBLE PRECISION   DLAMC3
721      EXTERNAL           DLAMC3
722*     ..
723*     .. Intrinsic Functions ..
724      INTRINSIC          MOD
725*     ..
726*     .. Executable Statements ..
727*
728*     First compute LEXP and UEXP, two powers of 2 that bound
729*     abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
730*     approximately to the bound that is closest to abs(EMIN).
731*     (EMAX is the exponent of the required number RMAX).
732*
733      LEXP = 1
734      EXBITS = 1
735   10 CONTINUE
736      TRY = LEXP*2
737      IF( TRY.LE.( -EMIN ) ) THEN
738         LEXP = TRY
739         EXBITS = EXBITS + 1
740         GO TO 10
741      END IF
742      IF( LEXP.EQ.-EMIN ) THEN
743         UEXP = LEXP
744      ELSE
745         UEXP = TRY
746         EXBITS = EXBITS + 1
747      END IF
748*
749*     Now -LEXP is less than or equal to EMIN, and -UEXP is greater
750*     than or equal to EMIN. EXBITS is the number of bits needed to
751*     store the exponent.
752*
753      IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
754         EXPSUM = 2*LEXP
755      ELSE
756         EXPSUM = 2*UEXP
757      END IF
758*
759*     EXPSUM is the exponent range, approximately equal to
760*     EMAX - EMIN + 1 .
761*
762      EMAX = EXPSUM + EMIN - 1
763      NBITS = 1 + EXBITS + P
764*
765*     NBITS is the total number of bits needed to store a
766*     floating-point number.
767*
768      IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
769*
770*        Either there are an odd number of bits used to store a
771*        floating-point number, which is unlikely, or some bits are
772*        not used in the representation of numbers, which is possible,
773*        (e.g. Cray machines) or the mantissa has an implicit bit,
774*        (e.g. IEEE machines, Dec Vax machines), which is perhaps the
775*        most likely. We have to assume the last alternative.
776*        If this is true, then we need to reduce EMAX by one because
777*        there must be some way of representing zero in an implicit-bit
778*        system. On machines like Cray, we are reducing EMAX by one
779*        unnecessarily.
780*
781         EMAX = EMAX - 1
782      END IF
783*
784      IF( IEEE ) THEN
785*
786*        Assume we are on an IEEE machine which reserves one exponent
787*        for infinity and NaN.
788*
789         EMAX = EMAX - 1
790      END IF
791*
792*     Now create RMAX, the largest machine number, which should
793*     be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
794*
795*     First compute 1.0 - BETA**(-P), being careful that the
796*     result is less than 1.0 .
797*
798      RECBAS = ONE / BETA
799      Z = BETA - ONE
800      Y = ZERO
801      DO 20 I = 1, P
802         Z = Z*RECBAS
803         IF( Y.LT.ONE )
804     $      OLDY = Y
805         Y = DLAMC3( Y, Z )
806   20 CONTINUE
807      IF( Y.GE.ONE )
808     $   Y = OLDY
809*
810*     Now multiply by BETA**EMAX to get RMAX.
811*
812      DO 30 I = 1, EMAX
813         Y = DLAMC3( Y*BETA, ZERO )
814   30 CONTINUE
815*
816      RMAX = Y
817      RETURN
818*
819*     End of DLAMC5
820*
821      END
822