ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
tools.f
Go to the documentation of this file.
1 * ================================================================
2 * This file contains the following LAPACK routines, for use by the
3 * BLACS tester: LSAME, SLAMCH, DLAMCH, DLARND, ZLARND, DLARAN,
4 * and ZLARAN. If you have ScaLAPACK or LAPACK, all of these files
5 * are present in your library, and you may discard this file and
6 * point to the appropriate archive instead.
7 * ================================================================
8 
9  DOUBLE PRECISION FUNCTION dlamch( CMACH )
10 *
11 * -- LAPACK auxiliary routine (version 2.0) --
12 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
13 * Courant Institute, Argonne National Lab, and Rice University
14 * October 31, 1992
15 *
16 * .. Scalar Arguments ..
17  CHARACTER cmach
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * DLAMCH determines double precision machine parameters.
24 *
25 * Arguments
26 * =========
27 *
28 * CMACH (input) CHARACTER*1
29 * Specifies the value to be returned by DLAMCH:
30 * = 'E' or 'e', DLAMCH := eps
31 * = 'S' or 's , DLAMCH := sfmin
32 * = 'B' or 'b', DLAMCH := base
33 * = 'P' or 'p', DLAMCH := eps*base
34 * = 'N' or 'n', DLAMCH := t
35 * = 'R' or 'r', DLAMCH := rnd
36 * = 'M' or 'm', DLAMCH := emin
37 * = 'U' or 'u', DLAMCH := rmin
38 * = 'L' or 'l', DLAMCH := emax
39 * = 'O' or 'o', DLAMCH := rmax
40 *
41 * where
42 *
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 *
54 * =====================================================================
55 *
56 * .. Parameters ..
57  DOUBLE PRECISION one, zero
58  parameter( one = 1.0d+0, zero = 0.0d+0 )
59 * ..
60 * .. Local Scalars ..
61  LOGICAL first, lrnd
62  INTEGER beta, imax, imin, it
63  DOUBLE PRECISION base, emax, emin, eps, prec, rmach, rmax, rmin,
64  $ rnd, sfmin, small, t
65 * ..
66 * .. External Functions ..
67  LOGICAL lsame
68  EXTERNAL lsame
69 * ..
70 * .. External Subroutines ..
71  EXTERNAL dlamc2
72 * ..
73 * .. Save statement ..
74  SAVE first, eps, sfmin, base, t, rnd, emin, rmin,
75  $ emax, rmax, prec
76 * ..
77 * .. Data statements ..
78  DATA first / .true. /
79 * ..
80 * .. Executable Statements ..
81 *
82  IF( first ) THEN
83  first = .false.
84  CALL dlamc2( beta, it, lrnd, eps, imin, rmin, imax, rmax )
85  base = beta
86  t = it
87  IF( lrnd ) THEN
88  rnd = one
89  eps = ( base**( 1-it ) ) / 2
90  ELSE
91  rnd = zero
92  eps = base**( 1-it )
93  END IF
94  prec = eps*base
95  emin = imin
96  emax = imax
97  sfmin = rmin
98  small = one / rmax
99  IF( small.GE.sfmin ) THEN
100 *
101 * Use SMALL plus a bit, to avoid the possibility of rounding
102 * causing overflow when computing 1/sfmin.
103 *
104  sfmin = small*( one+eps )
105  END IF
106  END IF
107 *
108  IF( lsame( cmach, 'E' ) ) THEN
109  rmach = eps
110  ELSE IF( lsame( cmach, 'S' ) ) THEN
111  rmach = sfmin
112  ELSE IF( lsame( cmach, 'B' ) ) THEN
113  rmach = base
114  ELSE IF( lsame( cmach, 'P' ) ) THEN
115  rmach = prec
116  ELSE IF( lsame( cmach, 'N' ) ) THEN
117  rmach = t
118  ELSE IF( lsame( cmach, 'R' ) ) THEN
119  rmach = rnd
120  ELSE IF( lsame( cmach, 'M' ) ) THEN
121  rmach = emin
122  ELSE IF( lsame( cmach, 'U' ) ) THEN
123  rmach = rmin
124  ELSE IF( lsame( cmach, 'L' ) ) THEN
125  rmach = emax
126  ELSE IF( lsame( cmach, 'O' ) ) THEN
127  rmach = rmax
128  END IF
129 *
130  dlamch = rmach
131  RETURN
132 *
133 * End of DLAMCH
134 *
135  END
136 *
137 ************************************************************************
138 *
139  SUBROUTINE dlamc1( BETA, T, RND, IEEE1 )
140 *
141 * -- LAPACK auxiliary routine (version 2.0) --
142 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
143 * Courant Institute, Argonne National Lab, and Rice University
144 * October 31, 1992
145 *
146 * .. Scalar Arguments ..
147  LOGICAL ieee1, rnd
148  INTEGER beta, t
149 * ..
150 *
151 * Purpose
152 * =======
153 *
154 * DLAMC1 determines the machine parameters given by BETA, T, RND, and
155 * IEEE1.
156 *
157 * Arguments
158 * =========
159 *
160 * BETA (output) INTEGER
161 * The base of the machine.
162 *
163 * T (output) INTEGER
164 * The number of ( BETA ) digits in the mantissa.
165 *
166 * RND (output) LOGICAL
167 * Specifies whether proper rounding ( RND = .TRUE. ) or
168 * chopping ( RND = .FALSE. ) occurs in addition. This may not
169 * be a reliable guide to the way in which the machine performs
170 * its arithmetic.
171 *
172 * IEEE1 (output) LOGICAL
173 * Specifies whether rounding appears to be done in the IEEE
174 * 'round to nearest' style.
175 *
176 * Further Details
177 * ===============
178 *
179 * The routine is based on the routine ENVRON by Malcolm and
180 * incorporates suggestions by Gentleman and Marovich. See
181 *
182 * Malcolm M. A. (1972) Algorithms to reveal properties of
183 * floating-point arithmetic. Comms. of the ACM, 15, 949-951.
184 *
185 * Gentleman W. M. and Marovich S. B. (1974) More on algorithms
186 * that reveal properties of floating point arithmetic units.
187 * Comms. of the ACM, 17, 276-277.
188 *
189 * =====================================================================
190 *
191 * .. Local Scalars ..
192  LOGICAL first, lieee1, lrnd
193  INTEGER lbeta, lt
194  DOUBLE PRECISION a, b, c, f, one, qtr, savec, t1, t2
195 * ..
196 * .. External Functions ..
197  DOUBLE PRECISION dlamc3
198  EXTERNAL dlamc3
199 * ..
200 * .. Save statement ..
201  SAVE first, lieee1, lbeta, lrnd, lt
202 * ..
203 * .. Data statements ..
204  DATA first / .true. /
205 * ..
206 * .. Executable Statements ..
207 *
208  IF( first ) THEN
209  first = .false.
210  one = 1
211 *
212 * LBETA, LIEEE1, LT and LRND are the local values of BETA,
213 * IEEE1, T and RND.
214 *
215 * Throughout this routine we use the function DLAMC3 to ensure
216 * that relevant values are stored and not held in registers, or
217 * are not affected by optimizers.
218 *
219 * Compute a = 2.0**m with the smallest positive integer m such
220 * that
221 *
222 * fl( a + 1.0 ) = a.
223 *
224  a = 1
225  c = 1
226 *
227 *+ WHILE( C.EQ.ONE )LOOP
228  10 CONTINUE
229  IF( c.EQ.one ) THEN
230  a = 2*a
231  c = dlamc3( a, one )
232  c = dlamc3( c, -a )
233  GO TO 10
234  END IF
235 *+ END WHILE
236 *
237 * Now compute b = 2.0**m with the smallest positive integer m
238 * such that
239 *
240 * fl( a + b ) .gt. a.
241 *
242  b = 1
243  c = dlamc3( a, b )
244 *
245 *+ WHILE( C.EQ.A )LOOP
246  20 CONTINUE
247  IF( c.EQ.a ) THEN
248  b = 2*b
249  c = dlamc3( a, b )
250  GO TO 20
251  END IF
252 *+ END WHILE
253 *
254 * Now compute the base. a and c are neighbouring floating point
255 * numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
256 * their difference is beta. Adding 0.25 to c is to ensure that it
257 * is truncated to beta and not ( beta - 1 ).
258 *
259  qtr = one / 4
260  savec = c
261  c = dlamc3( c, -a )
262  lbeta = c + qtr
263 *
264 * Now determine whether rounding or chopping occurs, by adding a
265 * bit less than beta/2 and a bit more than beta/2 to a.
266 *
267  b = lbeta
268  f = dlamc3( b / 2, -b / 100 )
269  c = dlamc3( f, a )
270  IF( c.EQ.a ) THEN
271  lrnd = .true.
272  ELSE
273  lrnd = .false.
274  END IF
275  f = dlamc3( b / 2, b / 100 )
276  c = dlamc3( f, a )
277  IF( ( lrnd ) .AND. ( c.EQ.a ) )
278  $ lrnd = .false.
279 *
280 * Try and decide whether rounding is done in the IEEE 'round to
281 * nearest' style. B/2 is half a unit in the last place of the two
282 * numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
283 * zero, and SAVEC is odd. Thus adding B/2 to A should not change
284 * A, but adding B/2 to SAVEC should change SAVEC.
285 *
286  t1 = dlamc3( b / 2, a )
287  t2 = dlamc3( b / 2, savec )
288  lieee1 = ( t1.EQ.a ) .AND. ( t2.GT.savec ) .AND. lrnd
289 *
290 * Now find the mantissa, t. It should be the integer part of
291 * log to the base beta of a, however it is safer to determine t
292 * by powering. So we find t as the smallest positive integer for
293 * which
294 *
295 * fl( beta**t + 1.0 ) = 1.0.
296 *
297  lt = 0
298  a = 1
299  c = 1
300 *
301 *+ WHILE( C.EQ.ONE )LOOP
302  30 CONTINUE
303  IF( c.EQ.one ) THEN
304  lt = lt + 1
305  a = a*lbeta
306  c = dlamc3( a, one )
307  c = dlamc3( c, -a )
308  GO TO 30
309  END IF
310 *+ END WHILE
311 *
312  END IF
313 *
314  beta = lbeta
315  t = lt
316  rnd = lrnd
317  ieee1 = lieee1
318  RETURN
319 *
320 * End of DLAMC1
321 *
322  END
323 *
324 ************************************************************************
325 *
326  SUBROUTINE dlamc2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
327 *
328 * -- LAPACK auxiliary routine (version 2.0) --
329 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
330 * Courant Institute, Argonne National Lab, and Rice University
331 * October 31, 1992
332 *
333 * .. Scalar Arguments ..
334  LOGICAL rnd
335  INTEGER beta, emax, emin, t
336  DOUBLE PRECISION eps, rmax, rmin
337 * ..
338 *
339 * Purpose
340 * =======
341 *
342 * DLAMC2 determines the machine parameters specified in its argument
343 * list.
344 *
345 * Arguments
346 * =========
347 *
348 * BETA (output) INTEGER
349 * The base of the machine.
350 *
351 * T (output) INTEGER
352 * The number of ( BETA ) digits in the mantissa.
353 *
354 * RND (output) LOGICAL
355 * Specifies whether proper rounding ( RND = .TRUE. ) or
356 * chopping ( RND = .FALSE. ) occurs in addition. This may not
357 * be a reliable guide to the way in which the machine performs
358 * its arithmetic.
359 *
360 * EPS (output) DOUBLE PRECISION
361 * The smallest positive number such that
362 *
363 * fl( 1.0 - EPS ) .LT. 1.0,
364 *
365 * where fl denotes the computed value.
366 *
367 * EMIN (output) INTEGER
368 * The minimum exponent before (gradual) underflow occurs.
369 *
370 * RMIN (output) DOUBLE PRECISION
371 * The smallest normalized number for the machine, given by
372 * BASE**( EMIN - 1 ), where BASE is the floating point value
373 * of BETA.
374 *
375 * EMAX (output) INTEGER
376 * The maximum exponent before overflow occurs.
377 *
378 * RMAX (output) DOUBLE PRECISION
379 * The largest positive number for the machine, given by
380 * BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
381 * value of BETA.
382 *
383 * Further Details
384 * ===============
385 *
386 * The computation of EPS is based on a routine PARANOIA by
387 * W. Kahan of the University of California at Berkeley.
388 *
389 * =====================================================================
390 *
391 * .. Local Scalars ..
392  LOGICAL first, ieee, iwarn, lieee1, lrnd
393  INTEGER gnmin, gpmin, i, lbeta, lemax, lemin, lt,
394  $ ngnmin, ngpmin
395  DOUBLE PRECISION a, b, c, half, leps, lrmax, lrmin, one, rbase,
396  $ sixth, small, third, two, zero
397 * ..
398 * .. External Functions ..
399  DOUBLE PRECISION dlamc3
400  EXTERNAL dlamc3
401 * ..
402 * .. External Subroutines ..
403  EXTERNAL dlamc1, dlamc4, dlamc5
404 * ..
405 * .. Intrinsic Functions ..
406  INTRINSIC abs, max, min
407 * ..
408 * .. Save statement ..
409  SAVE first, iwarn, lbeta, lemax, lemin, leps, lrmax,
410  $ lrmin, lt
411 * ..
412 * .. Data statements ..
413  DATA first / .true. / , iwarn / .false. /
414 * ..
415 * .. Executable Statements ..
416 *
417  IF( first ) THEN
418  first = .false.
419  zero = 0
420  one = 1
421  two = 2
422 *
423 * LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
424 * BETA, T, RND, EPS, EMIN and RMIN.
425 *
426 * Throughout this routine we use the function DLAMC3 to ensure
427 * that relevant values are stored and not held in registers, or
428 * are not affected by optimizers.
429 *
430 * DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
431 *
432  CALL dlamc1( lbeta, lt, lrnd, lieee1 )
433 *
434 * Start to find EPS.
435 *
436  b = lbeta
437  a = b**( -lt )
438  leps = a
439 *
440 * Try some tricks to see whether or not this is the correct EPS.
441 *
442  b = two / 3
443  half = one / 2
444  sixth = dlamc3( b, -half )
445  third = dlamc3( sixth, sixth )
446  b = dlamc3( third, -half )
447  b = dlamc3( b, sixth )
448  b = abs( b )
449  IF( b.LT.leps )
450  $ b = leps
451 *
452  leps = 1
453 *
454 *+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
455  10 CONTINUE
456  IF( ( leps.GT.b ) .AND. ( b.GT.zero ) ) THEN
457  leps = b
458  c = dlamc3( half*leps, ( two**5 )*( leps**2 ) )
459  c = dlamc3( half, -c )
460  b = dlamc3( half, c )
461  c = dlamc3( half, -b )
462  b = dlamc3( half, c )
463  GO TO 10
464  END IF
465 *+ END WHILE
466 *
467  IF( a.LT.leps )
468  $ leps = a
469 *
470 * Computation of EPS complete.
471 *
472 * Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
473 * Keep dividing A by BETA until (gradual) underflow occurs. This
474 * is detected when we cannot recover the previous A.
475 *
476  rbase = one / lbeta
477  small = one
478  DO 20 i = 1, 3
479  small = dlamc3( small*rbase, zero )
480  20 CONTINUE
481  a = dlamc3( one, small )
482  CALL dlamc4( ngpmin, one, lbeta )
483  CALL dlamc4( ngnmin, -one, lbeta )
484  CALL dlamc4( gpmin, a, lbeta )
485  CALL dlamc4( gnmin, -a, lbeta )
486  ieee = .false.
487 *
488  IF( ( ngpmin.EQ.ngnmin ) .AND. ( gpmin.EQ.gnmin ) ) THEN
489  IF( ngpmin.EQ.gpmin ) THEN
490  lemin = ngpmin
491 * ( Non twos-complement machines, no gradual underflow;
492 * e.g., VAX )
493  ELSE IF( ( gpmin-ngpmin ).EQ.3 ) THEN
494  lemin = ngpmin - 1 + lt
495  ieee = .true.
496 * ( Non twos-complement machines, with gradual underflow;
497 * e.g., IEEE standard followers )
498  ELSE
499  lemin = min( ngpmin, gpmin )
500 * ( A guess; no known machine )
501  iwarn = .true.
502  END IF
503 *
504  ELSE IF( ( ngpmin.EQ.gpmin ) .AND. ( ngnmin.EQ.gnmin ) ) THEN
505  IF( abs( ngpmin-ngnmin ).EQ.1 ) THEN
506  lemin = max( ngpmin, ngnmin )
507 * ( Twos-complement machines, no gradual underflow;
508 * e.g., CYBER 205 )
509  ELSE
510  lemin = min( ngpmin, ngnmin )
511 * ( A guess; no known machine )
512  iwarn = .true.
513  END IF
514 *
515  ELSE IF( ( abs( ngpmin-ngnmin ).EQ.1 ) .AND.
516  $ ( gpmin.EQ.gnmin ) ) THEN
517  IF( ( gpmin-min( ngpmin, ngnmin ) ).EQ.3 ) THEN
518  lemin = max( ngpmin, ngnmin ) - 1 + lt
519 * ( Twos-complement machines with gradual underflow;
520 * no known machine )
521  ELSE
522  lemin = min( ngpmin, ngnmin )
523 * ( A guess; no known machine )
524  iwarn = .true.
525  END IF
526 *
527  ELSE
528  lemin = min( ngpmin, ngnmin, gpmin, gnmin )
529 * ( A guess; no known machine )
530  iwarn = .true.
531  END IF
532 ***
533 * Comment out this if block if EMIN is ok
534  IF( iwarn ) THEN
535  first = .true.
536  WRITE( 6, fmt = 9999 )lemin
537  END IF
538 ***
539 *
540 * Assume IEEE arithmetic if we found denormalised numbers above,
541 * or if arithmetic seems to round in the IEEE style, determined
542 * in routine DLAMC1. A true IEEE machine should have both things
543 * true; however, faulty machines may have one or the other.
544 *
545  ieee = ieee .OR. lieee1
546 *
547 * Compute RMIN by successive division by BETA. We could compute
548 * RMIN as BASE**( EMIN - 1 ), but some machines underflow during
549 * this computation.
550 *
551  lrmin = 1
552  DO 30 i = 1, 1 - lemin
553  lrmin = dlamc3( lrmin*rbase, zero )
554  30 CONTINUE
555 *
556 * Finally, call DLAMC5 to compute EMAX and RMAX.
557 *
558  CALL dlamc5( lbeta, lt, lemin, ieee, lemax, lrmax )
559  END IF
560 *
561  beta = lbeta
562  t = lt
563  rnd = lrnd
564  eps = leps
565  emin = lemin
566  rmin = lrmin
567  emax = lemax
568  rmax = lrmax
569 *
570  RETURN
571 *
572  9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
573  $ ' EMIN = ', i8, /
574  $ ' If, after inspection, the value EMIN looks',
575  $ ' acceptable please comment out ',
576  $ / ' the IF block as marked within the code of routine',
577  $ ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / )
578 *
579 * End of DLAMC2
580 *
581  END
582 *
583 ************************************************************************
584 *
585  DOUBLE PRECISION FUNCTION dlamc3( A, B )
586 *
587 * -- LAPACK auxiliary routine (version 2.0) --
588 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
589 * Courant Institute, Argonne National Lab, and Rice University
590 * October 31, 1992
591 *
592 * .. Scalar Arguments ..
593  DOUBLE PRECISION a, b
594 * ..
595 *
596 * Purpose
597 * =======
598 *
599 * DLAMC3 is intended to force A and B to be stored prior to doing
600 * the addition of A and B , for use in situations where optimizers
601 * might hold one of these in a register.
602 *
603 * Arguments
604 * =========
605 *
606 * A, B (input) DOUBLE PRECISION
607 * The values A and B.
608 *
609 * =====================================================================
610 *
611 * .. Executable Statements ..
612 *
613  dlamc3 = a + b
614 *
615  RETURN
616 *
617 * End of DLAMC3
618 *
619  END
620 *
621 ************************************************************************
622 *
623  SUBROUTINE dlamc4( EMIN, START, BASE )
624 *
625 * -- LAPACK auxiliary routine (version 2.0) --
626 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
627 * Courant Institute, Argonne National Lab, and Rice University
628 * October 31, 1992
629 *
630 * .. Scalar Arguments ..
631  INTEGER base, emin
632  DOUBLE PRECISION start
633 * ..
634 *
635 * Purpose
636 * =======
637 *
638 * DLAMC4 is a service routine for DLAMC2.
639 *
640 * Arguments
641 * =========
642 *
643 * EMIN (output) EMIN
644 * The minimum exponent before (gradual) underflow, computed by
645 * setting A = START and dividing by BASE until the previous A
646 * can not be recovered.
647 *
648 * START (input) DOUBLE PRECISION
649 * The starting point for determining EMIN.
650 *
651 * BASE (input) INTEGER
652 * The base of the machine.
653 *
654 * =====================================================================
655 *
656 * .. Local Scalars ..
657  INTEGER i
658  DOUBLE PRECISION a, b1, b2, c1, c2, d1, d2, one, rbase, zero
659 * ..
660 * .. External Functions ..
661  DOUBLE PRECISION dlamc3
662  EXTERNAL dlamc3
663 * ..
664 * .. Executable Statements ..
665 *
666  a = start
667  one = 1
668  rbase = one / base
669  zero = 0
670  emin = 1
671  b1 = dlamc3( a*rbase, zero )
672  c1 = a
673  c2 = a
674  d1 = a
675  d2 = a
676 *+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
677 * $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP
678  10 CONTINUE
679  IF( ( c1.EQ.a ) .AND. ( c2.EQ.a ) .AND. ( d1.EQ.a ) .AND.
680  $ ( d2.EQ.a ) ) THEN
681  emin = emin - 1
682  a = b1
683  b1 = dlamc3( a / base, zero )
684  c1 = dlamc3( b1*base, zero )
685  d1 = zero
686  DO 20 i = 1, base
687  d1 = d1 + b1
688  20 CONTINUE
689  b2 = dlamc3( a*rbase, zero )
690  c2 = dlamc3( b2 / rbase, zero )
691  d2 = zero
692  DO 30 i = 1, base
693  d2 = d2 + b2
694  30 CONTINUE
695  GO TO 10
696  END IF
697 *+ END WHILE
698 *
699  RETURN
700 *
701 * End of DLAMC4
702 *
703  END
704 *
705 ************************************************************************
706 *
707  SUBROUTINE dlamc5( BETA, P, EMIN, IEEE, EMAX, RMAX )
708 *
709 * -- LAPACK auxiliary routine (version 2.0) --
710 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
711 * Courant Institute, Argonne National Lab, and Rice University
712 * October 31, 1992
713 *
714 * .. Scalar Arguments ..
715  LOGICAL ieee
716  INTEGER beta, emax, emin, p
717  DOUBLE PRECISION rmax
718 * ..
719 *
720 * Purpose
721 * =======
722 *
723 * DLAMC5 attempts to compute RMAX, the largest machine floating-point
724 * number, without overflow. It assumes that EMAX + abs(EMIN) sum
725 * approximately to a power of 2. It will fail on machines where this
726 * assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
727 * EMAX = 28718). It will also fail if the value supplied for EMIN is
728 * too large (i.e. too close to zero), probably with overflow.
729 *
730 * Arguments
731 * =========
732 *
733 * BETA (input) INTEGER
734 * The base of floating-point arithmetic.
735 *
736 * P (input) INTEGER
737 * The number of base BETA digits in the mantissa of a
738 * floating-point value.
739 *
740 * EMIN (input) INTEGER
741 * The minimum exponent before (gradual) underflow.
742 *
743 * IEEE (input) LOGICAL
744 * A logical flag specifying whether or not the arithmetic
745 * system is thought to comply with the IEEE standard.
746 *
747 * EMAX (output) INTEGER
748 * The largest exponent before overflow
749 *
750 * RMAX (output) DOUBLE PRECISION
751 * The largest machine floating-point number.
752 *
753 * =====================================================================
754 *
755 * .. Parameters ..
756  DOUBLE PRECISION zero, one
757  parameter( zero = 0.0d0, one = 1.0d0 )
758 * ..
759 * .. Local Scalars ..
760  INTEGER exbits, expsum, i, lexp, nbits, try, uexp
761  DOUBLE PRECISION oldy, recbas, y, z
762 * ..
763 * .. External Functions ..
764  DOUBLE PRECISION dlamc3
765  EXTERNAL dlamc3
766 * ..
767 * .. Intrinsic Functions ..
768  INTRINSIC mod
769 * ..
770 * .. Executable Statements ..
771 *
772 * First compute LEXP and UEXP, two powers of 2 that bound
773 * abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
774 * approximately to the bound that is closest to abs(EMIN).
775 * (EMAX is the exponent of the required number RMAX).
776 *
777  lexp = 1
778  exbits = 1
779  10 CONTINUE
780  try = lexp*2
781  IF( try.LE.( -emin ) ) THEN
782  lexp = try
783  exbits = exbits + 1
784  GO TO 10
785  END IF
786  IF( lexp.EQ.-emin ) THEN
787  uexp = lexp
788  ELSE
789  uexp = try
790  exbits = exbits + 1
791  END IF
792 *
793 * Now -LEXP is less than or equal to EMIN, and -UEXP is greater
794 * than or equal to EMIN. EXBITS is the number of bits needed to
795 * store the exponent.
796 *
797  IF( ( uexp+emin ).GT.( -lexp-emin ) ) THEN
798  expsum = 2*lexp
799  ELSE
800  expsum = 2*uexp
801  END IF
802 *
803 * EXPSUM is the exponent range, approximately equal to
804 * EMAX - EMIN + 1 .
805 *
806  emax = expsum + emin - 1
807  nbits = 1 + exbits + p
808 *
809 * NBITS is the total number of bits needed to store a
810 * floating-point number.
811 *
812  IF( ( mod( nbits, 2 ).EQ.1 ) .AND. ( beta.EQ.2 ) ) THEN
813 *
814 * Either there are an odd number of bits used to store a
815 * floating-point number, which is unlikely, or some bits are
816 * not used in the representation of numbers, which is possible,
817 * (e.g. Cray machines) or the mantissa has an implicit bit,
818 * (e.g. IEEE machines, Dec Vax machines), which is perhaps the
819 * most likely. We have to assume the last alternative.
820 * If this is true, then we need to reduce EMAX by one because
821 * there must be some way of representing zero in an implicit-bit
822 * system. On machines like Cray, we are reducing EMAX by one
823 * unnecessarily.
824 *
825  emax = emax - 1
826  END IF
827 *
828  IF( ieee ) THEN
829 *
830 * Assume we are on an IEEE machine which reserves one exponent
831 * for infinity and NaN.
832 *
833  emax = emax - 1
834  END IF
835 *
836 * Now create RMAX, the largest machine number, which should
837 * be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
838 *
839 * First compute 1.0 - BETA**(-P), being careful that the
840 * result is less than 1.0 .
841 *
842  recbas = one / beta
843  z = beta - one
844  y = zero
845  DO 20 i = 1, p
846  z = z*recbas
847  IF( y.LT.one )
848  $ oldy = y
849  y = dlamc3( y, z )
850  20 CONTINUE
851  IF( y.GE.one )
852  $ y = oldy
853 *
854 * Now multiply by BETA**EMAX to get RMAX.
855 *
856  DO 30 i = 1, emax
857  y = dlamc3( y*beta, zero )
858  30 CONTINUE
859 *
860  rmax = y
861  RETURN
862 *
863 * End of DLAMC5
864 *
865  END
866  REAL function slamch( cmach )
867 *
868 * -- LAPACK auxiliary routine (version 2.0) --
869 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
870 * Courant Institute, Argonne National Lab, and Rice University
871 * October 31, 1992
872 *
873 * .. Scalar Arguments ..
874  CHARACTER cmach
875 * ..
876 *
877 * Purpose
878 * =======
879 *
880 * SLAMCH determines single precision machine parameters.
881 *
882 * Arguments
883 * =========
884 *
885 * CMACH (input) CHARACTER*1
886 * Specifies the value to be returned by SLAMCH:
887 * = 'E' or 'e', SLAMCH := eps
888 * = 'S' or 's , SLAMCH := sfmin
889 * = 'B' or 'b', SLAMCH := base
890 * = 'P' or 'p', SLAMCH := eps*base
891 * = 'N' or 'n', SLAMCH := t
892 * = 'R' or 'r', SLAMCH := rnd
893 * = 'M' or 'm', SLAMCH := emin
894 * = 'U' or 'u', SLAMCH := rmin
895 * = 'L' or 'l', SLAMCH := emax
896 * = 'O' or 'o', SLAMCH := rmax
897 *
898 * where
899 *
900 * eps = relative machine precision
901 * sfmin = safe minimum, such that 1/sfmin does not overflow
902 * base = base of the machine
903 * prec = eps*base
904 * t = number of (base) digits in the mantissa
905 * rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
906 * emin = minimum exponent before (gradual) underflow
907 * rmin = underflow threshold - base**(emin-1)
908 * emax = largest exponent before overflow
909 * rmax = overflow threshold - (base**emax)*(1-eps)
910 *
911 * =====================================================================
912 *
913 * .. Parameters ..
914  REAL one, zero
915  parameter( one = 1.0e+0, zero = 0.0e+0 )
916 * ..
917 * .. Local Scalars ..
918  LOGICAL first, lrnd
919  INTEGER beta, imax, imin, it
920  REAL base, emax, emin, eps, prec, rmach, rmax, rmin,
921  $ rnd, sfmin, small, t
922 * ..
923 * .. External Functions ..
924  LOGICAL lsame
925  EXTERNAL lsame
926 * ..
927 * .. External Subroutines ..
928  EXTERNAL slamc2
929 * ..
930 * .. Save statement ..
931  SAVE first, eps, sfmin, base, t, rnd, emin, rmin,
932  $ emax, rmax, prec
933 * ..
934 * .. Data statements ..
935  DATA first / .true. /
936 * ..
937 * .. Executable Statements ..
938 *
939  IF( first ) THEN
940  first = .false.
941  CALL slamc2( beta, it, lrnd, eps, imin, rmin, imax, rmax )
942  base = beta
943  t = it
944  IF( lrnd ) THEN
945  rnd = one
946  eps = ( base**( 1-it ) ) / 2
947  ELSE
948  rnd = zero
949  eps = base**( 1-it )
950  END IF
951  prec = eps*base
952  emin = imin
953  emax = imax
954  sfmin = rmin
955  small = one / rmax
956  IF( small.GE.sfmin ) THEN
957 *
958 * Use SMALL plus a bit, to avoid the possibility of rounding
959 * causing overflow when computing 1/sfmin.
960 *
961  sfmin = small*( one+eps )
962  END IF
963  END IF
964 *
965  IF( lsame( cmach, 'E' ) ) THEN
966  rmach = eps
967  ELSE IF( lsame( cmach, 'S' ) ) THEN
968  rmach = sfmin
969  ELSE IF( lsame( cmach, 'B' ) ) THEN
970  rmach = base
971  ELSE IF( lsame( cmach, 'P' ) ) THEN
972  rmach = prec
973  ELSE IF( lsame( cmach, 'N' ) ) THEN
974  rmach = t
975  ELSE IF( lsame( cmach, 'R' ) ) THEN
976  rmach = rnd
977  ELSE IF( lsame( cmach, 'M' ) ) THEN
978  rmach = emin
979  ELSE IF( lsame( cmach, 'U' ) ) THEN
980  rmach = rmin
981  ELSE IF( lsame( cmach, 'L' ) ) THEN
982  rmach = emax
983  ELSE IF( lsame( cmach, 'O' ) ) THEN
984  rmach = rmax
985  END IF
986 *
987  slamch = rmach
988  RETURN
989 *
990 * End of SLAMCH
991 *
992  END
993 *
994 ************************************************************************
995 *
996  SUBROUTINE slamc1( BETA, T, RND, IEEE1 )
997 *
998 * -- LAPACK auxiliary routine (version 2.0) --
999 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1000 * Courant Institute, Argonne National Lab, and Rice University
1001 * October 31, 1992
1002 *
1003 * .. Scalar Arguments ..
1004  LOGICAL IEEE1, RND
1005  INTEGER BETA, T
1006 * ..
1007 *
1008 * Purpose
1009 * =======
1010 *
1011 * SLAMC1 determines the machine parameters given by BETA, T, RND, and
1012 * IEEE1.
1013 *
1014 * Arguments
1015 * =========
1016 *
1017 * BETA (output) INTEGER
1018 * The base of the machine.
1019 *
1020 * T (output) INTEGER
1021 * The number of ( BETA ) digits in the mantissa.
1022 *
1023 * RND (output) LOGICAL
1024 * Specifies whether proper rounding ( RND = .TRUE. ) or
1025 * chopping ( RND = .FALSE. ) occurs in addition. This may not
1026 * be a reliable guide to the way in which the machine performs
1027 * its arithmetic.
1028 *
1029 * IEEE1 (output) LOGICAL
1030 * Specifies whether rounding appears to be done in the IEEE
1031 * 'round to nearest' style.
1032 *
1033 * Further Details
1034 * ===============
1035 *
1036 * The routine is based on the routine ENVRON by Malcolm and
1037 * incorporates suggestions by Gentleman and Marovich. See
1038 *
1039 * Malcolm M. A. (1972) Algorithms to reveal properties of
1040 * floating-point arithmetic. Comms. of the ACM, 15, 949-951.
1041 *
1042 * Gentleman W. M. and Marovich S. B. (1974) More on algorithms
1043 * that reveal properties of floating point arithmetic units.
1044 * Comms. of the ACM, 17, 276-277.
1045 *
1046 * =====================================================================
1047 *
1048 * .. Local Scalars ..
1049  LOGICAL FIRST, LIEEE1, LRND
1050  INTEGER LBETA, LT
1051  REAL A, B, C, F, ONE, QTR, SAVEC, T1, T2
1052 * ..
1053 * .. External Functions ..
1054  REAL SLAMC3
1055  EXTERNAL slamc3
1056 * ..
1057 * .. Save statement ..
1058  SAVE first, lieee1, lbeta, lrnd, lt
1059 * ..
1060 * .. Data statements ..
1061  DATA first / .true. /
1062 * ..
1063 * .. Executable Statements ..
1064 *
1065  IF( first ) THEN
1066  first = .false.
1067  one = 1
1068 *
1069 * LBETA, LIEEE1, LT and LRND are the local values of BETA,
1070 * IEEE1, T and RND.
1071 *
1072 * Throughout this routine we use the function SLAMC3 to ensure
1073 * that relevant values are stored and not held in registers, or
1074 * are not affected by optimizers.
1075 *
1076 * Compute a = 2.0**m with the smallest positive integer m such
1077 * that
1078 *
1079 * fl( a + 1.0 ) = a.
1080 *
1081  a = 1
1082  c = 1
1083 *
1084 *+ WHILE( C.EQ.ONE )LOOP
1085  10 CONTINUE
1086  IF( c.EQ.one ) THEN
1087  a = 2*a
1088  c = slamc3( a, one )
1089  c = slamc3( c, -a )
1090  GO TO 10
1091  END IF
1092 *+ END WHILE
1093 *
1094 * Now compute b = 2.0**m with the smallest positive integer m
1095 * such that
1096 *
1097 * fl( a + b ) .gt. a.
1098 *
1099  b = 1
1100  c = slamc3( a, b )
1101 *
1102 *+ WHILE( C.EQ.A )LOOP
1103  20 CONTINUE
1104  IF( c.EQ.a ) THEN
1105  b = 2*b
1106  c = slamc3( a, b )
1107  GO TO 20
1108  END IF
1109 *+ END WHILE
1110 *
1111 * Now compute the base. a and c are neighbouring floating point
1112 * numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
1113 * their difference is beta. Adding 0.25 to c is to ensure that it
1114 * is truncated to beta and not ( beta - 1 ).
1115 *
1116  qtr = one / 4
1117  savec = c
1118  c = slamc3( c, -a )
1119  lbeta = c + qtr
1120 *
1121 * Now determine whether rounding or chopping occurs, by adding a
1122 * bit less than beta/2 and a bit more than beta/2 to a.
1123 *
1124  b = lbeta
1125  f = slamc3( b / 2, -b / 100 )
1126  c = slamc3( f, a )
1127  IF( c.EQ.a ) THEN
1128  lrnd = .true.
1129  ELSE
1130  lrnd = .false.
1131  END IF
1132  f = slamc3( b / 2, b / 100 )
1133  c = slamc3( f, a )
1134  IF( ( lrnd ) .AND. ( c.EQ.a ) )
1135  $ lrnd = .false.
1136 *
1137 * Try and decide whether rounding is done in the IEEE 'round to
1138 * nearest' style. B/2 is half a unit in the last place of the two
1139 * numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
1140 * zero, and SAVEC is odd. Thus adding B/2 to A should not change
1141 * A, but adding B/2 to SAVEC should change SAVEC.
1142 *
1143  t1 = slamc3( b / 2, a )
1144  t2 = slamc3( b / 2, savec )
1145  lieee1 = ( t1.EQ.a ) .AND. ( t2.GT.savec ) .AND. lrnd
1146 *
1147 * Now find the mantissa, t. It should be the integer part of
1148 * log to the base beta of a, however it is safer to determine t
1149 * by powering. So we find t as the smallest positive integer for
1150 * which
1151 *
1152 * fl( beta**t + 1.0 ) = 1.0.
1153 *
1154  lt = 0
1155  a = 1
1156  c = 1
1157 *
1158 *+ WHILE( C.EQ.ONE )LOOP
1159  30 CONTINUE
1160  IF( c.EQ.one ) THEN
1161  lt = lt + 1
1162  a = a*lbeta
1163  c = slamc3( a, one )
1164  c = slamc3( c, -a )
1165  GO TO 30
1166  END IF
1167 *+ END WHILE
1168 *
1169  END IF
1170 *
1171  beta = lbeta
1172  t = lt
1173  rnd = lrnd
1174  ieee1 = lieee1
1175  RETURN
1176 *
1177 * End of SLAMC1
1178 *
1179  END
1180 *
1181 ************************************************************************
1182 *
1183  SUBROUTINE slamc2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
1184 *
1185 * -- LAPACK auxiliary routine (version 2.0) --
1186 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1187 * Courant Institute, Argonne National Lab, and Rice University
1188 * October 31, 1992
1189 *
1190 * .. Scalar Arguments ..
1191  LOGICAL RND
1192  INTEGER BETA, EMAX, EMIN, T
1193  REAL EPS, RMAX, RMIN
1194 * ..
1195 *
1196 * Purpose
1197 * =======
1198 *
1199 * SLAMC2 determines the machine parameters specified in its argument
1200 * list.
1201 *
1202 * Arguments
1203 * =========
1204 *
1205 * BETA (output) INTEGER
1206 * The base of the machine.
1207 *
1208 * T (output) INTEGER
1209 * The number of ( BETA ) digits in the mantissa.
1210 *
1211 * RND (output) LOGICAL
1212 * Specifies whether proper rounding ( RND = .TRUE. ) or
1213 * chopping ( RND = .FALSE. ) occurs in addition. This may not
1214 * be a reliable guide to the way in which the machine performs
1215 * its arithmetic.
1216 *
1217 * EPS (output) REAL
1218 * The smallest positive number such that
1219 *
1220 * fl( 1.0 - EPS ) .LT. 1.0,
1221 *
1222 * where fl denotes the computed value.
1223 *
1224 * EMIN (output) INTEGER
1225 * The minimum exponent before (gradual) underflow occurs.
1226 *
1227 * RMIN (output) REAL
1228 * The smallest normalized number for the machine, given by
1229 * BASE**( EMIN - 1 ), where BASE is the floating point value
1230 * of BETA.
1231 *
1232 * EMAX (output) INTEGER
1233 * The maximum exponent before overflow occurs.
1234 *
1235 * RMAX (output) REAL
1236 * The largest positive number for the machine, given by
1237 * BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
1238 * value of BETA.
1239 *
1240 * Further Details
1241 * ===============
1242 *
1243 * The computation of EPS is based on a routine PARANOIA by
1244 * W. Kahan of the University of California at Berkeley.
1245 *
1246 * =====================================================================
1247 *
1248 * .. Local Scalars ..
1249  LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND
1250  INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
1251  $ NGNMIN, NGPMIN
1252  REAL A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
1253  $ SIXTH, SMALL, THIRD, TWO, ZERO
1254 * ..
1255 * .. External Functions ..
1256  REAL SLAMC3
1257  EXTERNAL slamc3
1258 * ..
1259 * .. External Subroutines ..
1260  EXTERNAL slamc1, slamc4, slamc5
1261 * ..
1262 * .. Intrinsic Functions ..
1263  INTRINSIC abs, max, min
1264 * ..
1265 * .. Save statement ..
1266  SAVE first, iwarn, lbeta, lemax, lemin, leps, lrmax,
1267  $ lrmin, lt
1268 * ..
1269 * .. Data statements ..
1270  DATA first / .true. / , iwarn / .false. /
1271 * ..
1272 * .. Executable Statements ..
1273 *
1274  IF( first ) THEN
1275  first = .false.
1276  zero = 0
1277  one = 1
1278  two = 2
1279 *
1280 * LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
1281 * BETA, T, RND, EPS, EMIN and RMIN.
1282 *
1283 * Throughout this routine we use the function SLAMC3 to ensure
1284 * that relevant values are stored and not held in registers, or
1285 * are not affected by optimizers.
1286 *
1287 * SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
1288 *
1289  CALL slamc1( lbeta, lt, lrnd, lieee1 )
1290 *
1291 * Start to find EPS.
1292 *
1293  b = lbeta
1294  a = b**( -lt )
1295  leps = a
1296 *
1297 * Try some tricks to see whether or not this is the correct EPS.
1298 *
1299  b = two / 3
1300  half = one / 2
1301  sixth = slamc3( b, -half )
1302  third = slamc3( sixth, sixth )
1303  b = slamc3( third, -half )
1304  b = slamc3( b, sixth )
1305  b = abs( b )
1306  IF( b.LT.leps )
1307  $ b = leps
1308 *
1309  leps = 1
1310 *
1311 *+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
1312  10 CONTINUE
1313  IF( ( leps.GT.b ) .AND. ( b.GT.zero ) ) THEN
1314  leps = b
1315  c = slamc3( half*leps, ( two**5 )*( leps**2 ) )
1316  c = slamc3( half, -c )
1317  b = slamc3( half, c )
1318  c = slamc3( half, -b )
1319  b = slamc3( half, c )
1320  GO TO 10
1321  END IF
1322 *+ END WHILE
1323 *
1324  IF( a.LT.leps )
1325  $ leps = a
1326 *
1327 * Computation of EPS complete.
1328 *
1329 * Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
1330 * Keep dividing A by BETA until (gradual) underflow occurs. This
1331 * is detected when we cannot recover the previous A.
1332 *
1333  rbase = one / lbeta
1334  small = one
1335  DO 20 i = 1, 3
1336  small = slamc3( small*rbase, zero )
1337  20 CONTINUE
1338  a = slamc3( one, small )
1339  CALL slamc4( ngpmin, one, lbeta )
1340  CALL slamc4( ngnmin, -one, lbeta )
1341  CALL slamc4( gpmin, a, lbeta )
1342  CALL slamc4( gnmin, -a, lbeta )
1343  ieee = .false.
1344 *
1345  IF( ( ngpmin.EQ.ngnmin ) .AND. ( gpmin.EQ.gnmin ) ) THEN
1346  IF( ngpmin.EQ.gpmin ) THEN
1347  lemin = ngpmin
1348 * ( Non twos-complement machines, no gradual underflow;
1349 * e.g., VAX )
1350  ELSE IF( ( gpmin-ngpmin ).EQ.3 ) THEN
1351  lemin = ngpmin - 1 + lt
1352  ieee = .true.
1353 * ( Non twos-complement machines, with gradual underflow;
1354 * e.g., IEEE standard followers )
1355  ELSE
1356  lemin = min( ngpmin, gpmin )
1357 * ( A guess; no known machine )
1358  iwarn = .true.
1359  END IF
1360 *
1361  ELSE IF( ( ngpmin.EQ.gpmin ) .AND. ( ngnmin.EQ.gnmin ) ) THEN
1362  IF( abs( ngpmin-ngnmin ).EQ.1 ) THEN
1363  lemin = max( ngpmin, ngnmin )
1364 * ( Twos-complement machines, no gradual underflow;
1365 * e.g., CYBER 205 )
1366  ELSE
1367  lemin = min( ngpmin, ngnmin )
1368 * ( A guess; no known machine )
1369  iwarn = .true.
1370  END IF
1371 *
1372  ELSE IF( ( abs( ngpmin-ngnmin ).EQ.1 ) .AND.
1373  $ ( gpmin.EQ.gnmin ) ) THEN
1374  IF( ( gpmin-min( ngpmin, ngnmin ) ).EQ.3 ) THEN
1375  lemin = max( ngpmin, ngnmin ) - 1 + lt
1376 * ( Twos-complement machines with gradual underflow;
1377 * no known machine )
1378  ELSE
1379  lemin = min( ngpmin, ngnmin )
1380 * ( A guess; no known machine )
1381  iwarn = .true.
1382  END IF
1383 *
1384  ELSE
1385  lemin = min( ngpmin, ngnmin, gpmin, gnmin )
1386 * ( A guess; no known machine )
1387  iwarn = .true.
1388  END IF
1389 ***
1390 * Comment out this if block if EMIN is ok
1391  IF( iwarn ) THEN
1392  first = .true.
1393  WRITE( 6, fmt = 9999 )lemin
1394  END IF
1395 ***
1396 *
1397 * Assume IEEE arithmetic if we found denormalised numbers above,
1398 * or if arithmetic seems to round in the IEEE style, determined
1399 * in routine SLAMC1. A true IEEE machine should have both things
1400 * true; however, faulty machines may have one or the other.
1401 *
1402  ieee = ieee .OR. lieee1
1403 *
1404 * Compute RMIN by successive division by BETA. We could compute
1405 * RMIN as BASE**( EMIN - 1 ), but some machines underflow during
1406 * this computation.
1407 *
1408  lrmin = 1
1409  DO 30 i = 1, 1 - lemin
1410  lrmin = slamc3( lrmin*rbase, zero )
1411  30 CONTINUE
1412 *
1413 * Finally, call SLAMC5 to compute EMAX and RMAX.
1414 *
1415  CALL slamc5( lbeta, lt, lemin, ieee, lemax, lrmax )
1416  END IF
1417 *
1418  beta = lbeta
1419  t = lt
1420  rnd = lrnd
1421  eps = leps
1422  emin = lemin
1423  rmin = lrmin
1424  emax = lemax
1425  rmax = lrmax
1426 *
1427  RETURN
1428 *
1429  9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
1430  $ ' EMIN = ', i8, /
1431  $ ' If, after inspection, the value EMIN looks',
1432  $ ' acceptable please comment out ',
1433  $ / ' the IF block as marked within the code of routine',
1434  $ ' SLAMC2,', / ' otherwise supply EMIN explicitly.', / )
1435 *
1436 * End of SLAMC2
1437 *
1438  END
1439 *
1440 ************************************************************************
1441 *
1442  REAL FUNCTION SLAMC3( A, B )
1443 *
1444 * -- LAPACK auxiliary routine (version 2.0) --
1445 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1446 * Courant Institute, Argonne National Lab, and Rice University
1447 * October 31, 1992
1448 *
1449 * .. Scalar Arguments ..
1450  REAL a, b
1451 * ..
1452 *
1453 * Purpose
1454 * =======
1455 *
1456 * SLAMC3 is intended to force A and B to be stored prior to doing
1457 * the addition of A and B , for use in situations where optimizers
1458 * might hold one of these in a register.
1459 *
1460 * Arguments
1461 * =========
1462 *
1463 * A, B (input) REAL
1464 * The values A and B.
1465 *
1466 * =====================================================================
1467 *
1468 * .. Executable Statements ..
1469 *
1470  slamc3 = a + b
1471 *
1472  RETURN
1473 *
1474 * End of SLAMC3
1475 *
1476  END
1477 *
1478 ************************************************************************
1479 *
1480  SUBROUTINE slamc4( EMIN, START, BASE )
1481 *
1482 * -- LAPACK auxiliary routine (version 2.0) --
1483 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1484 * Courant Institute, Argonne National Lab, and Rice University
1485 * October 31, 1992
1486 *
1487 * .. Scalar Arguments ..
1488  INTEGER BASE, EMIN
1489  REAL START
1490 * ..
1491 *
1492 * Purpose
1493 * =======
1494 *
1495 * SLAMC4 is a service routine for SLAMC2.
1496 *
1497 * Arguments
1498 * =========
1499 *
1500 * EMIN (output) EMIN
1501 * The minimum exponent before (gradual) underflow, computed by
1502 * setting A = START and dividing by BASE until the previous A
1503 * can not be recovered.
1504 *
1505 * START (input) REAL
1506 * The starting point for determining EMIN.
1507 *
1508 * BASE (input) INTEGER
1509 * The base of the machine.
1510 *
1511 * =====================================================================
1512 *
1513 * .. Local Scalars ..
1514  INTEGER I
1515  REAL A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
1516 * ..
1517 * .. External Functions ..
1518  REAL SLAMC3
1519  EXTERNAL slamc3
1520 * ..
1521 * .. Executable Statements ..
1522 *
1523  a = start
1524  one = 1
1525  rbase = one / base
1526  zero = 0
1527  emin = 1
1528  b1 = slamc3( a*rbase, zero )
1529  c1 = a
1530  c2 = a
1531  d1 = a
1532  d2 = a
1533 *+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
1534 * $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP
1535  10 CONTINUE
1536  IF( ( c1.EQ.a ) .AND. ( c2.EQ.a ) .AND. ( d1.EQ.a ) .AND.
1537  $ ( d2.EQ.a ) ) THEN
1538  emin = emin - 1
1539  a = b1
1540  b1 = slamc3( a / base, zero )
1541  c1 = slamc3( b1*base, zero )
1542  d1 = zero
1543  DO 20 i = 1, base
1544  d1 = d1 + b1
1545  20 CONTINUE
1546  b2 = slamc3( a*rbase, zero )
1547  c2 = slamc3( b2 / rbase, zero )
1548  d2 = zero
1549  DO 30 i = 1, base
1550  d2 = d2 + b2
1551  30 CONTINUE
1552  GO TO 10
1553  END IF
1554 *+ END WHILE
1555 *
1556  RETURN
1557 *
1558 * End of SLAMC4
1559 *
1560  END
1561 *
1562 ************************************************************************
1563 *
1564  SUBROUTINE slamc5( BETA, P, EMIN, IEEE, EMAX, RMAX )
1565 *
1566 * -- LAPACK auxiliary routine (version 2.0) --
1567 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1568 * Courant Institute, Argonne National Lab, and Rice University
1569 * October 31, 1992
1570 *
1571 * .. Scalar Arguments ..
1572  LOGICAL IEEE
1573  INTEGER BETA, EMAX, EMIN, P
1574  REAL RMAX
1575 * ..
1576 *
1577 * Purpose
1578 * =======
1579 *
1580 * SLAMC5 attempts to compute RMAX, the largest machine floating-point
1581 * number, without overflow. It assumes that EMAX + abs(EMIN) sum
1582 * approximately to a power of 2. It will fail on machines where this
1583 * assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
1584 * EMAX = 28718). It will also fail if the value supplied for EMIN is
1585 * too large (i.e. too close to zero), probably with overflow.
1586 *
1587 * Arguments
1588 * =========
1589 *
1590 * BETA (input) INTEGER
1591 * The base of floating-point arithmetic.
1592 *
1593 * P (input) INTEGER
1594 * The number of base BETA digits in the mantissa of a
1595 * floating-point value.
1596 *
1597 * EMIN (input) INTEGER
1598 * The minimum exponent before (gradual) underflow.
1599 *
1600 * IEEE (input) LOGICAL
1601 * A logical flag specifying whether or not the arithmetic
1602 * system is thought to comply with the IEEE standard.
1603 *
1604 * EMAX (output) INTEGER
1605 * The largest exponent before overflow
1606 *
1607 * RMAX (output) REAL
1608 * The largest machine floating-point number.
1609 *
1610 * =====================================================================
1611 *
1612 * .. Parameters ..
1613  REAL ZERO, ONE
1614  parameter( zero = 0.0e0, one = 1.0e0 )
1615 * ..
1616 * .. Local Scalars ..
1617  INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
1618  REAL OLDY, RECBAS, Y, Z
1619 * ..
1620 * .. External Functions ..
1621  REAL SLAMC3
1622  EXTERNAL slamc3
1623 * ..
1624 * .. Intrinsic Functions ..
1625  INTRINSIC mod
1626 * ..
1627 * .. Executable Statements ..
1628 *
1629 * First compute LEXP and UEXP, two powers of 2 that bound
1630 * abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
1631 * approximately to the bound that is closest to abs(EMIN).
1632 * (EMAX is the exponent of the required number RMAX).
1633 *
1634  lexp = 1
1635  exbits = 1
1636  10 CONTINUE
1637  try = lexp*2
1638  IF( try.LE.( -emin ) ) THEN
1639  lexp = try
1640  exbits = exbits + 1
1641  GO TO 10
1642  END IF
1643  IF( lexp.EQ.-emin ) THEN
1644  uexp = lexp
1645  ELSE
1646  uexp = try
1647  exbits = exbits + 1
1648  END IF
1649 *
1650 * Now -LEXP is less than or equal to EMIN, and -UEXP is greater
1651 * than or equal to EMIN. EXBITS is the number of bits needed to
1652 * store the exponent.
1653 *
1654  IF( ( uexp+emin ).GT.( -lexp-emin ) ) THEN
1655  expsum = 2*lexp
1656  ELSE
1657  expsum = 2*uexp
1658  END IF
1659 *
1660 * EXPSUM is the exponent range, approximately equal to
1661 * EMAX - EMIN + 1 .
1662 *
1663  emax = expsum + emin - 1
1664  nbits = 1 + exbits + p
1665 *
1666 * NBITS is the total number of bits needed to store a
1667 * floating-point number.
1668 *
1669  IF( ( mod( nbits, 2 ).EQ.1 ) .AND. ( beta.EQ.2 ) ) THEN
1670 *
1671 * Either there are an odd number of bits used to store a
1672 * floating-point number, which is unlikely, or some bits are
1673 * not used in the representation of numbers, which is possible,
1674 * (e.g. Cray machines) or the mantissa has an implicit bit,
1675 * (e.g. IEEE machines, Dec Vax machines), which is perhaps the
1676 * most likely. We have to assume the last alternative.
1677 * If this is true, then we need to reduce EMAX by one because
1678 * there must be some way of representing zero in an implicit-bit
1679 * system. On machines like Cray, we are reducing EMAX by one
1680 * unnecessarily.
1681 *
1682  emax = emax - 1
1683  END IF
1684 *
1685  IF( ieee ) THEN
1686 *
1687 * Assume we are on an IEEE machine which reserves one exponent
1688 * for infinity and NaN.
1689 *
1690  emax = emax - 1
1691  END IF
1692 *
1693 * Now create RMAX, the largest machine number, which should
1694 * be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
1695 *
1696 * First compute 1.0 - BETA**(-P), being careful that the
1697 * result is less than 1.0 .
1698 *
1699  recbas = one / beta
1700  z = beta - one
1701  y = zero
1702  DO 20 i = 1, p
1703  z = z*recbas
1704  IF( y.LT.one )
1705  $ oldy = y
1706  y = slamc3( y, z )
1707  20 CONTINUE
1708  IF( y.GE.one )
1709  $ y = oldy
1710 *
1711 * Now multiply by BETA**EMAX to get RMAX.
1712 *
1713  DO 30 i = 1, emax
1714  y = slamc3( y*beta, zero )
1715  30 CONTINUE
1716 *
1717  rmax = y
1718  RETURN
1719 *
1720 * End of SLAMC5
1721 *
1722  END
1723  LOGICAL FUNCTION lsame( CA, CB )
1724 *
1725 * -- LAPACK auxiliary routine (version 2.0) --
1726 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1727 * Courant Institute, Argonne National Lab, and Rice University
1728 * June 30, 1994
1729 *
1730 * .. Scalar Arguments ..
1731  CHARACTER ca, cb
1732 * ..
1733 *
1734 * Purpose
1735 * =======
1736 *
1737 * LSAME returns .TRUE. if CA is the same letter as CB regardless of
1738 * case.
1739 *
1740 * Arguments
1741 * =========
1742 *
1743 * CA (input) CHARACTER*1
1744 * CB (input) CHARACTER*1
1745 * CA and CB specify the single characters to be compared.
1746 *
1747 * =====================================================================
1748 *
1749 * .. Intrinsic Functions ..
1750  INTRINSIC ichar
1751 * ..
1752 * .. Local Scalars ..
1753  INTEGER inta, intb, zcode
1754 * ..
1755 * .. Executable Statements ..
1756 *
1757 * Test if the characters are equal
1758 *
1759  lsame = ca.EQ.cb
1760  IF( lsame )
1761  $ RETURN
1762 *
1763 * Now test for equivalence if both characters are alphabetic.
1764 *
1765  zcode = ichar( 'Z' )
1766 *
1767 * Use 'Z' rather than 'A' so that ASCII can be detected on Prime
1768 * machines, on which ICHAR returns a value with bit 8 set.
1769 * ICHAR('A') on Prime machines returns 193 which is the same as
1770 * ICHAR('A') on an EBCDIC machine.
1771 *
1772  inta = ichar( ca )
1773  intb = ichar( cb )
1774 *
1775  IF( zcode.EQ.90 .OR. zcode.EQ.122 ) THEN
1776 *
1777 * ASCII is assumed - ZCODE is the ASCII code of either lower or
1778 * upper case 'Z'.
1779 *
1780  IF( inta.GE.97 .AND. inta.LE.122 ) inta = inta - 32
1781  IF( intb.GE.97 .AND. intb.LE.122 ) intb = intb - 32
1782 *
1783  ELSE IF( zcode.EQ.233 .OR. zcode.EQ.169 ) THEN
1784 *
1785 * EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
1786 * upper case 'Z'.
1787 *
1788  IF( inta.GE.129 .AND. inta.LE.137 .OR.
1789  $ inta.GE.145 .AND. inta.LE.153 .OR.
1790  $ inta.GE.162 .AND. inta.LE.169 ) inta = inta + 64
1791  IF( intb.GE.129 .AND. intb.LE.137 .OR.
1792  $ intb.GE.145 .AND. intb.LE.153 .OR.
1793  $ intb.GE.162 .AND. intb.LE.169 ) intb = intb + 64
1794 *
1795  ELSE IF( zcode.EQ.218 .OR. zcode.EQ.250 ) THEN
1796 *
1797 * ASCII is assumed, on Prime machines - ZCODE is the ASCII code
1798 * plus 128 of either lower or upper case 'Z'.
1799 *
1800  IF( inta.GE.225 .AND. inta.LE.250 ) inta = inta - 32
1801  IF( intb.GE.225 .AND. intb.LE.250 ) intb = intb - 32
1802  END IF
1803  lsame = inta.EQ.intb
1804 *
1805 * RETURN
1806 *
1807 * End of LSAME
1808 *
1809  END
1810  DOUBLE PRECISION FUNCTION dlarnd( IDIST, ISEED )
1811 *
1812 * -- LAPACK auxiliary routine (version 2.0) --
1813 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1814 * Courant Institute, Argonne National Lab, and Rice University
1815 * June 30, 1994
1816 *
1817 * .. Scalar Arguments ..
1818  INTEGER idist
1819 * ..
1820 * .. Array Arguments ..
1821  INTEGER iseed( 4 )
1822 * ..
1823 *
1824 * Purpose
1825 * =======
1826 *
1827 * DLARND returns a random real number from a uniform or normal
1828 * distribution.
1829 *
1830 * Arguments
1831 * =========
1832 *
1833 * IDIST (input) INTEGER
1834 * Specifies the distribution of the random numbers:
1835 * = 1: uniform (0,1)
1836 * = 2: uniform (-1,1)
1837 * = 3: normal (0,1)
1838 *
1839 * ISEED (input/output) INTEGER array, dimension (4)
1840 * On entry, the seed of the random number generator; the array
1841 * elements must be between 0 and 4095, and ISEED(4) must be
1842 * odd.
1843 * On exit, the seed is updated.
1844 *
1845 * Further Details
1846 * ===============
1847 *
1848 * This routine calls the auxiliary routine DLARAN to generate a random
1849 * real number from a uniform (0,1) distribution. The Box-Muller method
1850 * is used to transform numbers from a uniform to a normal distribution.
1851 *
1852 * =====================================================================
1853 *
1854 * .. Parameters ..
1855  DOUBLE PRECISION one, two
1856  parameter( one = 1.0d+0, two = 2.0d+0 )
1857  DOUBLE PRECISION twopi
1858  parameter( twopi = 6.2831853071795864769252867663d+0 )
1859 * ..
1860 * .. Local Scalars ..
1861  DOUBLE PRECISION t1, t2
1862 * ..
1863 * .. External Functions ..
1864  DOUBLE PRECISION dlaran
1865  EXTERNAL dlaran
1866 * ..
1867 * .. Intrinsic Functions ..
1868  INTRINSIC cos, log, sqrt
1869 * ..
1870 * .. Executable Statements ..
1871 *
1872 * Generate a real random number from a uniform (0,1) distribution
1873 *
1874  t1 = dlaran( iseed )
1875 *
1876  IF( idist.EQ.1 ) THEN
1877 *
1878 * uniform (0,1)
1879 *
1880  dlarnd = t1
1881  ELSE IF( idist.EQ.2 ) THEN
1882 *
1883 * uniform (-1,1)
1884 *
1885  dlarnd = two*t1 - one
1886  ELSE IF( idist.EQ.3 ) THEN
1887 *
1888 * normal (0,1)
1889 *
1890  t2 = dlaran( iseed )
1891  dlarnd = sqrt( -two*log( t1 ) )*cos( twopi*t2 )
1892  END IF
1893  RETURN
1894 *
1895 * End of DLARND
1896 *
1897  END
1898  DOUBLE COMPLEX FUNCTION zlarnd( IDIST, ISEED )
1899 *
1900 * -- LAPACK auxiliary routine (version 2.0) --
1901 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1902 * Courant Institute, Argonne National Lab, and Rice University
1903 * June 30, 1994
1904 *
1905 * .. Scalar Arguments ..
1906  INTEGER idist
1907 * ..
1908 * .. Array Arguments ..
1909  INTEGER iseed( 4 )
1910 * ..
1911 *
1912 * Purpose
1913 * =======
1914 *
1915 * ZLARND returns a random complex number from a uniform or normal
1916 * distribution.
1917 *
1918 * Arguments
1919 * =========
1920 *
1921 * IDIST (input) INTEGER
1922 * Specifies the distribution of the random numbers:
1923 * = 1: real and imaginary parts each uniform (0,1)
1924 * = 2: real and imaginary parts each uniform (-1,1)
1925 * = 3: real and imaginary parts each normal (0,1)
1926 * = 4: uniformly distributed on the disc abs(z) <= 1
1927 * = 5: uniformly distributed on the circle abs(z) = 1
1928 *
1929 * ISEED (input/output) INTEGER array, dimension (4)
1930 * On entry, the seed of the random number generator; the array
1931 * elements must be between 0 and 4095, and ISEED(4) must be
1932 * odd.
1933 * On exit, the seed is updated.
1934 *
1935 * Further Details
1936 * ===============
1937 *
1938 * This routine calls the auxiliary routine DLARAN to generate a random
1939 * real number from a uniform (0,1) distribution. The Box-Muller method
1940 * is used to transform numbers from a uniform to a normal distribution.
1941 *
1942 * =====================================================================
1943 *
1944 * .. Parameters ..
1945  DOUBLE PRECISION zero, one, two
1946  parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
1947  DOUBLE PRECISION twopi
1948  parameter( twopi = 6.2831853071795864769252867663d+0 )
1949 * ..
1950 * .. Local Scalars ..
1951  DOUBLE PRECISION t1, t2
1952 * ..
1953 * .. External Functions ..
1954  DOUBLE PRECISION dlaran
1955  EXTERNAL dlaran
1956 * ..
1957 * .. Intrinsic Functions ..
1958  INTRINSIC dcmplx, exp, log, sqrt
1959 * ..
1960 * .. Executable Statements ..
1961 *
1962 * Generate a pair of real random numbers from a uniform (0,1)
1963 * distribution
1964 *
1965  t1 = dlaran( iseed )
1966  t2 = dlaran( iseed )
1967 *
1968  IF( idist.EQ.1 ) THEN
1969 *
1970 * real and imaginary parts each uniform (0,1)
1971 *
1972  zlarnd = dcmplx( t1, t2 )
1973  ELSE IF( idist.EQ.2 ) THEN
1974 *
1975 * real and imaginary parts each uniform (-1,1)
1976 *
1977  zlarnd = dcmplx( two*t1-one, two*t2-one )
1978  ELSE IF( idist.EQ.3 ) THEN
1979 *
1980 * real and imaginary parts each normal (0,1)
1981 *
1982  zlarnd = sqrt( -two*log( t1 ) )*exp( dcmplx( zero, twopi*t2 ) )
1983  ELSE IF( idist.EQ.4 ) THEN
1984 *
1985 * uniform distribution on the unit disc abs(z) <= 1
1986 *
1987  zlarnd = sqrt( t1 )*exp( dcmplx( zero, twopi*t2 ) )
1988  ELSE IF( idist.EQ.5 ) THEN
1989 *
1990 * uniform distribution on the unit circle abs(z) = 1
1991 *
1992  zlarnd = exp( dcmplx( zero, twopi*t2 ) )
1993  END IF
1994  RETURN
1995 *
1996 * End of ZLARND
1997 *
1998  END
1999  DOUBLE PRECISION FUNCTION dlaran( ISEED )
2000 *
2001 * -- LAPACK auxiliary routine (version 2.0) --
2002 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
2003 * Courant Institute, Argonne National Lab, and Rice University
2004 * February 29, 1992
2005 *
2006 * .. Array Arguments ..
2007  INTEGER iseed( 4 )
2008 * ..
2009 *
2010 * Purpose
2011 * =======
2012 *
2013 * DLARAN returns a random real number from a uniform (0,1)
2014 * distribution.
2015 *
2016 * Arguments
2017 * =========
2018 *
2019 * ISEED (input/output) INTEGER array, dimension (4)
2020 * On entry, the seed of the random number generator; the array
2021 * elements must be between 0 and 4095, and ISEED(4) must be
2022 * odd.
2023 * On exit, the seed is updated.
2024 *
2025 * Further Details
2026 * ===============
2027 *
2028 * This routine uses a multiplicative congruential method with modulus
2029 * 2**48 and multiplier 33952834046453 (see G.S.Fishman,
2030 * 'Multiplicative congruential random number generators with modulus
2031 * 2**b: an exhaustive analysis for b = 32 and a partial analysis for
2032 * b = 48', Math. Comp. 189, pp 331-344, 1990).
2033 *
2034 * 48-bit integers are stored in 4 integer array elements with 12 bits
2035 * per element. Hence the routine is portable across machines with
2036 * integers of 32 bits or more.
2037 *
2038 * =====================================================================
2039 *
2040 * .. Parameters ..
2041  INTEGER m1, m2, m3, m4
2042  parameter( m1 = 494, m2 = 322, m3 = 2508, m4 = 2549 )
2043  DOUBLE PRECISION one
2044  parameter( one = 1.0d+0 )
2045  INTEGER ipw2
2046  DOUBLE PRECISION r
2047  parameter( ipw2 = 4096, r = one / ipw2 )
2048 * ..
2049 * .. Local Scalars ..
2050  INTEGER it1, it2, it3, it4
2051 * ..
2052 * .. Intrinsic Functions ..
2053  INTRINSIC dble, mod
2054 * ..
2055 * .. Executable Statements ..
2056 *
2057 * multiply the seed by the multiplier modulo 2**48
2058 *
2059  it4 = iseed( 4 )*m4
2060  it3 = it4 / ipw2
2061  it4 = it4 - ipw2*it3
2062  it3 = it3 + iseed( 3 )*m4 + iseed( 4 )*m3
2063  it2 = it3 / ipw2
2064  it3 = it3 - ipw2*it2
2065  it2 = it2 + iseed( 2 )*m4 + iseed( 3 )*m3 + iseed( 4 )*m2
2066  it1 = it2 / ipw2
2067  it2 = it2 - ipw2*it1
2068  it1 = it1 + iseed( 1 )*m4 + iseed( 2 )*m3 + iseed( 3 )*m2 +
2069  $ iseed( 4 )*m1
2070  it1 = mod( it1, ipw2 )
2071 *
2072 * return updated seed
2073 *
2074  iseed( 1 ) = it1
2075  iseed( 2 ) = it2
2076  iseed( 3 ) = it3
2077  iseed( 4 ) = it4
2078 *
2079 * convert 48-bit integer to a real number in the interval (0,1)
2080 *
2081  dlaran = r*( dble( it1 )+r*( dble( it2 )+r*( dble( it3 )+r*
2082  $ ( dble( it4 ) ) ) ) )
2083  RETURN
2084 *
2085 * End of DLARAN
2086 *
2087  END
max
#define max(A, B)
Definition: pcgemr.c:180
slamc4
subroutine slamc4(EMIN, START, BASE)
Definition: tools.f:1481
slamc3
real function slamc3(A, B)
Definition: tools.f:1443
lsame
logical function lsame(CA, CB)
Definition: tools.f:1724
dlamc2
subroutine dlamc2(BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX)
Definition: tools.f:327
dlamc5
subroutine dlamc5(BETA, P, EMIN, IEEE, EMAX, RMAX)
Definition: tools.f:708
zlarnd
double complex function zlarnd(IDIST, ISEED)
Definition: tools.f:1899
dlamch
double precision function dlamch(CMACH)
Definition: tools.f:10
slamc5
subroutine slamc5(BETA, P, EMIN, IEEE, EMAX, RMAX)
Definition: tools.f:1565
slamch
real function slamch(CMACH)
Definition: tools.f:867
slamc2
subroutine slamc2(BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX)
Definition: tools.f:1184
slamc1
subroutine slamc1(BETA, T, RND, IEEE1)
Definition: tools.f:997
dlarnd
double precision function dlarnd(IDIST, ISEED)
Definition: tools.f:1811
dlaran
double precision function dlaran(ISEED)
Definition: tools.f:2000
dlamc1
subroutine dlamc1(BETA, T, RND, IEEE1)
Definition: tools.f:140
dlamc4
subroutine dlamc4(EMIN, START, BASE)
Definition: tools.f:624
dlamc3
double precision function dlamc3(A, B)
Definition: tools.f:586
min
#define min(A, B)
Definition: pcgemr.c:181