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