LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
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 *
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