5 /* *********************************************************************** */
7 doublereal dlamc3_(doublereal *a, doublereal *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 /* DLAMC3 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. */
30 /* A (input) DOUBLE PRECISION */
31 /* B (input) DOUBLE PRECISION */
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 dlamch_(char *cmach) {
59 double eps = DBL_EPSILON, 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 doublereal c_b32 = 0.;
98 doublereal dlamch_(char *cmach)
100 /* Initialized data */
102 static logical first = TRUE_;
104 /* System generated locals */
108 /* Builtin functions */
109 double pow_di(doublereal *, integer *);
111 /* Local variables */
114 static doublereal rnd, eps, base;
116 static doublereal emin, prec, emax;
119 static doublereal rmin, rmax;
121 extern logical lsame_(char *, char *);
123 static doublereal sfmin;
124 extern /* Subroutine */ int dlamc2_(integer *, integer *, logical *,
125 doublereal *, integer *, doublereal *, integer *, doublereal *);
128 /* -- LAPACK auxiliary routine (version 3.1) -- */
129 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
132 /* .. Scalar Arguments .. */
138 /* DLAMCH determines double precision machine parameters. */
143 /* CMACH (input) CHARACTER*1 */
144 /* Specifies the value to be returned by DLAMCH: */
145 /* = 'E' or 'e', DLAMCH := eps */
146 /* = 'S' or 's , DLAMCH := sfmin */
147 /* = 'B' or 'b', DLAMCH := base */
148 /* = 'P' or 'p', DLAMCH := eps*base */
149 /* = 'N' or 'n', DLAMCH := t */
150 /* = 'R' or 'r', DLAMCH := rnd */
151 /* = 'M' or 'm', DLAMCH := emin */
152 /* = 'U' or 'u', DLAMCH := rmin */
153 /* = 'L' or 'l', DLAMCH := emax */
154 /* = 'O' or 'o', DLAMCH := 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 dlamc2_(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
187 base = (doublereal) beta;
192 eps = pow_di(&base, &i__1) / 2;
196 eps = pow_di(&base, &i__1);
199 emin = (doublereal) imin;
200 emax = (doublereal) imax;
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.);
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 dlamc1_(integer *beta, integer *t, logical *rnd, logical
248 /* Initialized data */
250 static logical first = TRUE_;
252 /* System generated locals */
253 doublereal d__1, d__2;
255 /* Local variables */
256 doublereal a, b, c__, f, t1, t2;
260 static integer lbeta;
262 extern doublereal dlamc3_(doublereal *, doublereal *);
263 static logical lieee1;
266 /* -- LAPACK auxiliary routine (version 3.1) -- */
267 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
270 /* .. Scalar Arguments .. */
276 /* DLAMC1 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 DLAMC3 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__ = dlamc3_(&a, &one);
347 c__ = dlamc3_(&c__, &d__1);
352 /* Now compute b = 2.0**m with the smallest positive integer m */
355 /* fl( a + b ) .gt. a. */
358 c__ = dlamc3_(&a, &b);
360 /* + WHILE( C.EQ.A )LOOP */
364 c__ = dlamc3_(&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__ = dlamc3_(&c__, &d__1);
378 lbeta = (integer) (c__ + qtr);
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. */
383 b = (doublereal) lbeta;
386 f = dlamc3_(&d__1, &d__2);
387 c__ = dlamc3_(&f, &a);
395 f = dlamc3_(&d__1, &d__2);
396 c__ = dlamc3_(&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 = dlamc3_(&d__1, &a);
410 t2 = dlamc3_(&d__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__ = dlamc3_(&a, &one);
431 c__ = dlamc3_(&c__, &d__1);
450 /* *********************************************************************** */
452 /* Subroutine */ int dlamc2_(integer *beta, integer *t, logical *rnd,
453 doublereal *eps, integer *emin, doublereal *rmin, integer *emax,
456 /* Initialized data */
458 static logical first = TRUE_;
459 static logical iwarn = FALSE_;
462 static char fmt_9999[] = "(//\002 WARNING. The value EMIN may be incorre"
463 "ct:-\002,\002 EMIN = \002,i8,/\002 If, after inspection, the va"
464 "lue EMIN looks\002,\002 acceptable please comment out \002,/\002"
465 " the IF block as marked within the code of routine\002,\002 DLAM"
466 "C2,\002,/\002 otherwise supply EMIN explicitly.\002,/)";
468 /* System generated locals */
470 doublereal d__1, d__2, d__3, d__4, d__5;
472 /* Builtin functions */
473 double pow_di(doublereal *, integer *);
474 //integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
476 /* Local variables */
477 doublereal a, b, c__;
484 static doublereal leps;
486 static integer lbeta;
488 static integer lemin, lemax;
493 static doublereal lrmin, lrmax;
495 extern /* Subroutine */ int dlamc1_(integer *, integer *, logical *,
497 extern doublereal dlamc3_(doublereal *, doublereal *);
499 extern /* Subroutine */ int dlamc4_(integer *, doublereal *, integer *),
500 dlamc5_(integer *, integer *, integer *, logical *, integer *,
502 integer ngnmin, ngpmin;
504 /* Fortran I/O blocks */
505 static cilist io___58 = { 0, 6, 0, fmt_9999, 0 };
509 /* -- LAPACK auxiliary routine (version 3.1) -- */
510 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
513 /* .. Scalar Arguments .. */
519 /* DLAMC2 determines the machine parameters specified in its argument */
525 /* BETA (output) INTEGER */
526 /* The base of the machine. */
528 /* T (output) INTEGER */
529 /* The number of ( BETA ) digits in the mantissa. */
531 /* RND (output) LOGICAL */
532 /* Specifies whether proper rounding ( RND = .TRUE. ) or */
533 /* chopping ( RND = .FALSE. ) occurs in addition. This may not */
534 /* be a reliable guide to the way in which the machine performs */
535 /* its arithmetic. */
537 /* EPS (output) DOUBLE PRECISION */
538 /* The smallest positive number such that */
540 /* fl( 1.0 - EPS ) .LT. 1.0, */
542 /* where fl denotes the computed value. */
544 /* EMIN (output) INTEGER */
545 /* The minimum exponent before (gradual) underflow occurs. */
547 /* RMIN (output) DOUBLE PRECISION */
548 /* The smallest normalized number for the machine, given by */
549 /* BASE**( EMIN - 1 ), where BASE is the floating point value */
552 /* EMAX (output) INTEGER */
553 /* The maximum exponent before overflow occurs. */
555 /* RMAX (output) DOUBLE PRECISION */
556 /* The largest positive number for the machine, given by */
557 /* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point */
560 /* Further Details */
561 /* =============== */
563 /* The computation of EPS is based on a routine PARANOIA by */
564 /* W. Kahan of the University of California at Berkeley. */
566 /* ===================================================================== */
568 /* .. Local Scalars .. */
570 /* .. External Functions .. */
572 /* .. External Subroutines .. */
574 /* .. Intrinsic Functions .. */
576 /* .. Save statement .. */
578 /* .. Data statements .. */
580 /* .. Executable Statements .. */
587 /* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of */
588 /* BETA, T, RND, EPS, EMIN and RMIN. */
590 /* Throughout this routine we use the function DLAMC3 to ensure */
591 /* that relevant values are stored and not held in registers, or */
592 /* are not affected by optimizers. */
594 /* DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. */
596 dlamc1_(&lbeta, <, &lrnd, &lieee1);
598 /* Start to find EPS. */
600 b = (doublereal) lbeta;
602 a = pow_di(&b, &i__1);
605 /* Try some tricks to see whether or not this is the correct EPS. */
610 sixth = dlamc3_(&b, &d__1);
611 third = dlamc3_(&sixth, &sixth);
613 b = dlamc3_(&third, &d__1);
614 b = dlamc3_(&b, &sixth);
622 /* + WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */
624 if (leps > b && b > zero) {
627 /* Computing 5th power */
628 d__3 = two, d__4 = d__3, d__3 *= d__3;
629 /* Computing 2nd power */
631 d__2 = d__4 * (d__3 * d__3) * (d__5 * d__5);
632 c__ = dlamc3_(&d__1, &d__2);
634 c__ = dlamc3_(&half, &d__1);
635 b = dlamc3_(&half, &c__);
637 c__ = dlamc3_(&half, &d__1);
638 b = dlamc3_(&half, &c__);
647 /* Computation of EPS complete. */
649 /* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). */
650 /* Keep dividing A by BETA until (gradual) underflow occurs. This */
651 /* is detected when we cannot recover the previous A. */
655 for (i__ = 1; i__ <= 3; ++i__) {
656 d__1 = small * rbase;
657 small = dlamc3_(&d__1, &zero);
660 a = dlamc3_(&one, &small);
661 dlamc4_(&ngpmin, &one, &lbeta);
663 dlamc4_(&ngnmin, &d__1, &lbeta);
664 dlamc4_(&gpmin, &a, &lbeta);
666 dlamc4_(&gnmin, &d__1, &lbeta);
669 if (ngpmin == ngnmin && gpmin == gnmin) {
670 if (ngpmin == gpmin) {
672 /* ( Non twos-complement machines, no gradual underflow; */
674 } else if (gpmin - ngpmin == 3) {
675 lemin = ngpmin - 1 + lt;
677 /* ( Non twos-complement machines, with gradual underflow; */
678 /* e.g., IEEE standard followers ) */
680 lemin = min(ngpmin,gpmin);
681 /* ( A guess; no known machine ) */
685 } else if (ngpmin == gpmin && ngnmin == gnmin) {
686 if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1) {
687 lemin = max(ngpmin,ngnmin);
688 /* ( Twos-complement machines, no gradual underflow; */
689 /* e.g., CYBER 205 ) */
691 lemin = min(ngpmin,ngnmin);
692 /* ( A guess; no known machine ) */
696 } else if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1 && gpmin == gnmin)
698 if (gpmin - min(ngpmin,ngnmin) == 3) {
699 lemin = max(ngpmin,ngnmin) - 1 + lt;
700 /* ( Twos-complement machines with gradual underflow; */
701 /* no known machine ) */
703 lemin = min(ngpmin,ngnmin);
704 /* ( A guess; no known machine ) */
710 i__1 = min(ngpmin,ngnmin), i__1 = min(i__1,gpmin);
711 lemin = min(i__1,gnmin);
712 /* ( A guess; no known machine ) */
717 /* Comment out this if block if EMIN is ok */
720 printf("\n\n WARNING. The value EMIN may be incorrect:- ");
721 printf("EMIN = %8i\n",lemin);
722 printf("If, after inspection, the value EMIN looks acceptable");
723 printf("please comment out \n the IF block as marked within the");
724 printf("code of routine DLAMC2, \n otherwise supply EMIN");
725 printf("explicitly.\n");
728 do_fio(&c__1, (char *)&lemin, (ftnlen)sizeof(integer));
734 /* Assume IEEE arithmetic if we found denormalised numbers above, */
735 /* or if arithmetic seems to round in the IEEE style, determined */
736 /* in routine DLAMC1. A true IEEE machine should have both things */
737 /* true; however, faulty machines may have one or the other. */
739 ieee = ieee || lieee1;
741 /* Compute RMIN by successive division by BETA. We could compute */
742 /* RMIN as BASE**( EMIN - 1 ), but some machines underflow during */
743 /* this computation. */
747 for (i__ = 1; i__ <= i__1; ++i__) {
748 d__1 = lrmin * rbase;
749 lrmin = dlamc3_(&d__1, &zero);
753 /* Finally, call DLAMC5 to compute EMAX and RMAX. */
755 dlamc5_(&lbeta, <, &lemin, &ieee, &lemax, &lrmax);
776 /* *********************************************************************** */
778 /* Subroutine */ int dlamc4_(integer *emin, doublereal *start, integer *base)
780 /* System generated locals */
784 /* Local variables */
787 doublereal b1, b2, c1, c2, d1, d2, one, zero, rbase;
788 extern doublereal dlamc3_(doublereal *, doublereal *);
791 /* -- LAPACK auxiliary routine (version 3.1) -- */
792 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
795 /* .. Scalar Arguments .. */
801 /* DLAMC4 is a service routine for DLAMC2. */
806 /* EMIN (output) INTEGER */
807 /* The minimum exponent before (gradual) underflow, computed by */
808 /* setting A = START and dividing by BASE until the previous A */
809 /* can not be recovered. */
811 /* START (input) DOUBLE PRECISION */
812 /* The starting point for determining EMIN. */
814 /* BASE (input) INTEGER */
815 /* The base of the machine. */
817 /* ===================================================================== */
819 /* .. Local Scalars .. */
821 /* .. External Functions .. */
823 /* .. Executable Statements .. */
831 b1 = dlamc3_(&d__1, &zero);
836 /* + WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. */
837 /* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP */
839 if (c1 == a && c2 == a && d1 == a && d2 == a) {
843 b1 = dlamc3_(&d__1, &zero);
845 c1 = dlamc3_(&d__1, &zero);
848 for (i__ = 1; i__ <= i__1; ++i__) {
853 b2 = dlamc3_(&d__1, &zero);
855 c2 = dlamc3_(&d__1, &zero);
858 for (i__ = 1; i__ <= i__1; ++i__) {
873 /* *********************************************************************** */
875 /* Subroutine */ int dlamc5_(integer *beta, integer *p, integer *emin,
876 logical *ieee, integer *emax, doublereal *rmax)
878 /* System generated locals */
882 /* Local variables */
888 extern doublereal dlamc3_(doublereal *, doublereal *);
890 integer exbits, expsum;
893 /* -- LAPACK auxiliary routine (version 3.1) -- */
894 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
897 /* .. Scalar Arguments .. */
903 /* DLAMC5 attempts to compute RMAX, the largest machine floating-point */
904 /* number, without overflow. It assumes that EMAX + abs(EMIN) sum */
905 /* approximately to a power of 2. It will fail on machines where this */
906 /* assumption does not hold, for example, the Cyber 205 (EMIN = -28625, */
907 /* EMAX = 28718). It will also fail if the value supplied for EMIN is */
908 /* too large (i.e. too close to zero), probably with overflow. */
913 /* BETA (input) INTEGER */
914 /* The base of floating-point arithmetic. */
916 /* P (input) INTEGER */
917 /* The number of base BETA digits in the mantissa of a */
918 /* floating-point value. */
920 /* EMIN (input) INTEGER */
921 /* The minimum exponent before (gradual) underflow. */
923 /* IEEE (input) LOGICAL */
924 /* A logical flag specifying whether or not the arithmetic */
925 /* system is thought to comply with the IEEE standard. */
927 /* EMAX (output) INTEGER */
928 /* The largest exponent before overflow */
930 /* RMAX (output) DOUBLE PRECISION */
931 /* The largest machine floating-point number. */
933 /* ===================================================================== */
935 /* .. Parameters .. */
937 /* .. Local Scalars .. */
939 /* .. External Functions .. */
941 /* .. Intrinsic Functions .. */
943 /* .. Executable Statements .. */
945 /* First compute LEXP and UEXP, two powers of 2 that bound */
946 /* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum */
947 /* approximately to the bound that is closest to abs(EMIN). */
948 /* (EMAX is the exponent of the required number RMAX). */
954 if (try__ <= -(*emin)) {
959 if (lexp == -(*emin)) {
966 /* Now -LEXP is less than or equal to EMIN, and -UEXP is greater */
967 /* than or equal to EMIN. EXBITS is the number of bits needed to */
968 /* store the exponent. */
970 if (uexp + *emin > -lexp - *emin) {
976 /* EXPSUM is the exponent range, approximately equal to */
977 /* EMAX - EMIN + 1 . */
979 *emax = expsum + *emin - 1;
980 nbits = exbits + 1 + *p;
982 /* NBITS is the total number of bits needed to store a */
983 /* floating-point number. */
985 if (nbits % 2 == 1 && *beta == 2) {
987 /* Either there are an odd number of bits used to store a */
988 /* floating-point number, which is unlikely, or some bits are */
989 /* not used in the representation of numbers, which is possible, */
990 /* (e.g. Cray machines) or the mantissa has an implicit bit, */
991 /* (e.g. IEEE machines, Dec Vax machines), which is perhaps the */
992 /* most likely. We have to assume the last alternative. */
993 /* If this is true, then we need to reduce EMAX by one because */
994 /* there must be some way of representing zero in an implicit-bit */
995 /* system. On machines like Cray, we are reducing EMAX by one */
1003 /* Assume we are on an IEEE machine which reserves one exponent */
1004 /* for infinity and NaN. */
1009 /* Now create RMAX, the largest machine number, which should */
1010 /* be equal to (1.0 - BETA**(-P)) * BETA**EMAX . */
1012 /* First compute 1.0 - BETA**(-P), being careful that the */
1013 /* result is less than 1.0 . */
1015 recbas = 1. / *beta;
1019 for (i__ = 1; i__ <= i__1; ++i__) {
1024 y = dlamc3_(&y, &z__);
1031 /* Now multiply by BETA**EMAX to get RMAX. */
1034 for (i__ = 1; i__ <= i__1; ++i__) {
1036 y = dlamc3_(&d__1, &c_b32);