--- /dev/null
+#include "clapack.h"
+#include <float.h>
+#include <stdio.h>
+
+/* *********************************************************************** */
+
+doublereal slamc3_(real *a, real *b)
+{
+ /* System generated locals */
+ real ret_val;
+
+
+ /* -- LAPACK auxiliary routine (version 3.1) -- */
+ /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+ /* November 2006 */
+
+ /* .. Scalar Arguments .. */
+ /* .. */
+
+ /* Purpose */
+ /* ======= */
+
+ /* SLAMC3 is intended to force A and B to be stored prior to doing */
+ /* the addition of A and B , for use in situations where optimizers */
+ /* might hold one of these in a register. */
+
+ /* Arguments */
+ /* ========= */
+
+ /* A (input) REAL */
+ /* B (input) REAL */
+ /* The values A and B. */
+
+ /* ===================================================================== */
+
+ /* .. Executable Statements .. */
+
+ ret_val = *a + *b;
+
+ return ret_val;
+
+ /* End of SLAMC3 */
+
+} /* slamc3_ */
+
+
+#if 1
+
+/* simpler version of dlamch for the case of IEEE754-compliant FPU module by Piotr Luszczek S.
+ taken from http://www.mail-archive.com/numpy-discussion@lists.sourceforge.net/msg02448.html */
+
+#ifndef FLT_DIGITS
+#define FLT_DIGITS 24
+#endif
+
+doublereal
+slamch_(char *cmach) {
+ char ch = cmach[0];
+ float eps=FLT_EPSILON*0.5f, sfmin, small;
+
+ if ('B' == ch || 'b' == ch) {
+ return FLT_RADIX;
+ } else if ('E' == ch || 'e' == ch) {
+ return eps;
+ } else if ('L' == ch || 'l' == ch) {
+ return FLT_MAX_EXP;
+ } else if ('M' == ch || 'm' == ch) {
+ return FLT_MIN_EXP;
+ } else if ('N' == ch || 'n' == ch) {
+ return FLT_DIGITS;
+ } else if ('O' == ch || 'o' == ch) {
+ return FLT_MAX;
+ } else if ('P' == ch || 'p' == ch) {
+ return eps * FLT_RADIX;
+ } else if ('R' == ch || 'r' == ch) {
+ return FLT_ROUNDS < 2;
+ } else if ('S' == ch || 's' == ch) {
+ /* Use SMALL plus a bit, to avoid the possibility of rounding causing overflow
+ when computing 1/sfmin. */
+ sfmin = FLT_MIN;
+ small = 2. / FLT_MAX;
+ if (small <= sfmin) small = sfmin * (1 + eps);
+ return small;
+ } else if ('U' == ch || 'u' == ch) {
+ return FLT_MIN;
+ }
+
+ return 0;
+}
+
+#else
+
+/* Table of constant values */
+
+static integer c__1 = 1;
+static real c_b32 = 0.f;
+
+doublereal slamch_(char *cmach)
+{
+ /* Initialized data */
+
+ static logical first = TRUE_;
+
+ /* System generated locals */
+ integer i__1;
+ real ret_val;
+
+ /* Builtin functions */
+ double pow_ri(real *, integer *);
+
+ /* Local variables */
+ static real t;
+ integer it;
+ static real rnd, eps, base;
+ integer beta;
+ static real emin, prec, emax;
+ integer imin, imax;
+ logical lrnd;
+ static real rmin, rmax;
+ real rmach;
+ extern logical lsame_(char *, char *);
+ real small;
+ static real sfmin;
+ extern /* Subroutine */ int slamc2_(integer *, integer *, logical *, real
+ *, integer *, real *, integer *, real *);
+
+
+/* -- LAPACK auxiliary routine (version 3.1) -- */
+/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+/* November 2006 */
+
+/* .. Scalar Arguments .. */
+/* .. */
+
+/* Purpose */
+/* ======= */
+
+/* SLAMCH determines single precision machine parameters. */
+
+/* Arguments */
+/* ========= */
+
+/* CMACH (input) CHARACTER*1 */
+/* Specifies the value to be returned by SLAMCH: */
+/* = 'E' or 'e', SLAMCH := eps */
+/* = 'S' or 's , SLAMCH := sfmin */
+/* = 'B' or 'b', SLAMCH := base */
+/* = 'P' or 'p', SLAMCH := eps*base */
+/* = 'N' or 'n', SLAMCH := t */
+/* = 'R' or 'r', SLAMCH := rnd */
+/* = 'M' or 'm', SLAMCH := emin */
+/* = 'U' or 'u', SLAMCH := rmin */
+/* = 'L' or 'l', SLAMCH := emax */
+/* = 'O' or 'o', SLAMCH := rmax */
+
+/* where */
+
+/* eps = relative machine precision */
+/* sfmin = safe minimum, such that 1/sfmin does not overflow */
+/* base = base of the machine */
+/* prec = eps*base */
+/* t = number of (base) digits in the mantissa */
+/* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise */
+/* emin = minimum exponent before (gradual) underflow */
+/* rmin = underflow threshold - base**(emin-1) */
+/* emax = largest exponent before overflow */
+/* rmax = overflow threshold - (base**emax)*(1-eps) */
+
+/* ===================================================================== */
+
+/* .. Parameters .. */
+/* .. */
+/* .. Local Scalars .. */
+/* .. */
+/* .. External Functions .. */
+/* .. */
+/* .. External Subroutines .. */
+/* .. */
+/* .. Save statement .. */
+/* .. */
+/* .. Data statements .. */
+/* .. */
+/* .. Executable Statements .. */
+
+ if (first) {
+ slamc2_(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
+ base = (real) beta;
+ t = (real) it;
+ if (lrnd) {
+ rnd = 1.f;
+ i__1 = 1 - it;
+ eps = pow_ri(&base, &i__1) / 2;
+ } else {
+ rnd = 0.f;
+ i__1 = 1 - it;
+ eps = pow_ri(&base, &i__1);
+ }
+ prec = eps * base;
+ emin = (real) imin;
+ emax = (real) imax;
+ sfmin = rmin;
+ small = 1.f / rmax;
+ if (small >= sfmin) {
+
+/* Use SMALL plus a bit, to avoid the possibility of rounding */
+/* causing overflow when computing 1/sfmin. */
+
+ sfmin = small * (eps + 1.f);
+ }
+ }
+
+ if (lsame_(cmach, "E")) {
+ rmach = eps;
+ } else if (lsame_(cmach, "S")) {
+ rmach = sfmin;
+ } else if (lsame_(cmach, "B")) {
+ rmach = base;
+ } else if (lsame_(cmach, "P")) {
+ rmach = prec;
+ } else if (lsame_(cmach, "N")) {
+ rmach = t;
+ } else if (lsame_(cmach, "R")) {
+ rmach = rnd;
+ } else if (lsame_(cmach, "M")) {
+ rmach = emin;
+ } else if (lsame_(cmach, "U")) {
+ rmach = rmin;
+ } else if (lsame_(cmach, "L")) {
+ rmach = emax;
+ } else if (lsame_(cmach, "O")) {
+ rmach = rmax;
+ }
+
+ ret_val = rmach;
+ first = FALSE_;
+ return ret_val;
+
+/* End of SLAMCH */
+
+} /* slamch_ */
+
+
+/* *********************************************************************** */
+
+/* Subroutine */ int slamc1_(integer *beta, integer *t, logical *rnd, logical
+ *ieee1)
+{
+ /* Initialized data */
+
+ static logical first = TRUE_;
+
+ /* System generated locals */
+ real r__1, r__2;
+
+ /* Local variables */
+ real a, b, c__, f, t1, t2;
+ static integer lt;
+ real one, qtr;
+ static logical lrnd;
+ static integer lbeta;
+ real savec;
+ static logical lieee1;
+ extern doublereal slamc3_(real *, real *);
+
+
+/* -- LAPACK auxiliary routine (version 3.1) -- */
+/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+/* November 2006 */
+
+/* .. Scalar Arguments .. */
+/* .. */
+
+/* Purpose */
+/* ======= */
+
+/* SLAMC1 determines the machine parameters given by BETA, T, RND, and */
+/* IEEE1. */
+
+/* Arguments */
+/* ========= */
+
+/* BETA (output) INTEGER */
+/* The base of the machine. */
+
+/* T (output) INTEGER */
+/* The number of ( BETA ) digits in the mantissa. */
+
+/* RND (output) LOGICAL */
+/* Specifies whether proper rounding ( RND = .TRUE. ) or */
+/* chopping ( RND = .FALSE. ) occurs in addition. This may not */
+/* be a reliable guide to the way in which the machine performs */
+/* its arithmetic. */
+
+/* IEEE1 (output) LOGICAL */
+/* Specifies whether rounding appears to be done in the IEEE */
+/* 'round to nearest' style. */
+
+/* Further Details */
+/* =============== */
+
+/* The routine is based on the routine ENVRON by Malcolm and */
+/* incorporates suggestions by Gentleman and Marovich. See */
+
+/* Malcolm M. A. (1972) Algorithms to reveal properties of */
+/* floating-point arithmetic. Comms. of the ACM, 15, 949-951. */
+
+/* Gentleman W. M. and Marovich S. B. (1974) More on algorithms */
+/* that reveal properties of floating point arithmetic units. */
+/* Comms. of the ACM, 17, 276-277. */
+
+/* ===================================================================== */
+
+/* .. Local Scalars .. */
+/* .. */
+/* .. External Functions .. */
+/* .. */
+/* .. Save statement .. */
+/* .. */
+/* .. Data statements .. */
+/* .. */
+/* .. Executable Statements .. */
+
+ if (first) {
+ one = 1.f;
+
+/* LBETA, LIEEE1, LT and LRND are the local values of BETA, */
+/* IEEE1, T and RND. */
+
+/* Throughout this routine we use the function SLAMC3 to ensure */
+/* that relevant values are stored and not held in registers, or */
+/* are not affected by optimizers. */
+
+/* Compute a = 2.0**m with the smallest positive integer m such */
+/* that */
+
+/* fl( a + 1.0 ) = a. */
+
+ a = 1.f;
+ c__ = 1.f;
+
+/* + WHILE( C.EQ.ONE )LOOP */
+L10:
+ if (c__ == one) {
+ a *= 2;
+ c__ = slamc3_(&a, &one);
+ r__1 = -a;
+ c__ = slamc3_(&c__, &r__1);
+ goto L10;
+ }
+/* + END WHILE */
+
+/* Now compute b = 2.0**m with the smallest positive integer m */
+/* such that */
+
+/* fl( a + b ) .gt. a. */
+
+ b = 1.f;
+ c__ = slamc3_(&a, &b);
+
+/* + WHILE( C.EQ.A )LOOP */
+L20:
+ if (c__ == a) {
+ b *= 2;
+ c__ = slamc3_(&a, &b);
+ goto L20;
+ }
+/* + END WHILE */
+
+/* Now compute the base. a and c are neighbouring floating point */
+/* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so */
+/* their difference is beta. Adding 0.25 to c is to ensure that it */
+/* is truncated to beta and not ( beta - 1 ). */
+
+ qtr = one / 4;
+ savec = c__;
+ r__1 = -a;
+ c__ = slamc3_(&c__, &r__1);
+ lbeta = c__ + qtr;
+
+/* Now determine whether rounding or chopping occurs, by adding a */
+/* bit less than beta/2 and a bit more than beta/2 to a. */
+
+ b = (real) lbeta;
+ r__1 = b / 2;
+ r__2 = -b / 100;
+ f = slamc3_(&r__1, &r__2);
+ c__ = slamc3_(&f, &a);
+ if (c__ == a) {
+ lrnd = TRUE_;
+ } else {
+ lrnd = FALSE_;
+ }
+ r__1 = b / 2;
+ r__2 = b / 100;
+ f = slamc3_(&r__1, &r__2);
+ c__ = slamc3_(&f, &a);
+ if (lrnd && c__ == a) {
+ lrnd = FALSE_;
+ }
+
+/* Try and decide whether rounding is done in the IEEE 'round to */
+/* nearest' style. B/2 is half a unit in the last place of the two */
+/* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit */
+/* zero, and SAVEC is odd. Thus adding B/2 to A should not change */
+/* A, but adding B/2 to SAVEC should change SAVEC. */
+
+ r__1 = b / 2;
+ t1 = slamc3_(&r__1, &a);
+ r__1 = b / 2;
+ t2 = slamc3_(&r__1, &savec);
+ lieee1 = t1 == a && t2 > savec && lrnd;
+
+/* Now find the mantissa, t. It should be the integer part of */
+/* log to the base beta of a, however it is safer to determine t */
+/* by powering. So we find t as the smallest positive integer for */
+/* which */
+
+/* fl( beta**t + 1.0 ) = 1.0. */
+
+ lt = 0;
+ a = 1.f;
+ c__ = 1.f;
+
+/* + WHILE( C.EQ.ONE )LOOP */
+L30:
+ if (c__ == one) {
+ ++lt;
+ a *= lbeta;
+ c__ = slamc3_(&a, &one);
+ r__1 = -a;
+ c__ = slamc3_(&c__, &r__1);
+ goto L30;
+ }
+/* + END WHILE */
+
+ }
+
+ *beta = lbeta;
+ *t = lt;
+ *rnd = lrnd;
+ *ieee1 = lieee1;
+ first = FALSE_;
+ return 0;
+
+/* End of SLAMC1 */
+
+} /* slamc1_ */
+
+
+/* *********************************************************************** */
+
+/* Subroutine */ int slamc2_(integer *beta, integer *t, logical *rnd, real *
+ eps, integer *emin, real *rmin, integer *emax, real *rmax)
+{
+ /* Initialized data */
+
+ static logical first = TRUE_;
+ static logical iwarn = FALSE_;
+
+ /* Format strings */
+ static char fmt_9999[] = "(//\002 WARNING. The value EMIN may be incorre"
+ "ct:-\002,\002 EMIN = \002,i8,/\002 If, after inspection, the va"
+ "lue EMIN looks\002,\002 acceptable please comment out \002,/\002"
+ " the IF block as marked within the code of routine\002,\002 SLAM"
+ "C2,\002,/\002 otherwise supply EMIN explicitly.\002,/)";
+
+ /* System generated locals */
+ integer i__1;
+ real r__1, r__2, r__3, r__4, r__5;
+
+ /* Builtin functions */
+ double pow_ri(real *, integer *);
+ //integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
+
+ /* Local variables */
+ real a, b, c__;
+ integer i__;
+ static integer lt;
+ real one, two;
+ logical ieee;
+ real half;
+ logical lrnd;
+ static real leps;
+ real zero;
+ static integer lbeta;
+ real rbase;
+ static integer lemin, lemax;
+ integer gnmin;
+ real small;
+ integer gpmin;
+ real third;
+ static real lrmin, lrmax;
+ real sixth;
+ logical lieee1;
+ extern /* Subroutine */ int slamc1_(integer *, integer *, logical *,
+ logical *);
+ extern doublereal slamc3_(real *, real *);
+ extern /* Subroutine */ int slamc4_(integer *, real *, integer *),
+ slamc5_(integer *, integer *, integer *, logical *, integer *,
+ real *);
+ integer ngnmin, ngpmin;
+
+ /* Fortran I/O blocks */
+ static cilist io___58 = { 0, 6, 0, fmt_9999, 0 };
+
+
+
+/* -- LAPACK auxiliary routine (version 3.1) -- */
+/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+/* November 2006 */
+
+/* .. Scalar Arguments .. */
+/* .. */
+
+/* Purpose */
+/* ======= */
+
+/* SLAMC2 determines the machine parameters specified in its argument */
+/* list. */
+
+/* Arguments */
+/* ========= */
+
+/* BETA (output) INTEGER */
+/* The base of the machine. */
+
+/* T (output) INTEGER */
+/* The number of ( BETA ) digits in the mantissa. */
+
+/* RND (output) LOGICAL */
+/* Specifies whether proper rounding ( RND = .TRUE. ) or */
+/* chopping ( RND = .FALSE. ) occurs in addition. This may not */
+/* be a reliable guide to the way in which the machine performs */
+/* its arithmetic. */
+
+/* EPS (output) REAL */
+/* The smallest positive number such that */
+
+/* fl( 1.0 - EPS ) .LT. 1.0, */
+
+/* where fl denotes the computed value. */
+
+/* EMIN (output) INTEGER */
+/* The minimum exponent before (gradual) underflow occurs. */
+
+/* RMIN (output) REAL */
+/* The smallest normalized number for the machine, given by */
+/* BASE**( EMIN - 1 ), where BASE is the floating point value */
+/* of BETA. */
+
+/* EMAX (output) INTEGER */
+/* The maximum exponent before overflow occurs. */
+
+/* RMAX (output) REAL */
+/* The largest positive number for the machine, given by */
+/* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point */
+/* value of BETA. */
+
+/* Further Details */
+/* =============== */
+
+/* The computation of EPS is based on a routine PARANOIA by */
+/* W. Kahan of the University of California at Berkeley. */
+
+/* ===================================================================== */
+
+/* .. Local Scalars .. */
+/* .. */
+/* .. External Functions .. */
+/* .. */
+/* .. External Subroutines .. */
+/* .. */
+/* .. Intrinsic Functions .. */
+/* .. */
+/* .. Save statement .. */
+/* .. */
+/* .. Data statements .. */
+/* .. */
+/* .. Executable Statements .. */
+
+ if (first) {
+ zero = 0.f;
+ one = 1.f;
+ two = 2.f;
+
+/* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of */
+/* BETA, T, RND, EPS, EMIN and RMIN. */
+
+/* Throughout this routine we use the function SLAMC3 to ensure */
+/* that relevant values are stored and not held in registers, or */
+/* are not affected by optimizers. */
+
+/* SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. */
+
+ slamc1_(&lbeta, <, &lrnd, &lieee1);
+
+/* Start to find EPS. */
+
+ b = (real) lbeta;
+ i__1 = -lt;
+ a = pow_ri(&b, &i__1);
+ leps = a;
+
+/* Try some tricks to see whether or not this is the correct EPS. */
+
+ b = two / 3;
+ half = one / 2;
+ r__1 = -half;
+ sixth = slamc3_(&b, &r__1);
+ third = slamc3_(&sixth, &sixth);
+ r__1 = -half;
+ b = slamc3_(&third, &r__1);
+ b = slamc3_(&b, &sixth);
+ b = dabs(b);
+ if (b < leps) {
+ b = leps;
+ }
+
+ leps = 1.f;
+
+/* + WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */
+L10:
+ if (leps > b && b > zero) {
+ leps = b;
+ r__1 = half * leps;
+/* Computing 5th power */
+ r__3 = two, r__4 = r__3, r__3 *= r__3;
+/* Computing 2nd power */
+ r__5 = leps;
+ r__2 = r__4 * (r__3 * r__3) * (r__5 * r__5);
+ c__ = slamc3_(&r__1, &r__2);
+ r__1 = -c__;
+ c__ = slamc3_(&half, &r__1);
+ b = slamc3_(&half, &c__);
+ r__1 = -b;
+ c__ = slamc3_(&half, &r__1);
+ b = slamc3_(&half, &c__);
+ goto L10;
+ }
+/* + END WHILE */
+
+ if (a < leps) {
+ leps = a;
+ }
+
+/* Computation of EPS complete. */
+
+/* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). */
+/* Keep dividing A by BETA until (gradual) underflow occurs. This */
+/* is detected when we cannot recover the previous A. */
+
+ rbase = one / lbeta;
+ small = one;
+ for (i__ = 1; i__ <= 3; ++i__) {
+ r__1 = small * rbase;
+ small = slamc3_(&r__1, &zero);
+/* L20: */
+ }
+ a = slamc3_(&one, &small);
+ slamc4_(&ngpmin, &one, &lbeta);
+ r__1 = -one;
+ slamc4_(&ngnmin, &r__1, &lbeta);
+ slamc4_(&gpmin, &a, &lbeta);
+ r__1 = -a;
+ slamc4_(&gnmin, &r__1, &lbeta);
+ ieee = FALSE_;
+
+ if (ngpmin == ngnmin && gpmin == gnmin) {
+ if (ngpmin == gpmin) {
+ lemin = ngpmin;
+/* ( Non twos-complement machines, no gradual underflow; */
+/* e.g., VAX ) */
+ } else if (gpmin - ngpmin == 3) {
+ lemin = ngpmin - 1 + lt;
+ ieee = TRUE_;
+/* ( Non twos-complement machines, with gradual underflow; */
+/* e.g., IEEE standard followers ) */
+ } else {
+ lemin = min(ngpmin,gpmin);
+/* ( A guess; no known machine ) */
+ iwarn = TRUE_;
+ }
+
+ } else if (ngpmin == gpmin && ngnmin == gnmin) {
+ if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1) {
+ lemin = max(ngpmin,ngnmin);
+/* ( Twos-complement machines, no gradual underflow; */
+/* e.g., CYBER 205 ) */
+ } else {
+ lemin = min(ngpmin,ngnmin);
+/* ( A guess; no known machine ) */
+ iwarn = TRUE_;
+ }
+
+ } else if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1 && gpmin == gnmin)
+ {
+ if (gpmin - min(ngpmin,ngnmin) == 3) {
+ lemin = max(ngpmin,ngnmin) - 1 + lt;
+/* ( Twos-complement machines with gradual underflow; */
+/* no known machine ) */
+ } else {
+ lemin = min(ngpmin,ngnmin);
+/* ( A guess; no known machine ) */
+ iwarn = TRUE_;
+ }
+
+ } else {
+/* Computing MIN */
+ i__1 = min(ngpmin,ngnmin), i__1 = min(i__1,gpmin);
+ lemin = min(i__1,gnmin);
+/* ( A guess; no known machine ) */
+ iwarn = TRUE_;
+ }
+ first = FALSE_;
+/* ** */
+/* Comment out this if block if EMIN is ok */
+ if (iwarn) {
+ first = TRUE_;
+ printf("\n\n WARNING. The value EMIN may be incorrect:- ");
+ printf("EMIN = %8i\n",lemin);
+ printf("If, after inspection, the value EMIN looks acceptable");
+ printf("please comment out \n the IF block as marked within the");
+ printf("code of routine SLAMC2, \n otherwise supply EMIN");
+ printf("explicitly.\n");
+ /*
+ s_wsfe(&io___58);
+ do_fio(&c__1, (char *)&lemin, (ftnlen)sizeof(integer));
+ e_wsfe();
+ */
+ }
+/* ** */
+
+/* Assume IEEE arithmetic if we found denormalised numbers above, */
+/* or if arithmetic seems to round in the IEEE style, determined */
+/* in routine SLAMC1. A true IEEE machine should have both things */
+/* true; however, faulty machines may have one or the other. */
+
+ ieee = ieee || lieee1;
+
+/* Compute RMIN by successive division by BETA. We could compute */
+/* RMIN as BASE**( EMIN - 1 ), but some machines underflow during */
+/* this computation. */
+
+ lrmin = 1.f;
+ i__1 = 1 - lemin;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ r__1 = lrmin * rbase;
+ lrmin = slamc3_(&r__1, &zero);
+/* L30: */
+ }
+
+/* Finally, call SLAMC5 to compute EMAX and RMAX. */
+
+ slamc5_(&lbeta, <, &lemin, &ieee, &lemax, &lrmax);
+ }
+
+ *beta = lbeta;
+ *t = lt;
+ *rnd = lrnd;
+ *eps = leps;
+ *emin = lemin;
+ *rmin = lrmin;
+ *emax = lemax;
+ *rmax = lrmax;
+
+ return 0;
+
+
+/* End of SLAMC2 */
+
+} /* slamc2_ */
+
+
+/* *********************************************************************** */
+
+/* Subroutine */ int slamc4_(integer *emin, real *start, integer *base)
+{
+ /* System generated locals */
+ integer i__1;
+ real r__1;
+
+ /* Local variables */
+ real a;
+ integer i__;
+ real b1, b2, c1, c2, d1, d2, one, zero, rbase;
+ extern doublereal slamc3_(real *, real *);
+
+
+/* -- LAPACK auxiliary routine (version 3.1) -- */
+/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+/* November 2006 */
+
+/* .. Scalar Arguments .. */
+/* .. */
+
+/* Purpose */
+/* ======= */
+
+/* SLAMC4 is a service routine for SLAMC2. */
+
+/* Arguments */
+/* ========= */
+
+/* EMIN (output) INTEGER */
+/* The minimum exponent before (gradual) underflow, computed by */
+/* setting A = START and dividing by BASE until the previous A */
+/* can not be recovered. */
+
+/* START (input) REAL */
+/* The starting point for determining EMIN. */
+
+/* BASE (input) INTEGER */
+/* The base of the machine. */
+
+/* ===================================================================== */
+
+/* .. Local Scalars .. */
+/* .. */
+/* .. External Functions .. */
+/* .. */
+/* .. Executable Statements .. */
+
+ a = *start;
+ one = 1.f;
+ rbase = one / *base;
+ zero = 0.f;
+ *emin = 1;
+ r__1 = a * rbase;
+ b1 = slamc3_(&r__1, &zero);
+ c1 = a;
+ c2 = a;
+ d1 = a;
+ d2 = a;
+/* + WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. */
+/* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP */
+L10:
+ if (c1 == a && c2 == a && d1 == a && d2 == a) {
+ --(*emin);
+ a = b1;
+ r__1 = a / *base;
+ b1 = slamc3_(&r__1, &zero);
+ r__1 = b1 * *base;
+ c1 = slamc3_(&r__1, &zero);
+ d1 = zero;
+ i__1 = *base;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ d1 += b1;
+/* L20: */
+ }
+ r__1 = a * rbase;
+ b2 = slamc3_(&r__1, &zero);
+ r__1 = b2 / rbase;
+ c2 = slamc3_(&r__1, &zero);
+ d2 = zero;
+ i__1 = *base;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ d2 += b2;
+/* L30: */
+ }
+ goto L10;
+ }
+/* + END WHILE */
+
+ return 0;
+
+/* End of SLAMC4 */
+
+} /* slamc4_ */
+
+
+/* *********************************************************************** */
+
+/* Subroutine */ int slamc5_(integer *beta, integer *p, integer *emin,
+ logical *ieee, integer *emax, real *rmax)
+{
+ /* System generated locals */
+ integer i__1;
+ real r__1;
+
+ /* Local variables */
+ integer i__;
+ real y, z__;
+ integer try__, lexp;
+ real oldy;
+ integer uexp, nbits;
+ extern doublereal slamc3_(real *, real *);
+ real recbas;
+ integer exbits, expsum;
+
+
+/* -- LAPACK auxiliary routine (version 3.1) -- */
+/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+/* November 2006 */
+
+/* .. Scalar Arguments .. */
+/* .. */
+
+/* Purpose */
+/* ======= */
+
+/* SLAMC5 attempts to compute RMAX, the largest machine floating-point */
+/* number, without overflow. It assumes that EMAX + abs(EMIN) sum */
+/* approximately to a power of 2. It will fail on machines where this */
+/* assumption does not hold, for example, the Cyber 205 (EMIN = -28625, */
+/* EMAX = 28718). It will also fail if the value supplied for EMIN is */
+/* too large (i.e. too close to zero), probably with overflow. */
+
+/* Arguments */
+/* ========= */
+
+/* BETA (input) INTEGER */
+/* The base of floating-point arithmetic. */
+
+/* P (input) INTEGER */
+/* The number of base BETA digits in the mantissa of a */
+/* floating-point value. */
+
+/* EMIN (input) INTEGER */
+/* The minimum exponent before (gradual) underflow. */
+
+/* IEEE (input) LOGICAL */
+/* A logical flag specifying whether or not the arithmetic */
+/* system is thought to comply with the IEEE standard. */
+
+/* EMAX (output) INTEGER */
+/* The largest exponent before overflow */
+
+/* RMAX (output) REAL */
+/* The largest machine floating-point number. */
+
+/* ===================================================================== */
+
+/* .. Parameters .. */
+/* .. */
+/* .. Local Scalars .. */
+/* .. */
+/* .. External Functions .. */
+/* .. */
+/* .. Intrinsic Functions .. */
+/* .. */
+/* .. Executable Statements .. */
+
+/* First compute LEXP and UEXP, two powers of 2 that bound */
+/* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum */
+/* approximately to the bound that is closest to abs(EMIN). */
+/* (EMAX is the exponent of the required number RMAX). */
+
+ lexp = 1;
+ exbits = 1;
+L10:
+ try__ = lexp << 1;
+ if (try__ <= -(*emin)) {
+ lexp = try__;
+ ++exbits;
+ goto L10;
+ }
+ if (lexp == -(*emin)) {
+ uexp = lexp;
+ } else {
+ uexp = try__;
+ ++exbits;
+ }
+
+/* Now -LEXP is less than or equal to EMIN, and -UEXP is greater */
+/* than or equal to EMIN. EXBITS is the number of bits needed to */
+/* store the exponent. */
+
+ if (uexp + *emin > -lexp - *emin) {
+ expsum = lexp << 1;
+ } else {
+ expsum = uexp << 1;
+ }
+
+/* EXPSUM is the exponent range, approximately equal to */
+/* EMAX - EMIN + 1 . */
+
+ *emax = expsum + *emin - 1;
+ nbits = exbits + 1 + *p;
+
+/* NBITS is the total number of bits needed to store a */
+/* floating-point number. */
+
+ if (nbits % 2 == 1 && *beta == 2) {
+
+/* Either there are an odd number of bits used to store a */
+/* floating-point number, which is unlikely, or some bits are */
+/* not used in the representation of numbers, which is possible, */
+/* (e.g. Cray machines) or the mantissa has an implicit bit, */
+/* (e.g. IEEE machines, Dec Vax machines), which is perhaps the */
+/* most likely. We have to assume the last alternative. */
+/* If this is true, then we need to reduce EMAX by one because */
+/* there must be some way of representing zero in an implicit-bit */
+/* system. On machines like Cray, we are reducing EMAX by one */
+/* unnecessarily. */
+
+ --(*emax);
+ }
+
+ if (*ieee) {
+
+/* Assume we are on an IEEE machine which reserves one exponent */
+/* for infinity and NaN. */
+
+ --(*emax);
+ }
+
+/* Now create RMAX, the largest machine number, which should */
+/* be equal to (1.0 - BETA**(-P)) * BETA**EMAX . */
+
+/* First compute 1.0 - BETA**(-P), being careful that the */
+/* result is less than 1.0 . */
+
+ recbas = 1.f / *beta;
+ z__ = *beta - 1.f;
+ y = 0.f;
+ i__1 = *p;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ z__ *= recbas;
+ if (y < 1.f) {
+ oldy = y;
+ }
+ y = slamc3_(&y, &z__);
+/* L20: */
+ }
+ if (y >= 1.f) {
+ y = oldy;
+ }
+
+/* Now multiply by BETA**EMAX to get RMAX. */
+
+ i__1 = *emax;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ r__1 = y * *beta;
+ y = slamc3_(&r__1, &c_b32);
+/* L30: */
+ }
+
+ *rmax = y;
+ return 0;
+
+/* End of SLAMC5 */
+
+} /* slamc5_ */
+
+#endif
projects / opencv / blobdiff