text:
#ifndef MKL_TEMPLATE
#define MKL_TEMPLATE

//#include <cblas.h>
#include <stddef.h>
//#include <blas.h>
#ifdef small
#undef small
#endif
//#include <lapack.h>
#include <cblas_defvar.h>

#ifdef NEW_MATLAB
   typedef ptrdiff_t INTT;
#else
   typedef int INTT;
#endif

/// a few static variables for lapack
static char low='l';
static char lower='L';
static char nonUnit='n';
static char upper='u';
static INTT info=0;
static char incr='I';
static char decr='D';
static char all='A';
static char no='N';
static char reduced='S';
static char allV='V';

#ifdef REMOVE_
#define dnrm2_ dnrm2
#define snrm2_ snrm2
#define dcopy_ dcopy
#define scopy_ scopy
#define daxpy_ daxpy
#define saxpy_ saxpy
#define dscal_ dscal
#define sscal_ sscal
#define dasum_ dasum
#define sasum_ sasum
#define ddot_ ddot
#define sdot_ sdot
#define dgemv_ dgemv
#define sgemv_ sgemv
#define dger_ dger
#define sger_ sger
#define dtrmv_ dtrmv
#define strmv_ strmv
#define dsyr_ dsyr
#define ssyr_ ssyr
#define dsymv_ dsymv
#define ssymv_ ssymv
#define dgemm_ dgemm
#define sgemm_ sgemm
#define dsyrk_ dsyrk
#define ssyrk_ ssyrk
#define dtrmm_ dtrmm
#define strmm_ strmm
#define dtrtri_ dtrtri
#define strtri_ strtri
#define idamax_ idamax
#define isamax_ isamax
#define dsytrf_ dsytrf
#define ssytrf_ ssytrf
#define dsytri_ dsytri
#define ssytri_ ssytri
#define dlasrt_ dlasrt
#define slasrt_ slasrt
#define dgesvd_ dgesvd
#define sgesvd_ sgesvd
#define dsyev_ dsyev
#define ssyev_ ssyev
#endif

/// external functions
#ifdef HAVE_MKL   // obsolete
extern "C" {
#endif
   size_t cblas_idamin( int n,  double* X,  int incX);
   size_t cblas_isamin( int n,  float* X,  int incX);
#ifdef HAVE_MKL
};
#endif

#ifdef HAVE_MKL
extern "C" {
   void vdSqr( int n,  double* vecIn, double* vecOut);
   void vsSqr( int n,  float* vecIn, float* vecOut);
   void vdSqrt( int n,  double* vecIn, double* vecOut);
   void vsSqrt( int n,  float* vecIn, float* vecOut);
   void vdInvSqrt( int n,  double* vecIn, double* vecOut);
   void vsInvSqrt( int n,  float* vecIn, float* vecOut);
   void vdSub( int n,  double* vecIn,  double* vecIn2, double* vecOut);
   void vsSub( int n,  float* vecIn,  float* vecIn2, float* vecOut);
   void vdDiv( int n,  double* vecIn,  double* vecIn2, double* vecOut);
   void vsDiv( int n,  float* vecIn,  float* vecIn2, float* vecOut);
   void vdExp( int n,  double* vecIn, double* vecOut);
   void vsExp( int n,  float* vecIn, float* vecOut);
   void vdInv( int n,  double* vecIn, double* vecOut);
   void vsInv( int n,  float* vecIn, float* vecOut);
   void vdAdd( int n,  double* vecIn,  double* vecIn2, double* vecOut);
   void vsAdd( int n,  float* vecIn,  float* vecIn2, float* vecOut);
   void vdMul( int n,  double* vecIn,  double* vecIn2, double* vecOut);
   void vsMul( int n,  float* vecIn,  float* vecIn2, float* vecOut);
   void vdAbs( int n,  double* vecIn, double* vecOut);
   void vsAbs( int n,  float* vecIn, float* vecOut);
}
#endif

// INTTerfaces to a few BLAS function, Level 1
/// INTTerface to cblas_*nrm2
template <typename T> T cblas_nrm2( INTT n,  T* X,  INTT incX);
/// INTTerface to cblas_*copy
template <typename T> void cblas_copy( INTT n,  T* X,  INTT incX,
      T* Y,  INTT incY);
/// INTTerface to cblas_*axpy
template <typename T> void cblas_axpy( INTT n,  T a,  T* X,
       INTT incX, T* Y,  INTT incY);
/// INTTerface to cblas_*scal
template <typename T> void cblas_scal( INTT n,  T a, T* X,
       INTT incX);
/// INTTerface to cblas_*asum
template <typename T> T cblas_asum( INTT n,  T* X,  INTT incX);
/// INTTerface to cblas_*adot
template <typename T> T cblas_dot( INTT n,  T* X,  INTT incX,
       T* Y, INTT incY);
/// interface to cblas_i*amin
template <typename T> int cblas_iamin( INTT n,  T* X,  INTT incX);
/// interface to cblas_i*amax
template <typename T> int cblas_iamax( INTT n,  T* X,  INTT incX);

// INTTerfaces to a few BLAS function, Level 2

/// INTTerface to cblas_*gemv
template <typename T> void cblas_gemv( CBLAS_ORDER order,
       CBLAS_TRANSPOSE TransA,  INTT M,
       INTT N,  T alpha,  T *A,  INTT lda,  T *X,
       INTT incX,  T beta,T *Y,   INTT incY);
/// INTTerface to cblas_*trmv
template <typename T> void inline cblas_trmv( CBLAS_ORDER order,  CBLAS_UPLO Uplo,
       CBLAS_TRANSPOSE TransA,  CBLAS_DIAG Diag,  INTT N,
       T *A,  INTT lda, T *X,  INTT incX);
/// INTTerface to cblas_*syr
template <typename T> void inline cblas_syr( CBLAS_ORDER order,
        CBLAS_UPLO Uplo,  INTT N,  T alpha,
       T *X,  INTT incX, T *A,  INTT lda);

/// INTTerface to cblas_*symv
template <typename T> inline void cblas_symv( CBLAS_ORDER order,
       CBLAS_UPLO Uplo,  INTT N,
       T alpha,  T *A,  INTT lda,  T *X,
       INTT incX,  T beta,T *Y,   INTT incY);

// INTTerfaces to a few BLAS function, Level 3
/// INTTerface to cblas_*gemm
template <typename T> void cblas_gemm( CBLAS_ORDER order,
       CBLAS_TRANSPOSE TransA,  CBLAS_TRANSPOSE TransB,
       INTT M,  INTT N,  INTT K,  T alpha,
       T *A,  INTT lda,  T *B,  INTT ldb,
       T beta, T *C,  INTT ldc);
/// INTTerface to cblas_*syrk
template <typename T> void cblas_syrk( CBLAS_ORDER order,
       CBLAS_UPLO Uplo,  CBLAS_TRANSPOSE Trans,  INTT N,  INTT K,
       T alpha,  T *A,  INTT lda,
       T beta, T*C,  INTT ldc);
/// INTTerface to cblas_*ger
template <typename T> void cblas_ger( CBLAS_ORDER order,
       INTT M,  INTT N,  T alpha,  T *X,  INTT incX,
       T* Y,  INTT incY, T*A,  INTT lda);
/// INTTerface to cblas_*trmm
template <typename T> void cblas_trmm( CBLAS_ORDER order,
       CBLAS_SIDE Side,  CBLAS_UPLO Uplo,
       CBLAS_TRANSPOSE TransA,  CBLAS_DIAG Diag,
       INTT M,  INTT N,  T alpha,
       T*A,  INTT lda,T *B,  INTT ldb);

// interfaces to a few functions from the intel Vector Mathematical Library
/// interface to v*Sqr
template <typename T> void vSqrt( int n,  T* vecIn, T* vecOut);
/// interface to v*Sqr
template <typename T> void vInvSqrt( int n,  T* vecIn, T* vecOut);
/// interface to v*Sqr
template <typename T> void vSqr( int n,  T* vecIn, T* vecOut);
/// interface to v*Sub
template <typename T> void vSub( int n,  T* vecIn,  T* vecIn2, T* vecOut);
/// interface to v*Div
template <typename T> void vDiv( int n,  T* vecIn,  T* vecIn2, T* vecOut);
/// interface to v*Exp
template <typename T> void vExp( int n,  T* vecIn, T* vecOut);
/// interface to v*Inv
template <typename T> void vInv( int n,  T* vecIn, T* vecOut);
/// interface to v*Add
template <typename T> void vAdd( int n,  T* vecIn,  T* vecIn2, T* vecOut);
/// interface to v*Mul
template <typename T> void vMul( int n,  T* vecIn,  T* vecIn2, T* vecOut);
/// interface to v*Abs
template <typename T> void vAbs( int n,  T* vecIn, T* vecOut);

// interfaces to a few LAPACK functions
/// interface to *trtri
template <typename T> void trtri(char& uplo, char& diag,
      INTT n, T * a, INTT lda);
/// interface to *sytri  // call sytrf
template <typename T> void sytri(char& uplo, INTT n, T* a, INTT lda);
//, INTT* ipiv,
//      T* work);
/// interaface to *lasrt
template <typename T> void lasrt(char& id,  INTT n, T *d);
//template <typename T> void lasrt2(char& id,  INTT& n, T *d, int* key);
template <typename T> void gesvd( char& jobu, char& jobvt, INTT m,
      INTT n, T* a, INTT lda, T* s,
      T* u, INTT ldu, T* vt, INTT ldvt);
template <typename T> void syev( char& jobz, char& uplo, INTT n,
         T* a, INTT lda, T* w);

/* ******************
 * Implementations
 * *****************/

extern "C" {
   double dnrm2_(INTT *n,double *x,INTT *incX);
   float snrm2_(INTT *n,float *x,INTT *incX);
   void dcopy_(INTT *n,double *x,INTT *incX, double *y,INTT *incY);
   void scopy_(INTT *n,float *x,INTT *incX, float *y,INTT *incY);
   void daxpy_(INTT *n,double* a, double *x,INTT *incX, double *y,INTT *incY);
   void saxpy_(INTT *n,float* a, float *x,INTT *incX, float *y,INTT *incY);
   void dscal_(INTT *n,double* a, double *x,INTT *incX);
   void sscal_(INTT *n,float* a, float *x,INTT *incX);
   double dasum_(INTT *n,double *x,INTT *incX);
   float sasum_(INTT *n,float *x,INTT *incX);
   double ddot_(INTT *n,double *x,INTT *incX, double *y,INTT *incY);
   float sdot_(INTT *n,float *x,INTT *incX, float *y,INTT *incY);
   void dgemv_(char   *trans, INTT *m, INTT *n, double *alpha, double *a,
         INTT *lda, double *x, INTT *incx, double *beta, double *y,INTT *incy);
   void sgemv_(char   *trans, INTT *m, INTT *n, float *alpha, float *a,
         INTT *lda, float *x, INTT *incx, float *beta, float *y,INTT *incy);
   void dger_(INTT *m, INTT *n, double *alpha, double *x, INTT *incx,
         double *y, INTT *incy, double *a, INTT *lda);
   void sger_(INTT *m, INTT *n, float *alpha, float *x, INTT *incx,
         float *y, INTT *incy, float *a, INTT *lda);
   void dtrmv_(char *uplo, char *trans, char   *diag, INTT *n, double *a,
         INTT *lda, double *x, INTT *incx);
   void strmv_(char *uplo, char *trans, char   *diag, INTT *n, float *a,
         INTT *lda, float *x, INTT *incx);
   void dsyr_(char *uplo, INTT *n, double *alpha, double *x, INTT *incx,
         double *a, INTT *lda);
   void ssyr_(char *uplo, INTT *n, float *alpha, float *x, INTT *incx,
         float *a, INTT *lda);
   void dsymv_(char   *uplo, INTT *n, double *alpha, double *a, INTT *lda,
         double *x, INTT *incx, double *beta, double *y, INTT *incy);
   void ssymv_(char   *uplo, INTT *n, float *alpha, float *a, INTT *lda,
         float *x, INTT *incx, float *beta, float *y, INTT *incy);
   void dgemm_(char *transa, char *transb, INTT *m, INTT *n, INTT *k,
         double *alpha, double *a, INTT *lda, double *b, INTT *ldb, double *beta,
         double *c, INTT *ldc);
   void sgemm_(char *transa, char *transb, INTT *m, INTT *n, INTT *k,
         float *alpha, float *a, INTT *lda, float *b, INTT *ldb, float *beta,
         float *c, INTT *ldc);
   void dsyrk_(char *uplo, char *trans, INTT *n, INTT *k, double *alpha,
         double *a, INTT *lda, double *beta, double *c, INTT *ldc);
   void ssyrk_(char *uplo, char *trans, INTT *n, INTT *k, float *alpha,
         float *a, INTT *lda, float *beta, float *c, INTT *ldc);
   void dtrmm_(char *side,char *uplo,char *transa, char *diag, INTT *m,
         INTT *n, double *alpha, double *a, INTT *lda, double *b,
         INTT *ldb);
   void strmm_(char *side,char *uplo,char *transa, char *diag, INTT *m,
         INTT *n, float *alpha, float *a, INTT *lda, float *b,
         INTT *ldb);
   INTT idamax_(INTT *n, double *dx, INTT *incx);
   INTT isamax_(INTT *n, float *dx, INTT *incx);
   void dtrtri_(char* uplo, char* diag, INTT* n, double * a, INTT* lda,
         INTT* info);
   void strtri_(char* uplo, char* diag, INTT* n, float * a, INTT* lda,
         INTT* info);
   void dsytrf_(char* uplo, INTT* n, double* a, INTT* lda, INTT* ipiv,
         double* work, INTT* lwork, INTT* info);
   void ssytrf_(char* uplo, INTT* n, float* a, INTT* lda, INTT* ipiv,
         float* work, INTT* lwork, INTT* info);
   void dsytri_(char* uplo, INTT* n, double* a, INTT* lda, INTT* ipiv,
         double* work, INTT* info);
   void ssytri_(char* uplo, INTT* n, float* a, INTT* lda, INTT* ipiv,
         float* work, INTT* info);
   void dlasrt_(char* id,  INTT* n, double *d, INTT* info);
   void slasrt_(char* id,  INTT* n, float*d, INTT* info);
   void dgesvd_(char*jobu, char *jobvt, INTT *m, INTT *n, double *a,
         INTT *lda, double *s, double *u, INTT *ldu, double *vt,
         INTT *ldvt, double *work, INTT *lwork, INTT *info);
   void sgesvd_(char*jobu, char *jobvt, INTT *m, INTT *n, float *a,
         INTT *lda, float *s, float *u, INTT *ldu, float *vt,
         INTT *ldvt, float *work, INTT *lwork, INTT *info);
   void dsyev_(char *jobz, char *uplo, INTT *n, double *a, INTT *lda,
         double *w, double *work, INTT *lwork, INTT *info);
   void ssyev_(char *jobz, char *uplo, INTT *n, float *a, INTT *lda,
         float *w, float *work, INTT *lwork, INTT *info);
}

// Implementations of the INTTerfaces, BLAS Level 1
/// Implementation of the INTTerface for cblas_dnrm2
template <> inline double cblas_nrm2<double>( INTT n,  double* X,
       INTT incX) {
   //return cblas_dnrm2(n,X,incX);
   return dnrm2_(&n,X,&incX);
};
/// Implementation of the INTTerface for cblas_snrm2
template <> inline float cblas_nrm2<float>( INTT n,  float* X,
       INTT incX) {
   //return cblas_snrm2(n,X,incX);
   return snrm2_(&n,X,&incX);
};
/// Implementation of the INTTerface for cblas_dcopy
template <> inline void cblas_copy<double>( INTT n,  double* X,
       INTT incX, double* Y,  INTT incY) {
   //cblas_dcopy(n,X,incX,Y,incY);
   dcopy_(&n,X,&incX,Y,&incY);
};
/// Implementation of the INTTerface for cblas_scopy
template <> inline void cblas_copy<float>( INTT n,  float* X,  INTT incX,
      float* Y,  INTT incY) {
   //cblas_scopy(n,X,incX,Y,incY);
   scopy_(&n,X,&incX,Y,&incY);
};
/// Implementation of the INTTerface for cblas_scopy
template <> inline void cblas_copy<int>( INTT n,  int* X,  INTT incX,
      int* Y,  INTT incY) {
   for (int i = 0; i<n; ++i)
      Y[incY*i]=X[incX*i];
};
/// Implementation of the INTTerface for cblas_scopy
template <> inline void cblas_copy<bool>( INTT n,  bool* X,  INTT incX,
      bool* Y,  INTT incY) {
   for (int i = 0; i<n; ++i)
      Y[incY*i]=X[incX*i];
};

/// Implementation of the INTTerface for cblas_daxpy
template <> inline void cblas_axpy<double>( INTT n,  double a,  double* X,
       INTT incX, double* Y,  INTT incY) {
   //cblas_daxpy(n,a,X,incX,Y,incY);
   daxpy_(&n,&a,X,&incX,Y,&incY);
};
/// Implementation of the INTTerface for cblas_saxpy
template <> inline void cblas_axpy<float>( INTT n,  float a,  float* X,
       INTT incX, float* Y,  INTT incY) {
   //cblas_saxpy(n,a,X,incX,Y,incY);
   saxpy_(&n,&a,X,&incX,Y,&incY);
};

/// Implementation of the INTTerface for cblas_saxpy
template <> inline void cblas_axpy<int>( INTT n,  int a,  int* X,
       INTT incX, int* Y,  INTT incY) {
   for (int i = 0; i<n; ++i)
      Y[i] += a*X[i];
};

/// Implementation of the INTTerface for cblas_saxpy
template <> inline void cblas_axpy<bool>( INTT n,  bool a,  bool* X,
       INTT incX, bool* Y,  INTT incY) {
   for (int i = 0; i<n; ++i)
      Y[i] = a*X[i];
};

/// Implementation of the INTTerface for cblas_dscal
template <> inline void cblas_scal<double>( INTT n,  double a, double* X,
       INTT incX) {
   //cblas_dscal(n,a,X,incX);
   dscal_(&n,&a,X,&incX);
};
/// Implementation of the INTTerface for cblas_sscal
template <> inline void cblas_scal<float>( INTT n,  float a, float* X,
       INTT incX) {
   //cblas_sscal(n,a,X,incX);
   sscal_(&n,&a,X,&incX);
};
/// Implementation of the INTTerface for cblas_sscal
template <> inline void cblas_scal<int>( INTT n,  int a, int* X,
       INTT incX) {
   for (int i = 0; i<n; ++i) X[i*incX]*=a;
};
/// Implementation of the INTTerface for cblas_sscal
template <> inline void cblas_scal<bool>( INTT n,  bool a, bool* X,
       INTT incX) {
   /// not implemented
};

/// Implementation of the INTTerface for cblas_dasum
template <> inline double cblas_asum<double>( INTT n,  double* X,  INTT incX) {
   //return cblas_dasum(n,X,1);
   return dasum_(&n,X,&incX);
};
/// Implementation of the INTTerface for cblas_sasum
template <> inline float cblas_asum<float>( INTT n,  float* X,  INTT incX) {
   //return cblas_sasum(n,X,1);
   return sasum_(&n,X,&incX);
};
/// Implementation of the INTTerface for cblas_ddot
template <> inline double cblas_dot<double>( INTT n,  double* X,
       INTT incX,  double* Y, INTT incY) {
   //return cblas_ddot(n,X,incX,Y,incY);
   return ddot_(&n,X,&incX,Y,&incY);
};
/// Implementation of the INTTerface for cblas_sdot
template <> inline float cblas_dot<float>( INTT n,  float* X,
       INTT incX,  float* Y, INTT incY) {
   //return cblas_sdot(n,X,incX,Y,incY);
   return sdot_(&n,X,&incX,Y,&incY);
};
template <> inline int cblas_dot<int>( INTT n,  int* X,
       INTT incX,  int* Y, INTT incY) {
   int total=0;
   int i,j;
   j=0;
   for (i = 0; i<n; ++i) {
      total+=X[i*incX]*Y[j];
      //j+=incY;
      j+=(int)incY;
   }
   return total;
};
/// Implementation of the INTTerface for cblas_sdot
template <> inline bool cblas_dot<bool>( INTT n,  bool* X,
       INTT incX,  bool* Y, INTT incY) {
   /// not implemented
   return true;
};

// Implementations of the INTTerfaces, BLAS Level 2
///  Implementation of the INTTerface for cblas_dgemv
template <> inline void cblas_gemv<double>( CBLAS_ORDER order,
       CBLAS_TRANSPOSE TransA,  INTT M,  INTT N,
       double alpha,  double *A,  INTT lda,
       double *X,  INTT incX,  double beta,
      double *Y,  INTT incY) {
   //cblas_dgemv(order,TransA,M,N,alpha,A,lda,X,incX,beta,Y,incY);
   dgemv_(cblas_transpose(TransA),&M,&N,&alpha,A,&lda,X,&incX,&beta,Y,&incY);
};
///  Implementation of the INTTerface for cblas_sgemv
template <> inline void cblas_gemv<float>( CBLAS_ORDER order,
       CBLAS_TRANSPOSE TransA,  INTT M,  INTT N,
       float alpha,  float *A,  INTT lda,
       float *X,  INTT incX,  float beta,
      float *Y,  INTT incY) {
   //cblas_sgemv(order,TransA,M,N,alpha,A,lda,X,incX,beta,Y,incY);
   sgemv_(cblas_transpose(TransA),&M,&N,&alpha,A,&lda,X,&incX,&beta,Y,&incY);
};
///  Implementation of the INTTerface for cblas_sgemv
template <> inline void cblas_gemv<int>( CBLAS_ORDER order,
       CBLAS_TRANSPOSE TransA,  INTT M,  INTT N,
       int alpha,  int *A,  INTT lda,
       int *X,  INTT incX,  int beta,
      int *Y,  INTT incY) {
   ///  not implemented
};
///  Implementation of the INTTerface for cblas_sgemv
template <> inline void cblas_gemv<bool>( CBLAS_ORDER order,
       CBLAS_TRANSPOSE TransA,  INTT M,  INTT N,
       bool alpha,  bool *A,  INTT lda,
       bool *X,  INTT incX,  bool beta,
      bool *Y,  INTT incY) {
   /// not implemented
};

///  Implementation of the INTTerface for cblas_dger
template <> inline void cblas_ger<double>( CBLAS_ORDER order,
       INTT M,  INTT N,  double alpha,  double *X,  INTT incX,
       double* Y,  INTT incY, double *A,  INTT lda) {
   //cblas_dger(order,M,N,alpha,X,incX,Y,incY,A,lda);
   dger_(&M,&N,&alpha,X,&incX,Y,&incY,A,&lda);
};
///  Implementation of the INTTerface for cblas_sger
template <> inline void cblas_ger<float>( CBLAS_ORDER order,
       INTT M,  INTT N,  float alpha,  float *X,  INTT incX,
       float* Y,  INTT incY, float *A,  INTT lda) {
   //cblas_sger(order,M,N,alpha,X,incX,Y,incY,A,lda);
   sger_(&M,&N,&alpha,X,&incX,Y,&incY,A,&lda);
};
///  Implementation of the INTTerface for cblas_dtrmv
template <> inline void cblas_trmv<double>( CBLAS_ORDER order,  CBLAS_UPLO Uplo,
       CBLAS_TRANSPOSE TransA,  CBLAS_DIAG Diag,  INTT N,
       double *A,  INTT lda, double *X,  INTT incX) {
   //cblas_dtrmv(order,Uplo,TransA,Diag,N,A,lda,X,incX);
   dtrmv_(cblas_uplo(Uplo),cblas_transpose(TransA),cblas_diag(Diag),&N,A,&lda,X,&incX);
};
///  Implementation of the INTTerface for cblas_strmv
template <> inline void cblas_trmv<float>( CBLAS_ORDER order,  CBLAS_UPLO Uplo,
       CBLAS_TRANSPOSE TransA,  CBLAS_DIAG Diag,  INTT N,
       float *A,  INTT lda, float *X,  INTT incX) {
   //cblas_strmv(order,Uplo,TransA,Diag,N,A,lda,X,incX);
   strmv_(cblas_uplo(Uplo),cblas_transpose(TransA),cblas_diag(Diag),&N,A,&lda,X,&incX);
};
/// Implementation of cblas_dsyr
template <> inline void cblas_syr( CBLAS_ORDER order,
        CBLAS_UPLO Uplo,
       INTT N,  double alpha,  double*X,
       INTT incX, double *A,  INTT lda) {
   //cblas_dsyr(order,Uplo,N,alpha,X,incX,A,lda);
   dsyr_(cblas_uplo(Uplo),&N,&alpha,X,&incX,A,&lda);
};
/// Implementation of cblas_ssyr
template <> inline void cblas_syr( CBLAS_ORDER order,
        CBLAS_UPLO Uplo,
       INTT N,  float alpha,  float*X,
       INTT incX, float *A,  INTT lda) {
   //cblas_ssyr(order,Uplo,N,alpha,X,incX,A,lda);
   ssyr_(cblas_uplo(Uplo),&N,&alpha,X,&incX,A,&lda);
};
/// Implementation of cblas_ssymv
template <> inline void cblas_symv( CBLAS_ORDER order,
       CBLAS_UPLO Uplo,  INTT N,
       float alpha,  float *A,  INTT lda,  float *X,
       INTT incX,  float beta,float *Y,   INTT incY) {
   //cblas_ssymv(order,Uplo,N,alpha,A,lda,X,incX,beta,Y,incY);
   ssymv_(cblas_uplo(Uplo),&N,&alpha,A,&lda,X,&incX,&beta,Y,&incY);
}
/// Implementation of cblas_dsymv
template <> inline void cblas_symv( CBLAS_ORDER order,
       CBLAS_UPLO Uplo,  INTT N,
       double alpha,  double *A,  INTT lda,  double *X,
       INTT incX,  double beta,double *Y,   INTT incY) {
   //cblas_dsymv(order,Uplo,N,alpha,A,lda,X,incX,beta,Y,incY);
   dsymv_(cblas_uplo(Uplo),&N,&alpha,A,&lda,X,&incX,&beta,Y,&incY);
}

// Implementations of the INTTerfaces, BLAS Level 3
///  Implementation of the INTTerface for cblas_dgemm
template <> inline void cblas_gemm<double>( CBLAS_ORDER order,
       CBLAS_TRANSPOSE TransA,  CBLAS_TRANSPOSE TransB,
       INTT M,  INTT N,  INTT K,  double alpha,
       double *A,  INTT lda,  double *B,  INTT ldb,
       double beta, double *C,  INTT ldc) {
   //cblas_dgemm(Order,TransA,TransB,M,N,K,alpha,A,lda,B,ldb,beta,C,ldc);
   dgemm_(cblas_transpose(TransA),cblas_transpose(TransB),&M,&N,&K,&alpha,A,&lda,B,&ldb,&beta,C,&ldc);
};
///  Implementation of the INTTerface for cblas_sgemm
template <> inline void cblas_gemm<float>( CBLAS_ORDER order,
       CBLAS_TRANSPOSE TransA,  CBLAS_TRANSPOSE TransB,
       INTT M,  INTT N,  INTT K,  float alpha,
       float *A,  INTT lda,  float *B,  INTT ldb,
       float beta, float *C,  INTT ldc) {
   //cblas_sgemm(Order,TransA,TransB,M,N,K,alpha,A,lda,B,ldb,beta,C,ldc);
   sgemm_(cblas_transpose(TransA),cblas_transpose(TransB),&M,&N,&K,&alpha,A,&lda,B,&ldb,&beta,C,&ldc);
};
template <> inline void cblas_gemm<int>( CBLAS_ORDER order,
       CBLAS_TRANSPOSE TransA,  CBLAS_TRANSPOSE TransB,
       INTT M,  INTT N,  INTT K,  int alpha,
       int *A,  INTT lda,  int *B,  INTT ldb,
       int beta, int *C,  INTT ldc) {
   /// not implemented
};
///  Implementation of the INTTerface for cblas_sgemm
template <> inline void cblas_gemm<bool>( CBLAS_ORDER order,
       CBLAS_TRANSPOSE TransA,  CBLAS_TRANSPOSE TransB,
       INTT M,  INTT N,  INTT K,  bool alpha,
       bool *A,  INTT lda,  bool *B,  INTT ldb,
       bool beta, bool *C,  INTT ldc) {
   /// not implemented
};

///  Implementation of the INTTerface for cblas_dsyrk
template <> inline void cblas_syrk<double>( CBLAS_ORDER order,
       CBLAS_UPLO Uplo,  CBLAS_TRANSPOSE Trans,  INTT N,  INTT K,
       double alpha,  double *A,  INTT lda,
       double beta, double *C,  INTT ldc) {
   //cblas_dsyrk(Order,Uplo,Trans,N,K,alpha,A,lda,beta,C,ldc);
   dsyrk_(cblas_uplo(Uplo),cblas_transpose(Trans),&N,&K,&alpha,A,&lda,&beta,C,&ldc);
};
///  Implementation of the INTTerface for cblas_ssyrk
template <> inline void cblas_syrk<float>( CBLAS_ORDER order,
       CBLAS_UPLO Uplo,  CBLAS_TRANSPOSE Trans,  INTT N,  INTT K,
       float alpha,  float *A,  INTT lda,
       float beta, float *C,  INTT ldc) {
   //cblas_ssyrk(Order,Uplo,Trans,N,K,alpha,A,lda,beta,C,ldc);
   ssyrk_(cblas_uplo(Uplo),cblas_transpose(Trans),&N,&K,&alpha,A,&lda,&beta,C,&ldc);
};
///  Implementation of the INTTerface for cblas_ssyrk
template <> inline void cblas_syrk<int>( CBLAS_ORDER order,
       CBLAS_UPLO Uplo,  CBLAS_TRANSPOSE Trans,  INTT N,  INTT K,
       int alpha,  int *A,  INTT lda,
       int beta, int *C,  INTT ldc) {
   /// not implemented
};
///  Implementation of the INTTerface for cblas_ssyrk
template <> inline void cblas_syrk<bool>( CBLAS_ORDER order,
       CBLAS_UPLO Uplo,  CBLAS_TRANSPOSE Trans,  INTT N,  INTT K,
       bool alpha,  bool *A,  INTT lda,
       bool beta, bool *C,  INTT ldc) {
   /// not implemented
};

///  Implementation of the INTTerface for cblas_dtrmm
template <> inline void cblas_trmm<double>( CBLAS_ORDER order,
       CBLAS_SIDE Side,  CBLAS_UPLO Uplo,
       CBLAS_TRANSPOSE TransA,  CBLAS_DIAG Diag,
       INTT M,  INTT N,  double alpha,
       double *A,  INTT lda,double *B,  INTT ldb) {
   //cblas_dtrmm(Order,Side,Uplo,TransA,Diag,M,N,alpha,A,lda,B,ldb);
   dtrmm_(cblas_side(Side),cblas_uplo(Uplo),cblas_transpose(TransA),cblas_diag(Diag),&M,&N,&alpha,A,&lda,B,&ldb);
};
///  Implementation of the INTTerface for cblas_strmm
template <> inline void cblas_trmm<float>( CBLAS_ORDER order,
       CBLAS_SIDE Side,  CBLAS_UPLO Uplo,
       CBLAS_TRANSPOSE TransA,  CBLAS_DIAG Diag,
       INTT M,  INTT N,  float alpha,
       float *A,  INTT lda,float *B,  INTT ldb) {
   //cblas_strmm(Order,Side,Uplo,TransA,Diag,M,N,alpha,A,lda,B,ldb);
   strmm_(cblas_side(Side),cblas_uplo(Uplo),cblas_transpose(TransA),cblas_diag(Diag),&M,&N,&alpha,A,&lda,B,&ldb);
};
///  Implementation of the interface for cblas_idamax
template <> inline int cblas_iamax<double>( INTT n,  double* X,
       INTT incX) {
   //return cblas_idamax(n,X,incX);
   return static_cast<int>(idamax_(&n,X,&incX)-1);
};
///  Implementation of the interface for cblas_isamax
template <> inline int cblas_iamax<float>( INTT n,  float* X,
       INTT incX) {
   //return cblas_isamax(n,X,incX);
   return static_cast<int>(isamax_(&n,X,&incX)-1);
};

// Implementations of the interfaces, LAPACK
/// Implemenation of the interface for dtrtri
template <> inline void trtri<double>(char& uplo, char& diag,
      INTT n, double * a, INTT lda) {
   //dtrtri_(&uplo,&diag,&n,a,&lda,&info);
   dtrtri_(&uplo,&diag,&n,a,&lda,&info);
};
/// Implemenation of the interface for strtri
template <> inline void trtri<float>(char& uplo, char& diag,
      INTT n, float* a, INTT lda) {
   //strtri_(&uplo,&diag,&n,a,&lda,&info);
   strtri_(&uplo,&diag,&n,a,&lda,&info);
};

/// Implemenation of the interface for dsytri
template <> inline void sytri<double>(char& uplo, INTT n, double* a, INTT lda) {
//, INTT* ipiv, double* work) {
   //dsytri_(&uplo,&n,a,&lda,ipiv,work,&info);
   INTT lwork=-1;
   INTT* ipiv= new INTT[n];
   double* query, *work;
   query = new double[1];
   dsytrf_(&uplo,&n,a,&lda,ipiv,query,&lwork,&info);
   lwork=static_cast<INTT>(*query);
   delete[](query);
   work = new double[static_cast<int>(lwork)];
   dsytrf_(&uplo,&n,a,&lda,ipiv,work,&lwork,&info);
   delete[](work);
   work = new double[static_cast<int>(2*n)];
   dsytri_(&uplo,&n,a,&lda,ipiv,work,&info);
   delete[](work);
   delete[](ipiv);
};
/// Implemenation of the interface for ssytri
template <> inline void sytri<float>(char& uplo, INTT n, float* a, INTT lda) {
   INTT lwork=-1;
   INTT* ipiv= new INTT[n];
   float* query, *work;
   query = new float[1];
   ssytrf_(&uplo,&n,a,&lda,ipiv,query,&lwork,&info);
   lwork=static_cast<INTT>(*query);
   delete[](query);
   work = new float[static_cast<int>(lwork)];
   ssytrf_(&uplo,&n,a,&lda,ipiv,work,&lwork,&info);
   delete[](work);
   work = new float[static_cast<int>(2*n)];
   ssytri_(&uplo,&n,a,&lda,ipiv,work,&info);
   delete[](work);
   delete[](ipiv);
};
/// interaface to *lasrt
template <> inline void lasrt(char& id, INTT n, double *d) {
   //dlasrt_(&id,const_cast<int*>(&n),d,&info);
   dlasrt_(&id,&n,d,&info);
};
/// interaface to *lasrt
template <> inline void lasrt(char& id, INTT n, float *d) {
   //slasrt_(&id,const_cast<int*>(&n),d,&info);
   slasrt_(&id,&n,d,&info);
};
//template <> inline void lasrt2(char& id, INTT& n, double *d,int* key) {
//   //dlasrt2_(&id,const_cast<int*>(&n),d,key,&info);
//   dlasrt2(&id,&n,d,key,&info);
//};
///// interaface to *lasrt
//template <> inline void lasrt2(char& id, INTT& n, float *d, int* key) {
//   //slasrt2_(&id,const_cast<int*>(&n),d,key,&info);
//   slasrt2(&id,&n,d,key,&info);
//};
template <> void inline gesvd( char& jobu, char& jobvt, INTT m,
      INTT n, double* a, INTT lda, double* s,
      double* u, INTT ldu, double* vt, INTT ldvt) {
   double* query = new double[1];
   INTT lwork=-1;
   dgesvd_(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt,
         query, &lwork, &info );
   lwork=static_cast<INTT>(*query);
   delete[](query);
   double* work = new double[static_cast<int>(lwork)];
   dgesvd_(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt,
         work, &lwork, &info );
   delete[](work);
}
template <> void inline gesvd( char& jobu, char& jobvt, INTT m,
      INTT n, float* a, INTT lda, float* s,
      float* u, INTT ldu, float* vt, INTT ldvt) {
   float* query = new float[1];
   INTT lwork=-1;
   sgesvd_(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt,
         query, &lwork, &info );
   lwork=static_cast<INTT>(*query);
   delete[](query);
   float* work = new float[static_cast<int>(lwork)];
   sgesvd_(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt,
         work, &lwork, &info );
   delete[](work);
}

template <> void inline syev( char& jobz, char& uplo, INTT n,
         float* a, INTT lda, float* w) {
   float* query = new float[1];
   INTT lwork=-1;
   ssyev_(&jobz,&uplo,&n,a,&lda,w,query,&lwork,&info);
   lwork=static_cast<INTT>(*query);
   delete[](query);
   float* work = new float[static_cast<int>(lwork)];
   ssyev_(&jobz,&uplo,&n,a,&lda,w,work,&lwork,&info);
   delete[](work);
};

template <> void inline syev( char& jobz, char& uplo, INTT n,
         double* a, INTT lda, double* w) {
   double* query = new double[1];
   INTT lwork=-1;
   dsyev_(&jobz,&uplo,&n,a,&lda,w,query,&lwork,&info);
   lwork=static_cast<INTT>(*query);
   delete[](query);
   double* work = new double[static_cast<int>(lwork)];
   dsyev_(&jobz,&uplo,&n,a,&lda,w,work,&lwork,&info);
   delete[](work);
};

/// If the MKL is not present, a slow implementation is used instead.
#ifdef HAVE_MKL
/// Implemenation of the interface for vdSqr
template <> inline void vSqr<double>( int n,  double* vecIn,
      double* vecOut) {
   vdSqr(n,vecIn,vecOut);
};
/// Implemenation of the interface for vsSqr
template <> inline void vSqr<float>( int n,  float* vecIn,
      float* vecOut) {
   vsSqr(n,vecIn,vecOut);
};
template <> inline void vSqrt<double>( int n,  double* vecIn,
      double* vecOut) {
   vdSqrt(n,vecIn,vecOut);
};
/// Implemenation of the interface for vsSqr
template <> inline void vSqrt<float>( int n,  float* vecIn,
      float* vecOut) {
   vsSqrt(n,vecIn,vecOut);
};
template <> inline void vInvSqrt<double>( int n,  double* vecIn,
      double* vecOut) {
   vdInvSqrt(n,vecIn,vecOut);
};
/// Implemenation of the interface for vsSqr
template <> inline void vInvSqrt<float>( int n,  float* vecIn,
      float* vecOut) {
   vsInvSqrt(n,vecIn,vecOut);
};

/// Implemenation of the interface for vdSub
template <> inline void vSub<double>( int n,  double* vecIn,
       double* vecIn2, double* vecOut) {
   vdSub(n,vecIn,vecIn2,vecOut);
};
/// Implemenation of the interface for vsSub
template <> inline void vSub<float>( int n,  float* vecIn,
       float* vecIn2, float* vecOut) {
   vsSub(n,vecIn,vecIn2,vecOut);
};
/// Implemenation of the interface for vdDiv
template <> inline void vDiv<double>( int n,  double* vecIn,
       double* vecIn2, double* vecOut) {
   vdDiv(n,vecIn,vecIn2,vecOut);
};
/// Implemenation of the interface for vsDiv
template <> inline void vDiv<float>( int n,  float* vecIn,
       float* vecIn2, float* vecOut) {
   vsDiv(n,vecIn,vecIn2,vecOut);
};
/// Implemenation of the interface for vdExp
template <> inline void vExp<double>( int n,  double* vecIn,
      double* vecOut) {
   vdExp(n,vecIn,vecOut);
};
/// Implemenation of the interface for vsExp
template <> inline void vExp<float>( int n,  float* vecIn,
      float* vecOut) {
   vsExp(n,vecIn,vecOut);
};
/// Implemenation of the interface for vdInv
template <> inline void vInv<double>( int n,  double* vecIn,
      double* vecOut) {
   vdInv(n,vecIn,vecOut);
};
/// Implemenation of the interface for vsInv
template <> inline void vInv<float>( int n,  float* vecIn,
      float* vecOut) {
   vsInv(n,vecIn,vecOut);
};
/// Implemenation of the interface for vdAdd
template <> inline void vAdd<double>( int n,  double* vecIn,
       double* vecIn2, double* vecOut) {
   vdAdd(n,vecIn,vecIn2,vecOut);
};
/// Implemenation of the interface for vsAdd
template <> inline void vAdd<float>( int n,  float* vecIn,
       float* vecIn2, float* vecOut) {
   vsAdd(n,vecIn,vecIn2,vecOut);
};
/// Implemenation of the interface for vdMul
template <> inline void vMul<double>( int n,  double* vecIn,
       double* vecIn2, double* vecOut) {
   vdMul(n,vecIn,vecIn2,vecOut);
};
/// Implemenation of the interface for vsMul
template <> inline void vMul<float>( int n,  float* vecIn,
       float* vecIn2, float* vecOut) {
   vsMul(n,vecIn,vecIn2,vecOut);
};

/// interface to vdAbs
template <> inline void vAbs( int n,  double* vecIn,
      double* vecOut) {
   vdAbs(n,vecIn,vecOut);
};
/// interface to vdAbs
template <> inline void vAbs( int n,  float* vecIn,
      float* vecOut) {
   vsAbs(n,vecIn,vecOut);
};

/// implemenation of the interface of the non-offical blas, level 1 function
/// cblas_idamin
template <> inline int cblas_iamin<double>( int n,  double* x,
       int incx) {
   return (int) cblas_idamin(n,x,incx);
};
/// implemenation of the interface of the non-offical blas, level 1 function
/// cblas_isamin
template <> inline int cblas_iamin<float>( int n,  float* x,
       int incx) {
   return (int) cblas_isamin(n,x,incx);
};
/// slow alternative implementation of some MKL function
#else
/// Slow implementation of vdSqr and vsSqr
template <typename T> inline void vSqr( int n,  T* vecIn, T* vecOut) {
   for (int i = 0; i<n; ++i) vecOut[i]=vecIn[i]*vecIn[i];
};
template <typename T> inline void vSqrt( int n,  T* vecIn, T* vecOut) {
   for (int i = 0; i<n; ++i) vecOut[i]=sqr<T>(vecIn[i]);
};
template <typename T> inline void vInvSqrt( int n,  T* vecIn, T* vecOut) {
   for (int i = 0; i<n; ++i) vecOut[i]=T(1.0)/sqr<T>(vecIn[i]);
};

/// Slow implementation of vdSub and vsSub
template <typename T> inline void vSub( int n,  T* vecIn1,
       T* vecIn2, T* vecOut) {
   for (int i = 0; i<n; ++i) vecOut[i]=vecIn1[i]-vecIn2[i];
};
/// Slow implementation of vdInv and vsInv
template <typename T> inline void vInv( int n,  T* vecIn, T* vecOut) {
   for (int i = 0; i<n; ++i) vecOut[i]=1.0/vecIn[i];
};
/// Slow implementation of vdExp and vsExp
template <typename T> inline void vExp( int n,  T* vecIn, T* vecOut) {
   for (int i = 0; i<n; ++i) vecOut[i]=exp(vecIn[i]);
};
/// Slow implementation of vdAdd and vsAdd
template <typename T> inline void vAdd( int n,  T* vecIn1,
       T* vecIn2, T* vecOut) {
   for (int i = 0; i<n; ++i) vecOut[i]=vecIn1[i]+vecIn2[i];
};
/// Slow implementation of vdMul and vsMul
template <typename T> inline void vMul( int n,  T* vecIn1,
       T* vecIn2, T* vecOut) {
   for (int i = 0; i<n; ++i) vecOut[i]=vecIn1[i]*vecIn2[i];
};
/// Slow implementation of vdDiv and vsDiv
template <typename T> inline void vDiv( int n,  T* vecIn1,
       T* vecIn2, T* vecOut) {
   for (int i = 0; i<n; ++i) vecOut[i]=vecIn1[i]/vecIn2[i];
};
/// Slow implementation of vAbs
template <typename T> inline void vAbs( int n,  T* vecIn,
      T* vecOut) {
   for (int i = 0; i<n; ++i) vecOut[i]=abs<T>(vecIn[i]);
};

/// Slow implementation of cblas_idamin and cblas_isamin
template <typename T> int inline cblas_iamin(INTT n, T* X, INTT incX) {
   int imin=0;
   double min=fabs(X[0]);
   for (int j = 1; j<n; j+=incX) {
      double cur = fabs(X[j]);
      if (cur < min) {
         imin=j;
         min = cur;
      }
   }
   return imin;
}
#endif

#endif

max_stars_repo_name: ChrystleMyrnaLobo/person-reid-benchmark


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

#pragma once
#include "Data\Pose.h"
#include "Data\Data.h"

#include <gsl\gsl>

namespace mage
{
    /**
     * Converts the contents of a CVMat into the mage::Matrix data container by copying the values
     */
    Matrix ToMageMat(const cv::Matx34f& cvMat);
    Matrix ToMageMat(const cv::Matx44f& cvMat);

    Position ToMagePos(const cv::Point3f& cvPt);
    Direction ToMageDir(const cv::Point3f& cvDir);

    cv::Vec3f FromMageDir(const Direction& mageDir);

    Matrix ToMageMat(const std::array<float, 4 * 4>& f);

    Direction ToMageDir(const std::array<float, 3>& dir);

    mage::Pose MageMatrixToMagePose(const mage::Matrix& viewMatrix);

    void ConvertBGRToNV12(const cv::Mat& srcImg, gsl::span<uint8_t> destBuffer);

    bool IsIdentity(const Matrix& m);

    inline mage::Matrix CreateIdentityMageMatrix()
    {
        return{ 1.0f, 0.0f, 0.0f, 0.0f,
            0.0f, 1.0f, 0.0f, 0.0f,
            0.0f, 0.0f, 1.0f, 0.0f,
            0.0f, 0.0f, 0.0f, 1.0f };
    };

    // Interpolates values between min and max to the range defined from a to b,
    // [a:b] can be either an increasing or decreasing range, min is always mapped to a, and max to b.
    unsigned char LinearInterpolationToChar(float value, float min, float max, unsigned char a, unsigned char b);

    // Construct external normal (z positive) from the X,Y components of the internal normal (z negative).
    // Synthetic normals are facing the camera (external),
    // generated normals are not (internal), this function is used to make the conversion,
    // by changing the sign of the coordinates.
    cv::Vec3f ConstructExternalNormalFromXY(float normalX, float normalY);
}

ext: h

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
//
//  particle.h
//  EpiGenMCMC
//
//  Created by Lucy Li on 13/04/2016.
//  Copyright (c) 2016 Lucy Li, Imperial College London. All rights reserved.
//

#ifndef __EpiGenMCMC__particle__
#define __EpiGenMCMC__particle__

#include <vector>
#include <gsl/gsl_randist.h>
#include "trajectory.h"

class Particle {
    int num_particles;
    int num_groups;
    int num_time_steps;
    std::vector <Trajectory*> trajectories;
    std::vector <double> weights;
    std::vector <double> normalised_weights;
    std::vector <unsigned int> next_particles;
    std::vector <int> parents;
    std::vector <double> overall_traj; // Incidence: P1[Group1[T1, T2...], Group2[T1, T2...]...], P2[] ...
    std::vector <double> overall_traj2; // Prevalence: P1[Group1[T1, T2...], Group2[T1, T2...]...], P2[] ...
    std::vector <int> particle_ancestry;
    std::vector <int> non_zero_particles;
    std::vector <double> cumulative_weights;
    std::vector <int> empty;
    void replace(int, int);
    void normalise();
public:
    Particle();
    Particle(int, Trajectory);
    void start_particle_tracing(int, int);
    void save_traj_to_matrix(int, int);
    void save_traj_to_matrix(int, int, int);
    void save_ancestry(int, int);
    void save_ancestry(int, int, int);
    void reset_parents();
    int get_traj_random(gsl_rng *);
    void retrace_traj(Trajectory&, gsl_rng *);
    Trajectory* get_traj(int) const;
    void set_weight(double, int, bool);
    double get_weight(int);
    double get_total_weight();
    double get_ESS();
    void reset_weights();
    void resample(gsl_rng *);
    void clear();
};

#endif /* defined(__EpiGenMCMC__particle__) */

max_stars_repo_name: lucymli/EpiGenMCMC

max_issues_repo_name: lucymli/EpiGenMCMC

max_forks_repo_name: lucymli/EpiGenMCMC


text:
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <mpfr.h>
#include <petsc.h>

#include "../ellipsoid/ellipsoid.h"
#include "../ellipsoid/ellSolv.h"
#include "../constants.h"
#include "testFunctions.h"
#include "../tanhsinh.h"

/*
EX3 - PLOTS CONVERGENCE OF A SINGLE NORMALIZATION CONSTANT
normConv.py PLOTS RESULT
*/

#undef __FUNCT__
#define __FUNCT__ "NormConstantIntFixed"
PetscErrorCode NormConstantIntFixed(EllipsoidalSystem *e, PetscInt n, PetscInt p, PetscInt prec, PetscInt nPts, PetscReal *intVals, PetscReal *normConst)
{
  PetscErrorCode ierr;
  PetscInt flopCount;
  PetscFunctionBegin;
  flopCount = 0;

  FuncInfo2 ctx1 = { .e = e, .n = n, .p = p, .numeratorType = 0, .denomSign = 1 };
  FuncInfo2 ctx2 = { .e = e, .n = n, .p = p, .numeratorType = 1, .denomSign = 1 };
  FuncInfo2 ctx3 = { .e = e, .n = n, .p = p, .numeratorType = 0, .denomSign = -1 };
  FuncInfo2 ctx4 = { .e = e, .n = n, .p = p, .numeratorType = 1, .denomSign = -1 };

  PetscReal integrals[4];

  mpfr_t mpfrzero;
  mpfr_t mpfrone;
  mpfr_inits(mpfrzero, mpfrone, NULL);
  mpfr_set_d(mpfrzero, 0.0, MPFR_RNDN);
  mpfr_set_d(mpfrone, 1.0, MPFR_RNDN);

  ierr = DEQuad((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, e->hp_h, e->hp_k, prec, nPts, integrals+0, &ctx1);
  intVals[0] = integrals[0];
  ierr = DEQuad((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, e->hp_h, e->hp_k, prec, nPts, integrals+1, &ctx2);
  intVals[1] = integrals[1];
  ierr = DEQuad((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, mpfrzero, e->hp_h, prec, nPts, integrals+2, &ctx3);
  intVals[2] = integrals[2];
  ierr = DEQuad((PetscErrorCode (*)(mpfr_t*, mpfr_t*, void*)) normFunction1, mpfrzero, e->hp_h, prec, nPts, integrals+3, &ctx4);
  intVals[3] = integrals[3];
  *normConst = 8.0*(integrals[2]*integrals[1] - integrals[0]*integrals[3]); flopCount += 4;

  mpfr_clears(mpfrzero, mpfrone, NULL);

  ierr = PetscLogFlops(flopCount);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT_ "normConvergence"
PetscErrorCode normConvergence(PetscInt n, PetscInt p)
{
  PetscErrorCode ierr;

  const PetscInt MPFR_PREC = 64;
  const PetscInt NUM_SOLUTIONS = 80;
  const PetscInt POINTS_MIN = 4;
  const PetscInt POINTS_STEP = 2;
  const PetscInt prec = 16;
  const PetscReal a = 3.0;
  const PetscReal b = 2.0;
  const PetscReal c = 1.0;
  mpfr_t solExact;
  mpfr_t intExact[4];
  PetscReal solutions[NUM_SOLUTIONS];
  mpfr_t errors   [NUM_SOLUTIONS];
  PetscReal intSols[4*NUM_SOLUTIONS];
  mpfr_t intErrors[4*NUM_SOLUTIONS];
  PetscReal sol2, sol3, sol4;
  PetscLogEvent flopCounts[NUM_SOLUTIONS];
  EllipsoidalSystem e;
  PetscInt i;
  PetscEventPerfInfo info;
  //PetscInt n = 7;
  //PeatscInt p = 4;

  PetscFunctionBegin;

  mpfr_set_default_prec(4*MPFR_PREC);
  initEllipsoidalSystem(&e, a, b, c, MPFR_PREC);

  mpfr_init(solExact);
  for(PetscInt k=0; k < 4; ++k)
    mpfr_init(intExact[k]);
  for(PetscInt k=0; k < NUM_SOLUTIONS; ++k )
    mpfr_init(errors[k]);
  for(PetscInt k=0; k < 4*NUM_SOLUTIONS; ++k)
    mpfr_init(intErrors[k]);

  /* calculate "exact" solution */
  ierr = calcNormalization2MPFR(&e, n, p, intExact, &solExact);CHKERRQ(ierr);
  printf("old norm constant: %15.15f\n", solExact);

  /* calculate approximate solutions and record flops */
  char text[40] = "%d points";
  char sText[40];
  for(i=0; i<NUM_SOLUTIONS; ++i) {
    sprintf(sText, text, POINTS_MIN + POINTS_STEP*i);
    ierr = PetscLogEventRegister(sText, 0, flopCounts+i);CHKERRQ(ierr);
    ierr = PetscLogEventBegin(flopCounts[i], 0, 0, 0, 0);CHKERRQ(ierr);
    ierr = NormConstantIntFixed(&e, n, p, prec, POINTS_MIN + POINTS_STEP*i, intSols+4*i, solutions+i);CHKERRQ(ierr);
    ierr = PetscLogEventEnd(flopCounts[i], 0, 0, 0, 0);CHKERRQ(ierr);
    printf("new norm constant: %15.15f\n", solutions[i]);
  }
  /* calculate errors */
  for(i=0; i<NUM_SOLUTIONS; ++i) {
    //errors[i] = PetscAbsReal((solExact - solutions[i])/solExact);
    mpfr_sub_d(errors[i], solExact, solutions[i], MPFR_RNDN);
    mpfr_div(errors[i], errors[i], solExact, MPFR_RNDN);
    mpfr_abs(errors[i], errors[i], MPFR_RNDN);
    //intErrors[4*i+0] = PetscAbsReal((intExact[0] - intSols[4*i+0])/intExact[0]);
    mpfr_sub_d(intErrors[4*i+0], intExact[0], intSols[4*i+0], MPFR_RNDN);
    mpfr_div(intErrors[4*i+0], intErrors[4*i+0], intExact[0], MPFR_RNDN);
    mpfr_abs(intErrors[4*i+0], intErrors[4*i+0], MPFR_RNDN);
    //intErrors[4*i+1] = PetscAbsReal((intExact[1] - intSols[4*i+1])/intExact[1]);
    mpfr_sub_d(intErrors[4*i+1], intExact[1], intSols[4*i+1], MPFR_RNDN);
    mpfr_div(intErrors[4*i+1], intErrors[4*i+1], intExact[1], MPFR_RNDN);
    mpfr_abs(intErrors[4*i+1], intErrors[4*i+1], MPFR_RNDN);
    //intErrors[4*i+2] = PetscAbsReal((intExact[2] - intSols[4*i+2])/intExact[2]);
    mpfr_sub_d(intErrors[4*i+2], intExact[2], intSols[4*i+2], MPFR_RNDN);
    mpfr_div(intErrors[4*i+2], intErrors[4*i+2], intExact[2], MPFR_RNDN);
    mpfr_abs(intErrors[4*i+2], intErrors[4*i+2], MPFR_RNDN);
    //intErrors[4*i+3] = PetscAbsReal((intExact[3] - intSols[4*i+3])/intExact[3]);
    mpfr_sub_d(intErrors[4*i+3], intExact[3], intSols[4*i+3], MPFR_RNDN);
    mpfr_div(intErrors[4*i+3], intErrors[4*i+3], intExact[3], MPFR_RNDN);
    mpfr_abs(intErrors[4*i+3], intErrors[4*i+3], MPFR_RNDN);
    printf("errors[%d] = %15.15f\n", i, errors[i]);
  }

  FILE *fp1 = fopen("out/normInt1Prec.txt", "w");
  FILE *fp2 = fopen("out/normInt2Prec.txt", "w");
  FILE *fp3 = fopen("out/normInt3Prec.txt", "w");
  FILE *fp4 = fopen("out/normInt4Prec.txt", "w");
  FILE *fp = fopen("out/normWorkPrec.txt", "w");
  fprintf(fp, "points flops error\n");
  fprintf(fp1, "points error\n");
  fprintf(fp2, "points error\n");
  fprintf(fp3, "points error\n");
  fprintf(fp4, "points error\n");
  for(i=0; i<NUM_SOLUTIONS; ++i) {
    ierr = PetscLogEventGetPerfInfo(PETSC_DETERMINE, flopCounts[i], &info);CHKERRQ(ierr);
    fprintf(fp, "%d %4.4e %4.4e\n", POINTS_MIN + POINTS_STEP*i, info.flops, mpfr_get_d(errors[i], MPFR_RNDN));
    fprintf(fp1, "%d %4.4e\n", POINTS_MIN + POINTS_STEP*i, mpfr_get_d(intErrors[4*i+0], MPFR_RNDN));
    fprintf(fp2, "%d %4.4e\n", POINTS_MIN + POINTS_STEP*i, mpfr_get_d(intErrors[4*i+1], MPFR_RNDN));
    fprintf(fp3, "%d %4.4e\n", POINTS_MIN + POINTS_STEP*i, mpfr_get_d(intErrors[4*i+2], MPFR_RNDN));
    fprintf(fp4, "%d %4.4e\n", POINTS_MIN + POINTS_STEP*i, mpfr_get_d(intErrors[4*i+3], MPFR_RNDN));
  }
  fclose(fp);
  fclose(fp1);
  fclose(fp2);
  fclose(fp3);
  mpfr_clear(solExact);

  PetscFunctionReturn(0);

}

#undef __FUNCT__
#define __FUNCT__ "main"
PetscErrorCode main(int argc, char **argv)
{
  PetscErrorCode ierr;
  PetscFunctionBeginUser;

  ierr = PetscInitialize(&argc, &argv, NULL, NULL);CHKERRQ(ierr);
  ierr = PetscLogDefaultBegin();CHKERRQ(ierr);

  ierr = normConvergence(21,22);CHKERRQ(ierr);

  ierr = PetscFinalize();

  PetscFunctionReturn(0);
}

ext: c

max_stars_repo_name: tom-klotz/ellipsoid-solvation

max_issues_repo_name: tom-klotz/ellipsoid-solvation

max_forks_repo_name: tom-klotz/ellipsoid-solvation


text:
/*
This files contains some generic utilities
 */

/* ----- System includes -----*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/* ----- IFU and GSL includes -----*/
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_sort.h>
#include "IFU_io.h"
// #include "IFU_math.h"
/* ----- Specific IFU includes ----- */
#include "data_io.h"
#ifdef FITS
#include <fitsio.h>
#endif

/* Removing MIDAS by hand - not really clean, but should allow to compile with an IFU with --disable-midas enabled*/
#undef MIDAS

#ifdef MIDAS
#include "midas_defs.h"
#include "tbldef.h"
#include "fctext.h"
#include "ldbext.h"
#include "proto_st.h"

#endif
extern IO_Format InputIO, OutputIO;

/* ----- local includes ----- */
#include "utils.h"

/* ----- static functions ---------------------------------------- */
static int ascending(const void *a, const void*b) {
  if (*(double*)a==*(double*)b) return 0;
  return (*(double*)a > *(double*)b ) ? 1 : -1;
}

/* ----- ut_median ---------------------------------------- */
void ut_sort_ascending(double* values, int n)
     /* median of unsorted values */
{
  qsort(values,n,sizeof(double),&ascending);
}

/* ----- ut_median ---------------------------------------- */
double ut_median(double* values, int n)
     /* median of unsorted values */
{
  qsort(values,n,sizeof(double),&ascending);
  /* alternatively, I could have used  gsl_stats_median_from_sorted_data
     but it is fucking the flies ! */
  /* i saw also something usable in mathlib : the names might conflict !*/
  if (n%2)
    return values[n/2];
  else
    return (values[n/2-1]+values[n/2])/2;

}

/* ----- ut_min ---------------------------------------- */
double ut_min(double* values, int n)
     /* min of array */
{
  double min = values[0];
  int i;

  for (i=1;i<n;i++) {
    if (values[i]<min)
      min=values[i];
  }
  return min;
}

/* ----------------------------- ut_trunc_sigma_known -------------- */
void ut_trunc_sigma_known(double** values, int * nitems, double sigma, double sigcut){
     /* trunc the data set, removing any value outside +/- sigcut sigma
       if sigcut<0, does noting .
       it works iteratively by removing too far away values. => slow if
       the sigcut is too low !

       values : pointer to array of values - modified
       nitems : number of values -> modified
       sigcut : if <0, does nothing

       Beware : the number of remaining items might be 0 !

     */
  int ok=0;

  if (sigcut<0)
    return;

  qsort(*values,*nitems,sizeof(double),ascending);

  while (!ok && *nitems>1) {

    double mean=ut_mean(*values,*nitems);
    double valinf = (mean-(*values)[0])/sigma;
    double valsup = ((*values)[*nitems-1]-mean)/sigma;
    if (valsup>sigcut && valsup>=valinf) {
      double bound = MAX(sigcut,valinf);
      /* remove more than one for efficiency */
      while(valsup>=bound && *nitems>1) {
        (*nitems)--;
        valsup = ((*values)[*nitems-1]-mean)/sigma;
      }
    } else if (valinf>sigcut && valinf>=valsup) {
      double bound = MAX(sigcut,valinf);
      while(valinf>=bound && *nitems>1 ){
        (*nitems)--;
        (*values)++;
        valinf = (mean-(*values)[0])/sigma;
      }
    } else {
      ok=1;
    }
  }
}

/* ------------------------------- ut_trunc_sigma_unknown ------------------ */
void ut_trunc_sigma_unknown(double** values, int * nitems, double sigcut){
     /* trunc the data set, removing any value outside +/- sigcut sigma
       if sigcut<0, does noting .
       it works iteratively by removing too far away values. => slow if
       the sigcut is too low !
       Beware : sigma is computed from the distribution itself.
       Points far away will push the sigma towards greater values.

       values : pointer to array of values - modified
       nitems : number of values -> modified
       sigcut : if <0, does nothing

       Beware : the number of remaining items might be 0 !

     */
  int ok=0;

  if (sigcut<0)
    return;

  qsort(*values,*nitems,sizeof(double),ascending);

  while (!ok && *nitems>1) {
    double mean,sigma,valinf,valsup;

    mean = ut_mean((*values),*nitems);
    sigma = gsl_stats_sd_m(*values,1,*nitems,mean);
    valinf = (mean-(*values)[0])/sigma;
    valsup = ((*values)[*nitems-1]-mean)/sigma;
    if (valsup>sigcut && valsup>=valinf) {
      (*nitems)--;
    } else if (valinf>sigcut && valinf>=valsup) {
      (*nitems)--;
      (*values)++;
    } else {
      ok=1;
    }
  }
}

/* ------------------------------- ut_trunc_sigma_unknown ------------------ */
void ut_trunc_sigma_unknown_fast_sorted(double** values, int * nitems, double sigcut){
     /* trunc the data set, removing any value outside +/- sigcut sigma
       if sigcut<0, does noting .
       it works grossly, removing each time all what is outside the window,
       which is not accurate

       Beware : sigma is computed from the distribution itself.
       Points far away will push the sigma towards greater values.

       values : pointer to array of values - modified
       nitems : number of values -> modified
       sigcut : if <0, does nothing

       Beware : the number of remaining items might be 0 !

     */
  int ok=0;

  if (sigcut<0)
    return;

  /*  qsort(*values,*nitems,sizeof(double),ascending); */

  while (!ok && *nitems>1) {

    double mean,sigma,valinf,valsup;
    int nitems_old = *nitems;

    mean = ut_mean((*values),*nitems);
    sigma = gsl_stats_sd_m(*values,1,*nitems,mean);
    valinf = mean-sigma*sigcut;
    valsup = mean + sigma*sigcut;
    while ((*values)[0] < valinf && *nitems>1) {
      (*nitems)--;
      (*values)++;
    }
    while ((*values)[*nitems-1] > valsup && *nitems>1) {
      (*nitems)--;
    }
    if (nitems_old == (*nitems))
      ok=1;
  }
}

/*-------------------- ut_fraction_sigcut ----------------------*/
double ut_fraction_sigcut(double fraction) {
  /* given a distribution with fraction outliers, at infinite
   distance with respect to normal values,
  returns the standard deviation of the outliers, when the
  RMS is computed with it included.*/

  /* the computation is performed with fraction values at 1 and
   1-fraction at 0 */
  double sigcut,rms;
  /* mean = fraction */
  rms = sqrt(fraction*fraction*(1-fraction) + (1-fraction)*(1-fraction)*fraction);

  sigcut = (1-fraction)/rms;
  return sigcut;
}

/*-------------------- ut_mean ----------------------*/
double ut_mean(double* val,int n){

  int i;
  double sum = 0;
  for (i=0;i<n;i++)
    sum += val[i];
  return sum/n;
}

/*-------------------- ut_mean ----------------------*/
double ut_mode(double* val,int n){
  /* This function returns the highest probability density a set of
  datas. This only works reliabely if there is only 1 maximum,
  and in the limit where the data is sufficiently sampled.
  No error computation is performed ...*/
  /* The search is dyadic and iterative, each step looking for the part
  where the median is the narrower. The time is O(NlnN) (+ the time to run
  indexx)*/
  /* This is a fast and inaccurate algorithm. For a general algorithm, see the
  Parzenwindow method - the drawback is that it depends on a window size
  a priori*/

  size_t * index;
  double med;
  int nmin=0;
  int nmax=n; /*  out of the table */

  index = malloc(n*sizeof(size_t));

  /* test */

  gsl_sort_index(index,val,1,n);

  do {

    /* median computation */
    if ((nmax-nmin) % 2 == 0)
      med = (val[index[ (nmin + nmax )/2 ]] + val[index[ (nmin+nmax)/2 -1 ]])/2;
    else
      med = val[index[ (nmin + nmax)/2 ]];

    /* The density is higher for the upper part */
    if ( fabs(med - val[index[nmin]])  > fabs (med - val[index[nmax-1]]) )
      nmin =(nmin + nmax - 1 )/2;
    else if ( fabs(med - val[index[nmin]])  < fabs (med - val[index[nmax-1]]) )
      /* rounding in the good direction  in order to include the median place */
      nmax =  ( nmin + nmax)/2 + 1;
    else /* stop in case of equality */
      nmax = nmin;
  } while (nmax-nmin > 2);

  return med;
}

/*-------------------- ut_mean ----------------------*/
double ut_mode2(double* val,int n){
  /* This function returns the highest probability density a set of
  datas. This only works reliabely if there is only 1 maximum,
  and in the limit where the data is sufficiently sampled.
  No error computation is performed ...*/
  /* The search is dyadic and iterative, each step looking for the part
  which contains most of the data with respect to the highest value.
  The time is O(NlnN) (+ the time to run indexx)*/
  /* This is a fast and inaccurate algorithm. For a general algorithm, see the
  Parzenwindow method - the drawback is that it depends on a window size
  a priori*/

  size_t * index;
  double half;
  int nmin=0;
  int nmax=n; /*  out of the table */

  index = malloc(n*sizeof(size_t));

  /* test */

  gsl_sort_index(index,val,1,n);

  do {
    int i;
    /* half window computation */
    half=(val[index[nmin]]+val[index[nmax-1]])/2;
    for (i=nmin;i<nmax;i++) {
      if (val[index[i]]>half)
        break;
    }
    if (i==nmax)
      return half;
    else if (i<nmin+(nmax-nmin)/2)
      nmin=i;
    else if (i>nmin+(nmax-nmin)/2 || ((nmax-nmin)%2) )
      nmax=i;
    else {
      nmin++;
      nmax--;
    }
  } while (nmax-nmin>2);
  return half;
}

/*-------------------- ut_autocovariance ----------------------*/
void ut_autocovariance(double* data, double* autocorr, int n) {
/* data is of dim n, autocorr of dim n */
  int decal,i;
  double mean;
  /*  double var; */
  double *mdata = malloc(sizeof(double)*n);

  /* compute the mean */
  mean = ut_mean(data,n);

  /* substract it (n computations) */
  for (i=0;i<n;i++)
    mdata[i] = data[i]-mean;

  /* compute variance and what is needed for decal=0*/
  decal=0;
  autocorr[0]=0;
  for (i=0;i<n;i++) {
    autocorr[0]+=mdata[i]*mdata[i];
  }
  autocorr[0] /= n;
  /* var = autocorr[0]; */
  /* autocorr[0] = 1.0; */

  /* n^2 part */
  for (decal=1; decal<n; decal ++) {
    autocorr[decal]=0;
    for (i=0;i<n-decal;i++) {
      autocorr[decal]+=mdata[i]*mdata[i+decal];
    }
    autocorr[decal] /= (n-decal);
    /*    autocorr[decal] /= var; */
  }
  free(mdata);

}

/*-------------------- ut_varname_from_imname ----------------------*/
void ut_varname_from_imname(char* imname,char* varname){
  /* returns the name of the variance image from the name of
   the image */

  /* check standard name */
  char extname[lg_name+1];
  extname[0]='\0';
  varname[0]='\0';
  if (strstr(imname,"[chip00]"))
    strcpy(extname,"[var00]");
  if (strstr(imname,"[chip01]"))
    strcpy(extname,"[var01]");
  if (strstr(imname,"[chip02]"))
    strcpy(extname,"[var02]");
  if (strstr(imname,"[chip03]"))
    strcpy(extname,"[var03]");
  if (strstr(imname,"[image]"))
    strcpy(extname,"[variance]");

  if (extname[0]) {
    /* builds the variance name */
    char *ptName;
    strcpy(varname,imname);
    ptName = strchr(varname,'[');
    strcpy(ptName,extname);
  }
}

/*-------------------- ut_open_check_name ----------------------*/
char* ut_open_check_name(char* name){
  /* if name is a name without extension, tries to see if a
     [image] exists */

  if (!strchr(name,'[')) {
    char new_name[lg_name+1];
    sprintf(new_name,"%s[image]",name);
    if (exist(new_name))
      strcpy(name,new_name);
  }
  return name;
}

/*-------------------- ut_open_check_name ----------------------*/
char* ut_create_check_name(char* name){
  /* if name is a name without extension, appends an [image]
     by default */

  if (!strchr(name,'[')) {
    char new_name[lg_name+1];
    sprintf(new_name,"%s[image]",name);
    strcpy(name,new_name);
  }
  return name;
}

/*-------------------- ut_varname_from_imname ----------------------*/
int ut_is_bichip_detcom(char* filename){
  /* returns 1 if the name corresponds to a detcom image */

  char testName[lg_name+1];
  if (strstr(filename,"%d"))
    return 1;
  if (strrchr(filename,'['))
    return 0;
  // .fits file extension default
  if (strstr(filename,".fits"))
    sprintf(testName,"%s[chip00]",filename);
  else
    sprintf(testName,"%s.fits[chip00]",filename);
  if (exist(testName))
    return 1;
  else {
    if (strstr(filename,".fits"))
      sprintf(testName,"%s[CHIP00]",filename);
    else
      sprintf(testName,"%s.fits[CHIP00]",filename);
    if (exist(testName))
      return 1;
    else
      return 0;
  }
}

/*-------------------- ut_build_tmp_name ----------------------*/
void ut_build_tmp_name(char* filename, char* tmp_prefix){
  /* builds a temporary memory file name from the image name and a
  prefix */
  char imname[lg_name+1], *pt_name;

  if (!(pt_name = strrchr(filename,'/')))
    pt_name = filename;
  else
    pt_name = pt_name+1;
  sprintf(imname,"mem://%s_%s",tmp_prefix,pt_name);
  //  sprintf(imname,"%s_%s",tmp_prefix,pt_name);
  strcpy(filename,imname);
}

/*-------------------- ut_primary_header_name ----------------------*/
void ut_primary_header_name(char* full_name, char* primary_name){
  /* the primary header is the image name without any extension*/

  char* pt_name;
  strcpy(primary_name,full_name);

  if ((pt_name = strchr(primary_name,'[')))
    pt_name[0]='\0';

}

/* ===== IFU lib disagreement ======================================== */

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!.func                        open_frame_fast()
!
!.purp          opens a 2D frame and updates the image structure items
!.desc
! open_frame(frame,name,mode)
!
! IMAGE2D *frame;       image structure
! char *name;           frame name
! char *mode;           open mode (Input,Ouput,IO)
!.ed
-------------------------------------------------------------------- */
int
open_frame_fast(IMAGE2D *frame, char *name, char *mode)
{
  char errtext[132], filename[lg_name+1];
  int status, nbaxes, iomode, int_datatype;
  float cuts[4];
  int info[5];
#ifdef IRAF
  int two_dim=2;
  int len;
#endif
#ifdef FITS
  fitsfile *fptr;
  int nbread;
  int npix;
  int group = 0;
  double pixref;
#endif

  memset(frame->ident,' ',lg_ident);
  frame->ident[lg_ident] = '\0';
  memset(frame->cunit,' ',lg_unit);
  frame->cunit[lg_unit] = '\0';
  memset(frame->history,' ',lg_hist);
  frame->history[lg_hist] = '\0';
  frame->external_info = NULL;
  frame->file_type = T_IMA2D;
  frame->data_format = InputIO.basic_io;

  strcpy(filename,name);
  first_blk(filename);
  strcpy(frame->name,filename);
  append_ima_extension(frame->name,InputIO.basic_io);

  strcpy(filename,frame->name);

  if (!exist(filename)) { /* check if fil exists */
    status = ERR_OPEN;
    sprintf(errtext,"open_frame: frame %s",filename);
    Handle_Error(errtext,status);
    return(status);
  }

  switch(mode[0]) {
  case 'I' :
    if (mode[1] == 'O')
      frame->iomode = (int)IO_MODE;
    else
      frame->iomode = (int)I_MODE;
    break;
  case 'O' : frame->iomode = (int)O_MODE;
    break;
  default  : frame->iomode = (int)I_MODE;
    break;
  }

  iomode = get_iomode_code(InputIO.basic_io,frame->iomode);

  switch (InputIO.basic_io) {

#ifdef MIDAS
  case MIDAS_FORMAT :
    status = SCFINF(filename,2,info);
    if (status == 0) {
      status = SCIGET(filename, info[1], iomode, F_IMA_TYPE, 2,
                      &nbaxes, &(frame->nx), &(frame->startx), &(frame->stepx),
                      frame->ident, frame->cunit, (char **)(&(frame->data)),
                      &(frame->imno));
      frame->data_type = info[1];
      frame->data_type = decode_datatype(InputIO.basic_io,frame->data_type);

      if (nbaxes!=2) /* We open a spectrum like an image, and that's not good */
        status = ERR_OPEN;

    }
    break;
#endif
#ifdef IRAF
  case IRAF_FORMAT :
  case STSDAS_FORMAT :
    len = strlen(filename);
    uimopn(filename,&iomode,&(frame->imno),&status,len);
    if (status != 0)
      break;
    uimgid(&(frame->imno),&int_datatype,&two_dim,&(frame->nx),&status);
    frame->data_type = decode_datatype(InputIO.basic_io,(short)(int_datatype));
    if (status != 0)
      break;
    alloc_frame_mem(frame, datatype);
    switch(frame->data_type) {
    case SHORT :
      uigs2s(&(frame->imno),&one,&(frame->nx),&one,&(frame->ny),
             frame->data.s_data,&status);
      break;
    case INT :
    case LONG :
      uigs2l(&(frame->imno),&one,&(frame->nx),&one,&(frame->ny),
             frame->data.l_data,&status);
      break;
    case FLOAT :
      uigs2r(&(frame->imno),&one,&(frame->nx),&one,&(frame->ny),
             frame->data.f_data,&status);
      break;
    case DOUBLE :
      uigs2d(&(frame->imno),&one,&(frame->nx),&one,&(frame->ny),
             frame->data.d_data,&status);
      break;
    }
    disable_user_warnings();
    RD_desc(frame,"IDENT",CHAR,lg_ident,frame->ident);
    restore_user_warnings();
    break;
#endif
#ifdef FITS
  case FITS_A_FORMAT :
  case FITS_B_FORMAT :
    status =0;
    if (fits_open_file(&fptr,filename,iomode,&status)) {
      status = ERR_ACCESS; break;
    }
    frame->external_info = (void *)fptr;
    if (fits_read_key(fptr, TINT,"NAXIS", &nbaxes,NULL, &status)) {
      status = ERR_READ; break;
    }
    if (nbaxes != 2) {
      status = ERR_IMA_HEAD; break;
    }
    if (fits_read_key(fptr, TINT, "NAXIS1",
                      &(frame->nx), NULL, &status)) {
      status = ERR_READ; break;
    }
    if (fits_read_key(fptr, TINT, "NAXIS2",
                      &(frame->ny), NULL, &status)) {
      status = ERR_READ; break;
    }
    if (status == 0) {
      pixref = 1.0;
      fits_read_key(fptr, TDOUBLE, "CRPIX1", &pixref, NULL, &status);
      if (status) { status = 0; pixref = 1; }
      fits_read_key(fptr, TDOUBLE, "CRVAL1", &(frame->startx), NULL, &status);
      if (status) { status = 0; frame->startx = (double)1; }
      fits_read_key(fptr, TDOUBLE, "CDELT1", &(frame->stepx), NULL, &status);
      if (status) { status = 0; frame->stepx = (double)1; }
      frame->startx -= (pixref-1)*frame->stepx;
      pixref = 1.0;
      fits_read_key(fptr, TDOUBLE, "CRPIX2", &pixref, NULL, &status);
      if (status) { status = 0; pixref = 1; }
      fits_read_key(fptr, TDOUBLE, "CRVAL2", &(frame->starty), NULL, &status);
      if (status) { status = 0; frame->starty = (double)1; }
      fits_read_key(fptr, TDOUBLE, "CDELT2", &(frame->stepy), NULL, &status);
      if (status) { status = 0; frame->stepy = (double)1; }
      frame->starty -= (pixref-1)*frame->stepy;
    }
    else
      break;

    int_datatype = (fptr->Fptr->tableptr)->tdatatype;
    frame->data_type = decode_datatype(InputIO.basic_io,(short)int_datatype);
    if (frame->data_type == SHORT) {
      if (fptr->Fptr->tableptr[1].tscale == 1 && fptr->Fptr->tableptr[1].tzero == 32768)
        /* unsigned short !!! */
        frame->data_type = LONG;
    }

    if (alloc_frame_mem(frame, frame->data_type) < 0) {
      fits_close_file(fptr,&status);
      status = ERR_ALLOC;
      break;
    }

    npix = frame->nx*frame->ny;
    switch (frame->data_type) {
    case SHORT :
      if (fits_read_img_sht(fptr,group,1L,npix,(short)0,
                            frame->data.s_data,&nbread,&status)) {
        status = ERR_READ; break;
      }
      break;
    case LONG :
    case INT :
      if (fits_read_img_lng(fptr,group,1L,npix,(int)0,
                            frame->data.l_data,&nbread,&status)) {
        status = ERR_READ; break;
      }
      break;
    case FLOAT :
      if (fits_read_img_flt(fptr,group,1L,npix,(float)0,
                            frame->data.f_data,&nbread,&status)) {
        status = ERR_READ; break;
      }
      break;
    case DOUBLE :
      if (fits_read_img_dbl(fptr,group,1L,npix,(double)0,
                            frame->data.d_data,&nbread,&status)) {
        status = ERR_READ; break;
      }
      break;
    }
    break;
#endif
  }

  if (status) {
    sprintf(errtext,"open_frame: frame %s",filename);
    status = get_tiger_errcode(frame->data_format,status);
    Handle_Error(errtext,status);
  }
  else {
    disable_user_warnings();
    status = RD_desc(frame,"LHCUTS",FLOAT,4,cuts);
    RD_desc(frame,"HISTORY",CHAR,lg_hist,frame->history);
    restore_user_warnings();

    frame->endx = frame->startx + (frame->nx -1)*frame->stepx;
    frame->endy = frame->starty + (frame->ny -1)*frame->stepy;
    if (status <= 0) {
      /* image_minmax is a really slow routine, and most of the time useless */
      frame->min = -ut_big_value;
      frame->max = +ut_big_value;
      /*			image_minmax(frame); */
    }
    else {
      frame->min = cuts[2];
      frame->max = cuts[3];
    }
    status = 0;
    /* parse wcs if contained in file */
    status = parse_wcs(frame);
  }
  return(status);
}

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!.func                       close_frame_fast()
!
!.purp             closes a currently active 2D frame
!.desc
! close_frame(frame)
!
! IMAGE2D *frame;       image structure
!.ed
-------------------------------------------------------------------- */
int
close_frame_fast(IMAGE2D *frame)			/* close active frame */
{
  char  errtext[132], filename[lg_name+1];
  int   stat, int_datatype;
  float cuts[4];
#ifdef IRAF
  int one=1;
#endif
#ifdef FITS
  fitsfile *fptr;
  int npix;
#endif

  strcpy(filename,frame->name);

  if (frame->iomode == (int)I_MODE) {
    switch (frame->data_format) {
#ifdef MIDAS
    case MIDAS_FORMAT :
      stat = SCFCLO(frame->imno);
      break;
#endif
#ifdef IRAF
    case IRAF_FORMAT :
    case STSDAS_FORMAT :
      uimclo(&(frame->imno),&stat);
      break;
#endif
#ifdef FITS
    case FITS_A_FORMAT :
    case FITS_B_FORMAT :
      stat =0;
      fptr = (fitsfile *)frame->external_info;
      fits_close_file(fptr,&stat);
      free_frame_mem(frame);
      frame->external_info = NULL;
      break;
#endif
    }
    if (stat) {
      sprintf(errtext,"close_frame: frame %s",filename);
      stat = get_tiger_errcode(frame->data_format,stat);
      Handle_Error(errtext,stat);
    }
    return(stat);
  }

  /*
  if (frame->data.d_data != NULL) {
    image_minmax(frame);

    cuts[0]=(float)frame->min; cuts[2]=(float)frame->min;
    cuts[1]=(float)frame->max; cuts[3]=(float)frame->max;
    stat = WR_desc(frame,"LHCUTS",FLOAT,4,cuts);
  }
  */

  WR_history(frame, (Anyfile *)0);

  switch (frame->data_format) {
#ifdef MIDAS
  case MIDAS_FORMAT :
    stat = SCFCLO(frame->imno);
    break;
#endif
#ifdef IRAF
  case IRAF_FORMAT :
  case STSDAS_FORMAT :
    switch(frame->data_type) {
    case SHORT :
      uips2s(&(frame->imno),&one,&(frame->nx),&one,&(frame->ny),
             frame->data.s_data,&stat);
      break;
    case INT :
    case LONG :
      uips2l(&(frame->imno),&one,&(frame->nx),&one,&(frame->ny),
             frame->data.l_data,&stat);
      break;
    case FLOAT :
      uips2r(&(frame->imno),&one,&(frame->nx),&one,&(frame->ny),
             frame->data.f_data,&stat);
      break;
    case DOUBLE :
      uips2d(&(frame->imno),&one,&(frame->nx),&one,&(frame->ny),
             frame->data.d_data,&stat);
      break;
    }
    if (stat == 0)
      uimclo(&(frame->imno),&stat);
    free_frame_mem(frame);
    break;
#endif
#ifdef FITS
  case FITS_A_FORMAT :
  case FITS_B_FORMAT :
    stat = 0;
    fptr = (fitsfile *)frame->external_info;
    if (frame->iomode != (int)I_MODE) {
      if (frame->data.d_data != NULL) {
        int_datatype = get_datatype_code(OutputIO.basic_io,frame->data_type);
        npix = frame->nx*frame->ny;
        if (fits_write_img(fptr,int_datatype,1L,npix,
                           frame->data.s_data,&stat)) {
          stat = ERR_WRIT;
        }
      }
    }
    if (! stat) {
      fits_close_file(fptr,&stat);
      stat = wcs_free(frame);
    }
    free_frame_mem(frame);
    frame->external_info = NULL;
    break;
#endif
  }
  if (stat) {
    sprintf(errtext,"close_frame: frame %s",filename);
    stat = get_tiger_errcode(frame->data_format,stat);
    Handle_Error(errtext,stat);
  } else {
    if (TK && (frame->iomode == O_MODE || frame->iomode == IO_MODE))
      {
        printf("@ N {%s}\n",filename);
      }
  }

  return(stat);
}

/* ----- juldat --------------------------------------------------------- */
/*C
C       FUNCTION JULDAT RETURNS THE JULIAN DATE AS A DOUBLE
C       PRECISION NUMBER.
C       ALGORITHM FROM "ALAMANAC FOR COMPUTERS" p.B2
C*/
double juldat( int year, int month, int day, double ut)
     /* year, month, day, ut are the UTC ones
      ut is the fractional hour */
{
        double tmp;
        tmp = 367.0*year + 0.5 + ut/24.0;
        tmp = tmp - ((7*(year+((month+9)/12)))/4) + ((275*month)/9) + day + 1721013;
        return tmp;
}

/* =====  string utilities ===== */

int ut_parse_line(char* line, char** tokens) {
  /* the token strings will contain the line elements separated by spaces */
  /* the original line will be modified (some '/0' char added)*/
  char* pt_line = line;
  int itoken=0;

  while (pt_line && pt_line[0]!='\0') {
    tokens[itoken++] = pt_line;
    pt_line = strchr(pt_line,' ');
    if (pt_line)
      while((pt_line)[0]==' ') {
        pt_line[0]='\0';
        pt_line++;
      }
  }
  return itoken;
}

max_stars_repo_name: snfactory/preprocess

max_issues_repo_name: snfactory/preprocess

max_forks_repo_name: snfactory/preprocess


text:
/* Copyright 2018,2019, Alistair Boyle, 3-clause BSD License */
#include "config.h"
#include <string.h> /* memset */
#include <stdlib.h> /* free */
#include <lapacke.h> /* inv: dgetrf, dgetri */
#include <math.h> /* fabs, sqrt */
#include <limits.h> /* MAX_INT */
#include <assert.h> /* assert */
#include <stdio.h> /* printf */

#include "matrix.h"
#include "model.h"

/* references:
 * [1] A. Boyle, PhD Thesis, 2016, Geophysical Applications of Electrical Impedance Tomography, Carleton University
 */

#ifdef UNIT_TESTING
extern void * _mock_test_malloc(const size_t size, const char * file, const int line);
extern void _test_free(void * const ptr, const char * file, const int line);
#define malloc(size) _mock_test_malloc(size, __FILE__, __LINE__)
#define free(ptr) _test_free(ptr, __FILE__, __LINE__)
#endif

#define SUCCESS 1
#define FAILURE 0

static void init_mesh(mesh * m)
{
    memset(m, 0, sizeof(mesh));
}

static void init_model(model * m)
{
    memset(m, 0, sizeof(model));
}

model * malloc_model()
{
    model * m = malloc(sizeof(model));
    if(m) {
        init_model(m);
    }
    return m;
}

static void free_mesh(mesh * m)
{
    free(m->nodes);
    free(m->elems);
    free(m->matidx);
    free(m->surfaceelems);
    free(m->bc);
    free(m->pmap_param);
    free(m->pmap_elem);
    free(m->pmap_frac);
}

void reset_model(model * m)
{
    assert(m != NULL);
    free(m->data);
    free(m->params);
    free(m->stimmeas);
    free(m->measgain);
    free(m->elec_to_sys);
    free(m->zcbc);
    free(m->zc);
    free_mesh(&(m->fwd));
    free_mesh(&(m->rec));
    init_model(m);
}

void * free_model(model * m)
{
    if(m == NULL) {
        return NULL;
    }
    reset_model(m);
    free(m);
    return NULL;
}

int set_model_data(model * m, int rows, int cols)
{
    assert(rows > 0);
    assert(cols > 0);
    assert(m != NULL);
    m->n_data[0] = rows;
    m->n_data[1] = cols;
    free(m->data);
    m->data = malloc(sizeof(double) * rows * cols);
    if(m->data == NULL) {
        return FAILURE;
    }
    return SUCCESS;
}

int set_model_params(model * m, int rows, int cols)
{
    assert(rows > 0);
    assert(cols > 0);
    assert(m != NULL);
    m->n_params[0] = rows;
    m->n_params[1] = cols;
    free(m->params);
    m->params = malloc(sizeof(double) * rows * cols);
    if(m->params == NULL) {
        return FAILURE;
    }
    return SUCCESS;
}

int set_model_stimmeas(model * m, int rows)
{
    assert(rows > 0);
    assert(m != NULL);
    m->n_stimmeas = rows;
    m->n_elec = 0;
    free(m->elec_to_sys);
    m->elec_to_sys = NULL;
    free(m->stimmeas);
    m->stimmeas = malloc(sizeof(int) * rows * 4);
    if(m->stimmeas == NULL) {
        return FAILURE;
    }
    free(m->measgain);
    m->measgain = malloc(sizeof(double) * rows);
    if(m->measgain == NULL) {
        return FAILURE;
    }
    return SUCCESS;
}

void set_mesh_dim(mesh * m, int dim)
{
    assert(dim == 2 || dim == 3);
    assert(m != NULL);
    free_mesh(m);
    init_mesh(m);
    m->dim = dim;
}

int set_mesh_nodes(mesh * m, int n_nodes)
{
    assert(n_nodes >= 0);
    assert(m != NULL);
    assert(m->dim == 2 || m->dim == 3);
    m->n_nodes = n_nodes;
    free(m->nodes);
    m->nodes = malloc(sizeof(double) * n_nodes * m->dim);
    if(m->nodes == NULL) {
        return FAILURE;
    }
    return SUCCESS;
}

int set_mesh_elems(mesh * m, int n_elems)
{
    assert(n_elems >= 0);
    assert(m != NULL);
    assert(m->dim == 2 || m->dim == 3);
    m->n_elems = n_elems;
    free(m->elems);
    m->elems = malloc(sizeof(int) * n_elems * (m->dim + 1));
    if(m->elems == NULL) {
        return FAILURE;
    }
    free(m->matidx);
    m->matidx = malloc(sizeof(int) * n_elems);
    if(m->matidx == NULL) {
        return FAILURE;
    }
    return SUCCESS;
}

int set_mesh_surfaceelems(mesh * m, int n_se)
{
    assert(n_se >= 0);
    assert(m != NULL);
    assert(m->dim == 2 || m->dim == 3);
    m->n_se = n_se;
    free(m->surfaceelems);
    m->surfaceelems = malloc(sizeof(int) * n_se * m->dim);
    if(m->surfaceelems == NULL) {
        return FAILURE;
    }
    free(m->bc);
    m->bc = malloc(sizeof(int) * n_se);
    if(m->bc == NULL) {
        return FAILURE;
    }
    return SUCCESS;
}

int set_mesh_pmap(mesh * m, int n_pmap)
{
    assert(n_pmap >= 0);
    assert(m != NULL);
    m->n_pmap = n_pmap;
    free(m->pmap_param);
    m->pmap_param = malloc(sizeof(int) * n_pmap);
    if(m->pmap_param == NULL) {
        return FAILURE;
    }
    free(m->pmap_elem);
    m->pmap_elem = malloc(sizeof(int) * n_pmap);
    if(m->pmap_elem == NULL) {
        return FAILURE;
    }
    free(m->pmap_frac);
    m->pmap_frac = malloc(sizeof(double) * n_pmap);
    if(m->pmap_frac == NULL) {
        return FAILURE;
    }
    return SUCCESS;
}

int set_model_zc(model * m, int rows)
{
    assert(rows > 0);
    assert(m != NULL);
    m->n_zc = rows;
    free(m->zc);
    m->zc = malloc(sizeof(double) * rows);
    if(m->zc == NULL) {
        return FAILURE;
    }
    free(m->zcbc);
    m->zcbc = malloc(sizeof(int) * rows);
    if(m->zcbc == NULL) {
        return FAILURE;
    }
    return SUCCESS;
}

int calc_elec_zc_map(model * m)
{
    assert(m != NULL);
    assert(m->n_elec > 0);
    double * tmp = malloc(sizeof(double) * m->n_elec);
    if(tmp == NULL) {
        return FAILURE;
    }
    if((m->n_zc == m->n_elec) && (m->zcbc == NULL)) {
        free(tmp);
        return SUCCESS;
    }
    for(int i = 0; i < m->n_elec; i++) {
        tmp[i] = 1e-2; /* default zc = 0.01 Ω/m² */
    }
    for(int i = 0; i < m->n_zc; i++) {
        /* TODO for now we assume that boundary condition number = electrode# */
        const int elec = m->zcbc[i];
        assert(elec > 0);
        assert(elec <= m->n_elec);
        tmp[elec - 1] = m->zc[i];
    }
    free(m->zcbc);
    free(m->zc);
    m->zcbc = NULL;
    m->zc = tmp;
    m->n_zc = m->n_elec;
    return SUCCESS;
}

double * inv(int n, double A[n][n])
{
    double * Ap = &(A[0][0]);
    int ipiv[n + 1];
    memset(ipiv, 0, (n + 1)*sizeof(int));
    if(LAPACKE_dgetrf(LAPACK_ROW_MAJOR, n, n, Ap, n, ipiv) ||
       LAPACKE_dgetri(LAPACK_ROW_MAJOR, n, Ap, n, ipiv)) {
        return NULL;
    }
    return Ap;
}

double det(int n, double A[n][n])
{
    if(n == 2) {
        return A[0][0] * A[1][1] - A[0][1] * A[1][0];
    }
    double d = 0;
    double sign = 1.0;
    for(int i = 0; i < n; i++, sign *= -1.0) {
        /* construct sub-matrix B */
        double B[n - 1][n - 1];
        for(int j = 1; j < n; j++) {
            for(int k = 0, s = 0; k < n; k++) {
                if(k != i) {
                    B[j - 1][s] = A[j][k];
                    s++;
                }
            }
        }
        d += sign * A[0][i] * det(n - 1, B); /* (-1)^n * a1i * |B| */
    }
    return d;
}

/* shape function
 * Se from [1] eq (C.31) -- first-order tetrahedral/triangular elements
 * Se is symmetric, only the lower triangular portion is calculated (E1^T E1 is symmetric)
 * inputs: nd = # of dimensions (2=2D, 3=3D),
 *         N = ptrs to node coordinates (nd+1 nodes with nd double precision values)
 *             ex: {&n[120], &n[332], &n[556], &n[733]} for 3D, where n[120] = {1.1, 0.1, 3.2}, etc.
 *             NOTE: nodes should be in increasing order so that calc_sys_elem_ij
 *                   gives an upper triangular matrix in the global node numbering
 * output: Se = replaces Se[0..t-1] for t = n(n+1)/2 entries with n=nd+1 elements as
 *         upper triangular symmetric matrix, row-major order
 * returns: NULL if matrix inverse fails, Se otherwise
 */
static double calc_elem_vol(const int n, double const * const * const nodes, double E[n][n])
{
    for(int i = 0; i < n; i++) {
        E[i][0] = 1;
        for(int j = 1; j < n; j++) {
            E[i][j] = nodes[i][j - 1];
        }
    }
    if(inv(n, E) == NULL) {
        return -1.0;
    }
    return fabs(det(n, E));
}

static double * calc_sys_elem_v(const int nd, double const * const * const nodes, double * Se)
{
    const int n = nd + 1;
    double E[n][n]; /* from [1] eq (C.33) */
    const double detE = calc_elem_vol(n, nodes, E);
    if(detE < 0) {
        return NULL;
    }
    int ndf = 1; /* = !nd = factorial(number of dimensions) */
    for(int i = 2; i <= nd; i++ ) {
        ndf *= i;
    }
    const double c = 1 / (ndf * detE);
    /* Se = c * E1^T E1, where E1 is E without the top row */
    int m = 0;
    for(int i = 0; i < n; i++ ) {
        for(int j = i; j < n; j++ ) {
            Se[m] = 0;
            for(int k = 1; k < n; k++) {
                Se[m] += E[k][i] * E[k][j];
            }
            /* TODO pull this out and multiply by conductivity D, to save some multiply operations */
            Se[m] *= c;
            m++;
        }
    }
    return Se;
}

static void calc_sys_elem_ij(const int nd, int const * const elem, int * ii, int * jj)
{
    const int n = nd + 1;
    int idx = 0;
    for(int i = 0; i < n; i++ ) {
        for(int j = i; j < n; j++ ) {
            int ei = elem[i] - 1;
            int ej = elem[j] - 1;
            /* swap indices (symmetric matrix) if this entry would be in the
             * upper triangular portion of the matrix */
            if(ei < ej) {
                int tmp = ei;
                ei = ej;
                ej = tmp;
            }
            ii[idx] = ei;
            jj[idx] = ej;
            idx++;
        }
    }
}

int calc_sys_elem_n(const int nd)
{
    /* nd=2: 6, nd=3: 10 ... for symmetric (isotropic conductivity)*/
    const int n = nd + 1;
    return n * (n + 1) / 2;
}

/* returns 0 on success, > 0 on failure as the singular element number e+1 */
int calc_sys_elem(mesh const * const m, int * ii, int * jj, double * Se)
{
    const int nd = m->dim;
    const int n_elems = m->n_elems;
    const int * elems = m->elems;
    const double * nodes = m->nodes;
    const int n = calc_sys_elem_n(nd);
    for(int i = 0; i < n_elems; i++) {
        calc_sys_elem_ij(nd, elems, ii, jj);
        double const * node_list [4];
        for(int j = 0; j < nd + 1; j ++) {
            int idx = elems[j] - 1;
            if ( idx < 0 ) {
                idx = 0;
            }
            node_list[j] = &(nodes[idx * nd]);
        }
        if(calc_sys_elem_v(nd, node_list, Se) == NULL) {
            return i + 1;
        }
        elems += (nd + 1);
        ii += n;
        jj += n;
        Se += n;
    }
    return 0;
}

void calc_sys_gnd(const int gnd, size_t nnz, int * ii, int * jj, double * Se)
{
    const int gndidx = gnd - 1;
    int ret = 0;
    for (int i = 0; i < nnz; i++ ) {
        const int is_row = (ii[i] == gndidx);
        const int is_col = (jj[i] == gndidx);
        if ( is_row || is_col) {
            if(ret == 0) {
                /* replace with (gndidx,gndidx)=1.0 */
                ii[i] = gndidx;
                jj[i] = gndidx;
                Se[i] = 1.0;
            }
            else {
                /* drop this entry */
                Se[i] = 0.0;
            }
            ret++;
        }
    }
    assert(ret != 0);
}

#define dot2(a,b) (((a[0])*(b[0])) + ((a[1])*(b[1])))
#define dot3(a,b) (((a[0])*(b[0])) + ((a[1])*(b[1])) + ((a[2])*(b[2])))
#define sum_sq2(a) dot2(a,a)
#define sum_sq3(a) dot3(a,a)
#define cross3i(a,b) (((a[1])*(b[2])) - ((a[2])*(b[1])))
#define cross3j(a,b) (((a[2])*(b[0])) - ((a[0])*(b[2])))
#define cross3k(a,b) (((a[0])*(b[1])) - ((a[1])*(b[0])))

// static void printf_vec(const int nd, double x[nd])
// {
//     printf(" (");
//     int i;
//     for(i = 0; i < nd - 1; i++) {
//         printf("%g, ", x[i]);
//     }
//     printf("%g", x[i]);
//     printf(")\n");
// }

static double calc_elem_area(const int nd, double nodes[nd][nd])
{
    /* len = diff(nodes) .. point to point distances */
    // for(int i = 0; i < nd; i++) {
    //     printf("[%d] ", i); printf_vec(nd, nodes[i]);
    // }
    for(int i = 0; i < nd - 1; i++) {
        for(int j = 0; j < nd; j++) {
            nodes[i][j] -= nodes[i + 1][j]; /* subtract the row below */
        }
        // printf("[%d]-[%d] ", i, i + 1); printf_vec(nd, nodes[i]);
    }
    if(nd == 3) {
        /* cross product of lengths c = a x b*/
        // printf("%gi\n",nodes[0][1]*nodes[1][2] - nodes[0][2]*nodes[1][1]);
        // printf("%gj\n",nodes[0][2]*nodes[1][0] - nodes[0][0]*nodes[1][2]);
        // printf("%gk\n",nodes[0][0]*nodes[1][1] - nodes[0][1]*nodes[1][0]);
        const double tmp[3] = {
            cross3i(nodes[0], nodes[1]),
            cross3j(nodes[0], nodes[1]),
            cross3k(nodes[0], nodes[1]),
        };
        // printf("[0] x [1]: "); printf_vec(nd, tmp);
        assert(dot3(tmp, nodes[0]) == 0); /* confirm c . a = 0 */
        assert(dot3(tmp, nodes[1]) == 0); /* confirm c . b = 0 */
        nodes[0][0] = sum_sq3(tmp);
        // printf("∑[0]²: "); printf_vec(1, nodes[0]);
    }
    else {
        nodes[0][0] = sum_sq2(nodes[0]);
        // printf("∑[0]²: "); printf_vec(1, nodes[0]);
    }
    const double area = sqrt(nodes[0][0]) / (double)(nd - 1);
    return area;
}

static void copy_nodes_from_surfaceelems(mesh const * const m, const int se_idx, const int nd, double nodes[nd][nd])
{
    const int dim = m->dim;
    for(int i = 0; i < dim; i++) {
        const int nidx = m->surfaceelems[se_idx * dim + i] - 1; /* node index */
        for(int j = 0; j < dim; j++) {
            nodes[i][j] = m->nodes[nidx * dim + j];
        }
    }
}

static int max_bc(mesh const * const m)
{
    int max = -1;
    for(int i = 0; i < m->n_se; i++ ) {
        const int bci = m->bc[i];
        if (bci > max) {
            max = bci;
        }
    }
    return max;
}

int calc_elec_to_sys_map(model * m)
{
    assert(m != NULL);
    free(m->elec_to_sys);
    m->n_elec = 0;
    m->elec_to_sys = NULL;
    const int dim = m->fwd.dim;
    const int n_bc = max_bc(&(m->fwd)) + 1;
    if(n_bc < 1) {
        return SUCCESS;
    }
    int cnt[n_bc];
    int nodes[n_bc];
    int pem_nodes[n_bc];
    for(int i = 0; i < n_bc; i++) {
        cnt[i] = 0;
        nodes[i] = 0;
        pem_nodes[i] = 0;
    }
    for(int i = 0; i < m->fwd.n_se; i ++ ) {
        const int bci = m->fwd.bc[i];
        cnt[bci]++;
        for(int j = 0; j < dim; j++) {
            const int nidx = m->fwd.surfaceelems[i * dim + j];
            if(nidx > 0) {
                nodes[bci]++;
                pem_nodes[bci] = nidx - 1;
            }
        }
    }
    /* we assume the BC=0 is the generic boundary with no electrodes, so start at i=1 */
    int n_elec = 0;
    for(int i = 1; i < n_bc; i++) {
        if(cnt[i] > 0) {
            n_elec++;
        }
    }
    m->elec_to_sys = malloc(sizeof(int) * n_elec);
    if(m->elec_to_sys == NULL) {
        return FAILURE;
    }
    m->n_elec = n_elec;
    int n_cem = 0;
    int electrode = 0;
    for(int i = 1; i < n_bc; i++) {
        if((cnt[i] == 1) && (nodes[i] == 1)) { /* PEM */
            m->elec_to_sys[electrode++] = pem_nodes[i];
        }
        else if(cnt[i] > 0) { /* CEM */
            m->elec_to_sys[electrode++] = n_cem + m->fwd.n_nodes;
            n_cem++;
        }
    }
    assert(electrode == n_elec);
    return SUCCESS;
}

static int count_se_for_this_bc(mesh const * const m, int bc)
{
    int n = 0;
    for(int i = 0; i < m->n_se; i ++ ) {
        const int bci = m->bc[i];
        if (bci == bc) {
            n++;
        }
    }
    return n;
}

/* Construct a dense column vector for Neumann (current) stimulus 'amp' on
 * boundary m->bc = bc.
 * Returns number of nodes perturbed: 0 = failure. */
int calc_stim_neumann(mesh const * const m, double amp, int bc, double * b)
{
    assert(m != NULL);
    assert(b != NULL);
    const int dim = m->dim;
    const int n = count_se_for_this_bc(m, bc);
    int cnt_bc_local_nodes = n * dim;
    int list[n];
    double area[n];
    double total_area = 0;
    for(int i = 0, j = 0; i < m->n_se; i ++ ) {
        const int bci = m->bc[i];
        if (bci != bc) {
            continue;
        }
        double tmp[dim][dim];
        copy_nodes_from_surfaceelems(m, i, dim, tmp);
        area[j] = calc_elem_area(dim, tmp);
        total_area += area[j];
        // printf("bc#%d: area[%d]=%g\n", bc, j, area[j]);
        list[j] = i;
        j++;
    }
    /* set boundary condition */
    for(int j = 0; j < n; j++ ) {
        const int i = list[j];
        for(int k = 0; k < dim; k++) {
            const int nidx = m->surfaceelems[i * dim + k] - 1;
            // printf("#%d nidx=%d %g*%g/%g %d\n",j,nidx,amp, area[j],total_area,dim);
            b[ nidx ] += amp * (area[j] / total_area) / (double) dim;
        }
    }
    return cnt_bc_local_nodes;
}

int calc_sys_cem_n(const int nd)
{
    /* nn = 2D: 3x3 sym, 3D: 4x4 sym
     * n_{bc=i} = ∑(n_se_{bc=i})x(nn) for n_se_{bc=i} = number of surface elements with BC=i
     * total CEM system matrix elements is then
     *   n = ∑_{i=1..ne} n_{bc=i} */
    return (nd == 2) ? 6 : 10;
}

int calc_sys_cem(mesh const * const m, int bc, double zc, int cem_node, int Se_idx, int * ii, int * jj, double * Se)
{
    /* For each segment of the CEM add
     * 2D:    cem   --   n1    n2
     * A_2d = 6     0    -3    -3   * 1/6 * length / zc
     *        0     0     0     0
     *       -3     0     2     1
     *       -3     0     1     2
     * 3D:    cem   --   n1    n2    n3
     * A_3d = 12    0    -4    -4    -4 * 1/12 * area / zc
     *         0    0     0     0     0
     *        -4    0     2     1     1
     *        -4    0     1     2     1
     *        -4    0     1     1     2
     */
    assert(m != NULL);
    const int Se_idx_start = Se_idx;
    const int dim = m->dim;
    const int n = count_se_for_this_bc(m, bc);
    int list[n];
    double area[n];
    double total_area = 0;
    for(int i = 0, j = 0; i < m->n_se; i ++ ) {
        const int bci = m->bc[i];
        if (bci != bc) {
            continue;
        }
        double tmp[dim][dim];
        copy_nodes_from_surfaceelems(m, i, dim, tmp);
        area[j] = calc_elem_area(dim, tmp);
        total_area += area[j];
        // printf("bc#%d: area[%d]=%g\n", bc, j, area[j]);
        list[j] = i;
        j++;
    }
    for(int j = 0; j < n; j++ ) {
        const int i = list[j];
        const double scale = (dim == 2) ? 6.0 : 12.0;
        const double area_zc = area[j] / zc;
        Se[Se_idx] = area_zc;
        ii[Se_idx] = cem_node;
        jj[Se_idx] = cem_node;
        Se_idx++;
        for(int k = 0; k < dim; k++) {
            const int nidx = m->surfaceelems[i * dim + k] - 1;
            Se[Se_idx] = ((dim == 2) ? -3.0 : -4.0) / scale * area_zc;
            ii[Se_idx] = cem_node;
            jj[Se_idx] = nidx;
            Se_idx++;
            Se[Se_idx] = +2.0 / scale * area_zc;
            ii[Se_idx] = nidx;
            jj[Se_idx] = nidx;
            Se_idx++;
            for(int p = 0; p < k; p++) {
                const int pidx = m->surfaceelems[i * dim + p] - 1;
                Se[Se_idx] = +1.0 / scale * area_zc;
                ii[Se_idx] = pidx;
                jj[Se_idx] = nidx;
                Se_idx++;
            }
        }
    }
    const int nnz_added = calc_sys_cem_n(dim) * n;
    assert(Se_idx - Se_idx_start == nnz_added);
    return nnz_added;
}

enum se_type {UNKNOWN, NONE, PEM, CEM};
static int check_bc(mesh const * const msh)
{
    assert(msh != NULL);
    const int dim = msh->dim;
    const int n_bc = max_bc(msh) + 1;
    enum se_type type[n_bc];
    for(int i = 0; i < n_bc; i++) {
        type[i] = UNKNOWN;
    }
    for(int i = 0; i < msh->n_se; i ++ ) {
        const int bci = msh->bc[i];
        /* a PEM is recorded in surfaceelems with a single node
         * as the first entry and the remaining entries as zero;
         * so we test the second entry to see if its zero */
        int cnt = 0;
        for(int j = 0; j < msh->dim; j++) {
            cnt += (msh->surfaceelems[i * dim + j] > 0) ? 1 : 0;
        }
        enum se_type newtype;
        if(cnt == msh->dim) {
            newtype = CEM;
        }
        else if(cnt == 1) {
            newtype = PEM;
        }
        else {
            printf("error: bad number of nodes on boundary %d\n", bci);
            return FAILURE;    /* bad node#s */
        }
        if(type[bci] == UNKNOWN) {
            type[bci] = newtype;
        }
        else if(type[bci] != newtype) {
            printf("error: mixed PEM and CEM on boundary %d\n", bci);
            return FAILURE;    /* mixed PEM and CEM syntax in same  BC */
        }
        else if(newtype == PEM) {
            printf("error: multiple PEM surface elements have the same boundary#%d, see surface element %d\n", bci, i + 1);
            return FAILURE;    /* duplicate PEM for same BC */
        }
        if((newtype == PEM) && (msh->surfaceelems[i * dim + 0] <= 0)) {
            printf("error: for PEM on boundary %d, must have node# on first entry\n", bci);
            return FAILURE;    /* must have se[0] be the node for a PEM */
        }
    }
    /* expect all BC to be used */
    for(int bc = 0; bc < n_bc; bc++) {
        if(type[bc] == UNKNOWN) {
            printf("error: unused boundary %d\n", bc);
            return FAILURE;
        }
    }
    return SUCCESS;
}

int check_model(model const * const mdl)
{
    return check_bc(&(mdl->fwd));
}

static void printf_mesh(mesh const * const msh, int verbosity)
{
    assert(msh != NULL);
    if(msh->n_elems == 0) {
        printf(" (none)\n");
        return;
    }
    printf(" %dD, %d nodes, %d elements, %d surfaces, %d entry parameter map\n",
           msh->dim, msh->n_nodes, msh->n_elems, msh->n_se, msh->n_pmap);
    if(verbosity <= 1) {
        return;
    }
    if(msh->n_elems > 0) {
        printf("    elements\n");
    }
    const int nd = msh->dim;
    const int nd1 = nd + 1;
    for(int i = 0; i < msh->n_elems; i++) {
        printf("     ");
        for(int j = 0; j < nd1; j++) {
            printf(" %-5d", msh->elems[i * nd1 + j]);
        }
        printf("\n");
    }
    if(msh->n_nodes > 0) {
        printf("    nodes\n");
    }
    for(int i = 0; i < msh->n_nodes; i++) {
        printf("     ");
        for(int j = 0; j < nd; j++) {
            printf(" %-+7.3f", msh->nodes[i * nd + j]);
        }
        printf("\n");
    }
    if(msh->n_se > 0) {
        printf("    surface elements\n");
    }
    for(int i = 0; i < msh->n_se; i++) {
        printf("     ");
        for(int j = 0; j < nd; j++) {
            printf(" %-5d", msh->surfaceelems[i * nd + j]);
        }
        if(msh->bc[i] != 0) {
            printf(" (bc#%d)", msh->bc[i]);
        }
        printf("\n");
    }
    if(msh->n_pmap > 0) {
        printf("    parameter map\n");
    }
    for(int i = 0; i < msh->n_pmap; i++) {
        printf("     (%5d,%5d)  = %0.3f\n",
               msh->pmap_elem[i], msh->pmap_param[i],
               msh->pmap_frac[i]);
    }
}

#define plural(a) ((a!=1)?"s":"")
void printf_model(model const * const mdl, int verbosity)
{
    printf("model: %d electrodes, %d measurements,", mdl->n_elec, mdl->n_stimmeas);
    printf(" λ = %g\n", mdl->hp);
    printf("  %d parameter frame%s", mdl->n_params[1], plural(mdl->n_params[1]));
    if(mdl->n_params[1] > 0) {
        printf(" (%d parameter%s/frame)\n", mdl->n_params[0], plural(mdl->n_params[0]));
    }
    else {
        printf("\n");
    }
    printf("  %d measurement frame%s", mdl->n_data[1], plural(mdl->n_data[1]));
    if(mdl->n_data[1] > 0) {
        printf(" (%d measurement%s/frame)\n", mdl->n_data[0], plural(mdl->n_data[0]));
    }
    else {
        printf("\n");
    }
    if((verbosity > 1) && (mdl->n_stimmeas > 0)) {
        printf("  %d stimulus and measurement pairs: A-B M-N (gain)\n", mdl->n_stimmeas);
        for(int i = 0; i < mdl->n_stimmeas; i++) {
            printf("    %d-%d %d-%d (%5.2e)\n",
                   mdl->stimmeas[i * 4 + 0],
                   mdl->stimmeas[i * 4 + 1],
                   mdl->stimmeas[i * 4 + 2],
                   mdl->stimmeas[i * 4 + 3],
                   mdl->measgain[i]);
        }
    }
    if((verbosity > 1) && (mdl->n_zc > 0)) {
        printf("  %d contact impedances\n", mdl->n_zc);
        if(mdl->zcbc == NULL) {
            for(int i = 0; i < mdl->n_zc; i++) {
                printf("    %5.2e (electrode#%d)\n", mdl->zc[i], i + 1);
            }
        }
        else {
            for(int i = 0; i < mdl->n_zc; i++) {
                printf("    %5.2e (bc#%d)\n", mdl->zc[i], mdl->zcbc[i]);
            }
        }
    }
    printf("  forward model:");
    printf_mesh(&(mdl->fwd), verbosity);
    printf("  reconstruction model:");
    printf_mesh(&(mdl->rec), verbosity);
    if((verbosity > 1) && (mdl->n_params[0] > 0)) {
        const int rows = mdl->n_params[0];
        printf("  parameters\n");
        for(int i = 0; i < mdl->n_params[0]; i++) {
            printf("   ");
            for(int j = 0; j < mdl->n_params[1]; j++) {
                printf(" %18.8e", mdl->params[i + rows * j]);
            }
            printf("\n");
        }
    }
    if((verbosity > 1) && (mdl->n_data[0] > 0)) {
        const int rows = mdl->n_data[0];
        printf("  data\n");
        for(int i = 0; i < mdl->n_data[0]; i++) {
            printf("   ");
            for(int j = 0; j < mdl->n_data[1]; j++) {
                printf(" %18.8e", mdl->data[i + rows * j]);
            }
            printf("\n");
        }
    }
}

int calc_sys_size(model const * const m)
{
    const int dim = m->fwd.dim;
    const int n_bc = max_bc(&(m->fwd)) + 1;
    enum se_type type[n_bc];
    for(int i = 0; i < n_bc; i++) {
        type[i] = UNKNOWN;
    }
    type[0] = NONE; /* skip bc=0, assuming its the exterior */
    for(int i = 0; i < m->fwd.n_se; i ++ ) {
        const int bci = m->fwd.bc[i];
        if (type[bci] != UNKNOWN) {
            continue;
        }
        type[bci] = (m->fwd.surfaceelems[i * dim + 1] == 0) ? PEM : CEM;
    }
    int n_cem = 0;
    for(int bc = 0; bc < n_bc; bc++) {
        n_cem += (type[bc] == CEM) ? 1 : 0;
    }
    return (m->fwd.n_nodes + n_cem);
}

size_t calc_sys_nnz(model const * const m)
{
    const int dim = m->fwd.dim;
    const int n_bc = max_bc(&(m->fwd)) + 1;
    enum se_type type[n_bc];
    for(int i = 0; i < n_bc; i++) {
        type[i] = UNKNOWN;
    }
    int cnt[n_bc];
    for(int i = 0; i < n_bc; i++) {
        cnt[i] = 0;
    }
    for(int i = 0; i < m->fwd.n_se; i ++ ) {
        const int bci = m->fwd.bc[i];
        type[bci] = (m->fwd.surfaceelems[i * dim + 1] == 0) ? PEM : CEM;
        cnt[bci]++;
    }
    type[0] = NONE; /* force skipping bc=0, assuming its the exterior */
    /* count surfaceelems that are CEM */
    int se_cem = 0;
    for(int bc = 0; bc < n_bc; bc++) {
        se_cem += (type[bc] == CEM) ? cnt[bc] : 0;
    }
    const int ne = m->fwd.n_elems;
    const int se_n = calc_sys_elem_n(dim); /* sparse matrix entries per mesh element */
    const int nnz_per_se = calc_sys_cem_n(dim);
    return (se_n * ne) + (se_cem * nnz_per_se);
}

max_stars_repo_name: boyle/zedhat

max_issues_repo_name: boyle/zedhat

max_forks_repo_name: boyle/zedhat