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