5 /* *********************************************************************** */
7 doublereal slamc3_(real *a, real *b)
9 /* System generated locals */
13 /* -- LAPACK auxiliary routine (version 3.1) -- */
14 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
17 /* .. Scalar Arguments .. */
23 /* SLAMC3 is intended to force A and B to be stored prior to doing */
24 /* the addition of A and B , for use in situations where optimizers */
25 /* might hold one of these in a register. */
32 /* The values A and B. */
34 /* ===================================================================== */
36 /* .. Executable Statements .. */
49 /* simpler version of dlamch for the case of IEEE754-compliant FPU module by Piotr Luszczek S.
50 taken from http://www.mail-archive.com/numpy-discussion@lists.sourceforge.net/msg02448.html */
57 slamch_(char *cmach) {
59 float eps=FLT_EPSILON*0.5f, sfmin, small;
61 if ('B' == ch || 'b' == ch) {
63 } else if ('E' == ch || 'e' == ch) {
65 } else if ('L' == ch || 'l' == ch) {
67 } else if ('M' == ch || 'm' == ch) {
69 } else if ('N' == ch || 'n' == ch) {
71 } else if ('O' == ch || 'o' == ch) {
73 } else if ('P' == ch || 'p' == ch) {
74 return eps * FLT_RADIX;
75 } else if ('R' == ch || 'r' == ch) {
76 return FLT_ROUNDS < 2;
77 } else if ('S' == ch || 's' == ch) {
78 /* Use SMALL plus a bit, to avoid the possibility of rounding causing overflow
79 when computing 1/sfmin. */
82 if (small <= sfmin) small = sfmin * (1 + eps);
84 } else if ('U' == ch || 'u' == ch) {
93 /* Table of constant values */
95 static integer c__1 = 1;
96 static real c_b32 = 0.f;
98 doublereal slamch_(char *cmach)
100 /* Initialized data */
102 static logical first = TRUE_;
104 /* System generated locals */
108 /* Builtin functions */
109 double pow_ri(real *, integer *);
111 /* Local variables */
114 static real rnd, eps, base;
116 static real emin, prec, emax;
119 static real rmin, rmax;
121 extern logical lsame_(char *, char *);
124 extern /* Subroutine */ int slamc2_(integer *, integer *, logical *, real
125 *, integer *, real *, integer *, real *);
128 /* -- LAPACK auxiliary routine (version 3.1) -- */
129 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
132 /* .. Scalar Arguments .. */
138 /* SLAMCH determines single precision machine parameters. */
143 /* CMACH (input) CHARACTER*1 */
144 /* Specifies the value to be returned by SLAMCH: */
145 /* = 'E' or 'e', SLAMCH := eps */
146 /* = 'S' or 's , SLAMCH := sfmin */
147 /* = 'B' or 'b', SLAMCH := base */
148 /* = 'P' or 'p', SLAMCH := eps*base */
149 /* = 'N' or 'n', SLAMCH := t */
150 /* = 'R' or 'r', SLAMCH := rnd */
151 /* = 'M' or 'm', SLAMCH := emin */
152 /* = 'U' or 'u', SLAMCH := rmin */
153 /* = 'L' or 'l', SLAMCH := emax */
154 /* = 'O' or 'o', SLAMCH := rmax */
158 /* eps = relative machine precision */
159 /* sfmin = safe minimum, such that 1/sfmin does not overflow */
160 /* base = base of the machine */
161 /* prec = eps*base */
162 /* t = number of (base) digits in the mantissa */
163 /* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise */
164 /* emin = minimum exponent before (gradual) underflow */
165 /* rmin = underflow threshold - base**(emin-1) */
166 /* emax = largest exponent before overflow */
167 /* rmax = overflow threshold - (base**emax)*(1-eps) */
169 /* ===================================================================== */
171 /* .. Parameters .. */
173 /* .. Local Scalars .. */
175 /* .. External Functions .. */
177 /* .. External Subroutines .. */
179 /* .. Save statement .. */
181 /* .. Data statements .. */
183 /* .. Executable Statements .. */
186 slamc2_(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
192 eps = pow_ri(&base, &i__1) / 2;
196 eps = pow_ri(&base, &i__1);
203 if (small >= sfmin) {
205 /* Use SMALL plus a bit, to avoid the possibility of rounding */
206 /* causing overflow when computing 1/sfmin. */
208 sfmin = small * (eps + 1.f);
212 if (lsame_(cmach, "E")) {
214 } else if (lsame_(cmach, "S")) {
216 } else if (lsame_(cmach, "B")) {
218 } else if (lsame_(cmach, "P")) {
220 } else if (lsame_(cmach, "N")) {
222 } else if (lsame_(cmach, "R")) {
224 } else if (lsame_(cmach, "M")) {
226 } else if (lsame_(cmach, "U")) {
228 } else if (lsame_(cmach, "L")) {
230 } else if (lsame_(cmach, "O")) {
243 /* *********************************************************************** */
245 /* Subroutine */ int slamc1_(integer *beta, integer *t, logical *rnd, logical
248 /* Initialized data */
250 static logical first = TRUE_;
252 /* System generated locals */
255 /* Local variables */
256 real a, b, c__, f, t1, t2;
260 static integer lbeta;
262 static logical lieee1;
263 extern doublereal slamc3_(real *, real *);
266 /* -- LAPACK auxiliary routine (version 3.1) -- */
267 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
270 /* .. Scalar Arguments .. */
276 /* SLAMC1 determines the machine parameters given by BETA, T, RND, and */
282 /* BETA (output) INTEGER */
283 /* The base of the machine. */
285 /* T (output) INTEGER */
286 /* The number of ( BETA ) digits in the mantissa. */
288 /* RND (output) LOGICAL */
289 /* Specifies whether proper rounding ( RND = .TRUE. ) or */
290 /* chopping ( RND = .FALSE. ) occurs in addition. This may not */
291 /* be a reliable guide to the way in which the machine performs */
292 /* its arithmetic. */
294 /* IEEE1 (output) LOGICAL */
295 /* Specifies whether rounding appears to be done in the IEEE */
296 /* 'round to nearest' style. */
298 /* Further Details */
299 /* =============== */
301 /* The routine is based on the routine ENVRON by Malcolm and */
302 /* incorporates suggestions by Gentleman and Marovich. See */
304 /* Malcolm M. A. (1972) Algorithms to reveal properties of */
305 /* floating-point arithmetic. Comms. of the ACM, 15, 949-951. */
307 /* Gentleman W. M. and Marovich S. B. (1974) More on algorithms */
308 /* that reveal properties of floating point arithmetic units. */
309 /* Comms. of the ACM, 17, 276-277. */
311 /* ===================================================================== */
313 /* .. Local Scalars .. */
315 /* .. External Functions .. */
317 /* .. Save statement .. */
319 /* .. Data statements .. */
321 /* .. Executable Statements .. */
326 /* LBETA, LIEEE1, LT and LRND are the local values of BETA, */
327 /* IEEE1, T and RND. */
329 /* Throughout this routine we use the function SLAMC3 to ensure */
330 /* that relevant values are stored and not held in registers, or */
331 /* are not affected by optimizers. */
333 /* Compute a = 2.0**m with the smallest positive integer m such */
336 /* fl( a + 1.0 ) = a. */
341 /* + WHILE( C.EQ.ONE )LOOP */
345 c__ = slamc3_(&a, &one);
347 c__ = slamc3_(&c__, &r__1);
352 /* Now compute b = 2.0**m with the smallest positive integer m */
355 /* fl( a + b ) .gt. a. */
358 c__ = slamc3_(&a, &b);
360 /* + WHILE( C.EQ.A )LOOP */
364 c__ = slamc3_(&a, &b);
369 /* Now compute the base. a and c are neighbouring floating point */
370 /* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so */
371 /* their difference is beta. Adding 0.25 to c is to ensure that it */
372 /* is truncated to beta and not ( beta - 1 ). */
377 c__ = slamc3_(&c__, &r__1);
380 /* Now determine whether rounding or chopping occurs, by adding a */
381 /* bit less than beta/2 and a bit more than beta/2 to a. */
386 f = slamc3_(&r__1, &r__2);
387 c__ = slamc3_(&f, &a);
395 f = slamc3_(&r__1, &r__2);
396 c__ = slamc3_(&f, &a);
397 if (lrnd && c__ == a) {
401 /* Try and decide whether rounding is done in the IEEE 'round to */
402 /* nearest' style. B/2 is half a unit in the last place of the two */
403 /* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit */
404 /* zero, and SAVEC is odd. Thus adding B/2 to A should not change */
405 /* A, but adding B/2 to SAVEC should change SAVEC. */
408 t1 = slamc3_(&r__1, &a);
410 t2 = slamc3_(&r__1, &savec);
411 lieee1 = t1 == a && t2 > savec && lrnd;
413 /* Now find the mantissa, t. It should be the integer part of */
414 /* log to the base beta of a, however it is safer to determine t */
415 /* by powering. So we find t as the smallest positive integer for */
418 /* fl( beta**t + 1.0 ) = 1.0. */
424 /* + WHILE( C.EQ.ONE )LOOP */
429 c__ = slamc3_(&a, &one);
431 c__ = slamc3_(&c__, &r__1);
450 /* *********************************************************************** */
452 /* Subroutine */ int slamc2_(integer *beta, integer *t, logical *rnd, real *
453 eps, integer *emin, real *rmin, integer *emax, real *rmax)
455 /* Initialized data */
457 static logical first = TRUE_;
458 static logical iwarn = FALSE_;
461 static char fmt_9999[] = "(//\002 WARNING. The value EMIN may be incorre"
462 "ct:-\002,\002 EMIN = \002,i8,/\002 If, after inspection, the va"
463 "lue EMIN looks\002,\002 acceptable please comment out \002,/\002"
464 " the IF block as marked within the code of routine\002,\002 SLAM"
465 "C2,\002,/\002 otherwise supply EMIN explicitly.\002,/)";
467 /* System generated locals */
469 real r__1, r__2, r__3, r__4, r__5;
471 /* Builtin functions */
472 double pow_ri(real *, integer *);
473 //integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
475 /* Local variables */
485 static integer lbeta;
487 static integer lemin, lemax;
492 static real lrmin, lrmax;
495 extern /* Subroutine */ int slamc1_(integer *, integer *, logical *,
497 extern doublereal slamc3_(real *, real *);
498 extern /* Subroutine */ int slamc4_(integer *, real *, integer *),
499 slamc5_(integer *, integer *, integer *, logical *, integer *,
501 integer ngnmin, ngpmin;
503 /* Fortran I/O blocks */
504 static cilist io___58 = { 0, 6, 0, fmt_9999, 0 };
508 /* -- LAPACK auxiliary routine (version 3.1) -- */
509 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
512 /* .. Scalar Arguments .. */
518 /* SLAMC2 determines the machine parameters specified in its argument */
524 /* BETA (output) INTEGER */
525 /* The base of the machine. */
527 /* T (output) INTEGER */
528 /* The number of ( BETA ) digits in the mantissa. */
530 /* RND (output) LOGICAL */
531 /* Specifies whether proper rounding ( RND = .TRUE. ) or */
532 /* chopping ( RND = .FALSE. ) occurs in addition. This may not */
533 /* be a reliable guide to the way in which the machine performs */
534 /* its arithmetic. */
536 /* EPS (output) REAL */
537 /* The smallest positive number such that */
539 /* fl( 1.0 - EPS ) .LT. 1.0, */
541 /* where fl denotes the computed value. */
543 /* EMIN (output) INTEGER */
544 /* The minimum exponent before (gradual) underflow occurs. */
546 /* RMIN (output) REAL */
547 /* The smallest normalized number for the machine, given by */
548 /* BASE**( EMIN - 1 ), where BASE is the floating point value */
551 /* EMAX (output) INTEGER */
552 /* The maximum exponent before overflow occurs. */
554 /* RMAX (output) REAL */
555 /* The largest positive number for the machine, given by */
556 /* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point */
559 /* Further Details */
560 /* =============== */
562 /* The computation of EPS is based on a routine PARANOIA by */
563 /* W. Kahan of the University of California at Berkeley. */
565 /* ===================================================================== */
567 /* .. Local Scalars .. */
569 /* .. External Functions .. */
571 /* .. External Subroutines .. */
573 /* .. Intrinsic Functions .. */
575 /* .. Save statement .. */
577 /* .. Data statements .. */
579 /* .. Executable Statements .. */
586 /* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of */
587 /* BETA, T, RND, EPS, EMIN and RMIN. */
589 /* Throughout this routine we use the function SLAMC3 to ensure */
590 /* that relevant values are stored and not held in registers, or */
591 /* are not affected by optimizers. */
593 /* SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. */
595 slamc1_(&lbeta, <, &lrnd, &lieee1);
597 /* Start to find EPS. */
601 a = pow_ri(&b, &i__1);
604 /* Try some tricks to see whether or not this is the correct EPS. */
609 sixth = slamc3_(&b, &r__1);
610 third = slamc3_(&sixth, &sixth);
612 b = slamc3_(&third, &r__1);
613 b = slamc3_(&b, &sixth);
621 /* + WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */
623 if (leps > b && b > zero) {
626 /* Computing 5th power */
627 r__3 = two, r__4 = r__3, r__3 *= r__3;
628 /* Computing 2nd power */
630 r__2 = r__4 * (r__3 * r__3) * (r__5 * r__5);
631 c__ = slamc3_(&r__1, &r__2);
633 c__ = slamc3_(&half, &r__1);
634 b = slamc3_(&half, &c__);
636 c__ = slamc3_(&half, &r__1);
637 b = slamc3_(&half, &c__);
646 /* Computation of EPS complete. */
648 /* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). */
649 /* Keep dividing A by BETA until (gradual) underflow occurs. This */
650 /* is detected when we cannot recover the previous A. */
654 for (i__ = 1; i__ <= 3; ++i__) {
655 r__1 = small * rbase;
656 small = slamc3_(&r__1, &zero);
659 a = slamc3_(&one, &small);
660 slamc4_(&ngpmin, &one, &lbeta);
662 slamc4_(&ngnmin, &r__1, &lbeta);
663 slamc4_(&gpmin, &a, &lbeta);
665 slamc4_(&gnmin, &r__1, &lbeta);
668 if (ngpmin == ngnmin && gpmin == gnmin) {
669 if (ngpmin == gpmin) {
671 /* ( Non twos-complement machines, no gradual underflow; */
673 } else if (gpmin - ngpmin == 3) {
674 lemin = ngpmin - 1 + lt;
676 /* ( Non twos-complement machines, with gradual underflow; */
677 /* e.g., IEEE standard followers ) */
679 lemin = min(ngpmin,gpmin);
680 /* ( A guess; no known machine ) */
684 } else if (ngpmin == gpmin && ngnmin == gnmin) {
685 if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1) {
686 lemin = max(ngpmin,ngnmin);
687 /* ( Twos-complement machines, no gradual underflow; */
688 /* e.g., CYBER 205 ) */
690 lemin = min(ngpmin,ngnmin);
691 /* ( A guess; no known machine ) */
695 } else if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1 && gpmin == gnmin)
697 if (gpmin - min(ngpmin,ngnmin) == 3) {
698 lemin = max(ngpmin,ngnmin) - 1 + lt;
699 /* ( Twos-complement machines with gradual underflow; */
700 /* no known machine ) */
702 lemin = min(ngpmin,ngnmin);
703 /* ( A guess; no known machine ) */
709 i__1 = min(ngpmin,ngnmin), i__1 = min(i__1,gpmin);
710 lemin = min(i__1,gnmin);
711 /* ( A guess; no known machine ) */
716 /* Comment out this if block if EMIN is ok */
719 printf("\n\n WARNING. The value EMIN may be incorrect:- ");
720 printf("EMIN = %8i\n",lemin);
721 printf("If, after inspection, the value EMIN looks acceptable");
722 printf("please comment out \n the IF block as marked within the");
723 printf("code of routine SLAMC2, \n otherwise supply EMIN");
724 printf("explicitly.\n");
727 do_fio(&c__1, (char *)&lemin, (ftnlen)sizeof(integer));
733 /* Assume IEEE arithmetic if we found denormalised numbers above, */
734 /* or if arithmetic seems to round in the IEEE style, determined */
735 /* in routine SLAMC1. A true IEEE machine should have both things */
736 /* true; however, faulty machines may have one or the other. */
738 ieee = ieee || lieee1;
740 /* Compute RMIN by successive division by BETA. We could compute */
741 /* RMIN as BASE**( EMIN - 1 ), but some machines underflow during */
742 /* this computation. */
746 for (i__ = 1; i__ <= i__1; ++i__) {
747 r__1 = lrmin * rbase;
748 lrmin = slamc3_(&r__1, &zero);
752 /* Finally, call SLAMC5 to compute EMAX and RMAX. */
754 slamc5_(&lbeta, <, &lemin, &ieee, &lemax, &lrmax);
774 /* *********************************************************************** */
776 /* Subroutine */ int slamc4_(integer *emin, real *start, integer *base)
778 /* System generated locals */
782 /* Local variables */
785 real b1, b2, c1, c2, d1, d2, one, zero, rbase;
786 extern doublereal slamc3_(real *, real *);
789 /* -- LAPACK auxiliary routine (version 3.1) -- */
790 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
793 /* .. Scalar Arguments .. */
799 /* SLAMC4 is a service routine for SLAMC2. */
804 /* EMIN (output) INTEGER */
805 /* The minimum exponent before (gradual) underflow, computed by */
806 /* setting A = START and dividing by BASE until the previous A */
807 /* can not be recovered. */
809 /* START (input) REAL */
810 /* The starting point for determining EMIN. */
812 /* BASE (input) INTEGER */
813 /* The base of the machine. */
815 /* ===================================================================== */
817 /* .. Local Scalars .. */
819 /* .. External Functions .. */
821 /* .. Executable Statements .. */
829 b1 = slamc3_(&r__1, &zero);
834 /* + WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. */
835 /* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP */
837 if (c1 == a && c2 == a && d1 == a && d2 == a) {
841 b1 = slamc3_(&r__1, &zero);
843 c1 = slamc3_(&r__1, &zero);
846 for (i__ = 1; i__ <= i__1; ++i__) {
851 b2 = slamc3_(&r__1, &zero);
853 c2 = slamc3_(&r__1, &zero);
856 for (i__ = 1; i__ <= i__1; ++i__) {
871 /* *********************************************************************** */
873 /* Subroutine */ int slamc5_(integer *beta, integer *p, integer *emin,
874 logical *ieee, integer *emax, real *rmax)
876 /* System generated locals */
880 /* Local variables */
886 extern doublereal slamc3_(real *, real *);
888 integer exbits, expsum;
891 /* -- LAPACK auxiliary routine (version 3.1) -- */
892 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
895 /* .. Scalar Arguments .. */
901 /* SLAMC5 attempts to compute RMAX, the largest machine floating-point */
902 /* number, without overflow. It assumes that EMAX + abs(EMIN) sum */
903 /* approximately to a power of 2. It will fail on machines where this */
904 /* assumption does not hold, for example, the Cyber 205 (EMIN = -28625, */
905 /* EMAX = 28718). It will also fail if the value supplied for EMIN is */
906 /* too large (i.e. too close to zero), probably with overflow. */
911 /* BETA (input) INTEGER */
912 /* The base of floating-point arithmetic. */
914 /* P (input) INTEGER */
915 /* The number of base BETA digits in the mantissa of a */
916 /* floating-point value. */
918 /* EMIN (input) INTEGER */
919 /* The minimum exponent before (gradual) underflow. */
921 /* IEEE (input) LOGICAL */
922 /* A logical flag specifying whether or not the arithmetic */
923 /* system is thought to comply with the IEEE standard. */
925 /* EMAX (output) INTEGER */
926 /* The largest exponent before overflow */
928 /* RMAX (output) REAL */
929 /* The largest machine floating-point number. */
931 /* ===================================================================== */
933 /* .. Parameters .. */
935 /* .. Local Scalars .. */
937 /* .. External Functions .. */
939 /* .. Intrinsic Functions .. */
941 /* .. Executable Statements .. */
943 /* First compute LEXP and UEXP, two powers of 2 that bound */
944 /* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum */
945 /* approximately to the bound that is closest to abs(EMIN). */
946 /* (EMAX is the exponent of the required number RMAX). */
952 if (try__ <= -(*emin)) {
957 if (lexp == -(*emin)) {
964 /* Now -LEXP is less than or equal to EMIN, and -UEXP is greater */
965 /* than or equal to EMIN. EXBITS is the number of bits needed to */
966 /* store the exponent. */
968 if (uexp + *emin > -lexp - *emin) {
974 /* EXPSUM is the exponent range, approximately equal to */
975 /* EMAX - EMIN + 1 . */
977 *emax = expsum + *emin - 1;
978 nbits = exbits + 1 + *p;
980 /* NBITS is the total number of bits needed to store a */
981 /* floating-point number. */
983 if (nbits % 2 == 1 && *beta == 2) {
985 /* Either there are an odd number of bits used to store a */
986 /* floating-point number, which is unlikely, or some bits are */
987 /* not used in the representation of numbers, which is possible, */
988 /* (e.g. Cray machines) or the mantissa has an implicit bit, */
989 /* (e.g. IEEE machines, Dec Vax machines), which is perhaps the */
990 /* most likely. We have to assume the last alternative. */
991 /* If this is true, then we need to reduce EMAX by one because */
992 /* there must be some way of representing zero in an implicit-bit */
993 /* system. On machines like Cray, we are reducing EMAX by one */
1001 /* Assume we are on an IEEE machine which reserves one exponent */
1002 /* for infinity and NaN. */
1007 /* Now create RMAX, the largest machine number, which should */
1008 /* be equal to (1.0 - BETA**(-P)) * BETA**EMAX . */
1010 /* First compute 1.0 - BETA**(-P), being careful that the */
1011 /* result is less than 1.0 . */
1013 recbas = 1.f / *beta;
1017 for (i__ = 1; i__ <= i__1; ++i__) {
1022 y = slamc3_(&y, &z__);
1029 /* Now multiply by BETA**EMAX to get RMAX. */
1032 for (i__ = 1; i__ <= i__1; ++i__) {
1034 y = slamc3_(&r__1, &c_b32);