#include #include #include #include #include #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);} static char junk[] = "\n@(#)LIBF77 VERSION 19990503\n"; #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 #ifdef __cplusplus typedef logical (*L_fp)(...); #else typedef logical (*L_fp)(); #endif static float spow_ui(float x, integer n) { float 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; } 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; } #ifdef _MSC_VER static _Fcomplex cpow_ui(complex x, integer n) { complex pow={1.0,0.0}; unsigned long int u; if(n != 0) { if(n < 0) n = -n, x.r = 1/x.r, x.i=1/x.i; for(u = n; ; ) { if(u & 01) pow.r *= x.r, pow.i *= x.i; if(u >>= 1) x.r *= x.r, x.i *= x.i; else break; } } _Fcomplex p={pow.r, pow.i}; return p; } #else static _Complex float cpow_ui(_Complex float x, integer n) { _Complex float 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; } #endif #ifdef _MSC_VER static _Dcomplex zpow_ui(_Dcomplex x, integer n) { _Dcomplex pow={1.0,0.0}; unsigned long int u; if(n != 0) { if(n < 0) n = -n, x._Val[0] = 1/x._Val[0], x._Val[1] =1/x._Val[1]; for(u = n; ; ) { if(u & 01) pow._Val[0] *= x._Val[0], pow._Val[1] *= x._Val[1]; if(u >>= 1) x._Val[0] *= x._Val[0], x._Val[1] *= x._Val[1]; else break; } } _Dcomplex p = {pow._Val[0], pow._Val[1]}; return p; } #else static _Complex double zpow_ui(_Complex double x, integer n) { _Complex 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; } #endif static integer pow_ii(integer x, integer n) { integer pow; unsigned long int u; if (n <= 0) { if (n == 0 || x == 1) pow = 1; else if (x != -1) pow = x == 0 ? 1/x : 0; else n = -n; } if ((n > 0) || !(n == 0 || x == 1 || x != -1)) { u = n; for(pow = 1; ; ) { if(u & 01) pow *= x; if(u >>= 1) x *= x; else break; } } return pow; } static integer dmaxloc_(double *w, integer s, integer e, integer *n) { double m; integer i, mi; for(m=w[s-1], mi=s, i=s+1; i<=e; i++) if (w[i-1]>m) mi=i ,m=w[i-1]; return mi-s+1; } static integer smaxloc_(float *w, integer s, integer e, integer *n) { float m; integer i, mi; for(m=w[s-1], mi=s, i=s+1; i<=e; i++) if (w[i-1]>m) mi=i ,m=w[i-1]; return mi-s+1; } static inline void cdotc_(complex *z, integer *n_, complex *x, integer *incx_, complex *y, integer *incy_) { integer n = *n_, incx = *incx_, incy = *incy_, i; #ifdef _MSC_VER _Fcomplex zdotc = {0.0, 0.0}; if (incx == 1 && incy == 1) { for (i=0;i \brief \b SLAMCHF77 deprecated */ /* =========== DOCUMENTATION =========== */ /* Online html documentation available at */ /* http://www.netlib.org/lapack/explore-html/ */ /* Definition: */ /* =========== */ /* REAL FUNCTION SLAMCH( CMACH ) */ /* CHARACTER CMACH */ /* > \par Purpose: */ /* ============= */ /* > */ /* > \verbatim */ /* > */ /* > SLAMCH determines single precision machine parameters. */ /* > \endverbatim */ /* Arguments: */ /* ========== */ /* > \param[in] CMACH */ /* > \verbatim */ /* > 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) */ /* > \endverbatim */ /* Authors: */ /* ======== */ /* > \author Univ. of Tennessee */ /* > \author Univ. of California Berkeley */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ /* > \date April 2012 */ /* > \ingroup auxOTHERauxiliary */ /* ===================================================================== */ real slamch_(char *cmach) { /* Initialized data */ static logical first = TRUE_; /* System generated locals */ integer i__1; real ret_val; /* Local variables */ static real base; integer beta; static real emin, prec, emax; integer imin, imax; logical lrnd; static real rmin, rmax, t; real rmach; extern logical lsame_(char *, char *); real small; static real sfmin; extern /* Subroutine */ int slamc2_(integer *, integer *, logical *, real *, integer *, real *, integer *, real *); integer it; static real 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) { 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_ */ /* *********************************************************************** */ /* > \brief \b SLAMC1 */ /* > \details */ /* > \b Purpose: */ /* > \verbatim */ /* > SLAMC1 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 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 */ static logical lrnd; real a, b, c__, f; static integer lbeta; real savec; static logical lieee1; real t1, t2; extern real slamc3_(real *, real *); static integer lt; real 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.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_ */ /* *********************************************************************** */ /* > \brief \b SLAMC2 */ /* > \details */ /* > \b Purpose: */ /* > \verbatim */ /* > SLAMC2 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 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; /* Local variables */ logical ieee; real half; logical lrnd; static real leps; real zero, a, b, c__; integer i__; 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 real slamc3_(real *, real *); extern /* Subroutine */ int slamc4_(integer *, real *, integer *), slamc5_(integer *, integer *, integer *, logical *, integer *, real *); static integer lt; integer ngnmin, ngpmin; real 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.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 = abs(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 = 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 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_ */ /* *********************************************************************** */ /* > \brief \b SLAMC3 */ /* > \details */ /* > \b Purpose: */ /* > \verbatim */ /* > 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. */ /* > \endverbatim */ /* > */ /* > \param[in] A */ /* > */ /* > \param[in] B */ /* > \verbatim */ /* > The values A and B. */ /* > \endverbatim */ real slamc3_(real *a, real *b) { /* System generated locals */ real 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 SLAMC3 */ } /* slamc3_ */ /* *********************************************************************** */ /* > \brief \b SLAMC4 */ /* > \details */ /* > \b Purpose: */ /* > \verbatim */ /* > SLAMC4 is a service routine for SLAMC2. */ /* > \endverbatim */ /* > */ /* > \param[out] EMIN */ /* > \verbatim */ /* > The minimum exponent before (gradual) underflow, computed by */ /* > setting A = START and dividing by BASE until the previous A */ /* > can not be recovered. */ /* > \endverbatim */ /* > */ /* > \param[in] START */ /* > \verbatim */ /* > The starting point for determining EMIN. */ /* > \endverbatim */ /* > */ /* > \param[in] BASE */ /* > \verbatim */ /* > The base of the machine. */ /* > \endverbatim */ /* > */ /* Subroutine */ int slamc4_(integer *emin, real *start, integer *base) { /* System generated locals */ integer i__1; real r__1; /* Local variables */ real zero, a; integer i__; real rbase, b1, b2, c1, c2, d1, d2; extern real slamc3_(real *, real *); real one; /* -- LAPACK auxiliary routine (version 3.7.0) -- */ /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ /* November 2010 */ /* ===================================================================== */ 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_ */ /* *********************************************************************** */ /* > \brief \b SLAMC5 */ /* > \details */ /* > \b Purpose: */ /* > \verbatim */ /* > 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. */ /* > \endverbatim */ /* > */ /* > \param[in] BETA */ /* > \verbatim */ /* > The base of floating-point arithmetic. */ /* > \endverbatim */ /* > */ /* > \param[in] P */ /* > \verbatim */ /* > The number of base BETA digits in the mantissa of a */ /* > floating-point value. */ /* > \endverbatim */ /* > */ /* > \param[in] EMIN */ /* > \verbatim */ /* > The minimum exponent before (gradual) underflow. */ /* > \endverbatim */ /* > */ /* > \param[in] IEEE */ /* > \verbatim */ /* > A logical flag specifying whether or not the arithmetic */ /* > system is thought to comply with the IEEE standard. */ /* > \endverbatim */ /* > */ /* > \param[out] EMAX */ /* > \verbatim */ /* > The largest exponent before overflow */ /* > \endverbatim */ /* > */ /* > \param[out] RMAX */ /* > \verbatim */ /* > The largest machine floating-point number. */ /* > \endverbatim */ /* > */ /* 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 lexp; real oldy; integer uexp, i__; real y, z__; integer nbits; extern real slamc3_(real *, real *); real recbas; integer exbits, expsum, try__; /* -- LAPACK auxiliary routine (version 3.7.0) -- */ /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ /* November 2010 */ /* ===================================================================== */ /* 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_ */