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