From 1c81f50c5823983cdfba358341209be0154ce5ce Mon Sep 17 00:00:00 2001 From: Martin Kroeker Date: Tue, 22 Feb 2022 19:04:04 +0100 Subject: [PATCH] Add C equivalents as fallbacks --- lapack-netlib/INSTALL/Makefile | 31 + lapack-netlib/INSTALL/dlamch.c | 1379 ++++++++++++++++++++++ lapack-netlib/INSTALL/ilaver.c | 458 +++++++ lapack-netlib/INSTALL/lsame.c | 521 ++++++++ lapack-netlib/INSTALL/lsametst.c | 595 ++++++++++ lapack-netlib/INSTALL/second_INT_ETIME.c | 445 +++++++ lapack-netlib/INSTALL/second_NONE.c | 445 +++++++ lapack-netlib/INSTALL/secondtst.c | 566 +++++++++ lapack-netlib/INSTALL/slamch.c | 1378 +++++++++++++++++++++ 9 files changed, 5818 insertions(+) create mode 100644 lapack-netlib/INSTALL/dlamch.c create mode 100644 lapack-netlib/INSTALL/ilaver.c create mode 100644 lapack-netlib/INSTALL/lsame.c create mode 100644 lapack-netlib/INSTALL/lsametst.c create mode 100644 lapack-netlib/INSTALL/second_INT_ETIME.c create mode 100644 lapack-netlib/INSTALL/second_NONE.c create mode 100644 lapack-netlib/INSTALL/secondtst.c create mode 100644 lapack-netlib/INSTALL/slamch.c diff --git a/lapack-netlib/INSTALL/Makefile b/lapack-netlib/INSTALL/Makefile index 1007c1bca..1c86bd77d 100644 --- a/lapack-netlib/INSTALL/Makefile +++ b/lapack-netlib/INSTALL/Makefile @@ -4,6 +4,7 @@ include $(TOPSRCDIR)/make.inc .PHONY: all testlsame testslamch testdlamch testsecond testdsecnd testieee testversion all: testlsame testslamch testdlamch testsecond testdsecnd testieee testversion +ifneq ($(C_LAPACK), 1) testlsame: lsame.o lsametst.o $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^ @@ -27,6 +28,31 @@ testieee: tstiee.o testversion: ilaver.o LAPACK_version.o $(FC) $(FFLAGS) $(LDFLAGS) -o $@ $^ +else +testlsame: lsame.o lsametst.o + $(CC) -O2 -o $@ $^ + +testslamch: slamch.o lsame.o slamchtst.o + $(CC) $(CFLAGS) $(LDFLAGS) -o $@ $^ + +testdlamch: dlamch.o lsame.o dlamchtst.o + $(CC) $(CFLAGS) $(LDFLAGS) -o $@ $^ + +testsecond: second_$(TIMER).o secondtst.o + @echo "[INFO] : TIMER value: $(TIMER) (given by make.inc)" + $(CC) $(CFLAGS) $(LDFLAGS) -o $@ $^ + +testdsecnd: dsecnd_$(TIMER).o dsecndtst.o + @echo "[INFO] : TIMER value: $(TIMER) (given by make.inc)" + $(CC) $(CFLAGS) $(LDFLAGS) -o $@ $^ + +testieee: tstiee.o + $(CC) $(CFLAGS) $(LDFLAGS) -o $@ $^ + +testversion: ilaver.o LAPACK_version.o + $(CC) $(CFLAGS) $(LDFLAGS) -o $@ $^ +endif + .PHONY: run run: all ./testlsame @@ -46,5 +72,10 @@ cleanexe: cleantest: rm -f core +ifneq ($(C_LAPACK), 1) slamch.o: slamch.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $< dlamch.o: dlamch.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $< +else +slamch.o: slamch.c ; $(CC) $(CFLAGS) -c -o $@ $< +dlamch.o: dlamch.c ; $(CC) $(CFLAGS) -c -o $@ $< +endif diff --git a/lapack-netlib/INSTALL/dlamch.c b/lapack-netlib/INSTALL/dlamch.c new file mode 100644 index 000000000..e91a11e71 --- /dev/null +++ b/lapack-netlib/INSTALL/dlamch.c @@ -0,0 +1,1379 @@ +/* f2c.h -- Standard Fortran to C header file */ + +/** barf [ba:rf] 2. "He suggested using FORTRAN, and everybody barfed." + + - From The Shogakukan DICTIONARY OF NEW ENGLISH (Second edition) */ + +#ifndef F2C_INCLUDE +#define F2C_INCLUDE + +#include +#include +#include +#include +#include +#ifdef complex +#undef complex +#endif +#ifdef I +#undef I +#endif + +typedef int 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; +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;} +#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)); } +#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);} +#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);} +#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) = conj(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) (cimag(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; +} +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; +} +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; +} +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; + _Complex float zdotc = 0.0; + if (incx == 1 && incy == 1) { + for (i=0;i \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, <, &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, <, &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_ */ + + +/* *********************************************************************** */ + +/* > \brief \b DLAMC4 */ +/* > \details */ +/* > \b Purpose: */ +/* > \verbatim */ +/* > DLAMC4 is a service routine for DLAMC2. */ +/* > \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 dlamc4_(integer *emin, doublereal *start, integer *base) +{ + /* System generated locals */ + integer i__1; + doublereal d__1; + + /* Local variables */ + doublereal zero, a; + integer i__; + doublereal rbase, b1, b2, c1, c2, d1, d2; + extern doublereal dlamc3_(doublereal *, doublereal *); + doublereal one; + + +/* -- LAPACK auxiliary routine (version 3.7.0) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ +/* November 2010 */ + +/* ===================================================================== */ + + + a = *start; + one = 1.; + rbase = one / *base; + zero = 0.; + *emin = 1; + d__1 = a * rbase; + b1 = dlamc3_(&d__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; + d__1 = a / *base; + b1 = dlamc3_(&d__1, &zero); + d__1 = b1 * *base; + c1 = dlamc3_(&d__1, &zero); + d1 = zero; + i__1 = *base; + for (i__ = 1; i__ <= i__1; ++i__) { + d1 += b1; +/* L20: */ + } + d__1 = a * rbase; + b2 = dlamc3_(&d__1, &zero); + d__1 = b2 / rbase; + c2 = dlamc3_(&d__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 DLAMC4 */ + +} /* dlamc4_ */ + + +/* *********************************************************************** */ + +/* > \brief \b DLAMC5 */ +/* > \details */ +/* > \b Purpose: */ +/* > \verbatim */ +/* > DLAMC5 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 dlamc5_(integer *beta, integer *p, integer *emin, + logical *ieee, integer *emax, doublereal *rmax) +{ + /* System generated locals */ + integer i__1; + doublereal d__1; + + /* Local variables */ + integer lexp; + doublereal oldy; + integer uexp, i__; + doublereal y, z__; + integer nbits; + extern doublereal dlamc3_(doublereal *, doublereal *); + doublereal 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. / *beta; + z__ = *beta - 1.; + y = 0.; + i__1 = *p; + for (i__ = 1; i__ <= i__1; ++i__) { + z__ *= recbas; + if (y < 1.) { + oldy = y; + } + y = dlamc3_(&y, &z__); +/* L20: */ + } + if (y >= 1.) { + y = oldy; + } + +/* Now multiply by BETA**EMAX to get RMAX. */ + + i__1 = *emax; + for (i__ = 1; i__ <= i__1; ++i__) { + d__1 = y * *beta; + y = dlamc3_(&d__1, &c_b32); +/* L30: */ + } + + *rmax = y; + return 0; + +/* End of DLAMC5 */ + +} /* dlamc5_ */ + diff --git a/lapack-netlib/INSTALL/ilaver.c b/lapack-netlib/INSTALL/ilaver.c new file mode 100644 index 000000000..6f8e6a9f5 --- /dev/null +++ b/lapack-netlib/INSTALL/ilaver.c @@ -0,0 +1,458 @@ +/* f2c.h -- Standard Fortran to C header file */ + +/** barf [ba:rf] 2. "He suggested using FORTRAN, and everybody barfed." + + - From The Shogakukan DICTIONARY OF NEW ENGLISH (Second edition) */ + +#ifndef F2C_INCLUDE +#define F2C_INCLUDE + +#include +#include +#include +#include +#include +#ifdef complex +#undef complex +#endif +#ifdef I +#undef I +#endif + +typedef int 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; +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;} +#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)); } +#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);} +#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);} +#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) = conj(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) (cimag(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; +} +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; +} +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; +} +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; + _Complex float zdotc = 0.0; + if (incx == 1 && incy == 1) { + for (i=0;i \brief \b ILAVER returns the LAPACK version. */ +/* * */ +/* =========== DOCUMENTATION =========== */ + +/* Online html documentation available at */ +/* http://www.netlib.org/lapack/explore-html/ */ + +/* Definition: */ +/* =========== */ + +/* SUBROUTINE ILAVER( VERS_MAJOR, VERS_MINOR, VERS_PATCH ) */ + +/* INTEGER VERS_MAJOR, VERS_MINOR, VERS_PATCH */ + + +/* > \par Purpose: */ +/* ============= */ +/* > */ +/* > \verbatim */ +/* > */ +/* > This subroutine returns the LAPACK version. */ +/* > \endverbatim */ + +/* Arguments: */ +/* ========== */ + +/* > \param[out] VERS_MAJOR */ +/* > VERS_MAJOR is INTEGER */ +/* > return the lapack major version */ +/* > */ +/* > \param[out] VERS_MINOR */ +/* > VERS_MINOR is INTEGER */ +/* > return the lapack minor version from the major version */ +/* > */ +/* > \param[out] VERS_PATCH */ +/* > VERS_PATCH is INTEGER */ +/* > return the lapack patch version from the minor version */ + +/* Authors: */ +/* ======== */ + +/* > \author Univ. of Tennessee */ +/* > \author Univ. of California Berkeley */ +/* > \author Univ. of Colorado Denver */ +/* > \author NAG Ltd. */ + +/* > \date November 2019 */ + +/* > \ingroup auxOTHERauxiliary */ + +/* ===================================================================== */ +/* Subroutine */ int ilaver_(integer *vers_major__, integer *vers_minor__, + integer *vers_patch__) +{ + +/* -- LAPACK computational routine -- */ +/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ +/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ + +/* ===================================================================== */ + +/* ===================================================================== */ + *vers_major__ = 3; + *vers_minor__ = 9; + *vers_patch__ = 0; +/* ===================================================================== */ + + return 0; +} /* ilaver_ */ + diff --git a/lapack-netlib/INSTALL/lsame.c b/lapack-netlib/INSTALL/lsame.c new file mode 100644 index 000000000..b7ae8d447 --- /dev/null +++ b/lapack-netlib/INSTALL/lsame.c @@ -0,0 +1,521 @@ +/* f2c.h -- Standard Fortran to C header file */ + +/** barf [ba:rf] 2. "He suggested using FORTRAN, and everybody barfed." + + - From The Shogakukan DICTIONARY OF NEW ENGLISH (Second edition) */ + +#ifndef F2C_INCLUDE +#define F2C_INCLUDE + +#include +#include +#include +#include +#include +#ifdef complex +#undef complex +#endif +#ifdef I +#undef I +#endif + +typedef int 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; +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;} +#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)); } +#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);} +#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);} +#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) = conj(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) (cimag(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; +} +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; +} +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; +} +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; + _Complex float zdotc = 0.0; + if (incx == 1 && incy == 1) { + for (i=0;i \brief \b LSAME */ + +/* =========== DOCUMENTATION =========== */ + +/* Online html documentation available at */ +/* http://www.netlib.org/lapack/explore-html/ */ + +/* Definition: */ +/* =========== */ + +/* LOGICAL FUNCTION LSAME( CA, CB ) */ + +/* CHARACTER CA, CB */ + + +/* > \par Purpose: */ +/* ============= */ +/* > */ +/* > \verbatim */ +/* > */ +/* > LSAME returns .TRUE. if CA is the same letter as CB regardless of */ +/* > case. */ +/* > \endverbatim */ + +/* Arguments: */ +/* ========== */ + +/* > \param[in] CA */ +/* > \verbatim */ +/* > \endverbatim */ +/* > */ +/* > \param[in] CB */ +/* > \verbatim */ +/* > CA and CB specify the single characters to be compared. */ +/* > \endverbatim */ + +/* Authors: */ +/* ======== */ + +/* > \author Univ. of Tennessee */ +/* > \author Univ. of California Berkeley */ +/* > \author Univ. of Colorado Denver */ +/* > \author NAG Ltd. */ + +/* > \date December 2016 */ + +/* > \ingroup auxOTHERauxiliary */ + +/* ===================================================================== */ +logical lsame_(char *ca, char *cb) +{ + /* System generated locals */ + logical ret_val; + + /* Local variables */ + integer inta, intb, zcode; + + +/* -- 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..-- */ +/* December 2016 */ + + +/* ===================================================================== */ + + +/* Test if the characters are equal */ + + ret_val = *(unsigned char *)ca == *(unsigned char *)cb; + if (ret_val) { + return ret_val; + } + +/* Now test for equivalence if both characters are alphabetic. */ + + zcode = 'Z'; + +/* Use 'Z' rather than 'A' so that ASCII can be detected on Prime */ +/* machines, on which ICHAR returns a value with bit 8 set. */ +/* ICHAR('A') on Prime machines returns 193 which is the same as */ +/* ICHAR('A') on an EBCDIC machine. */ + + inta = *(unsigned char *)ca; + intb = *(unsigned char *)cb; + + if (zcode == 90 || zcode == 122) { + +/* ASCII is assumed - ZCODE is the ASCII code of either lower or */ +/* upper case 'Z'. */ + + if (inta >= 97 && inta <= 122) { + inta += -32; + } + if (intb >= 97 && intb <= 122) { + intb += -32; + } + + } else if (zcode == 233 || zcode == 169) { + +/* EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or */ +/* upper case 'Z'. */ + + if (inta >= 129 && inta <= 137 || inta >= 145 && inta <= 153 || inta + >= 162 && inta <= 169) { + inta += 64; + } + if (intb >= 129 && intb <= 137 || intb >= 145 && intb <= 153 || intb + >= 162 && intb <= 169) { + intb += 64; + } + + } else if (zcode == 218 || zcode == 250) { + +/* ASCII is assumed, on Prime machines - ZCODE is the ASCII code */ +/* plus 128 of either lower or upper case 'Z'. */ + + if (inta >= 225 && inta <= 250) { + inta += -32; + } + if (intb >= 225 && intb <= 250) { + intb += -32; + } + } + ret_val = inta == intb; + +/* RETURN */ + +/* End of LSAME */ + + return ret_val; +} /* lsame_ */ + diff --git a/lapack-netlib/INSTALL/lsametst.c b/lapack-netlib/INSTALL/lsametst.c new file mode 100644 index 000000000..4b46115fc --- /dev/null +++ b/lapack-netlib/INSTALL/lsametst.c @@ -0,0 +1,595 @@ +/* f2c.h -- Standard Fortran to C header file */ + +/** barf [ba:rf] 2. "He suggested using FORTRAN, and everybody barfed." + + - From The Shogakukan DICTIONARY OF NEW ENGLISH (Second edition) */ + +#ifndef F2C_INCLUDE +#define F2C_INCLUDE + +#include +#include +#include +#include +#include +#ifdef complex +#undef complex +#endif +#ifdef I +#undef I +#endif + +typedef int 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; +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;} +#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)); } +#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);} +#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);} +#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) = conj(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) (cimag(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; +} +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; +} +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; +} +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; + _Complex float zdotc = 0.0; + if (incx == 1 && incy == 1) { + for (i=0;i \brief \b LSAMETST */ + +/* =========== DOCUMENTATION =========== */ + +/* Online html documentation available at */ +/* http://www.netlib.org/lapack/explore-html/ */ + +/* Definition: */ +/* =========== */ + +/* PROGRAM LSAMETST */ + +/* Authors: */ +/* ======== */ + +/* > \author Univ. of Tennessee */ +/* > \author Univ. of California Berkeley */ +/* > \author Univ. of Colorado Denver */ +/* > \author NAG Ltd. */ + +/* > \date December 2016 */ + +/* > \ingroup auxOTHERauxiliary */ + +/* ===================================================================== PROGRAM LSAMETST */ + +/* -- LAPACK test routine (version 3.7.0) -- */ + +/* -- LAPACK computational 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..-- */ +/* December 2016 */ + +/* ===================================================================== */ +/* Main program */ main(void) +{ + /* Format strings */ + static char fmt_9999[] = "(\002 *** Error: LSAME( \002,a1,\002, \002," + "a1,\002) is .FALSE.\002)"; + static char fmt_9998[] = "(\002 *** Error: LSAME( \002,a1,\002, \002," + "a1,\002) is .TRUE.\002)"; + + /* System generated locals */ + integer i__1; + + /* Local variables */ + extern logical lsame_(char *, char *); + integer i1, i2; + + /* Fortran I/O blocks */ + static cilist io___3 = { 0, 6, 0, 0, 0 }; + static cilist io___4 = { 0, 6, 0, 0, 0 }; + static cilist io___5 = { 0, 6, 0, fmt_9999, 0 }; + static cilist io___6 = { 0, 6, 0, fmt_9999, 0 }; + static cilist io___7 = { 0, 6, 0, fmt_9999, 0 }; + static cilist io___8 = { 0, 6, 0, fmt_9999, 0 }; + static cilist io___9 = { 0, 6, 0, fmt_9998, 0 }; + static cilist io___10 = { 0, 6, 0, fmt_9998, 0 }; + static cilist io___11 = { 0, 6, 0, fmt_9998, 0 }; + static cilist io___12 = { 0, 6, 0, fmt_9998, 0 }; + static cilist io___13 = { 0, 6, 0, fmt_9998, 0 }; + static cilist io___14 = { 0, 6, 0, fmt_9998, 0 }; + static cilist io___15 = { 0, 6, 0, fmt_9998, 0 }; + static cilist io___16 = { 0, 6, 0, fmt_9998, 0 }; + static cilist io___17 = { 0, 6, 0, 0, 0 }; + + + + +/* Determine the character set. */ + + i1 = 'A'; + i2 = 'a'; + if (i2 - i1 == 32) { +/* + s_wsle(&io___3); + do_lio(&c__9, &c__1, " ASCII character set", (ftnlen)20); + e_wsle(); +*/ + printf(" ASCII character set"); +} else { + + printf(" Non-ASCII character set, IOFF should be %d",i2-i1); +/* + s_wsle(&io___4); + do_lio(&c__9, &c__1, " Non-ASCII character set, IOFF should be ", ( + ftnlen)41); + i__1 = i2 - i1; + do_lio(&c__3, &c__1, (char *)&i__1, (ftnlen)sizeof(integer)); + e_wsle(); +*/ + } + +/* Test LSAME. */ + + if (! lsame_("A", "A")) { +printf(" *** Error: LSAME(A,A) is .FALSE.\n"); +/* s_wsfe(&io___5); + do_fio(&c__1, "A", (ftnlen)1); + do_fio(&c__1, "A", (ftnlen)1); + e_wsfe(); +*/ + } + if (! lsame_("A", "a")) { +printf(" *** Error: LSAME(A,a) is .FALSE.\n"); +/* + s_wsfe(&io___6); + do_fio(&c__1, "A", (ftnlen)1); + do_fio(&c__1, "a", (ftnlen)1); + e_wsfe(); +*/ + } + if (! lsame_("a", "A")) { +printf(" *** Error: LSAME(a,A) is .FALSE.\n"); +/* s_wsfe(&io___7); + do_fio(&c__1, "a", (ftnlen)1); + do_fio(&c__1, "A", (ftnlen)1); + e_wsfe(); +*/ + } + if (! lsame_("a", "a")) { +printf(" *** Error: LSAME(a,a) is .FALSE.\n"); +/* s_wsfe(&io___8); + do_fio(&c__1, "a", (ftnlen)1); + do_fio(&c__1, "a", (ftnlen)1); + e_wsfe(); +*/ + } + if (lsame_("A", "B")) { +printf(" *** Error: LSAME(A,B) is .TRUE.\n"); +/* s_wsfe(&io___9); + do_fio(&c__1, "A", (ftnlen)1); + do_fio(&c__1, "B", (ftnlen)1); + e_wsfe(); +*/ + } + if (lsame_("A", "b")) { +printf(" *** Error: LSAME(A,b) is .TRUE.\n"); +/* s_wsfe(&io___10); + do_fio(&c__1, "A", (ftnlen)1); + do_fio(&c__1, "b", (ftnlen)1); + e_wsfe(); +*/ + } + if (lsame_("a", "B")) { +printf(" *** Error: LSAME(a,B) is .TRUE.\n"); +/* s_wsfe(&io___11); + do_fio(&c__1, "a", (ftnlen)1); + do_fio(&c__1, "B", (ftnlen)1); + e_wsfe(); +*/ + } + if (lsame_("a", "b")) { +printf(" *** Error: LSAME(a,b) is .TRUE.\n"); +/* s_wsfe(&io___12); + do_fio(&c__1, "a", (ftnlen)1); + do_fio(&c__1, "b", (ftnlen)1); + e_wsfe(); +*/ + } + if (lsame_("O", "/")) { +printf(" *** Error: LSAME(O,/) is .TRUE.\n"); +/* s_wsfe(&io___13); + do_fio(&c__1, "O", (ftnlen)1); + do_fio(&c__1, "/", (ftnlen)1); + e_wsfe(); +*/ + } + if (lsame_("/", "O")) { +printf(" *** Error: LSAME(/,O) is .TRUE.\n"); +/* s_wsfe(&io___14); + do_fio(&c__1, "/", (ftnlen)1); + do_fio(&c__1, "O", (ftnlen)1); + e_wsfe(); +*/ + } + if (lsame_("o", "/")) { +printf(" *** Error: LSAME(o,/) is .TRUE.\n"); +/* s_wsfe(&io___15); + do_fio(&c__1, "o", (ftnlen)1); + do_fio(&c__1, "/", (ftnlen)1); + e_wsfe(); +*/ + } + if (lsame_("/", "o")) { +printf(" *** Error: LSAME(/,o) is .TRUE.\n"); +/* s_wsfe(&io___16); + do_fio(&c__1, "/", (ftnlen)1); + do_fio(&c__1, "o", (ftnlen)1); + e_wsfe(); +*/ + } +printf(" Tests completed"); + +/* s_wsle(&io___17); + do_lio(&c__9, &c__1, " Tests completed", (ftnlen)16); + e_wsle(); +*/ + return 0; +} /* MAIN__ */ + diff --git a/lapack-netlib/INSTALL/second_INT_ETIME.c b/lapack-netlib/INSTALL/second_INT_ETIME.c new file mode 100644 index 000000000..51b7ddb54 --- /dev/null +++ b/lapack-netlib/INSTALL/second_INT_ETIME.c @@ -0,0 +1,445 @@ +/* f2c.h -- Standard Fortran to C header file */ + +/** barf [ba:rf] 2. "He suggested using FORTRAN, and everybody barfed." + + - From The Shogakukan DICTIONARY OF NEW ENGLISH (Second edition) */ + +#ifndef F2C_INCLUDE +#define F2C_INCLUDE + +#include +#include +#include +#include +#include +#ifdef complex +#undef complex +#endif +#ifdef I +#undef I +#endif + +typedef int 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; +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;} +#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)); } +#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);} +#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);} +#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) = conj(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) (cimag(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; +} +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; +} +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; +} +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; + _Complex float zdotc = 0.0; + if (incx == 1 && incy == 1) { + for (i=0;i \brief \b SECOND returns nothing */ + +/* =========== DOCUMENTATION =========== */ + +/* Online html documentation available at */ +/* http://www.netlib.org/lapack/explore-html/ */ + +/* Definition: */ +/* =========== */ + +/* REAL FUNCTION SECOND( ) */ + + +/* > \par Purpose: */ +/* ============= */ +/* > */ +/* > \verbatim */ +/* > */ +/* > SECOND returns nothing instead of returning the user time for a process in seconds. */ +/* > If you are using that routine, it means that neither EXTERNAL ETIME, */ +/* > EXTERNAL ETIME_, INTERNAL ETIME, INTERNAL CPU_TIME is available on */ +/* > your machine. */ +/* > \endverbatim */ + +/* Authors: */ +/* ======== */ + +/* > \author Univ. of Tennessee */ +/* > \author Univ. of California Berkeley */ +/* > \author Univ. of Colorado Denver */ +/* > \author NAG Ltd. */ + +/* > \date December 2016 */ + +/* > \ingroup auxOTHERauxiliary */ + +/* ===================================================================== */ +real second_(void) +{ + /* System generated locals */ + real ret_val; + + +/* -- 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..-- */ +/* December 2016 */ + +/* ===================================================================== */ + + ret_val = 0.f; + return ret_val; + +/* End of SECOND */ + +} /* second_ */ + diff --git a/lapack-netlib/INSTALL/second_NONE.c b/lapack-netlib/INSTALL/second_NONE.c new file mode 100644 index 000000000..51b7ddb54 --- /dev/null +++ b/lapack-netlib/INSTALL/second_NONE.c @@ -0,0 +1,445 @@ +/* f2c.h -- Standard Fortran to C header file */ + +/** barf [ba:rf] 2. "He suggested using FORTRAN, and everybody barfed." + + - From The Shogakukan DICTIONARY OF NEW ENGLISH (Second edition) */ + +#ifndef F2C_INCLUDE +#define F2C_INCLUDE + +#include +#include +#include +#include +#include +#ifdef complex +#undef complex +#endif +#ifdef I +#undef I +#endif + +typedef int 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; +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;} +#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)); } +#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);} +#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);} +#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) = conj(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) (cimag(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; +} +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; +} +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; +} +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; + _Complex float zdotc = 0.0; + if (incx == 1 && incy == 1) { + for (i=0;i \brief \b SECOND returns nothing */ + +/* =========== DOCUMENTATION =========== */ + +/* Online html documentation available at */ +/* http://www.netlib.org/lapack/explore-html/ */ + +/* Definition: */ +/* =========== */ + +/* REAL FUNCTION SECOND( ) */ + + +/* > \par Purpose: */ +/* ============= */ +/* > */ +/* > \verbatim */ +/* > */ +/* > SECOND returns nothing instead of returning the user time for a process in seconds. */ +/* > If you are using that routine, it means that neither EXTERNAL ETIME, */ +/* > EXTERNAL ETIME_, INTERNAL ETIME, INTERNAL CPU_TIME is available on */ +/* > your machine. */ +/* > \endverbatim */ + +/* Authors: */ +/* ======== */ + +/* > \author Univ. of Tennessee */ +/* > \author Univ. of California Berkeley */ +/* > \author Univ. of Colorado Denver */ +/* > \author NAG Ltd. */ + +/* > \date December 2016 */ + +/* > \ingroup auxOTHERauxiliary */ + +/* ===================================================================== */ +real second_(void) +{ + /* System generated locals */ + real ret_val; + + +/* -- 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..-- */ +/* December 2016 */ + +/* ===================================================================== */ + + ret_val = 0.f; + return ret_val; + +/* End of SECOND */ + +} /* second_ */ + diff --git a/lapack-netlib/INSTALL/secondtst.c b/lapack-netlib/INSTALL/secondtst.c new file mode 100644 index 000000000..694679bb5 --- /dev/null +++ b/lapack-netlib/INSTALL/secondtst.c @@ -0,0 +1,566 @@ +/* f2c.h -- Standard Fortran to C header file */ + +/** barf [ba:rf] 2. "He suggested using FORTRAN, and everybody barfed." + + - From The Shogakukan DICTIONARY OF NEW ENGLISH (Second edition) */ + +#ifndef F2C_INCLUDE +#define F2C_INCLUDE + +#include +#include +#include +#include +#include +#ifdef complex +#undef complex +#endif +#ifdef I +#undef I +#endif + +typedef int 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; +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;} +#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)); } +#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);} +#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);} +#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) = conj(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) (cimag(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; +} +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; +} +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; +} +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; + _Complex float zdotc = 0.0; + if (incx == 1 && incy == 1) { + for (i=0;i \brief \b SECONDTST */ + +/* =========== DOCUMENTATION =========== */ + +/* Online html documentation available at */ +/* http://www.netlib.org/lapack/explore-html/ */ + + +/* Authors: */ +/* ======== */ + +/* > \author Univ. of Tennessee */ +/* > \author Univ. of California Berkeley */ +/* > \author Univ. of Colorado Denver */ +/* > \author NAG Ltd. */ + +/* > \date November 2017 */ + +/* > \ingroup auxOTHERcomputational */ + +/* ===================================================================== PROGRAM SECONDTST */ + +/* -- LAPACK test routine (version 3.8.0) -- */ + +/* -- LAPACK computational routine (version 3.8.0) -- */ +/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ +/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ +/* November 2017 */ + +/* ===================================================================== */ + +/* Main program */ main(void) +{ + /* Format strings */ + static char fmt_9999[] = "(\002 Time for \002,g10.3,\002 SAXPY ops = " + "\002,g10.3,\002 seconds\002)"; + static char fmt_9998[] = "(\002 SAXPY performance rate = \002,g10" + ".3,\002 mflops \002)"; + static char fmt_9994[] = "(\002 *** Warning: Time for operations was le" + "ss or equal\002,\002 than zero => timing in TESTING might be dub" + "ious\002)"; + static char fmt_9997[] = "(\002 Including SECOND, time = \002,g10" + ".3,\002 seconds\002)"; + static char fmt_9996[] = "(\002 Average time for SECOND = \002,g10" + ".3,\002 milliseconds\002)"; + static char fmt_9995[] = "(\002 Equivalent floating point ops = \002,g10" + ".3,\002 ops\002)"; + + /* System generated locals */ + real r__1; + + /* Local variables */ + integer i__, j; + real alpha, x[1000], y[1000], total; + extern /* Subroutine */ int mysub_(integer *, real *, real *); + real t1, t2; + extern real second_(void); + real tnosec, avg; + + /* Fortran I/O blocks */ + static cilist io___10 = { 0, 6, 0, fmt_9999, 0 }; + static cilist io___11 = { 0, 6, 0, fmt_9998, 0 }; + static cilist io___12 = { 0, 6, 0, fmt_9994, 0 }; + static cilist io___13 = { 0, 6, 0, fmt_9997, 0 }; + static cilist io___15 = { 0, 6, 0, fmt_9996, 0 }; + static cilist io___16 = { 0, 6, 0, fmt_9995, 0 }; + + + + total = 1e8f; + +/* Initialize X and Y */ + + for (i__ = 1; i__ <= 1000; ++i__) { + x[i__ - 1] = 1.f / (real) i__; + y[i__ - 1] = (real) (1000 - i__) / 1e3f; +/* L10: */ + } + alpha = .315f; + +/* Time TOTAL SAXPY operations */ + + t1 = second_(); + for (j = 1; j <= 50000; ++j) { + for (i__ = 1; i__ <= 1000; ++i__) { + y[i__ - 1] += alpha * x[i__ - 1]; +/* L20: */ + } + alpha = -alpha; +/* L30: */ + } + t2 = second_(); + tnosec = t2 - t1; +/* + s_wsfe(&io___10); + do_fio(&c__1, (char *)&total, (ftnlen)sizeof(real)); + do_fio(&c__1, (char *)&tnosec, (ftnlen)sizeof(real)); + e_wsfe(); + if (tnosec > 0.f) { + s_wsfe(&io___11); + r__1 = total / 1e6f / tnosec; + do_fio(&c__1, (char *)&r__1, (ftnlen)sizeof(real)); + e_wsfe(); + } else { + s_wsfe(&io___12); + e_wsfe(); + } +*/ + printf("Time for %f10.3 SAXPY ops = %f10.3 seconds\n",total,tnosec); + if (tnosec > 0.f) { + printf("SAXPY performance rate = %f10.3 mflops\n",total/1.e6/tnosec ); + } else { + printf("*** Warning: Time for operations was less or equal than zero => timing in TESTING might be dubious\n" ); + } +/* Time TOTAL SAXPY operations with SECOND in the outer loop */ + + t1 = second_(); + for (j = 1; j <= 50000; ++j) { + for (i__ = 1; i__ <= 1000; ++i__) { + y[i__ - 1] += alpha * x[i__ - 1]; +/* L40: */ + } + alpha = -alpha; + t2 = second_(); +/* L50: */ + } + +/* Compute the time used in milliseconds used by an average call */ +/* to SECOND. */ +/* + s_wsfe(&io___13); + r__1 = t2 - t1; + do_fio(&c__1, (char *)&r__1, (ftnlen)sizeof(real)); + e_wsfe(); +*/ + printf("Including SECOND, time = %f10.3 seconds\n",t2-t1); + avg = (t2 - t1 - tnosec) * 1e3f / 5e4f; + if (avg > 0.f) { + printf("Average time for SECOND = %f10.3 milliseconds\n",avg ); + + /* + s_wsfe(&io___15); + do_fio(&c__1, (char *)&avg, (ftnlen)sizeof(real)); + e_wsfe(); +*/ + } + +/* Compute the equivalent number of floating point operations used */ +/* by an average call to SECOND. */ + + if (avg > 0.f && tnosec > 0.f) { +printf("Equivalent floating point ops = %f10.3 ops\n", avg/1000*total/tnosec); +/* s_wsfe(&io___16); + r__1 = avg / 1000 * total / tnosec; + do_fio(&c__1, (char *)&r__1, (ftnlen)sizeof(real)); + e_wsfe(); +*/ + } + + mysub_(&c__1000, x, y); + return 0; +} /* MAIN__ */ + +/* Subroutine */ int mysub_(integer *n, real *x, real *y) +{ + /* Parameter adjustments */ + --y; + --x; + + /* Function Body */ + return 0; +} /* mysub_ */ + diff --git a/lapack-netlib/INSTALL/slamch.c b/lapack-netlib/INSTALL/slamch.c new file mode 100644 index 000000000..734de03e5 --- /dev/null +++ b/lapack-netlib/INSTALL/slamch.c @@ -0,0 +1,1378 @@ +/* f2c.h -- Standard Fortran to C header file */ + +/** barf [ba:rf] 2. "He suggested using FORTRAN, and everybody barfed." + + - From The Shogakukan DICTIONARY OF NEW ENGLISH (Second edition) */ + +#ifndef F2C_INCLUDE +#define F2C_INCLUDE + +#include +#include +#include +#include +#include +#ifdef complex +#undef complex +#endif +#ifdef I +#undef I +#endif + +typedef int 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; +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;} +#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)); } +#define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);} +#define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);} +#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) = conj(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) (cimag(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; +} +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; +} +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; +} +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; + _Complex float zdotc = 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_ */ +