/* 

        Transverse lattice Global definitions

These are definitions relevant to diagonalizing the
Hamiltonian, fitting to "experimental" data, extrapolating
to the continuum, etc.  This header gives flags and definitions for
various external libraries.

The main function of this header is to provide an
interface between FORTRAN and C routines for various
platforms. 

*/

/* Number of digits of accuracy in Lanczos when doing 
   extraplations.  I have normally been using 8.           */

#define EDIGITS 9
#define PDIGITS "%.10g" /* number of digits to print out */

/* ---------Definitions to interface with FORTRAN routines-------*/
/*  Some of the following definitions are borrowed from the f2c  */
/*   header file (with various obscure changes.                  */

typedef double doublereal; 
typedef struct { doublereal r, i; } doublecomplex;

/*  Some of the definitions are Machine dependant               */
/*  Define a data type for fortran subroutines                  */
/*  The IBM SP2 and CRAY FORTRAN subroutines return void        */

#ifdef IBM_SP2  /* IBM SP2 at Argonne, for example */
#define SUBRETURN
typedef void fsub;
typedef int ftnlen;
typedef long logical;
typedef long integer; 

#elif defined(_CRAY) /* Cray machines */
#include <fortran.h> /* header file that includes conversion stuff */
#define SUBRETURN 
typedef void fsub;
typedef int logical;
typedef long integer; 

#elif  defined(__alpha) /* OSU DEC machines */
#define SUBRETURN 0
typedef int fsub;
typedef int ftnlen;
typedef int logical;
typedef int integer; 

#elif defined(BABBAGE)  /* Babbage at Cambridge (Hitachi) */
#define SUBRETURN
typedef void fsub;
typedef int ftnlen;
typedef int logical;
typedef int integer; 
extern int hf_fint(char *);  /* this must be called once at the start */

#else  /* everything else (especially GNU/LINUX) */
#define SUBRETURN 0
typedef int fsub;
typedef long ftnlen;
typedef long logical;
typedef long integer; 

#endif

/*  function names to export to fortran programs */

#ifdef __uxp__   /*  Fujitsu vpp */
#define MAIN MAIN__
typedef int /* Unknown procedure type (function name) */ (**U_fp)();
typedef /* Subroutine from FORTRAN*/ int (**S_fp)();
typedef doublereal (**D_fp)();

#else  /*      This is the more standard case                    */
#define MAIN main
typedef fsub (*U_fp)(); /* Unknown procedure type (function name) */
typedef /* Subroutine */ fsub (*S_fp)();
typedef doublereal (*D_fp)();
#endif


/* --------------Define some FORTRAN subroutine names ---------------*/
/* Since the conversion of FORTRAN subroutine names is somewhat
   machine dependant, the names are all first defined as macros      */

#ifdef _CRAY /* on CRAY, change library routines to "single precision" */
             /* Also, define macros to handle strings correctly */
#define DGEQRF SGEQRF
#define DORGQR SORGQR
#define DPOTRF(C1,A2,A3,A4,A5,L1) do{_fcd c1; c1 = _cptofcd(C1,L1); \
     SPOTRF(c1,A2,A3,A4,A5,L1);}while(0)
#define DSYGST(A1,C2,A3,A4,A5,A6,A7,A8,L2) do{_fcd c2; \
     c2 = _cptofcd(C2,L2); SSYGST(A1,c2,A3,A4,A5,A6,A7,A8);}while(0)
#define DSYGV(A1,C2,C3,A4,A5,A6,A7,A8,A9,A10,A11,A12,L2,L3) \
     do{_fcd c2,c3; c2 = _cptofcd(C2,L2); c3 = _cptofcd(C3,L3); \
     SSYGV(A1,c2,c3,A4,A5,A6,A7,A8,A9,A10,A11,A12);}while(0)
#define DSYEV(C1,C2,A3,A4,A5,A6,A7,A8,A9,L1,L2) \
     do{_fcd c1,c2; c1 = _cptofcd(C1,L1); c2 = _cptofcd(C2,L2); \
     SSYEV(c1,c2,A3,A4,A5,A6,A7,A8,A9);}while(0)
#define ZHEEV(C1,C2,A3,A4,A5,A6,A7,A8,A9,A10,L1,L2) \
     do{_fcd c1,c2; c1 = _cptofcd(C1,L1); c2 = _cptofcd(C2,L2); \
     CHEEV(c1,c2,A3,A4,A5,A6,A7,A8,A9,A10);}while(0)

#elif  defined(__alpha) /* on DEC machines string lengths are not passed*/

#define DPOTRF(C1,A2,A3,A4,A5,L1) \
       dpotrf_(C1,A2,A3,A4,A5,L1)
#define DSYGST(A1,C2,A3,A4,A5,A6,A7,A8,L2) \
       dsygst_(A1,C2,A3,A4,A5,A6,A7,A8)
#define DSYGV(A1,C2,C3,A4,A5,A6,A7,A8,A9,A10,A11,A12,L2,L3) \
       dsygv_(A1,C2,C3,A4,A5,A6,A7,A8,A9,A10,A11,A12)
#define DSYEV(C1,C2,A3,A4,A5,A6,A7,A8,A9,L1,L2) \
       dsyev_(C1,C2,A3,A4,A5,A6,A7,A8,A9)
#define ZHEEV(C1,C2,A3,A4,A5,A6,A7,A8,A9,A10,L1,L2) \
       zheev_(C1,C2,A3,A4,A5,A6,A7,A8,A9,A10)
#define PSEUDOINVERSE pseudoinverse_
#define FINDC  findc_ 
#define MYGAMMA mygamma_
#define MYBESSELJ mybesselj_
#define PLANDR2 plandr2_
#define DNLASO dnlaso_
#define LMDIF  lmdif_
#define F02FJZ f02fjz_
#define F02FJF f02fjf_
#define DGAUS8 dgaus8_
#define GAUSSQ gaussq_
#define DGEQRF dgeqrf_
#define DORGQR dorgqr_
#define DSAUPD dsaupd_
#define DSEUPD dseupd_
#define FORELOP forelop_
#define OP op_
#define OPM opm_


#elif defined(BABBAGE) /* don't have to do anything */

#elif defined(__hpux) || defined(IBM_SP2)
#define PSEUDOINVERSE pseudoinverse
#define FINDC  findc 
#define MYGAMMA mygamma
#define MYBESSELJ mybesselj
#define ZHEEV  zheev
#define DSYEV  dsyev
#define PLANDR2 plandr2
#define DNLASO dnlaso
#define LMDIF  lmdif
#define F02FJZ f02fjz
#define F02FJF f02fjf
#define DGAUS8 dgaus8
#define GAUSSQ gaussq
#define DPOTRF dpotrf
#define DSYGST dsygst
#define DGEQRF dgeqrf
#define DORGQR dorgqr
#define DSYGV  dsygv
#define DSAUPD dsaupd
#define DSEUPD dseupd
#define FORELOP forelop
#define OP op
#define OPM opm

#else  /* Everything else  */

#define PSEUDOINVERSE pseudoinverse_
#define FINDC  findc_ 
#define MYGAMMA mygamma_
#define MYBESSELJ mybesselj_
#define ZHEEV  zheev_
#define DSYEV  dsyev_
#define PLANDR2 plandr2_
#define DNLASO dnlaso_
#define LMDIF  lmdif_
#define F02FJZ f02fjz_
#define F02FJF f02fjf_
#define DGAUS8 dgaus8_
#define GAUSSQ gaussq_
#define DPOTRF dpotrf_
#define DSYGST dsygst_
#define DGEQRF dgeqrf_
#define DORGQR dorgqr_
#define DSYGV  dsygv_
#define DSAUPD dsaupd_
#define DSEUPD dseupd_
#define FORELOP forelop_
#define OP op_
#define OPM opm_

#endif 


/* ---------------FORTRAN Function prototypes -------------------------- */
/*     These were generated by f2c (fortran to c converter)              */
/*       f2c -!c -g -Wall  -P  *.f                                       */
/*     This is the content of the relevant *.P files except that         */
/*     the function names are now replaced by machine dependant macros.  */

/*  from weights.f   */
extern doublereal FINDC(doublereal *, doublereal *, doublereal *);
extern doublereal MYGAMMA(doublereal *);
extern doublereal MYBESSELJ(doublereal *,doublereal *);
extern fsub PSEUDOINVERSE(doublereal *, integer *, integer *);
extern fsub FORELOP(integer *, integer *, doublereal *, doublereal *, 
		    doublereal *, integer *, integer *,integer *);


/*  The following prototypes are from LAPACK. */

extern fsub DGEQRF(integer *,integer *,doublereal *,integer *,
		  doublereal *,doublereal *,integer *,integer *);
extern fsub DORGQR(integer *,integer *,integer *,doublereal *,integer *,
		  doublereal *,doublereal *,integer *,integer *);
#ifdef _CRAY /* on the CRAY, strings are handled differently */
extern fsub SPOTRF(_fcd,integer *,doublereal *,integer *,integer *);
extern fsub SSYGST(integer *,_fcd,integer *,doublereal *,integer *,
		  doublereal *,integer *,integer *);
extern fsub SSYGV(integer *,_fcd,_fcd,integer *,doublereal *,
		 integer *,doublereal *,integer *,doublereal *,
		 doublereal *,integer *,integer *); 
extern fsub SSYEV(_fcd, _fcd, integer *, doublereal *, 
     integer *, doublereal *, doublereal *, integer *, integer *);
extern fsub CHEEV(_fcd, _fcd, integer *, doublecomplex *, 
     integer *, doublereal *, doublecomplex *, integer *, 
     doublereal *, integer *);
#else /* not CRAY case */
extern fsub DPOTRF(char *,integer *,doublereal *,integer *,integer *, ftnlen);
extern fsub DSYGST(integer *,char *,integer *,doublereal *,integer *,
		  doublereal *,integer *,integer *, ftnlen);
extern fsub DSYGV(integer *,char *,char *,integer *,doublereal *,
		 integer *,doublereal *,integer *,doublereal *,
		 doublereal *,integer *,integer *, ftnlen, ftnlen); 
extern fsub DSYEV(char *, char *, integer *, doublereal *, 
     integer *, doublereal *, doublereal *, integer *, integer *, 
		 ftnlen, ftnlen);
extern fsub ZHEEV(char *, char *, integer *, doublecomplex *, 
     integer *, doublereal *, doublecomplex *, integer *, 
     doublereal *, integer *, ftnlen, ftnlen);
#endif /* different for CRAY */

/* Lanczos Routines from ARPACK */

#ifdef _CRAY /* on the CRAY, strings are handled differently */

extern fsub DSAUPD(integer *ido, _fcd, integer *n, _fcd, 
		  integer *nev, doublereal *tol, doublereal *resid, 
		  integer *ncv, doublereal *v, integer *ldv, 
		  integer *iparam, integer *ipntr, doublereal *workd, 
		  doublereal *workl, integer *lworkl, integer *info);

#define DSAUPD(A1,C2,A3,C4,A5,A6,A7,A8,A9,A10,A11,A12,A13,A14,A15,A16,\
     L2,L4) do{_fcd c2,c4; c2=_cptofcd(C2,L2); c4=_cptofcd(C4,L4); \
     DSAUPD(A1,c2,A3,c4,A5,A6,A7,A8,A9,A10,A11,A12,A13,A14,A15,A16);}\
     while(0)

extern fsub DSEUPD(logical *rvec, _fcd, logical *select, 
		  doublereal *d__, doublereal *z__, integer *ldz, 
		  doublereal *sigma, _fcd, integer *n, _fcd, 
		  integer *nev, doublereal *tol, doublereal *resid, 
		  integer *ncv, doublereal *v, integer *ldv, 
		  integer *iparam, integer *ipntr, doublereal *workd, 
		  doublereal *workl, integer *lworkl, integer *info);

#define DSEUPD(A1,C2,A3,A4,A5,A6,A7,C8,A9,C10,A11,A12,A13,A14,A15,A16,\
     A17,A18,A19,A20,A21,A22,L2,L8,L10) do{_fcd c2,c8,c10; int a1,a3; \
     a1 = _btol(A1); c2 = _cptofcd(C2,L2); a3 = _btol(A3); \
     c8 = _cptofcd(C8,L8); c10 = _cptofcd(C10,L10); \
     DSEUPD(a1,c2,a3,c4,A5,A6,A7,A8,A9,A10,A11,A12,A13,A14,A15,A16,\
     A17,A18,A19,A20,A21,A22);}while(0)

#else /* different for CRAY */

extern fsub DSAUPD(integer *ido, char *bmat, integer *n, char *which, 
		  integer *nev, doublereal *tol, doublereal *resid, 
		  integer *ncv, doublereal *v, integer *ldv, 
		  integer *iparam, integer *ipntr, doublereal *workd, 
		  doublereal *workl, integer *lworkl, integer *info, 
		  ftnlen, ftnlen);

extern fsub DSEUPD(logical *rvec, char *howmny, logical *select, 
		  doublereal *d__, doublereal *z__, integer *ldz, 
		  doublereal *sigma, char *bmat, integer *n, char *which, 
		  integer *nev, doublereal *tol, doublereal *resid, 
		  integer *ncv, doublereal *v, integer *ldv, 
		  integer *iparam, integer *ipntr, doublereal *workd, 
		  doublereal *workl, integer *lworkl, integer *info, 
		  ftnlen, ftnlen, ftnlen);

#endif /* different for CRAY */

/*  Lanczos program dnlaso.f  */

extern fsub DNLASO(S_fp, S_fp, integer *n, integer *nval, 
  integer *nfig, integer *nperm, integer *nmval, doublereal *val, 
  integer *nmvec, doublereal *vec, integer *nblock, integer *maxop, 
  integer *maxj, doublereal *work, integer *ind, integer *ierr);

/*  Minpack least squares lmdif */

extern fsub LMDIF(S_fp, integer *m, integer *n, doublereal *x, 
  doublereal *fvec, doublereal *ftol, doublereal *xtol, doublereal *gtol, 
  integer *maxfev, doublereal *epsfcn, doublereal *diag, integer *mode, 
  doublereal *factor, integer *nprint, integer *info, integer *nfev, 
  doublereal *fjac, integer *ldfjac, integer *ipvt, doublereal *qtf, 
  doublereal *wa1, doublereal *wa2, doublereal *wa3, doublereal *wa4);

/* call to PLANSO function PLANDR2 */
 
extern fsub PLANDR2(integer *n, integer *lanmax, integer *maxprs, integer *lohi, 
                    doublereal *condm, doublereal *kappa, integer *j, integer *neig, 
                    doublereal *ritz, doublereal *bnd, doublereal *w, integer *nw, 
		     integer *ierr, integer *msglvl, integer *mpicom);

/*  Integration subroutine */

extern fsub DGAUS8(D_fp,doublereal *,doublereal *,
			 doublereal *,doublereal *,integer *);

/*  Gauss-Jacobi weights */

extern fsub GAUSSQ(integer *, integer *n, doublereal *alpha, doublereal *, 
   integer *kpts, doublereal *endpts, doublereal *b, doublereal *t, 
		  doublereal *w);

/*  NAG library routines */
#if 0

extern fsub F02FJZ(void);
extern fsub F02FJF(integer *,integer *,integer *,integer *,doublereal *,
	     D_fp,U_fp,U_fp, integer *,doublereal *,integer *,
             doublereal *, doublereal *, integer *,
             doublereal *,integer *,integer *,integer *,integer *);

#endif

/* ----------------------- Function prototypes -------------------------- */

/* from file: extralaso.c */

void spectrum(sectors *s,         /* basis, filled by initbasis(...) */ 
	      sectors *s2,        /* second basis for nonzero momentum */
	      doublereal **vals,  /* list of eigenvalues */
	      doublereal **vecin, /* list of eigenvectors */
	      integer nval,       /* number of eigenvalues */
	      element *couplings, /* coupling constants */
	      element *angle,     /* angles and momenta, length HINDEX+2 */
	      int globalkt,       /* maximum value of kt, for mapcouplings */
	      integer nfig);      /* number of significant digits in answer */

void printhels(FILE *); /* print memory usage of structure hels */
void freespectrum(void); /* free memory used in spectrum */
size_t max_els; /* maximum length of els used used so far */
void printlaso(FILE *); /* Print out the memory used by qstore */

/* functions concerned with PLANDR2*/
extern fsub OP(integer *n, doublereal *P, doublereal *Q, doublereal *R, 
	       integer *our_comm_world);  /*R Answer, P and Q are Given*/
 
extern fsub OPM(integer *n, doublereal *x, doublereal *y, integer *our_comm_world);
 /*OPM receives X and returns M*X as Y, but M is the identity matrix so*/

extern fsub elop(integer *, integer *, doublereal *, doublereal *);
element *extrapolate1(integer,sectors *,element *,
		      integer *,int,element *);
element *extrapolate(integer, sectors *, sectors *, integer *,
		     integer, element *,element *,int *,element *);
void kpweights(element *,integer *,int,sectors *,element);
          /* function to extrapolate in K and p-truncation */
element pextrapolate(integer, sectors *,element *);
element pextrapolate1(integer, sectors *,element *);
extern fsub fcn2(integer *, integer *, doublereal *,doublereal *,integer *);
element *lextrapolate(integer, sectors *, sectors *s, integer *,
		      integer, element *, element *,element *,element);
element *lextrapolate1(integer, sectors *,element *, integer *,
		      integer, element *,element *);
                   /* extrapolation procedure for nonzero L  */
void kmaxweights(element *,integer *,int,sectors *,element *,element);
          /* function to extrapolate in K, p-truncation, and K_max. */
void wweights(element *,integer *,int,int (*)[HINDEX]); 
                                     /* fit winding modes */ 
void lweights(element *,integer *,int,element *); 
                                  /*fit longitudinal potential */
element fdot(const element *,const element *,int); /* dot product */
int dot(const int *,const int *,int); /* dot product */

/* storage array used by iovect.  Memory must be allocated 
   before accessing iovect.  */

#if 0 /* this is now disabled */
doublereal *qstore=NULL;
#endif

/* from file: improved.c */

#if IMPROVED==3
element inner13(int,int,int,int,element,integer *);
#endif 
             /* integration routines for 1->3 improved matrix elements*/   
/* the following are internal functions to be
   called by FORTRAN routines */
extern doublereal f1(doublereal *);
extern doublereal f2(doublereal *);
extern doublereal f3(doublereal *);
extern doublereal f11(doublereal *);
extern doublereal f22(doublereal *);
extern doublereal f33(doublereal *);

extern fsub iovect(integer *, integer *, doublereal *, integer *, integer *);
extern fsub vectop(integer *, integer *, doublereal *, doublereal *);

/* from file jessica.c  */

extern fsub fcn(integer *, integer *, doublereal *,doublereal *,integer *);

/* ---------- Pass function names to FORTRAN (Machine dependant)------*/
/*  Machine dependant functions used for passing subroutine names
    to and from FORTRAN routines.              */

#if defined(__uxp__)   /*  Fujitsu vpp */ 

static int (*elop_point)()=&elop;
#define ELOP_TO_FORTRAN &elop_point
static int (*fcn_point)()=&fcn;
#define FCN_TO_FORTRAN &fcn_point
static int (*fcn2_point)()=&fcn2;
#define FCN2_TO_FORTRAN &fcn2_point
static int (*iovect_point)()=&iovect;
#define IOVECT_TO_FORTRAN &iovect_point
static int (*vectop_point)()=&vectop;
#define VECTOP_TO_FORTRAN &vectop_point
static doublereal (*selfinertia_point)()=&selfinertia;
#define SELFINERTIA_TO_FORTRAN &selfinertia_point
static doublereal (*seenorm_point)()=&seenorm;
#define SEENORM_TO_FORTRAN &seenorm_point
static doublereal (*vinst_point)()=&vinst;
#define VINST_TO_FORTRAN &vinst_point
static doublereal (*vtau_point)()=&vtau;
#define VTAU_TO_FORTRAN &vtau_point
static doublereal (*vnewtau_point)()=&vnewtau;
#define VNEWTAU_TO_FORTRAN &vnewtau_point
static doublereal (*f1_point)()=&f1;
#define F1_TO_FORTRAN &f1_point
static doublereal (*f2_point)()=&f2;
#define F2_TO_FORTRAN &f2_point
static doublereal (*f3_point)()=&f3;
#define F3_TO_FORTRAN &f3_point
static doublereal (*f11_point)()=&f11;
#define F11_TO_FORTRAN &f11_point
static doublereal (*f22_point)()=&f22;
#define F22_TO_FORTRAN &f1_point
static doublereal (*f33_point)()=&f33;
#define F33_TO_FORTRAN &f33_point

#elif defined(__hpux)  /* Some HP machines are more picky   */

#define DOT_TO_FORTRAN (D_fp) dot
#define ELOP_TO_FORTRAN (S_fp) elop
#define FCN_TO_FORTRAN (S_fp) fcn
#define FCN2_TO_FORTRAN (S_fp) fcn2
#define ELOP2_TO_FORTRAN (S_fp) elop2
#define IOVECT_TO_FORTRAN (S_fp) iovect
#define VECTOP_TO_FORTRAN (S_fp)vectop
#define MONIT_TO_FORTRAN (S_fp) monit
#define SELFINERTIA_TO_FORTRAN (D_fp) selfinertia
#define SEENORM_TO_FORTRAN (D_fp) seenorm
#define VINST_TO_FORTRAN (D_fp) vinst
#define VTAU_TO_FORTRAN (D_fp) vtau
#define VNEWTAU_TO_FORTRAN (D_fp) vnewtau
#define F1_TO_FORTRAN (D_fp) f1
#define F2_TO_FORTRAN (D_fp) f2
#define F3_TO_FORTRAN  (D_fp) f3
#define F11_TO_FORTRAN (D_fp) f11
#define F22_TO_FORTRAN (D_fp) f22
#define F33_TO_FORTRAN (D_fp) f33

#else  /*  This is the more standard case with gcc    */

#define DOT_TO_FORTRAN dot
#define ELOP_TO_FORTRAN elop
#define FCN_TO_FORTRAN fcn
#define FCN2_TO_FORTRAN fcn2
#define ELOP2_TO_FORTRAN elop2
#define IOVECT_TO_FORTRAN iovect
#define VECTOP_TO_FORTRAN vectop
#define MONIT_TO_FORTRAN monit
#define SELFINERTIA_TO_FORTRAN selfinertia
#define SEENORM_TO_FORTRAN seenorm
#define VINST_TO_FORTRAN vinst
#define VTAU_TO_FORTRAN vtau
#define VNEWTAU_TO_FORTRAN vnewtau
#define F1_TO_FORTRAN f1
#define F2_TO_FORTRAN f2
#define F3_TO_FORTRAN f3
#define F11_TO_FORTRAN f11
#define F22_TO_FORTRAN f22
#define F33_TO_FORTRAN f33

#endif