LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
slamchf77.f
Go to the documentation of this file.
1 *> \brief \b SLAMCHF77 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 * REAL FUNCTION SLAMCH( CMACH )
12 *
13 * .. Scalar Arguments ..
14 * CHARACTER CMACH
15 * ..
16 *
17 *
18 *> \par Purpose:
19 * =============
20 *>
21 *> \verbatim
22 *>
23 *> SLAMCH determines single precision machine parameters.
24 *> \endverbatim
25 *
26 * Arguments:
27 * ==========
28 *
29 *> \param[in] CMACH
30 *> \verbatim
31 *> Specifies the value to be returned by SLAMCH:
32 *> = 'E' or 'e', SLAMCH := eps
33 *> = 'S' or 's , SLAMCH := sfmin
34 *> = 'B' or 'b', SLAMCH := base
35 *> = 'P' or 'p', SLAMCH := eps*base
36 *> = 'N' or 'n', SLAMCH := t
37 *> = 'R' or 'r', SLAMCH := rnd
38 *> = 'M' or 'm', SLAMCH := emin
39 *> = 'U' or 'u', SLAMCH := rmin
40 *> = 'L' or 'l', SLAMCH := emax
41 *> = 'O' or 'o', SLAMCH := 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  REAL FUNCTION slamch( CMACH )
69 *
70 * -- LAPACK auxiliary routine (version 3.4.1) --
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  REAL one, zero
80  parameter( one = 1.0e+0, zero = 0.0e+0 )
81 * ..
82 * .. Local Scalars ..
83  LOGICAL first, lrnd
84  INTEGER beta, imax, imin, it
85  REAL 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 slamc2
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 slamc2( 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  slamch = rmach
152  first = .false.
153  return
154 *
155 * End of SLAMCH
156 *
157  END
158 *
159 ************************************************************************
160 *
161 *> \brief \b SLAMC1
162 *> \details
163 *> \b Purpose:
164 *> \verbatim
165 *> SLAMC1 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 slamc1( BETA, T, RND, IEEE1 )
211 *
212 * -- LAPACK auxiliary routine (version 3.4.1) --
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  REAL a, b, c, f, one, qtr, savec, t1, t2
226 * ..
227 * .. External Functions ..
228  REAL slamc3
229  EXTERNAL slamc3
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 SLAMC3 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 = slamc3( a, one )
262  c = slamc3( 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 = slamc3( a, b )
274 *
275 *+ WHILE( C.EQ.A )LOOP
276  20 continue
277  IF( c.EQ.a ) THEN
278  b = 2*b
279  c = slamc3( 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 = slamc3( 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 = slamc3( b / 2, -b / 100 )
299  c = slamc3( f, a )
300  IF( c.EQ.a ) THEN
301  lrnd = .true.
302  ELSE
303  lrnd = .false.
304  END IF
305  f = slamc3( b / 2, b / 100 )
306  c = slamc3( 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 = slamc3( b / 2, a )
317  t2 = slamc3( 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 = slamc3( a, one )
337  c = slamc3( 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 SLAMC1
352 *
353  END
354 *
355 ************************************************************************
356 *
357 *> \brief \b SLAMC2
358 *> \details
359 *> \b Purpose:
360 *> \verbatim
361 *> SLAMC2 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 slamc2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
424 *
425 * -- LAPACK auxiliary routine (version 3.4.1) --
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  REAL 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  REAL a, b, c, half, leps, lrmax, lrmin, one, rbase,
441  $ sixth, small, third, two, zero
442 * ..
443 * .. External Functions ..
444  REAL slamc3
445  EXTERNAL slamc3
446 * ..
447 * .. External Subroutines ..
448  EXTERNAL slamc1, slamc4, slamc5
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 SLAMC3 to ensure
471 * that relevant values are stored and not held in registers, or
472 * are not affected by optimizers.
473 *
474 * SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
475 *
476  CALL slamc1( 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 = slamc3( b, -half )
489  third = slamc3( sixth, sixth )
490  b = slamc3( third, -half )
491  b = slamc3( 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 = slamc3( half*leps, ( two**5 )*( leps**2 ) )
503  c = slamc3( half, -c )
504  b = slamc3( half, c )
505  c = slamc3( half, -b )
506  b = slamc3( 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 = slamc3( small*rbase, zero )
524  20 continue
525  a = slamc3( one, small )
526  CALL slamc4( ngpmin, one, lbeta )
527  CALL slamc4( ngnmin, -one, lbeta )
528  CALL slamc4( gpmin, a, lbeta )
529  CALL slamc4( 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 SLAMC1. 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 = slamc3( lrmin*rbase, zero )
599  30 continue
600 *
601 * Finally, call SLAMC5 to compute EMAX and RMAX.
602 *
603  CALL slamc5( 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  $ ' SLAMC2,', / ' otherwise supply EMIN explicitly.', / )
623 *
624 * End of SLAMC2
625 *
626  END
627 *
628 ************************************************************************
629 *
630 *> \brief \b SLAMC3
631 *> \details
632 *> \b Purpose:
633 *> \verbatim
634 *> SLAMC3 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  REAL FUNCTION slamc3( A, B )
647 *
648 * -- LAPACK auxiliary routine (version 3.4.1) --
649 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
650 * November 2010
651 *
652 * .. Scalar Arguments ..
653  REAL a, b
654 * ..
655 * =====================================================================
656 *
657 * .. Executable Statements ..
658 *
659  slamc3 = a + b
660 *
661  return
662 *
663 * End of SLAMC3
664 *
665  END
666 *
667 ************************************************************************
668 *
669 *> \brief \b SLAMC4
670 *> \details
671 *> \b Purpose:
672 *> \verbatim
673 *> SLAMC4 is a service routine for SLAMC2.
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 slamc4( EMIN, START, BASE )
694 *
695 * -- LAPACK auxiliary routine (version 3.4.1) --
696 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
697 * November 2010
698 *
699 * .. Scalar Arguments ..
700  INTEGER base
701  INTEGER emin
702  REAL start
703 * ..
704 * =====================================================================
705 *
706 * .. Local Scalars ..
707  INTEGER i
708  REAL a, b1, b2, c1, c2, d1, d2, one, rbase, zero
709 * ..
710 * .. External Functions ..
711  REAL slamc3
712  EXTERNAL slamc3
713 * ..
714 * .. Executable Statements ..
715 *
716  a = start
717  one = 1
718  rbase = one / base
719  zero = 0
720  emin = 1
721  b1 = slamc3( a*rbase, zero )
722  c1 = a
723  c2 = a
724  d1 = a
725  d2 = a
726 *+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
727 * $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP
728  10 continue
729  IF( ( c1.EQ.a ) .AND. ( c2.EQ.a ) .AND. ( d1.EQ.a ) .AND.
730  $ ( d2.EQ.a ) ) THEN
731  emin = emin - 1
732  a = b1
733  b1 = slamc3( a / base, zero )
734  c1 = slamc3( b1*base, zero )
735  d1 = zero
736  DO 20 i = 1, base
737  d1 = d1 + b1
738  20 continue
739  b2 = slamc3( a*rbase, zero )
740  c2 = slamc3( b2 / rbase, zero )
741  d2 = zero
742  DO 30 i = 1, base
743  d2 = d2 + b2
744  30 continue
745  go to 10
746  END IF
747 *+ END WHILE
748 *
749  return
750 *
751 * End of SLAMC4
752 *
753  END
754 *
755 ************************************************************************
756 *
757 *> \brief \b SLAMC5
758 *> \details
759 *> \b Purpose:
760 *> \verbatim
761 *> SLAMC5 attempts to compute RMAX, the largest machine floating-point
762 *> number, without overflow. It assumes that EMAX + abs(EMIN) sum
763 *> approximately to a power of 2. It will fail on machines where this
764 *> assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
765 *> EMAX = 28718). It will also fail if the value supplied for EMIN is
766 *> too large (i.e. too close to zero), probably with overflow.
767 *> \endverbatim
768 *>
769 *> \param[in] BETA
770 *> \verbatim
771 *> The base of floating-point arithmetic.
772 *> \endverbatim
773 *>
774 *> \param[in] P
775 *> \verbatim
776 *> The number of base BETA digits in the mantissa of a
777 *> floating-point value.
778 *> \endverbatim
779 *>
780 *> \param[in] EMIN
781 *> \verbatim
782 *> The minimum exponent before (gradual) underflow.
783 *> \endverbatim
784 *>
785 *> \param[in] IEEE
786 *> \verbatim
787 *> A logical flag specifying whether or not the arithmetic
788 *> system is thought to comply with the IEEE standard.
789 *> \endverbatim
790 *>
791 *> \param[out] EMAX
792 *> \verbatim
793 *> The largest exponent before overflow
794 *> \endverbatim
795 *>
796 *> \param[out] RMAX
797 *> \verbatim
798 *> The largest machine floating-point number.
799 *> \endverbatim
800 *>
801  SUBROUTINE slamc5( BETA, P, EMIN, IEEE, EMAX, RMAX )
802 *
803 * -- LAPACK auxiliary routine (version 3.4.1) --
804 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
805 * November 2010
806 *
807 * .. Scalar Arguments ..
808  LOGICAL ieee
809  INTEGER beta, emax, emin, p
810  REAL rmax
811 * ..
812 * =====================================================================
813 *
814 * .. Parameters ..
815  REAL zero, one
816  parameter( zero = 0.0e0, one = 1.0e0 )
817 * ..
818 * .. Local Scalars ..
819  INTEGER exbits, expsum, i, lexp, nbits, try, uexp
820  REAL oldy, recbas, y, z
821 * ..
822 * .. External Functions ..
823  REAL slamc3
824  EXTERNAL slamc3
825 * ..
826 * .. Intrinsic Functions ..
827  INTRINSIC mod
828 * ..
829 * .. Executable Statements ..
830 *
831 * First compute LEXP and UEXP, two powers of 2 that bound
832 * abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
833 * approximately to the bound that is closest to abs(EMIN).
834 * (EMAX is the exponent of the required number RMAX).
835 *
836  lexp = 1
837  exbits = 1
838  10 continue
839  try = lexp*2
840  IF( try.LE.( -emin ) ) THEN
841  lexp = try
842  exbits = exbits + 1
843  go to 10
844  END IF
845  IF( lexp.EQ.-emin ) THEN
846  uexp = lexp
847  ELSE
848  uexp = try
849  exbits = exbits + 1
850  END IF
851 *
852 * Now -LEXP is less than or equal to EMIN, and -UEXP is greater
853 * than or equal to EMIN. EXBITS is the number of bits needed to
854 * store the exponent.
855 *
856  IF( ( uexp+emin ).GT.( -lexp-emin ) ) THEN
857  expsum = 2*lexp
858  ELSE
859  expsum = 2*uexp
860  END IF
861 *
862 * EXPSUM is the exponent range, approximately equal to
863 * EMAX - EMIN + 1 .
864 *
865  emax = expsum + emin - 1
866  nbits = 1 + exbits + p
867 *
868 * NBITS is the total number of bits needed to store a
869 * floating-point number.
870 *
871  IF( ( mod( nbits, 2 ).EQ.1 ) .AND. ( beta.EQ.2 ) ) THEN
872 *
873 * Either there are an odd number of bits used to store a
874 * floating-point number, which is unlikely, or some bits are
875 * not used in the representation of numbers, which is possible,
876 * (e.g. Cray machines) or the mantissa has an implicit bit,
877 * (e.g. IEEE machines, Dec Vax machines), which is perhaps the
878 * most likely. We have to assume the last alternative.
879 * If this is true, then we need to reduce EMAX by one because
880 * there must be some way of representing zero in an implicit-bit
881 * system. On machines like Cray, we are reducing EMAX by one
882 * unnecessarily.
883 *
884  emax = emax - 1
885  END IF
886 *
887  IF( ieee ) THEN
888 *
889 * Assume we are on an IEEE machine which reserves one exponent
890 * for infinity and NaN.
891 *
892  emax = emax - 1
893  END IF
894 *
895 * Now create RMAX, the largest machine number, which should
896 * be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
897 *
898 * First compute 1.0 - BETA**(-P), being careful that the
899 * result is less than 1.0 .
900 *
901  recbas = one / beta
902  z = beta - one
903  y = zero
904  DO 20 i = 1, p
905  z = z*recbas
906  IF( y.LT.one )
907  $ oldy = y
908  y = slamc3( y, z )
909  20 continue
910  IF( y.GE.one )
911  $ y = oldy
912 *
913 * Now multiply by BETA**EMAX to get RMAX.
914 *
915  DO 30 i = 1, emax
916  y = slamc3( y*beta, zero )
917  30 continue
918 *
919  rmax = y
920  return
921 *
922 * End of SLAMC5
923 *
924  END