#include <math.h>

#include <stdlib.h>

#include <string.h>

#include <stdio.h>

#include <complex.h>

#ifdef complex

#undef complex

#endif

#ifdef I

#undef I

#endif

#if defined(_WIN64)

typedef long long BLASLONG;

typedef unsigned long long BLASULONG;

#else

typedef long BLASLONG;

typedef unsigned long BLASULONG;

#endif

#ifdef LAPACK_ILP64

typedef BLASLONG blasint;

#if defined(_WIN64)

#define blasabs(x) llabs(x)

#else

#define blasabs(x) labs(x)

#endif

#else

typedef int blasint;

#define blasabs(x) abs(x)

#endif

typedef blasint integer;

typedef unsigned int uinteger;

typedef char *address;

typedef short int shortint;

typedef float real;

typedef double doublereal;

typedef struct { real r, i; } complex;

typedef struct { doublereal r, i; } doublecomplex;

#ifdef _MSC_VER

static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;}

static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;}

static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;}

static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;}

#else

static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}

static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}

static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}

static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}

#endif

#define pCf(z) (*_pCf(z))

#define pCd(z) (*_pCd(z))

typedef int logical;

typedef short int shortlogical;

typedef char logical1;

typedef char integer1;

#define TRUE_ (1)

#define FALSE_ (0)

/* Extern is for use with -E */

#ifndef Extern

#define Extern extern

#endif

/* I/O stuff */

typedef int flag;

typedef int ftnlen;

typedef int ftnint;

/*external read, write*/

typedef struct

{ flag cierr;

  ftnint ciunit;

  flag ciend;

  char *cifmt;

  ftnint cirec;

} cilist;

/*internal read, write*/

typedef struct

{ flag icierr;

  char *iciunit;

  flag iciend;

  char *icifmt;

  ftnint icirlen;

  ftnint icirnum;

} icilist;

/*open*/

typedef struct

{ flag oerr;

  ftnint ounit;

  char *ofnm;

  ftnlen ofnmlen;

  char *osta;

  char *oacc;

  char *ofm;

  ftnint orl;

  char *oblnk;

} olist;

/*close*/

typedef struct

{ flag cerr;

  ftnint cunit;

  char *csta;

} cllist;

/*rewind, backspace, endfile*/

typedef struct

{ flag aerr;

  ftnint aunit;

} alist;

/* inquire */

typedef struct

{ flag inerr;

  ftnint inunit;

  char *infile;

  ftnlen infilen;

  ftnint  *inex;  /*parameters in standard's order*/

  ftnint  *inopen;

  ftnint  *innum;

  ftnint  *innamed;

  char  *inname;

  ftnlen  innamlen;

  char  *inacc;

  ftnlen  inacclen;

  char  *inseq;

  ftnlen  inseqlen;

  char  *indir;

  ftnlen  indirlen;

  char  *infmt;

  ftnlen  infmtlen;

  char  *inform;

  ftnint  informlen;

  char  *inunf;

  ftnlen  inunflen;

  ftnint  *inrecl;

  ftnint  *innrec;

  char  *inblank;

  ftnlen  inblanklen;

} inlist;

#define VOID void

union Multitype { /* for multiple entry points */

  integer1 g;

  shortint h;

  integer i;

  /* longint j; */

  real r;

  doublereal d;

  complex c;

  doublecomplex z;

  };

typedef union Multitype Multitype;

struct Vardesc {  /* for Namelist */

  char *name;

  char *addr;

  ftnlen *dims;

  int  type;

  };

typedef struct Vardesc Vardesc;

struct Namelist {

  char *name;

  Vardesc **vars;

  int nvars;

  };

typedef struct Namelist Namelist;

#define abs(x) ((x) >= 0 ? (x) : -(x))

#define dabs(x) (fabs(x))

#define f2cmin(a,b) ((a) <= (b) ? (a) : (b))

#define f2cmax(a,b) ((a) >= (b) ? (a) : (b))

#define dmin(a,b) (f2cmin(a,b))

#define dmax(a,b) (f2cmax(a,b))

#define bit_test(a,b) ((a) >> (b) & 1)

#define bit_clear(a,b)  ((a) & ~((uinteger)1 << (b)))

#define bit_set(a,b)  ((a) |  ((uinteger)1 << (b)))

#define abort_() { sig_die("Fortran abort routine called", 1); }

#define c_abs(z) (cabsf(Cf(z)))

#define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }

#ifdef _MSC_VER

#define c_div(c, a, b) {Cf(c)._Val[0] = (Cf(a)._Val[0]/Cf(b)._Val[0]); Cf(c)._Val[1]=(Cf(a)._Val[1]/Cf(b)._Val[1]);}

#define z_div(c, a, b) {Cd(c)._Val[0] = (Cd(a)._Val[0]/Cd(b)._Val[0]); Cd(c)._Val[1]=(Cd(a)._Val[1]/df(b)._Val[1]);}

#else

#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}

#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}

#endif

#define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}

#define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}

#define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}

//#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}

#define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}

#define d_abs(x) (fabs(*(x)))

#define d_acos(x) (acos(*(x)))

#define d_asin(x) (asin(*(x)))

#define d_atan(x) (atan(*(x)))

#define d_atn2(x, y) (atan2(*(x),*(y)))

#define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }

#define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); }

#define d_cos(x) (cos(*(x)))

#define d_cosh(x) (cosh(*(x)))

#define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )

#define d_exp(x) (exp(*(x)))

#define d_imag(z) (cimag(Cd(z)))

#define r_imag(z) (cimagf(Cf(z)))

#define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))

#define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))

#define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )

#define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) )

#define d_log(x) (log(*(x)))

#define d_mod(x, y) (fmod(*(x), *(y)))

#define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))

#define d_nint(x) u_nint(*(x))

#define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))

#define d_sign(a,b) u_sign(*(a),*(b))

#define r_sign(a,b) u_sign(*(a),*(b))

#define d_sin(x) (sin(*(x)))

#define d_sinh(x) (sinh(*(x)))

#define d_sqrt(x) (sqrt(*(x)))

#define d_tan(x) (tan(*(x)))

#define d_tanh(x) (tanh(*(x)))

#define i_abs(x) abs(*(x))

#define i_dnnt(x) ((integer)u_nint(*(x)))

#define i_len(s, n) (n)

#define i_nint(x) ((integer)u_nint(*(x)))

#define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))

#define pow_dd(ap, bp) ( pow(*(ap), *(bp)))

#define pow_si(B,E) spow_ui(*(B),*(E))

#define pow_ri(B,E) spow_ui(*(B),*(E))

#define pow_di(B,E) dpow_ui(*(B),*(E))

#define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}

#define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));}

#define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}

#define s_cat(lpp, rpp, rnp, np, llp) {   ftnlen i, nc, ll; char *f__rp, *lp;   ll = (llp); lp = (lpp);   for(i=0; i < (int)*(np); ++i) {           nc = ll;          if((rnp)[i] < nc) nc = (rnp)[i];          ll -= nc;           f__rp = (rpp)[i];           while(--nc >= 0) *lp++ = *(f__rp)++;         }  while(--ll >= 0) *lp++ = ' '; }

#define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))

#define s_copy(A,B,C,D) { int __i,__m; for (__i=0, __m=f2cmin((C),(D)); __i<__m && (B)[__i] != 0; ++__i) (A)[__i] = (B)[__i]; }

#define sig_die(s, kill) { exit(1); }

#define s_stop(s, n) {exit(0);}

#define z_abs(z) (cabs(Cd(z)))

#define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}

#define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}

#define myexit_() break;

#define mycycle() continue;

#define myceiling(w) {ceil(w)}

#define myhuge(w) {HUGE_VAL}

//#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);}

#define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)}

/* procedure parameter types for -A and -C++ */

#define F2C_proc_par_types 1

static double dpow_ui(double x, integer n) {

  double pow=1.0; unsigned long int u;

  if(n != 0) {

    if(n < 0) n = -n, x = 1/x;

    for(u = n; ; ) {

      if(u & 01) pow *= x;

      if(u >>= 1) x *= x;

      else break;

    }

  }

  return pow;

}

/*  -- translated by f2c (version 20000121).

   You must link the resulting object file with the libraries:

  -lf2c -lm   (in that order)

*/

/* Table of constant values */

static integer c__1 = 1;

static doublereal c_b32 = 0.;

/* > \brief \b DLAMCHF77 deprecated */

/*  =========== DOCUMENTATION =========== */

/* Online html documentation available at */

/*            http://www.netlib.org/lapack/explore-html/ */

/*  Definition: */

/*  =========== */

/*      DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) */

/*     CHARACTER          CMACH */

/* > \par Purpose: */

/*  ============= */

/* > */

/* > \verbatim */

/* > */

/* > DLAMCHF77 determines double precision machine parameters. */

/* > \endverbatim */

/*  Arguments: */

/*  ========== */

/* > \param[in] CMACH */

/* > \verbatim */

/* >          Specifies the value to be returned by DLAMCH: */

/* >          = 'E' or 'e',   DLAMCH := eps */

/* >          = 'S' or 's ,   DLAMCH := sfmin */

/* >          = 'B' or 'b',   DLAMCH := base */

/* >          = 'P' or 'p',   DLAMCH := eps*base */

/* >          = 'N' or 'n',   DLAMCH := t */

/* >          = 'R' or 'r',   DLAMCH := rnd */

/* >          = 'M' or 'm',   DLAMCH := emin */

/* >          = 'U' or 'u',   DLAMCH := rmin */

/* >          = 'L' or 'l',   DLAMCH := emax */

/* >          = 'O' or 'o',   DLAMCH := 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) */

/* > \endverbatim */

/*  Authors: */

/*  ======== */

/* > \author Univ. of Tennessee */

/* > \author Univ. of California Berkeley */

/* > \author Univ. of Colorado Denver */

/* > \author NAG Ltd. */

/* > \date April 2012 */

/* > \ingroup auxOTHERauxiliary */

/*  ===================================================================== */

doublereal dlamch_(char *cmach)

{

    /* Initialized data */

    static logical first = TRUE_;

    /* System generated locals */

    integer i__1;

    doublereal ret_val;

    /* Local variables */

    static doublereal base;

    integer beta;

    static doublereal emin, prec, emax;

    integer imin, imax;

    logical lrnd;

    static doublereal rmin, rmax, t;

    doublereal rmach;

    extern logical lsame_(char *, char *);

    doublereal small;

    static doublereal sfmin;

    extern /* Subroutine */ int dlamc2_(integer *, integer *, logical *, 

      doublereal *, integer *, doublereal *, integer *, doublereal *);

    integer it;

    static doublereal rnd, eps;

/*  -- LAPACK auxiliary routine (version 3.7.0) -- */

/*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */

/*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */

/*     April 2012 */

    if (first) {

  dlamc2_(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);

  base = (doublereal) beta;

  t = (doublereal) it;

  if (lrnd) {

      rnd = 1.;

      i__1 = 1 - it;

      eps = pow_di(&base, &i__1) / 2;

  } else {

      rnd = 0.;

      i__1 = 1 - it;

      eps = pow_di(&base, &i__1);

  }

  prec = eps * base;

  emin = (doublereal) imin;

  emax = (doublereal) imax;

  sfmin = rmin;

  small = 1. / 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.);

  }

    }

    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 DLAMCH */

} /* dlamch_ */

/* *********************************************************************** */

/* > \brief \b DLAMC1 */

/* > \details */

/* > \b Purpose: */

/* > \verbatim */

/* > DLAMC1 determines the machine parameters given by BETA, T, RND, and */

/* > IEEE1. */

/* > \endverbatim */

/* > */

/* > \param[out] BETA */

/* > \verbatim */

/* >          The base of the machine. */

/* > \endverbatim */

/* > */

/* > \param[out] T */

/* > \verbatim */

/* >          The number of ( BETA ) digits in the mantissa. */

/* > \endverbatim */

/* > */

/* > \param[out] RND */

/* > \verbatim */

/* >          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. */

/* > \endverbatim */

/* > */

/* > \param[out] IEEE1 */

/* > \verbatim */

/* >          Specifies whether rounding appears to be done in the IEEE */

/* >          'round to nearest' style. */

/* > \endverbatim */

/* > \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. 

of Colorado Denver and NAG Ltd.. */

/* > \date April 2012 */

/* > \ingroup auxOTHERauxiliary */

/* > */

/* > \details \b Further \b Details */

/* > \verbatim */

/* > */

/* >  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. */

/* > \endverbatim */

/* > */

/* Subroutine */ int dlamc1_(integer *beta, integer *t, logical *rnd, logical 

  *ieee1)

{

    /* Initialized data */

    static logical first = TRUE_;

    /* System generated locals */

    doublereal d__1, d__2;

    /* Local variables */

    static logical lrnd;

    doublereal a, b, c__, f;

    static integer lbeta;

    doublereal savec;

    extern doublereal dlamc3_(doublereal *, doublereal *);

    static logical lieee1;

    doublereal t1, t2;

    static integer lt;

    doublereal one, qtr;

/*  -- LAPACK auxiliary routine (version 3.7.0) -- */

/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */

/*     November 2010 */

/* ===================================================================== */

    if (first) {

  one = 1.;

/*        LBETA,  LIEEE1,  LT and  LRND  are the  local values  of  BETA, */

/*        IEEE1, T and RND. */

/*        Throughout this routine  we use the function  DLAMC3  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.;

  c__ = 1.;

/* +       WHILE( C.EQ.ONE )LOOP */

L10:

  if (c__ == one) {

      a *= 2;

      c__ = dlamc3_(&a, &one);

      d__1 = -a;

      c__ = dlamc3_(&c__, &d__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.;

  c__ = dlamc3_(&a, &b);

/* +       WHILE( C.EQ.A )LOOP */

L20:

  if (c__ == a) {

      b *= 2;

      c__ = dlamc3_(&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__;

  d__1 = -a;

  c__ = dlamc3_(&c__, &d__1);

  lbeta = (integer) (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 = (doublereal) lbeta;

  d__1 = b / 2;

  d__2 = -b / 100;

  f = dlamc3_(&d__1, &d__2);

  c__ = dlamc3_(&f, &a);

  if (c__ == a) {

      lrnd = TRUE_;

  } else {

      lrnd = FALSE_;

  }

  d__1 = b / 2;

  d__2 = b / 100;

  f = dlamc3_(&d__1, &d__2);

  c__ = dlamc3_(&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. */

  d__1 = b / 2;

  t1 = dlamc3_(&d__1, &a);

  d__1 = b / 2;

  t2 = dlamc3_(&d__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.;

  c__ = 1.;

/* +       WHILE( C.EQ.ONE )LOOP */

L30:

  if (c__ == one) {

      ++lt;

      a *= lbeta;

      c__ = dlamc3_(&a, &one);

      d__1 = -a;

      c__ = dlamc3_(&c__, &d__1);

      goto L30;

  }

/* +       END WHILE */

    }

    *beta = lbeta;

    *t = lt;

    *rnd = lrnd;

    *ieee1 = lieee1;

    first = FALSE_;

    return 0;

/*     End of DLAMC1 */

} /* dlamc1_ */

/* *********************************************************************** */

/* > \brief \b DLAMC2 */

/* > \details */

/* > \b Purpose: */

/* > \verbatim */

/* > DLAMC2 determines the machine parameters specified in its argument */

/* > list. */

/* > \endverbatim */

/* > \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. 

of Colorado Denver and NAG Ltd.. */

/* > \date April 2012 */

/* > \ingroup auxOTHERauxiliary */

/* > */

/* > \param[out] BETA */

/* > \verbatim */

/* >          The base of the machine. */

/* > \endverbatim */

/* > */

/* > \param[out] T */

/* > \verbatim */

/* >          The number of ( BETA ) digits in the mantissa. */

/* > \endverbatim */

/* > */

/* > \param[out] RND */

/* > \verbatim */

/* >          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. */

/* > \endverbatim */

/* > */

/* > \param[out] EPS */

/* > \verbatim */

/* >          The smallest positive number such that */

/* >             fl( 1.0 - EPS ) .LT. 1.0, */

/* >          where fl denotes the computed value. */

/* > \endverbatim */

/* > */

/* > \param[out] EMIN */

/* > \verbatim */

/* >          The minimum exponent before (gradual) underflow occurs. */

/* > \endverbatim */

/* > */

/* > \param[out] RMIN */

/* > \verbatim */

/* >          The smallest normalized number for the machine, given by */

/* >          BASE**( EMIN - 1 ), where  BASE  is the floating point value */

/* >          of BETA. */

/* > \endverbatim */

/* > */

/* > \param[out] EMAX */

/* > \verbatim */

/* >          The maximum exponent before overflow occurs. */

/* > \endverbatim */

/* > */

/* > \param[out] RMAX */

/* > \verbatim */

/* >          The largest positive number for the machine, given by */

/* >          BASE**EMAX * ( 1 - EPS ), where  BASE  is the floating point */

/* >          value of BETA. */

/* > \endverbatim */

/* > */

/* > \details \b Further \b Details */

/* > \verbatim */

/* > */

/* >  The computation of  EPS  is based on a routine PARANOIA by */

/* >  W. Kahan of the University of California at Berkeley. */

/* > \endverbatim */

/* Subroutine */ int dlamc2_(integer *beta, integer *t, logical *rnd, 

  doublereal *eps, integer *emin, doublereal *rmin, integer *emax, 

  doublereal *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 DLAM"

      "C2,\002,/\002 otherwise supply EMIN explicitly.\002,/)";

    /* System generated locals */

    integer i__1;

    doublereal d__1, d__2, d__3, d__4, d__5;

    /* Local variables */

    logical ieee;

    doublereal half;

    logical lrnd;

    static doublereal leps;

    doublereal zero, a, b, c__;

    integer i__;

    static integer lbeta;

    doublereal rbase;

    static integer lemin, lemax;

    integer gnmin;

    doublereal small;

    integer gpmin;

    doublereal third;

    static doublereal lrmin, lrmax;

    doublereal sixth;

    extern /* Subroutine */ int dlamc1_(integer *, integer *, logical *, 

      logical *);

    extern doublereal dlamc3_(doublereal *, doublereal *);

    logical lieee1;

    extern /* Subroutine */ int dlamc4_(integer *, doublereal *, integer *), 

      dlamc5_(integer *, integer *, integer *, logical *, integer *, 

      doublereal *);

    static integer lt;

    integer ngnmin, ngpmin;

    doublereal one, two;

    /* Fortran I/O blocks */

    static cilist io___58 = { 0, 6, 0, fmt_9999, 0 };

/*  -- LAPACK auxiliary routine (version 3.7.0) -- */

/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */

/*     November 2010 */

/* ===================================================================== */

    if (first) {

  zero = 0.;

  one = 1.;

  two = 2.;

/*        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  DLAMC3  to ensure */

/*        that relevant values are stored  and not held in registers,  or */

/*        are not affected by optimizers. */

/*        DLAMC1 returns the parameters  LBETA, LT, LRND and LIEEE1. */

  dlamc1_(&lbeta, &lt, &lrnd, &lieee1);

/*        Start to find EPS. */

  b = (doublereal) lbeta;

  i__1 = -lt;

  a = pow_di(&b, &i__1);

  leps = a;

/*        Try some tricks to see whether or not this is the correct  EPS. */

  b = two / 3;

  half = one / 2;

  d__1 = -half;

  sixth = dlamc3_(&b, &d__1);

  third = dlamc3_(&sixth, &sixth);

  d__1 = -half;

  b = dlamc3_(&third, &d__1);

  b = dlamc3_(&b, &sixth);

  b = abs(b);

  if (b < leps) {

      b = leps;

  }

  leps = 1.;

/* +       WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */

L10:

  if (leps > b && b > zero) {

      leps = b;

      d__1 = half * leps;

/* Computing 5th power */

      d__3 = two, d__4 = d__3, d__3 *= d__3;

/* Computing 2nd power */

      d__5 = leps;

      d__2 = d__4 * (d__3 * d__3) * (d__5 * d__5);

      c__ = dlamc3_(&d__1, &d__2);

      d__1 = -c__;

      c__ = dlamc3_(&half, &d__1);

      b = dlamc3_(&half, &c__);

      d__1 = -b;

      c__ = dlamc3_(&half, &d__1);

      b = dlamc3_(&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__) {

      d__1 = small * rbase;

      small = dlamc3_(&d__1, &zero);

/* L20: */

  }

  a = dlamc3_(&one, &small);

  dlamc4_(&ngpmin, &one, &lbeta);

  d__1 = -one;

  dlamc4_(&ngnmin, &d__1, &lbeta);

  dlamc4_(&gpmin, &a, &lbeta);

  d__1 = -a;

  dlamc4_(&gnmin, &d__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 = f2cmin(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 = f2cmax(ngpmin,ngnmin);

/*            ( Twos-complement machines, no gradual underflow; */

/*              e.g., CYBER 205 ) */

      } else {

    lemin = f2cmin(ngpmin,ngnmin);

/*            ( A guess; no known machine ) */

    iwarn = TRUE_;

      }

  } else if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1 && gpmin == gnmin)

     {

      if (gpmin - f2cmin(ngpmin,ngnmin) == 3) {

    lemin = f2cmax(ngpmin,ngnmin) - 1 + lt;

/*            ( Twos-complement machines with gradual underflow; */

/*              no known machine ) */

      } else {

    lemin = f2cmin(ngpmin,ngnmin);

/*            ( A guess; no known machine ) */

    iwarn = TRUE_;

      }

  } else {

/* Computing MIN */

      i__1 = f2cmin(ngpmin,ngnmin), i__1 = f2cmin(i__1,gpmin);

      lemin = f2cmin(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_;

      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 DLAMC1. 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.;

  i__1 = 1 - lemin;

  for (i__ = 1; i__ <= i__1; ++i__) {

      d__1 = lrmin * rbase;

      lrmin = dlamc3_(&d__1, &zero);

/* L30: */

  }

/*        Finally, call DLAMC5 to compute EMAX and RMAX. */

  dlamc5_(&lbeta, &lt, &lemin, &ieee, &lemax, &lrmax);

    }

    *beta = lbeta;

    *t = lt;

    *rnd = lrnd;

    *eps = leps;

    *emin = lemin;

    *rmin = lrmin;

    *emax = lemax;

    *rmax = lrmax;

    return 0;

/*     End of DLAMC2 */

} /* dlamc2_ */

/* *********************************************************************** */

/* > \brief \b DLAMC3 */

/* > \details */

/* > \b Purpose: */

/* > \verbatim */

/* > DLAMC3  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. */

/* > \endverbatim */

/* > */

/* > \param[in] A */

/* > */

/* > \param[in] B */

/* > \verbatim */

/* >          The values A and B. */

/* > \endverbatim */

doublereal dlamc3_(doublereal *a, doublereal *b)

{

    /* System generated locals */

    doublereal ret_val;

/*  -- LAPACK auxiliary routine (version 3.7.0) -- */

/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */

/*     November 2010 */

/* ===================================================================== */

    ret_val = *a + *b;

    return ret_val;

/*     End of DLAMC3 */

} /* dlamc3_ */