LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
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 *> \ingroup auxOTHERauxiliary
64 *
65 * =====================================================================
66  REAL function slamch( cmach )
67 *
68 * -- LAPACK auxiliary routine --
69 * -- LAPACK is a software package provided by Univ. of Tennessee, --
70 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
71 *
72 * .. Scalar Arguments ..
73  CHARACTER cmach
74 * ..
75 * .. Parameters ..
76  REAL one, zero
77  parameter( one = 1.0e+0, zero = 0.0e+0 )
78 * ..
79 * .. Local Scalars ..
80  LOGICAL first, lrnd
81  INTEGER beta, imax, imin, it
82  REAL base, emax, emin, eps, prec, rmach, rmax, rmin,
83  $ rnd, sfmin, small, t
84 * ..
85 * .. External Functions ..
86  LOGICAL lsame
87  EXTERNAL lsame
88 * ..
89 * .. External Subroutines ..
90  EXTERNAL slamc2
91 * ..
92 * .. Save statement ..
93  SAVE first, eps, sfmin, base, t, rnd, emin, rmin,
94  $ emax, rmax, prec
95 * ..
96 * .. Data statements ..
97  DATA first / .true. /
98 * ..
99 * .. Executable Statements ..
100 *
101  IF( first ) THEN
102  CALL slamc2( beta, it, lrnd, eps, imin, rmin, imax, rmax )
103  base = beta
104  t = it
105  IF( lrnd ) THEN
106  rnd = one
107  eps = ( base**( 1-it ) ) / 2
108  ELSE
109  rnd = zero
110  eps = base**( 1-it )
111  END IF
112  prec = eps*base
113  emin = imin
114  emax = imax
115  sfmin = rmin
116  small = one / rmax
117  IF( small.GE.sfmin ) THEN
118 *
119 * Use SMALL plus a bit, to avoid the possibility of rounding
120 * causing overflow when computing 1/sfmin.
121 *
122  sfmin = small*( one+eps )
123  END IF
124  END IF
125 *
126  IF( lsame( cmach, 'E' ) ) THEN
127  rmach = eps
128  ELSE IF( lsame( cmach, 'S' ) ) THEN
129  rmach = sfmin
130  ELSE IF( lsame( cmach, 'B' ) ) THEN
131  rmach = base
132  ELSE IF( lsame( cmach, 'P' ) ) THEN
133  rmach = prec
134  ELSE IF( lsame( cmach, 'N' ) ) THEN
135  rmach = t
136  ELSE IF( lsame( cmach, 'R' ) ) THEN
137  rmach = rnd
138  ELSE IF( lsame( cmach, 'M' ) ) THEN
139  rmach = emin
140  ELSE IF( lsame( cmach, 'U' ) ) THEN
141  rmach = rmin
142  ELSE IF( lsame( cmach, 'L' ) ) THEN
143  rmach = emax
144  ELSE IF( lsame( cmach, 'O' ) ) THEN
145  rmach = rmax
146  END IF
147 *
148  slamch = rmach
149  first = .false.
150  RETURN
151 *
152 * End of SLAMCH
153 *
154  END
155 *
156 ************************************************************************
157 *
158 *> \brief \b SLAMC1
159 *> \details
160 *> \b Purpose:
161 *> \verbatim
162 *> SLAMC1 determines the machine parameters given by BETA, T, RND, and
163 *> IEEE1.
164 *> \endverbatim
165 *>
166 *> \param[out] BETA
167 *> \verbatim
168 *> The base of the machine.
169 *> \endverbatim
170 *>
171 *> \param[out] T
172 *> \verbatim
173 *> The number of ( BETA ) digits in the mantissa.
174 *> \endverbatim
175 *>
176 *> \param[out] RND
177 *> \verbatim
178 *> Specifies whether proper rounding ( RND = .TRUE. ) or
179 *> chopping ( RND = .FALSE. ) occurs in addition. This may not
180 *> be a reliable guide to the way in which the machine performs
181 *> its arithmetic.
182 *> \endverbatim
183 *>
184 *> \param[out] IEEE1
185 *> \verbatim
186 *> Specifies whether rounding appears to be done in the IEEE
187 *> 'round to nearest' style.
188 *> \endverbatim
189 *> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
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 slamc1( BETA, T, RND, IEEE1 )
207 *
208 * -- LAPACK auxiliary routine --
209 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
210 *
211 * .. Scalar Arguments ..
212  LOGICAL IEEE1, RND
213  INTEGER BETA, T
214 * ..
215 * =====================================================================
216 *
217 * .. Local Scalars ..
218  LOGICAL FIRST, LIEEE1, LRND
219  INTEGER LBETA, LT
220  REAL A, B, C, F, ONE, QTR, SAVEC, T1, T2
221 * ..
222 * .. External Functions ..
223  REAL SLAMC3
224  EXTERNAL slamc3
225 * ..
226 * .. Save statement ..
227  SAVE first, lieee1, lbeta, lrnd, lt
228 * ..
229 * .. Data statements ..
230  DATA first / .true. /
231 * ..
232 * .. Executable Statements ..
233 *
234  IF( first ) THEN
235  one = 1
236 *
237 * LBETA, LIEEE1, LT and LRND are the local values of BETA,
238 * IEEE1, T and RND.
239 *
240 * Throughout this routine we use the function SLAMC3 to ensure
241 * that relevant values are stored and not held in registers, or
242 * are not affected by optimizers.
243 *
244 * Compute a = 2.0**m with the smallest positive integer m such
245 * that
246 *
247 * fl( a + 1.0 ) = a.
248 *
249  a = 1
250  c = 1
251 *
252 *+ WHILE( C.EQ.ONE )LOOP
253  10 CONTINUE
254  IF( c.EQ.one ) THEN
255  a = 2*a
256  c = slamc3( a, one )
257  c = slamc3( c, -a )
258  GO TO 10
259  END IF
260 *+ END WHILE
261 *
262 * Now compute b = 2.0**m with the smallest positive integer m
263 * such that
264 *
265 * fl( a + b ) .gt. a.
266 *
267  b = 1
268  c = slamc3( a, b )
269 *
270 *+ WHILE( C.EQ.A )LOOP
271  20 CONTINUE
272  IF( c.EQ.a ) THEN
273  b = 2*b
274  c = slamc3( a, b )
275  GO TO 20
276  END IF
277 *+ END WHILE
278 *
279 * Now compute the base. a and c are neighbouring floating point
280 * numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
281 * their difference is beta. Adding 0.25 to c is to ensure that it
282 * is truncated to beta and not ( beta - 1 ).
283 *
284  qtr = one / 4
285  savec = c
286  c = slamc3( c, -a )
287  lbeta = c + qtr
288 *
289 * Now determine whether rounding or chopping occurs, by adding a
290 * bit less than beta/2 and a bit more than beta/2 to a.
291 *
292  b = lbeta
293  f = slamc3( b / 2, -b / 100 )
294  c = slamc3( f, a )
295  IF( c.EQ.a ) THEN
296  lrnd = .true.
297  ELSE
298  lrnd = .false.
299  END IF
300  f = slamc3( b / 2, b / 100 )
301  c = slamc3( f, a )
302  IF( ( lrnd ) .AND. ( c.EQ.a ) )
303  $ lrnd = .false.
304 *
305 * Try and decide whether rounding is done in the IEEE 'round to
306 * nearest' style. B/2 is half a unit in the last place of the two
307 * numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
308 * zero, and SAVEC is odd. Thus adding B/2 to A should not change
309 * A, but adding B/2 to SAVEC should change SAVEC.
310 *
311  t1 = slamc3( b / 2, a )
312  t2 = slamc3( b / 2, savec )
313  lieee1 = ( t1.EQ.a ) .AND. ( t2.GT.savec ) .AND. lrnd
314 *
315 * Now find the mantissa, t. It should be the integer part of
316 * log to the base beta of a, however it is safer to determine t
317 * by powering. So we find t as the smallest positive integer for
318 * which
319 *
320 * fl( beta**t + 1.0 ) = 1.0.
321 *
322  lt = 0
323  a = 1
324  c = 1
325 *
326 *+ WHILE( C.EQ.ONE )LOOP
327  30 CONTINUE
328  IF( c.EQ.one ) THEN
329  lt = lt + 1
330  a = a*lbeta
331  c = slamc3( a, one )
332  c = slamc3( c, -a )
333  GO TO 30
334  END IF
335 *+ END WHILE
336 *
337  END IF
338 *
339  beta = lbeta
340  t = lt
341  rnd = lrnd
342  ieee1 = lieee1
343  first = .false.
344  RETURN
345 *
346 * End of SLAMC1
347 *
348  END
349 *
350 ************************************************************************
351 *
352 *> \brief \b SLAMC2
353 *> \details
354 *> \b Purpose:
355 *> \verbatim
356 *> SLAMC2 determines the machine parameters specified in its argument
357 *> list.
358 *> \endverbatim
359 *> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
360 *> \ingroup auxOTHERauxiliary
361 *>
362 *> \param[out] BETA
363 *> \verbatim
364 *> The base of the machine.
365 *> \endverbatim
366 *>
367 *> \param[out] T
368 *> \verbatim
369 *> The number of ( BETA ) digits in the mantissa.
370 *> \endverbatim
371 *>
372 *> \param[out] RND
373 *> \verbatim
374 *> Specifies whether proper rounding ( RND = .TRUE. ) or
375 *> chopping ( RND = .FALSE. ) occurs in addition. This may not
376 *> be a reliable guide to the way in which the machine performs
377 *> its arithmetic.
378 *> \endverbatim
379 *>
380 *> \param[out] EPS
381 *> \verbatim
382 *> The smallest positive number such that
383 *> fl( 1.0 - EPS ) .LT. 1.0,
384 *> where fl denotes the computed value.
385 *> \endverbatim
386 *>
387 *> \param[out] EMIN
388 *> \verbatim
389 *> The minimum exponent before (gradual) underflow occurs.
390 *> \endverbatim
391 *>
392 *> \param[out] RMIN
393 *> \verbatim
394 *> The smallest normalized number for the machine, given by
395 *> BASE**( EMIN - 1 ), where BASE is the floating point value
396 *> of BETA.
397 *> \endverbatim
398 *>
399 *> \param[out] EMAX
400 *> \verbatim
401 *> The maximum exponent before overflow occurs.
402 *> \endverbatim
403 *>
404 *> \param[out] RMAX
405 *> \verbatim
406 *> The largest positive number for the machine, given by
407 *> BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
408 *> value of BETA.
409 *> \endverbatim
410 *>
411 *> \details \b Further \b Details
412 *> \verbatim
413 *>
414 *> The computation of EPS is based on a routine PARANOIA by
415 *> W. Kahan of the University of California at Berkeley.
416 *> \endverbatim
417  SUBROUTINE slamc2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
418 *
419 * -- LAPACK auxiliary routine --
420 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
421 *
422 * .. Scalar Arguments ..
423  LOGICAL RND
424  INTEGER BETA, EMAX, EMIN, T
425  REAL EPS, RMAX, RMIN
426 * ..
427 * =====================================================================
428 *
429 * .. Local Scalars ..
430  LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND
431  INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
432  $ NGNMIN, NGPMIN
433  REAL A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
434  $ SIXTH, SMALL, THIRD, TWO, ZERO
435 * ..
436 * .. External Functions ..
437  REAL SLAMC3
438  EXTERNAL slamc3
439 * ..
440 * .. External Subroutines ..
441  EXTERNAL slamc1, slamc4, slamc5
442 * ..
443 * .. Intrinsic Functions ..
444  INTRINSIC abs, max, min
445 * ..
446 * .. Save statement ..
447  SAVE first, iwarn, lbeta, lemax, lemin, leps, lrmax,
448  $ lrmin, lt
449 * ..
450 * .. Data statements ..
451  DATA first / .true. / , iwarn / .false. /
452 * ..
453 * .. Executable Statements ..
454 *
455  IF( first ) THEN
456  zero = 0
457  one = 1
458  two = 2
459 *
460 * LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
461 * BETA, T, RND, EPS, EMIN and RMIN.
462 *
463 * Throughout this routine we use the function SLAMC3 to ensure
464 * that relevant values are stored and not held in registers, or
465 * are not affected by optimizers.
466 *
467 * SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
468 *
469  CALL slamc1( lbeta, lt, lrnd, lieee1 )
470 *
471 * Start to find EPS.
472 *
473  b = lbeta
474  a = b**( -lt )
475  leps = a
476 *
477 * Try some tricks to see whether or not this is the correct EPS.
478 *
479  b = two / 3
480  half = one / 2
481  sixth = slamc3( b, -half )
482  third = slamc3( sixth, sixth )
483  b = slamc3( third, -half )
484  b = slamc3( b, sixth )
485  b = abs( b )
486  IF( b.LT.leps )
487  $ b = leps
488 *
489  leps = 1
490 *
491 *+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
492  10 CONTINUE
493  IF( ( leps.GT.b ) .AND. ( b.GT.zero ) ) THEN
494  leps = b
495  c = slamc3( half*leps, ( two**5 )*( leps**2 ) )
496  c = slamc3( half, -c )
497  b = slamc3( half, c )
498  c = slamc3( half, -b )
499  b = slamc3( half, c )
500  GO TO 10
501  END IF
502 *+ END WHILE
503 *
504  IF( a.LT.leps )
505  $ leps = a
506 *
507 * Computation of EPS complete.
508 *
509 * Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
510 * Keep dividing A by BETA until (gradual) underflow occurs. This
511 * is detected when we cannot recover the previous A.
512 *
513  rbase = one / lbeta
514  small = one
515  DO 20 i = 1, 3
516  small = slamc3( small*rbase, zero )
517  20 CONTINUE
518  a = slamc3( one, small )
519  CALL slamc4( ngpmin, one, lbeta )
520  CALL slamc4( ngnmin, -one, lbeta )
521  CALL slamc4( gpmin, a, lbeta )
522  CALL slamc4( gnmin, -a, lbeta )
523  ieee = .false.
524 *
525  IF( ( ngpmin.EQ.ngnmin ) .AND. ( gpmin.EQ.gnmin ) ) THEN
526  IF( ngpmin.EQ.gpmin ) THEN
527  lemin = ngpmin
528 * ( Non twos-complement machines, no gradual underflow;
529 * e.g., VAX )
530  ELSE IF( ( gpmin-ngpmin ).EQ.3 ) THEN
531  lemin = ngpmin - 1 + lt
532  ieee = .true.
533 * ( Non twos-complement machines, with gradual underflow;
534 * e.g., IEEE standard followers )
535  ELSE
536  lemin = min( ngpmin, gpmin )
537 * ( A guess; no known machine )
538  iwarn = .true.
539  END IF
540 *
541  ELSE IF( ( ngpmin.EQ.gpmin ) .AND. ( ngnmin.EQ.gnmin ) ) THEN
542  IF( abs( ngpmin-ngnmin ).EQ.1 ) THEN
543  lemin = max( ngpmin, ngnmin )
544 * ( Twos-complement machines, no gradual underflow;
545 * e.g., CYBER 205 )
546  ELSE
547  lemin = min( ngpmin, ngnmin )
548 * ( A guess; no known machine )
549  iwarn = .true.
550  END IF
551 *
552  ELSE IF( ( abs( ngpmin-ngnmin ).EQ.1 ) .AND.
553  $ ( gpmin.EQ.gnmin ) ) THEN
554  IF( ( gpmin-min( ngpmin, ngnmin ) ).EQ.3 ) THEN
555  lemin = max( ngpmin, ngnmin ) - 1 + lt
556 * ( Twos-complement machines with gradual underflow;
557 * no known machine )
558  ELSE
559  lemin = min( ngpmin, ngnmin )
560 * ( A guess; no known machine )
561  iwarn = .true.
562  END IF
563 *
564  ELSE
565  lemin = min( ngpmin, ngnmin, gpmin, gnmin )
566 * ( A guess; no known machine )
567  iwarn = .true.
568  END IF
569  first = .false.
570 ***
571 * Comment out this if block if EMIN is ok
572  IF( iwarn ) THEN
573  first = .true.
574  WRITE( 6, fmt = 9999 )lemin
575  END IF
576 ***
577 *
578 * Assume IEEE arithmetic if we found denormalised numbers above,
579 * or if arithmetic seems to round in the IEEE style, determined
580 * in routine SLAMC1. A true IEEE machine should have both things
581 * true; however, faulty machines may have one or the other.
582 *
583  ieee = ieee .OR. lieee1
584 *
585 * Compute RMIN by successive division by BETA. We could compute
586 * RMIN as BASE**( EMIN - 1 ), but some machines underflow during
587 * this computation.
588 *
589  lrmin = 1
590  DO 30 i = 1, 1 - lemin
591  lrmin = slamc3( lrmin*rbase, zero )
592  30 CONTINUE
593 *
594 * Finally, call SLAMC5 to compute EMAX and RMAX.
595 *
596  CALL slamc5( lbeta, lt, lemin, ieee, lemax, lrmax )
597  END IF
598 *
599  beta = lbeta
600  t = lt
601  rnd = lrnd
602  eps = leps
603  emin = lemin
604  rmin = lrmin
605  emax = lemax
606  rmax = lrmax
607 *
608  RETURN
609 *
610  9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
611  $ ' EMIN = ', i8, /
612  $ ' If, after inspection, the value EMIN looks',
613  $ ' acceptable please comment out ',
614  $ / ' the IF block as marked within the code of routine',
615  $ ' SLAMC2,', / ' otherwise supply EMIN explicitly.', / )
616 *
617 * End of SLAMC2
618 *
619  END
620 *
621 ************************************************************************
622 *
623 *> \brief \b SLAMC3
624 *> \details
625 *> \b Purpose:
626 *> \verbatim
627 *> SLAMC3 is intended to force A and B to be stored prior to doing
628 *> the addition of A and B , for use in situations where optimizers
629 *> might hold one of these in a register.
630 *> \endverbatim
631 *>
632 *> \param[in] A
633 *>
634 *> \param[in] B
635 *> \verbatim
636 *> The values A and B.
637 *> \endverbatim
638 
639  REAL function slamc3( a, b )
640 *
641 * -- LAPACK auxiliary routine --
642 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
643 *
644 * .. Scalar Arguments ..
645  REAL a, b
646 * ..
647 * =====================================================================
648 *
649 * .. Executable Statements ..
650 *
651  slamc3 = a + b
652 *
653  RETURN
654 *
655 * End of SLAMC3
656 *
657  END
658 *
659 ************************************************************************
660 *
661 *> \brief \b SLAMC4
662 *> \details
663 *> \b Purpose:
664 *> \verbatim
665 *> SLAMC4 is a service routine for SLAMC2.
666 *> \endverbatim
667 *>
668 *> \param[out] EMIN
669 *> \verbatim
670 *> The minimum exponent before (gradual) underflow, computed by
671 *> setting A = START and dividing by BASE until the previous A
672 *> can not be recovered.
673 *> \endverbatim
674 *>
675 *> \param[in] START
676 *> \verbatim
677 *> The starting point for determining EMIN.
678 *> \endverbatim
679 *>
680 *> \param[in] BASE
681 *> \verbatim
682 *> The base of the machine.
683 *> \endverbatim
684 *>
685  SUBROUTINE slamc4( EMIN, START, BASE )
686 *
687 * -- LAPACK auxiliary routine --
688 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
689 *
690 * .. Scalar Arguments ..
691  INTEGER BASE
692  INTEGER EMIN
693  REAL START
694 * ..
695 * =====================================================================
696 *
697 * .. Local Scalars ..
698  INTEGER I
699  REAL A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
700 * ..
701 * .. External Functions ..
702  REAL SLAMC3
703  EXTERNAL slamc3
704 * ..
705 * .. Executable Statements ..
706 *
707  a = start
708  one = 1
709  rbase = one / base
710  zero = 0
711  emin = 1
712  b1 = slamc3( 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 = slamc3( a / base, zero )
725  c1 = slamc3( b1*base, zero )
726  d1 = zero
727  DO 20 i = 1, base
728  d1 = d1 + b1
729  20 CONTINUE
730  b2 = slamc3( a*rbase, zero )
731  c2 = slamc3( 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 SLAMC4
743 *
744  END
745 *
746 ************************************************************************
747 *
748 *> \brief \b SLAMC5
749 *> \details
750 *> \b Purpose:
751 *> \verbatim
752 *> SLAMC5 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 slamc5( 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  REAL RMAX
801 * ..
802 * =====================================================================
803 *
804 * .. Parameters ..
805  REAL ZERO, ONE
806  parameter( zero = 0.0e0, one = 1.0e0 )
807 * ..
808 * .. Local Scalars ..
809  INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
810  REAL OLDY, RECBAS, Y, Z
811 * ..
812 * .. External Functions ..
813  REAL SLAMC3
814  EXTERNAL slamc3
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 = slamc3( 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 = slamc3( y*beta, zero )
907  30 CONTINUE
908 *
909  rmax = y
910  RETURN
911 *
912 * End of SLAMC5
913 *
914  END
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
real function slamc3(A, B)
SLAMC3
Definition: slamch.f:169
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
subroutine slamc1(BETA, T, RND, IEEE1)
SLAMC1
Definition: slamchf77.f:207
subroutine slamc5(BETA, P, EMIN, IEEE, EMAX, RMAX)
SLAMC5
Definition: slamchf77.f:793
subroutine slamc2(BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX)
SLAMC2
Definition: slamchf77.f:418
subroutine slamc4(EMIN, START, BASE)
SLAMC4
Definition: slamchf77.f:686