text:
//
// This file is subject to the terms and conditions defined in
// file 'LICENSE.txt', which is part of this source code package.
//
// Principal developer : Felipe Bordeu (Felipe.Bordeu@ec-nantes.fr)
//

#include <iostream>

#ifdef MATLAB_MEX_FILE
#include "mex.h"
#ifdef __cplusplus
    extern "C" bool utIsInterruptPending();
    extern "C" void utSetInterruptPending(bool);
#else
    extern bool utIsInterruptPending();
    extern void utSetInterruptPending(bool);
#endif
#endif /* MATLAB_MEX_FILE */

//if in matlab we use the matlab blas library
#ifdef MATLAB_MEX_FILE
    #include <blas.h>
    inline double cblas_ddot( ptrdiff_t& a, double*& b,const int& c, double* d,const int& e){
        ptrdiff_t cc = c;
        ptrdiff_t ee = e;
        return ddot(&a,b,&cc,d,&ee);
   }
   enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102 };
   enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113,
                         AtlasConj=114};

template<typename T5,typename T6,typename T7,typename T8,typename T9,typename T12>
void cblas_dgemv(const CBLAS_ORDER& a, const CBLAS_TRANSPOSE& b, const ptrdiff_t & c, const ptrdiff_t& d,const T5& e,T6 f,T7& g,T8 h,const T9& i,const double& j,double* k,const T12& l){
    char bb= 'n';
    if (b ==CblasTrans) bb = 't';
    ptrdiff_t cc = c;
    ptrdiff_t dd = d;
    double ee = e;
    ptrdiff_t ii = i;
    double jj = j;
    ptrdiff_t ll  =  l;
    dgemv( &bb,  &cc,    &dd,    &ee,    f,    &g,    h,    &ii,    &jj,    k,   &ll);
};

template<typename T5,typename T6,typename T7, typename T8,typename T9,typename T12>
void cblas_sgemv(const CBLAS_ORDER& a, const CBLAS_TRANSPOSE& b, const ptrdiff_t& c, const ptrdiff_t& d, const T5& e,T6 f,T7& g,T8 h,const  T9& i,const float& j,float* k,const T12& l){
    char bb= 'n';
    if (b ==CblasTrans) bb = 't';
    ptrdiff_t cc = c;
    ptrdiff_t dd = d;
    float ee = e;
    ptrdiff_t ii = i;
    float jj = j;
    ptrdiff_t ll  =  l;

    sgemv( &bb,  &cc,    &dd,    &ee,    f,    &g,    h,    &ii,    &jj,    k,   &ll);
};

 void cblas_dgemm(const CBLAS_ORDER& a, const CBLAS_TRANSPOSE& b, const CBLAS_TRANSPOSE& c, ptrdiff_t d, ptrdiff_t e, ptrdiff_t f ,double g, double* h, ptrdiff_t i,double* j , ptrdiff_t k,double l,double* m,ptrdiff_t n){
    char bb= 'n';
    if (b ==CblasTrans) bb = 't';
    char cc= 'n';
    if (c ==CblasTrans) cc = 't';

    dgemm( &bb, &cc, &d, &e, &f, &g, h, &i, j, &k, &l,m, &n);
}

template <typename pdiff>
void  cblas_dscal (pdiff& a,double& b, double* c,const pdiff& d){
    pdiff dd = d;
    dscal(&a,&b,c,&dd);
}

template <typename pdiff>
void  cblas_sscal (pdiff& a,float& b, float* c,const pdiff& d){
    pdiff dd = d;
    sscal(&a,&b,c,&dd);
}

template <typename pdiff>
inline float cblas_sdot( pdiff& a, float*& b,const pdiff& c, float* d,const pdiff& e){
    pdiff cc = c;
    pdiff ee = e;
    return sdot(&a,b,&cc,d,&ee);
}
#else
// if not matlab we use the 'system' blas

#ifdef __APPLE__
   #include   <Accelerate/Accelerate.h>
#else
#if defined(_WIN32)
   /// for now no compilation can be done in windows without matlab
#else
   /// linux
   #include <cblas.h>
   #include <vector>
#endif
#endif
#endif /* MATLAB_MEX_FILE */

template <typename ptrdiff, typename T>
T cblas_Xdot( ptrdiff&, T*&,const ptrdiff&, T*,const ptrdiff& );

template <typename ptrdiff>
double cblas_Xdot( ptrdiff& a, double*& b,const ptrdiff& c, double* d,const ptrdiff& e){
return cblas_ddot(a,b,c,d,e);
}

template <typename ptrdiff>
float cblas_Xdot( ptrdiff& a, float*& b,const ptrdiff& c, float* d,const ptrdiff& e){
return cblas_sdot(a,b,c,d,e);
}

//template<typename T>
//void cblas_Xgemv(const CBLAS_ORDER&  , const CBLAS_TRANSPOSE&  ,int&   ,int&  ,const T&       ,     T*  , int&  ,T*      , const int&   ,const T&      , T*       ,const int&);
//

template <typename ptrdiff>
void cblas_Xgemv(const CBLAS_ORDER& a, const CBLAS_TRANSPOSE& b,const int& c,const int& d,const double& e,double* f, ptrdiff& g,double* h,const int& i ,const double& j,double* k,const int& l){
            cblas_dgemv(a,b,c,d,e,f,g,h,i,j,k,l);
};

template <typename ptrdiff>
void cblas_Xgemv(const CBLAS_ORDER& a,const CBLAS_TRANSPOSE& b,const int& c,const int& d,const float& e,float* f,ptrdiff& g,float* h,const  int& i,const float& j,float* k,const int& l){
            cblas_sgemv(a,b,c,d,e,f,g,h,i,j,k,l);
};

//template<typename T1,typename T2, typename T3,typename T4>
//void cblas_Xscal(T1&,T2&,T3,const T4&);

//template<typename T1>
void  cblas_Xscal(ptrdiff_t& a,double& b,double* c,const ptrdiff_t& d){
        cblas_dscal(a,b,c ,d);
}

//template<typename T1>
void  cblas_Xscal(ptrdiff_t& a,float& b, float* c,const ptrdiff_t& d){
        cblas_sscal(a,b,c,d);
}

template <typename T>
void update_C(MatLabDataMatrix <T >& C,T*  r, MatLabDataMatrix <T >& MM,ptrdiff_t& N, const ptrdiff_t& dim);

template<>
void update_C(MatLabDataMatrix <double >& C,double*  r, MatLabDataMatrix <double >& MM,ptrdiff_t& NN, const ptrdiff_t& dim) {
    int N = NN;
    int M = MM.dsizes[0];
    int INCY = C.dsizes[0];
    cblas_dgemv(CblasColMajor,CblasTrans, M, N, 1.0, &MM.Data[0], MM.dsizes[0],      r, 1,      0, &C.Data[dim], C.dsizes[0] );
    #ifdef USE_MPI
       MPI_Allreduce(MPI_IN_PLACE, &C.Data[dim], NN, everydim, myOp, MPI_COMM_WORLD);
    #endif /* USE_MPI */
};

template<>
void update_C(MatLabDataMatrix <float >& C,float*  r, MatLabDataMatrix <float >& MM,ptrdiff_t& NN, const ptrdiff_t& dim) {
    int N = NN;
    int M = MM.dsizes[0];
    int INCY = C.dsizes[0];
    cblas_sgemv(CblasColMajor,CblasTrans, M, N, (float)1.0, &MM.Data[0], MM.dsizes[0],      r, 1,      0, &C.Data[dim], C.dsizes[0] );
    #ifdef USE_MPI
       MPI_Allreduce(MPI_IN_PLACE, &C.Data[dim], NN, everydim, myOp, MPI_COMM_WORLD);
    #endif /* USE_MPI */
};

template<typename T>
void update_D(MatLabDataMatrix <T >& D, MatLabDataMatrix <T >& FF, T* r){
  int cpt = -1;
  for(int j =0; j < FF.dsizes[1]; ++j){
    for(int i =0; i < FF.dsizes[0]; ++i){
        D.Data[++cpt] = FF(i,j)*r[i];
    }
  }
};

//template <typename T>
//void update_CWithFF2(MatLabDataMatrix <T >& C,T*  r, MatLabDataMatrix <T >& MM);
//,unsigned& N, const unsigned& dim, MatLabDataMatrix <T >& DD, std::vector<T> &workingplace);

//template<>
void update_CU(std::vector<double >& C,MatLabDataMatrix <double >&  sol, MatLabDataMatrix <double >& D) {

    if( sol.dsizes[1] == 0) return;

    int Sr,Sc,Dr,Dc;
    Sr = sol.dsizes[0];
    Sc = sol.dsizes[1];

    Dr = D.dsizes[0];
    Dc = D.dsizes[1];

    ///C(Cr,Cc) = sol(Sr,Sc)'*D(Dr,Dc);

    //printf (" Initializing data for matrix multiplication C=sol'*D for matrix \n"
    //        " sol(%ix%i) and matrix D(%ix%i)\n\n", Sr, Sc, Dr, Dc);

    //std::cout << "sol.size() : " << sol.fullsize << " =? " <<Sr*Sc << std::endl;
    //std::cout << "D.size() : " << D.fullsize << " =? " <<Dr*Dc << std::endl;

   //std::cout << "--------------------------------------" << std::endl;

    cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans,
                Sc,Dc, Dr, 1.0, sol.Data, sol.LD(), D.Data , D.LD(), 0, &C[0], sol.dsizes[1]);

    //std::cout << " ------------------------- DONE ------------------------" << std::endl;
    #ifdef USE_MPI
       MPI_Allreduce(MPI_IN_PLACE, &C[0], C.size(), everydim, myOp, MPI_COMM_WORLD);
    #endif /* USE_MPI */
};
//
void update_CU(std::vector<float >& C,MatLabDataMatrix <float >&  sol, MatLabDataMatrix <float >& D) {
//#warning to correct for the double implementation
    //cblas_sgemm(CblasColMajor, CblasTrans, CblasNoTrans, sol.dsizes[1], D.dsizes[1], sol.dsizes[0], 1.0, sol.Data, sol.dsizes[0], D.Data , D.dsizes[0], 0, &C[0],sol.dsizes[0]);
    #ifdef USE_MPI
       MPI_Allreduce(MPI_IN_PLACE, &C[0], C.size(), everydim, myOp, MPI_COMM_WORLD);
    #endif /* USE_MPI */
};

#define PrintM(X) std::cout << #X " : "  << std::endl << "  "; \
                  for(int ii =0; ii < X.dsizes[0] ; ++ii){\
                    for(int jj =0; jj < X.dsizes[1] ; ++jj){\
                      std::cout << X(ii,jj) << " ";\
                    }\
                    std::cout << std::endl << "  "; \
                  }\
                  std::cout << std::endl;

#define PrintV(X) std::cout << #X " : "  << std::endl << "  "; \
                  for(int ii=0; ii < X.size(); ++ii){\
                    std::cout << X[ii] << " "; \
                  } \
                  std::cout << std::endl;

// template<>
// void update_CWithFF2(MatLabDataMatrix <float >& C,float*  r, MatLabDataMatrix <float >& MM,unsigned& NN, const unsigned& dim, MatLabDataMatrix <float >& DD, std::vector<float> &workingplace) {    int N = NN;
//     ///TODO function to be coded
//     int M = MM.dsizes[0];
//     int LDA = MM.dsizes[0];
//     int INCY = C.dsizes[0];
//     cblas_sgemv(CblasColMajor,CblasTrans, M, N, 1.0, &MM.Data[0], LDA,      r, 1,      0, &C.Data[dim], C.dsizes[0] );
//     #ifdef USE_MPI
//        MPI_Allreduce(MPI_IN_PLACE, &C.Data[dim], NN, everydim, myOp, MPI_COMM_WORLD);
//     #endif /* USE_MPI */
//
// };

max_stars_repo_name: fbordeu/PGDforMatlab

max_issues_repo_name: fbordeu/PGDforMatlab

max_forks_repo_name: fbordeu/PGDforMatlab


text:
static char help[] =
"Load a matrix  A  and right-hand-side  b  from binary files (PETSc format).\n"
"Then solve the system  A x = b  using KSPSolve().\n"
"Example.  First save a system from tri.c:\n"
"  ./tri -ksp_view_mat binary:A.dat -ksp_view_rhs binary:b.dat\n"
"then load it and solve:\n"
"  ./loadsolve -fA A.dat -fb b.dat\n"
"To time the solution read the third printed number:\n"
"  ./loadsolve -fA A.dat -fb b.dat -log_view |grep KSPSolve\n"
"(This is a simpler code than src/ksp/ksp/examples/tutorials/ex10.c.)\n";

/*
small system example where RHS missing:
./tri -ksp_view_mat binary:A.dat
./loadsolve -fA A.dat -ksp_view_mat -ksp_view_rhs

large tridiagonal system (m=10^7) example:
./tri -tri_m 10000000 -ksp_view_mat binary:A.dat -ksp_view_rhs binary:b.dat
./loadsolve -fA A.dat -fb b.dat -log_view |grep KSPSolve
*/

#include <petsc.h>

int main(int argc,char **args) {
  PetscErrorCode ierr;
  Vec         x, b;
  Mat         A;
  KSP         ksp;
  int         m, n, mb;
  PetscBool   flg,
              verbose = PETSC_FALSE;
  char        nameA[PETSC_MAX_PATH_LEN] = "",
              nameb[PETSC_MAX_PATH_LEN] = "";
  PetscViewer fileA, fileb;

  PetscInitialize(&argc,&args,NULL,help);

  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","options for loadsolve",""); CHKERRQ(ierr);
  ierr = PetscOptionsString("-fA","input file containing matrix A",
                            "loadsolve.c",nameA,nameA,PETSC_MAX_PATH_LEN,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsString("-fb","input file containing vector b",
                            "loadsolve.c",nameb,nameb,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsBool("-verbose","say what is going on",
                          "loadsolve.c",verbose,&verbose,NULL); CHKERRQ(ierr);
  ierr = PetscOptionsEnd(); CHKERRQ(ierr);
  if (strlen(nameA) == 0) {
      SETERRQ(PETSC_COMM_WORLD,1,
              "no input matrix provided ... ending  (usage: loadsolve -fA A.dat)\n");
  }

  if (verbose) {
      ierr = PetscPrintf(PETSC_COMM_WORLD,
         "reading matrix from %s ...\n",nameA); CHKERRQ(ierr);
  }
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,nameA,FILE_MODE_READ,&fileA);CHKERRQ(ierr);
  ierr = MatLoad(A,fileA);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&fileA);CHKERRQ(ierr);
  ierr = MatGetSize(A,&m,&n); CHKERRQ(ierr);
  if (verbose) {
      ierr = PetscPrintf(PETSC_COMM_WORLD,
         "matrix has size m x n = %d x %d ...\n",m,n); CHKERRQ(ierr);
  }
  if (m != n) {
      SETERRQ(PETSC_COMM_WORLD,2,"only works for square matrices\n");
  }

  ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
  ierr = VecSetFromOptions(b);CHKERRQ(ierr);
  if (flg) {
      ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,nameb,FILE_MODE_READ,&fileb);CHKERRQ(ierr);
      if (verbose) {
          ierr = PetscPrintf(PETSC_COMM_WORLD,
             "reading vector from %s ...\n",nameb); CHKERRQ(ierr);
      }
      ierr = VecLoad(b,fileb);CHKERRQ(ierr);
      ierr = PetscViewerDestroy(&fileb);CHKERRQ(ierr);
      ierr = VecGetSize(b,&mb); CHKERRQ(ierr);
      if (mb != m) {
          SETERRQ(PETSC_COMM_WORLD,3,"size of matrix and vector do not match\n");
      }
  } else {
      if (verbose) {
          ierr = PetscPrintf(PETSC_COMM_WORLD,
             "right-hand-side vector b not provided ... using zero vector of length %d\n",m); CHKERRQ(ierr);
      }
      ierr = VecSetSizes(b,PETSC_DECIDE,m); CHKERRQ(ierr);
      ierr = VecSet(b,0.0); CHKERRQ(ierr);
  }

  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp); CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp,A,A); CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);

  ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
  ierr = KSPSolve(ksp,b,x); CHKERRQ(ierr);

  KSPDestroy(&ksp);  MatDestroy(&A);
  VecDestroy(&x);  VecDestroy(&b);
  return PetscFinalize();
}

ext: c


text:
#pragma once

#include <cstring>
#include <string>
#include <tuple>
#include <utility>

#include <gsl/assert>

namespace xmol::utils {

constexpr static int powers_of_10[] = {1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000};

template <int LEN> inline std::pair<bool, int> parse_uint_fixed_length(const std::string& line, int pos) noexcept {
  bool success = GSL_LIKELY(line.size() >= pos + LEN);
  if (GSL_UNLIKELY(!success)) {
    return {false, 0};
  }
  int number = 0;
#pragma unroll
  for (int i = 0; i < LEN; i++) {
    unsigned digit = static_cast<unsigned>(line[pos + i]) - '0';
    success &= digit <= 9;
    number += digit * powers_of_10[LEN - i - 1];
  }
  return std::make_pair(success, number);
};

inline std::pair<bool, int> parse_uint(const std::string& line, int pos, int LEN) noexcept {
  switch (LEN) {
  case (1):
    return parse_uint_fixed_length<1>(line, pos);
  case (2):
    return parse_uint_fixed_length<2>(line, pos);
  case (3):
    return parse_uint_fixed_length<3>(line, pos);
  case (4):
    return parse_uint_fixed_length<4>(line, pos);
  case (5):
    return parse_uint_fixed_length<5>(line, pos);
  case (6):
    return parse_uint_fixed_length<6>(line, pos);
  case (7):
    return parse_uint_fixed_length<7>(line, pos);
  case (8):
    return parse_uint_fixed_length<8>(line, pos);
  case (9):
    return parse_uint_fixed_length<9>(line, pos);
  default:
    return {false, 0};
  }
};

enum class SpaceStrip { NONE, LEFT, RIGHT, LEFT_AND_RIGHT };

template <SpaceStrip stripping = SpaceStrip::NONE>
std::pair<bool, int> parse_int(const std::string& line, int pos, int LEN) noexcept {
  // short-circuit for invalid
  if (GSL_UNLIKELY(line.size() < pos + LEN)) {
    return {false, 0};
  }

  if (stripping == SpaceStrip::RIGHT || stripping == SpaceStrip::LEFT_AND_RIGHT) {
    while (line[pos + LEN - 1] == ' ' && LEN > 0) {
      --LEN;
    }
  }

  if (stripping == SpaceStrip::LEFT || stripping == SpaceStrip::LEFT_AND_RIGHT) {
    while (line[pos] == ' ' && LEN > 0) {
      ++pos;
      --LEN;
    }
  }

  if (line[pos] == '-') {
    bool success;
    int value;
    std::tie(success, value) = parse_uint(line, pos + 1, LEN - 1);
    return {success, -value};
  } else {
    return parse_uint(line, pos, LEN);
  }
};

template <int WIDTH, int PRECISION, SpaceStrip STRIP> struct parse_fixed_precision_fn {
  static_assert(PRECISION >= 0);
  static_assert(WIDTH > 0);
  static_assert(PRECISION == 0 || WIDTH >= PRECISION + 2);

  inline std::pair<bool, double> operator()(const std::string& line, int pos) const noexcept {

    if (GSL_UNLIKELY(line.size() < pos + WIDTH || (PRECISION > 0 && line[pos + WIDTH - PRECISION - 1] != '.'))) {
      return {false, 0};
    }

    int whole = WIDTH - PRECISION;
    if (PRECISION > 0) {
      --whole;
    }
    int precision = PRECISION;

    if (STRIP == SpaceStrip::LEFT_AND_RIGHT || STRIP == SpaceStrip::LEFT) {
      while (line[pos] == ' ' && whole > 0) {
        --whole;
        pos++;
      }
    }

    if (STRIP == SpaceStrip::LEFT_AND_RIGHT || STRIP == SpaceStrip::RIGHT) {
      while (line[pos + whole + precision] == ' ' && precision > 0) {
        --precision;
      }
    }

    int sign = 1;

    if (line[pos] == '-') {
      sign = -1;
      ++pos;
      --whole;
    }
    bool success;
    int whole_part;
    std::tie(success, whole_part) = parse_uint(line, pos, whole);

    if (GSL_UNLIKELY(!success)) {
      return {false, 0};
    }
    if (PRECISION == 0) {
      return {true, sign * whole_part};
    }
    bool success2;
    int fraction_part;
    std::tie(success2, fraction_part) = parse_uint(line, pos + whole + 1, precision);
    if (GSL_UNLIKELY(!success2)) {
      return {false, 0};
    }

    return {true, sign * (whole_part + double(fraction_part) / powers_of_10[precision])};
  };
};

inline namespace functional_objects {
struct parse_fixed_precision_fn___ {
  template <int WIDTH, int PRECISION, SpaceStrip STRIP = SpaceStrip::LEFT_AND_RIGHT>
  inline std::pair<bool, double> parse(const std::string& line, int pos) const {
    return parse_fixed_precision_fn<WIDTH,PRECISION,STRIP>{}(line,pos);
  }
};
constexpr parse_fixed_precision_fn___ parse_fixed_precision{};
}

std::pair<bool, double> parse_fixed_precision_rt(const std::string& line, int pos, int width) noexcept;
}

ext: h


text:
/* cdf/gauss.c
 *
 * Copyright (C) 2002, 2004 Jason H. Stover.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 */

/*
 * Computes the cumulative distribution function for the Gaussian
 * distribution using a rational function approximation.  The
 * computation is for the standard Normal distribution, i.e., mean 0
 * and standard deviation 1. If you want to compute Pr(X < t) for a
 * Gaussian random variable X with non-zero mean m and standard
 * deviation sd not equal to 1, find gsl_cdf_ugaussian ((t-m)/sd).
 * This approximation is accurate to at least double precision. The
 * accuracy was verified with a pari-gp script.  The largest error
 * found was about 1.4E-20. The coefficients were derived by Cody.
 *
 * References:
 *
 * W.J. Cody. "Rational Chebyshev Approximations for the Error
 * Function," Mathematics of Computation, v23 n107 1969, 631-637.
 *
 * W. Fraser, J.F Hart. "On the Computation of Rational Approximations
 * to Continuous Functions," Communications of the ACM, v5 1962.
 *
 * W.J. Kennedy Jr., J.E. Gentle. "Statistical Computing." Marcel Dekker. 1980.
 *
 *
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cdf.h>

#ifndef M_1_SQRT2PI
#define M_1_SQRT2PI (M_2_SQRTPI * M_SQRT1_2 / 2.0)
#endif

#define SQRT32 (4.0 * M_SQRT2)

/*
 * IEEE double precision dependent constants.
 *
 * GAUSS_EPSILON: Smallest positive value such that
 *                gsl_cdf_gaussian(x) > 0.5.
 * GAUSS_XUPPER: Largest value x such that gsl_cdf_gaussian(x) < 1.0.
 * GAUSS_XLOWER: Smallest value x such that gsl_cdf_gaussian(x) > 0.0.
 */

#define GAUSS_EPSILON  (GSL_DBL_EPSILON / 2)
#define GAUSS_XUPPER (8.572)
#define GAUSS_XLOWER (-37.519)

#define GAUSS_SCALE (16.0)

static double
get_del (double x, double rational)
{
  double xsq = 0.0;
  double del = 0.0;
  double result = 0.0;

  xsq = floor (x * GAUSS_SCALE) / GAUSS_SCALE;
  del = (x - xsq) * (x + xsq);
  del *= 0.5;

  result = exp (-0.5 * xsq * xsq) * exp (-1.0 * del) * rational;

  return result;
}

/*
 * Normal cdf for fabs(x) < 0.66291
 */
static double
gauss_small (const double x)
{
  unsigned int i;
  double result = 0.0;
  double xsq;
  double xnum;
  double xden;

  const double a[5] = {
    2.2352520354606839287,
    161.02823106855587881,
    1067.6894854603709582,
    18154.981253343561249,
    0.065682337918207449113
  };
  const double b[4] = {
    47.20258190468824187,
    976.09855173777669322,
    10260.932208618978205,
    45507.789335026729956
  };

  xsq = x * x;
  xnum = a[4] * xsq;
  xden = xsq;

  for (i = 0; i < 3; i++)
    {
      xnum = (xnum + a[i]) * xsq;
      xden = (xden + b[i]) * xsq;
    }

  result = x * (xnum + a[3]) / (xden + b[3]);

  return result;
}

/*
 * Normal cdf for 0.66291 < fabs(x) < sqrt(32).
 */
static double
gauss_medium (const double x)
{
  unsigned int i;
  double temp = 0.0;
  double result = 0.0;
  double xnum;
  double xden;
  double absx;

  const double c[9] = {
    0.39894151208813466764,
    8.8831497943883759412,
    93.506656132177855979,
    597.27027639480026226,
    2494.5375852903726711,
    6848.1904505362823326,
    11602.651437647350124,
    9842.7148383839780218,
    1.0765576773720192317e-8
  };
  const double d[8] = {
    22.266688044328115691,
    235.38790178262499861,
    1519.377599407554805,
    6485.558298266760755,
    18615.571640885098091,
    34900.952721145977266,
    38912.003286093271411,
    19685.429676859990727
  };

  absx = fabs (x);

  xnum = c[8] * absx;
  xden = absx;

  for (i = 0; i < 7; i++)
    {
      xnum = (xnum + c[i]) * absx;
      xden = (xden + d[i]) * absx;
    }

  temp = (xnum + c[7]) / (xden + d[7]);

  result = get_del (x, temp);

  return result;
}

/*
 * Normal cdf for
 * {sqrt(32) < x < GAUSS_XUPPER} union { GAUSS_XLOWER < x < -sqrt(32) }.
 */
static double
gauss_large (const double x)
{
  int i;
  double result;
  double xsq;
  double temp;
  double xnum;
  double xden;
  double absx;

  const double p[6] = {
    0.21589853405795699,
    0.1274011611602473639,
    0.022235277870649807,
    0.001421619193227893466,
    2.9112874951168792e-5,
    0.02307344176494017303
  };
  const double q[5] = {
    1.28426009614491121,
    0.468238212480865118,
    0.0659881378689285515,
    0.00378239633202758244,
    7.29751555083966205e-5
  };

  absx = fabs (x);
  xsq = 1.0 / (x * x);
  xnum = p[5] * xsq;
  xden = xsq;

  for (i = 0; i < 4; i++)
    {
      xnum = (xnum + p[i]) * xsq;
      xden = (xden + q[i]) * xsq;
    }

  temp = xsq * (xnum + p[4]) / (xden + q[4]);
  temp = (M_1_SQRT2PI - temp) / absx;

  result = get_del (x, temp);

  return result;
}

double
gsl_cdf_ugaussian_P (const double x)
{
  double result;
  double absx = fabs (x);

  if (absx < GAUSS_EPSILON)
    {
      result = 0.5;
      return result;
    }
  else if (absx < 0.66291)
    {
      result = 0.5 + gauss_small (x);
      return result;
    }
  else if (absx < SQRT32)
    {
      result = gauss_medium (x);

      if (x > 0.0)
        {
          result = 1.0 - result;
        }

      return result;
    }
  else if (x > GAUSS_XUPPER)
    {
      result = 1.0;
      return result;
    }
  else if (x < GAUSS_XLOWER)
    {
      result = 0.0;
      return result;
    }
  else
    {
      result = gauss_large (x);

      if (x > 0.0)
        {
          result = 1.0 - result;
        }
    }

  return result;
}

double
gsl_cdf_ugaussian_Q (const double x)
{
  double result;
  double absx = fabs (x);

  if (absx < GAUSS_EPSILON)
    {
      result = 0.5;
      return result;
    }
  else if (absx < 0.66291)
    {
      result = gauss_small (x);

      if (x < 0.0)
        {
          result = fabs (result) + 0.5;
        }
      else
        {
          result = 0.5 - result;
        }

      return result;
    }
  else if (absx < SQRT32)
    {
      result = gauss_medium (x);

      if (x < 0.0)
        {
          result = 1.0 - result;
        }

      return result;
    }
  else if (x > -(GAUSS_XLOWER))
    {
      result = 0.0;
      return result;
    }
  else if (x < -(GAUSS_XUPPER))
    {
      result = 1.0;
      return result;
    }
  else
    {
      result = gauss_large (x);

      if (x < 0.0)
        {
          result = 1.0 - result;
        }

    }

  return result;
}

double
gsl_cdf_gaussian_P (const double x, const double sigma)
{
  return gsl_cdf_ugaussian_P (x / sigma);
}

double
gsl_cdf_gaussian_Q (const double x, const double sigma)
{
  return gsl_cdf_ugaussian_Q (x / sigma);
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#if !defined(fvSupport_h)
#define fvSupport_h

#include <petsc.h>

PETSC_EXTERN PetscErrorCode DMPlexReconstructGradientsFVM_MulfiField(DM dm, PetscFV fvm,  Vec locX, Vec grad);
PETSC_EXTERN PetscErrorCode DMPlexGetDataFVM_MulfiField(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM);

#endif

max_stars_repo_name: pakserep/ablate

max_issues_repo_name: pakserep/ablate

max_forks_repo_name: pakserep/ablate


text:
#ifndef _CSPLINE_H_
#define _CSPLINE_H_

#include <iostream>
#include <sstream>
#include <vector>
#include <cmath>
#include <vector_types.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>

template <typename T> class cspline{
    std::vector<T> y;
    std::vector<T> t;
    std::vector<T> h;
    std::vector<T> z;
    int N;

    void h_calc();
    void b_calc(std::vector<T> &b);
    void v_calc(std::vector<T> &v);
    void u_calc(std::vector<T> &u, std::vector<T> &b);
    void z_calc(std::vector<T> &v, std::vector<T> &u);

    public:
        cspline();

        cspline(std::vector<T> &x, std::vector<T> &fx);

        void initialize(std::vector<T> &x, std::vector<T> &fx);

        void set_pointer_for_device(std::vector<float4> &spline);

        void set_pointer_for_device(std::vector<double4> &spline);

        T evaluate(T x);

};

template<typename T> void cspline<T>::h_calc() {
    for (int i = 0; i < cspline<T>::N - 1; ++i)
        cspline<T>::h.push_back(t[i + 1] - t[i]);
}

template<typename T> void cspline<T>::b_calc(std::vector<T> &b) {
    for (int i = 0; i < cspline<T>::N - 1; ++i)
        b[i] = (cspline<T>::y[i + 1] - cspline<T>::y[i])/h[i];
}

template<typename T> void cspline<T>::v_calc(std::vector<T> &v) {
    for (int i = 1; i < cspline<T>::N - 1; ++i)
        v[i - 1] = 2.0*(cspline<T>::h[i - 1] + cspline<T>::h[i]);
}

template<typename T> void cspline<T>::u_calc(std::vector<T> &u, std::vector<T> &b) {
    for (int i = 1; i < cspline<T>::N - 1; ++i)
        u[i - 1] = 6.0*(b[i] - b[i - 1]);
}

template<typename T> void cspline<T>::z_calc(std::vector<T> &v, std::vector<T> &u) {
    gsl_matrix *A = gsl_matrix_alloc(cspline<T>::N - 2, cspline<T>::N - 2);
    gsl_vector *Z = gsl_vector_alloc(cspline<T>::N - 2);
    gsl_vector *U = gsl_vector_alloc(cspline<T>::N - 2);

    for (int i = 1; i < cspline<T>::N - 1; ++i) {
        gsl_vector_set(U, i - 1, u[i - 1]);
        for (int j = 1; j < cspline<T>::N - 1; ++j) {
            if (i == j) gsl_matrix_set(A, i - 1, j - 1, v[i - 1]);
            else if (j == i + 1) gsl_matrix_set(A, i - 1, j - 1, cspline<T>::h[i]);
            else if (j == i - 1) gsl_matrix_set(A, i - 1, j - 1, cspline<T>::h[j]);
            else gsl_matrix_set(A, i - 1, j - 1, 0.0);
        }
    }

    int s;
    gsl_permutation *p = gsl_permutation_alloc(cspline<T>::N - 2);
    gsl_linalg_LU_decomp(A, p, &s);
    gsl_linalg_LU_solve(A, p, U, Z);

    for (int i = 0; i < cspline<T>::N; ++i) {
        if (i == 0 || i == cspline<T>::N - 1) cspline<T>::z.push_back(0.0);
        else cspline<T>::z.push_back(gsl_vector_get(Z, i - 1));
    }

    gsl_vector_free(Z);
    gsl_vector_free(U);
    gsl_matrix_free(A);
    gsl_permutation_free(p);
}

template<typename T> cspline<T>::cspline() {
    cspline<T>::N = 1;
}

template<typename T> cspline<T>::cspline(std::vector<T> &x, std::vector<T> &fx) {
    cspline<T>::N = fx.size();

    for (int i = 0; i < cspline<T>::N; ++i) {
        cspline<T>::t.push_back(x[i]);
        cspline<T>::y.push_back(fx[i]);
    }

    std::vector<T> b(cspline<T>::N - 1);
    std::vector<T> v(cspline<T>::N - 2);
    std::vector<T> u(cspline<T>::N - 2);

    cspline<T>::h_calc();
    cspline<T>::b_calc(b);
    cspline<T>::v_calc(v);
    cspline<T>::u_calc(u, b);
    cspline<T>::z_calc(v, u);
}

template<typename T> void cspline<T>::initialize(std::vector<T> &x, std::vector<T> &fx) {
    cspline<T>::N = fx.size();

    for (int i = 0; i < cspline<T>::N; ++i) {
        cspline<T>::t.push_back(x[i]);
        cspline<T>::y.push_back(fx[i]);
    }

    std::vector<T> b(cspline<T>::N - 1);
    std::vector<T> v(cspline<T>::N - 2);
    std::vector<T> u(cspline<T>::N - 2);

    cspline<T>::h_calc();
    cspline<T>::b_calc(b);
    cspline<T>::v_calc(v);
    cspline<T>::u_calc(u, b);
    cspline<T>::z_calc(v, u);
}

template<typename T> void cspline<T>::set_pointer_for_device(std::vector<float4> &spline) {
    for (int i = 0; i < cspline<T>::N; ++i) {
        float4 temp = {cspline<T>::t[i], cspline<T>::y[i], cspline<T>::z[i], cspline<T>::h[i]};
        spline.push_back(temp);
    }
}

template<typename T> void cspline<T>::set_pointer_for_device(std::vector<double4> &spline) {
    for (int i = 0; i < cspline<T>::N; ++i) {
        double4 temp = {cspline<T>::t[i], cspline<T>::y[i], cspline<T>::z[i], cspline<T>::h[i]};
        spline.push_back(temp);
    }
}

template<typename T> T cspline<T>::evaluate(T x) {
    if (x < cspline<T>::t[0] || x > cspline<T>::t[cspline<T>::N - 1]) {
        std::stringstream message;
        message << "The requested x value is outside of the allowed range." << std::endl;
        throw std::runtime_error(message.str());
    }

    int i = 0;
    while (cspline<T>::t[i] < x) ++i;
    i--;

    T val = (cspline<T>::z[i + 1]*(x - cspline<T>::t[i])*(x - cspline<T>::t[i])*(x - cspline<T>::t[i]))/(6.0*cspline<T>::h[i])
             + (cspline<T>::z[i]*(cspline<T>::t[i + 1] - x)*(cspline<T>::t[i + 1] - x)*(cspline<T>::t[i + 1] - x))/(6.0*cspline<T>::h[i])
             + (cspline<T>::y[i + 1]/cspline<T>::h[i] - (cspline<T>::z[i + 1]*cspline<T>::h[i])/6.0)*(x - cspline<T>::t[i])
             + (cspline<T>::y[i]/cspline<T>::h[i] - (cspline<T>::h[i]*cspline<T>::z[i])/6.0)*(cspline<T>::t[i + 1] - x);

    return val;
}

#endif


text:
/* sum/levin_u.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007, 2009 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sum.h>

int
gsl_sum_levin_u_accel (const double *array, const size_t array_size,
                       gsl_sum_levin_u_workspace * w,
                       double *sum_accel, double *abserr)
{
  return gsl_sum_levin_u_minmax (array, array_size,
                                 0, array_size - 1, w, sum_accel, abserr);
}

int
gsl_sum_levin_u_minmax (const double *array, const size_t array_size,
                        const size_t min_terms, const size_t max_terms,
                        gsl_sum_levin_u_workspace * w,
                        double *sum_accel, double *abserr)
{
  /* Ignore any trailing zeros in the array */
  size_t size = array_size;

  while (size > 0 && array[size - 1] == 0) {
    size--;
  }

  if (size == 0)
    {
      *sum_accel = 0.0;
      *abserr = 0.0;
      w->sum_plain = 0.0;
      w->terms_used = 0;
      return GSL_SUCCESS;
    }
  else if (size == 1)
    {
      *sum_accel = array[0];
      *abserr = 0.0;
      w->sum_plain = array[0];
      w->terms_used = 1;
      return GSL_SUCCESS;
    }
  else
    {
      const double SMALL = 0.01;
      const size_t nmax = GSL_MAX (max_terms, array_size) - 1;
      double noise_n = 0.0, noise_nm1 = 0.0;
      double trunc_n = 0.0, trunc_nm1 = 0.0;
      double actual_trunc_n = 0.0, actual_trunc_nm1 = 0.0;
      double result_n = 0.0, result_nm1 = 0.0;
      double variance = 0;
      size_t n;
      unsigned int i;
      int better = 0;
      int before = 0;
      int converging = 0;
      double least_trunc = GSL_DBL_MAX;
      double least_trunc_noise = GSL_DBL_MAX;
      double least_trunc_result;

      /* Calculate specified minimum number of terms.  No convergence
         tests are made, and no truncation information is stored.  */

      for (n = 0; n < min_terms; n++)
        {
          const double t = array[n];
          result_nm1 = result_n;
          gsl_sum_levin_u_step (t, n, nmax, w, &result_n);
        }

      least_trunc_result = result_n;

      variance = 0;
      for (i = 0; i < n; i++)
        {
          double dn = w->dsum[i] * GSL_MACH_EPS * array[i];
          variance += dn * dn;
        }
      noise_n = sqrt (variance);

      /* Calculate up to maximum number of terms.  Check truncation
         condition.  */

      for (; n <= nmax; n++)
        {
          const double t = array[n];

          result_nm1 = result_n;
          gsl_sum_levin_u_step (t, n, nmax, w, &result_n);

          /* Compute the truncation error directly */

          actual_trunc_nm1 = actual_trunc_n;
          actual_trunc_n = fabs (result_n - result_nm1);

          /* Average results to make a more reliable estimate of the
             real truncation error */

          trunc_nm1 = trunc_n;
          trunc_n = 0.5 * (actual_trunc_n + actual_trunc_nm1);

          noise_nm1 = noise_n;
          variance = 0;

          for (i = 0; i <= n; i++)
            {
              double dn = w->dsum[i] * GSL_MACH_EPS * array[i];
              variance += dn * dn;
            }

          noise_n = sqrt (variance);

          /* Determine if we are in the convergence region.  */

          better = (trunc_n < trunc_nm1 || trunc_n < SMALL * fabs (result_n));
          converging = converging || (better && before);
          before = better;

          if (converging)
            {
              if (trunc_n < least_trunc)
                {
                  /* Found a low truncation point in the convergence
                     region. Save it. */

                  least_trunc_result = result_n;
                  least_trunc = trunc_n;
                  least_trunc_noise = noise_n;
                }

              if (noise_n > trunc_n / 3.0)
                break;

              if (trunc_n < 10.0 * GSL_MACH_EPS * fabs (result_n))
                break;
            }

        }

      if (converging)
        {
          /* Stopped in the convergence region.  Return result and
             error estimate.  */

          *sum_accel = least_trunc_result;
          *abserr = GSL_MAX_DBL (least_trunc, least_trunc_noise);
          w->terms_used = n;
          return GSL_SUCCESS;
        }
      else
        {
          /* Never reached the convergence region.  Use the last
             calculated values.  */

          *sum_accel = result_n;
          *abserr = GSL_MAX_DBL (trunc_n, noise_n);
          w->terms_used = n;
          return GSL_SUCCESS;
        }
    }
}

int
gsl_sum_levin_u_step (const double term, const size_t n, const size_t nmax,
                      gsl_sum_levin_u_workspace * w, double *sum_accel)
{

#define I(i,j) ((i)*(nmax+1) + (j))

  if (n == 0)
    {
      *sum_accel = term;
      w->sum_plain = term;

      w->q_den[0] = 1.0 / term;
      w->q_num[0] = 1.0;

      w->dq_den[I (0, 0)] = -1.0 / (term * term);
      w->dq_num[I (0, 0)] = 0.0;

      w->dsum[0] = 1.0;

      return GSL_SUCCESS;
    }
  else
    {
      double result;
      double factor = 1.0;
      double ratio = (double) n / (n + 1.0);
      unsigned int i;
      int j;

      w->sum_plain += term;

      w->q_den[n] = 1.0 / (term * (n + 1.0) * (n + 1.0));
      w->q_num[n] = w->sum_plain * w->q_den[n];

      for (i = 0; i < n; i++)
        {
          w->dq_den[I (i, n)] = 0;
          w->dq_num[I (i, n)] = w->q_den[n];
        }

      w->dq_den[I (n, n)] = -w->q_den[n] / term;
      w->dq_num[I (n, n)] =
        w->q_den[n] + w->sum_plain * (w->dq_den[I (n, n)]);

      for (j = n - 1; j >= 0; j--)
        {
          double c = factor * (j + 1) / (n + 1);
          factor *= ratio;
          w->q_den[j] = w->q_den[j + 1] - c * w->q_den[j];
          w->q_num[j] = w->q_num[j + 1] - c * w->q_num[j];

          for (i = 0; i < n; i++)
            {
              w->dq_den[I (i, j)] =
                w->dq_den[I (i, j + 1)] - c * w->dq_den[I (i, j)];
              w->dq_num[I (i, j)] =
                w->dq_num[I (i, j + 1)] - c * w->dq_num[I (i, j)];
            }

          w->dq_den[I (n, j)] = w->dq_den[I (n, j + 1)];
          w->dq_num[I (n, j)] = w->dq_num[I (n, j + 1)];
        }

      result = w->q_num[0] / w->q_den[0];

      *sum_accel = result;

      for (i = 0; i <= n; i++)
        {
          w->dsum[i] =
            (w->dq_num[I (i, 0)] -
             result * w->dq_den[I (i, 0)]) / w->q_den[0];
        }

      return GSL_SUCCESS;
    }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/*
 * testgen.c
 * Patrick Alken
 *
 * Compile: gcc -g -O2 -Wall -o testgen testgen.c -lm -lgsl -lcblas -latlas
 *
 * Usage: testgen [options]
 *
 * -i             : incremental matrices
 * -z             : compute Schur vectors and test them
 * -n size        : size of matrices
 * -l lower-bound : lower bound for matrix elements
 * -u upper-bound : upper bound for matrix elements
 * -c num         : number of matrices to solve
 */

#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_errno.h>

typedef struct
{
  gsl_eigen_gen_workspace *gen_p;
  gsl_matrix *A;
  gsl_matrix *B;
  gsl_vector_complex *alpha;
  gsl_vector *beta;
  gsl_vector_complex *evals;

  gsl_eigen_genv_workspace *genv_p;
  gsl_matrix *Av;
  gsl_matrix *Bv;
  gsl_vector_complex *alphav;
  gsl_vector *betav;
  gsl_matrix_complex *evec;

  gsl_matrix *Q;
  gsl_matrix *Z;
  int compute_schur;

  size_t n_evals;
} gen_workspace;

gen_workspace *gen_alloc(size_t n, int compute_schur);
void gen_free(gen_workspace *w);
int gen_proc(gen_workspace *w);

/*
 * Global variables
 */
unsigned long count = 0;

/*
 * Prototypes
 */

void make_random_matrix(gsl_matrix *m, gsl_rng *r, int lower, int upper);
void make_random_integer_matrix(gsl_matrix *m, gsl_rng *r, int lower,
                                int upper);
void make_start_matrix(gsl_matrix *m, int lower);
int inc_matrix (gsl_matrix *m, int lower, int upper);
void output_matrix(const gsl_matrix *m);
void print_matrix(const gsl_matrix *m, const char *str);
int test_evals(gsl_vector_complex *obs, gsl_vector_complex *expected,
               gsl_matrix *A, gsl_matrix *B,
               const char *obsname, const char *expname);
int test_alpha(gsl_vector_complex *obs, gsl_vector_complex *expected,
               gsl_matrix *A, gsl_matrix *B, const char *obsname,
               const char *expname);
int test_beta(gsl_vector *obs, gsl_vector *expected,
              gsl_matrix *A, gsl_matrix *B, const char *obsname,
              const char *expname);
void test_schur(gsl_matrix *A, gsl_matrix *S, gsl_matrix *Q, gsl_matrix *Z);
void print_vector(const gsl_vector_complex *eval, const char *str);
int cmp(double a, double b);
int compare(const void *a, const void *b);
void sort_complex_vector(gsl_vector_complex *v);
int test_eigenvectors(const gsl_matrix *A, const gsl_matrix *B,
                      const gsl_vector_complex *alpha,
                      const gsl_vector *beta,
                      const gsl_matrix_complex *evec);

gen_workspace *
gen_alloc(size_t n, int compute_schur)
{
  gen_workspace *w;

  w = (gen_workspace *) calloc(1, sizeof(gen_workspace));

  w->gen_p = gsl_eigen_gen_alloc(n);
  w->genv_p = gsl_eigen_genv_alloc(n);

  w->A = gsl_matrix_alloc(n, n);
  w->B = gsl_matrix_alloc(n, n);
  w->alpha = gsl_vector_complex_alloc(n);
  w->beta = gsl_vector_alloc(n);
  w->evals = gsl_vector_complex_alloc(n);
  w->Av = gsl_matrix_alloc(n, n);
  w->Bv = gsl_matrix_alloc(n, n);
  w->alphav = gsl_vector_complex_alloc(n);
  w->betav = gsl_vector_alloc(n);
  w->evec = gsl_matrix_complex_alloc(n, n);
  w->compute_schur = compute_schur;

  if (compute_schur)
    {
      w->Q = gsl_matrix_alloc(n, n);
      w->Z = gsl_matrix_alloc(n, n);
      gsl_eigen_gen_params(1, 1, 0, w->gen_p);
    }

  return (w);
} /* gen_alloc() */

void
gen_free(gen_workspace *w)
{
  if (w->gen_p)
    gsl_eigen_gen_free(w->gen_p);

  if (w->genv_p)
    gsl_eigen_genv_free(w->genv_p);

  if (w->A)
    gsl_matrix_free(w->A);

  if (w->B)
    gsl_matrix_free(w->B);

  if (w->alpha)
    gsl_vector_complex_free(w->alpha);

  if (w->beta)
    gsl_vector_free(w->beta);

  if (w->evals)
    gsl_vector_complex_free(w->evals);

  if (w->Av)
    gsl_matrix_free(w->Av);

  if (w->Bv)
    gsl_matrix_free(w->Bv);

  if (w->alphav)
    gsl_vector_complex_free(w->alphav);

  if (w->betav)
    gsl_vector_free(w->betav);

  if (w->evec)
    gsl_matrix_complex_free(w->evec);

  if (w->Q)
    gsl_matrix_free(w->Q);

  if (w->Z)
    gsl_matrix_free(w->Z);

  free(w);
}

int
gen_proc(gen_workspace *w)
{
  int s1, s2, s;

  s1 = gsl_eigen_gen_QZ(w->A, w->B, w->alpha, w->beta, w->Q, w->Z, w->gen_p);
  s2 = gsl_eigen_genv(w->Av, w->Bv, w->alphav, w->betav, w->evec, w->genv_p);

  w->n_evals = w->gen_p->n_evals;

  s = 0;
  if (s1)
    s = s1;
  else if (s2)
    s = s2;

  return s;
} /* gen_proc() */

/**********************************************
 * General routines
 **********************************************/

void
make_random_matrix(gsl_matrix *m, gsl_rng *r, int lower, int upper)
{
  size_t i, j;
  size_t N = m->size1;

  for (i = 0; i < N; ++i)
  {
    for (j = 0; j < N; ++j)
    {
      gsl_matrix_set(m,
                     i,
                     j,
                     gsl_rng_uniform(r) * (upper - lower) + lower);
    }
  }
} /* make_random_matrix() */

void
make_random_integer_matrix(gsl_matrix *m, gsl_rng *r, int lower, int upper)
{
  size_t i, j;
  size_t N = m->size1;

  for (i = 0; i < N; ++i)
  {
    for (j = 0; j < N; ++j)
    {
      double a = gsl_rng_uniform(r) * (upper - lower) + lower;
      gsl_matrix_set(m, i, j, floor(a));
    }
  }
} /* make_random_integer_matrix() */

void
make_start_matrix(gsl_matrix *m, int lower)

{
  size_t i, j;
  size_t N = m->size1;

  for (i = 0; i < N; ++i)
    for (j = 0; j < N; ++j)
      gsl_matrix_set(m, i, j, (double)lower);
} /* make_start_matrix() */

int
inc_matrix (gsl_matrix *m, int lower, int upper)
{
  size_t i = 0;
  size_t N = m->size1 * m->size2;
  int carry = 1;

  for (i = 0; carry > 0 && i < N; i++)
    {
      double v = m->data[i] + carry;
      carry = (v > upper) ? 1 : 0;
      m->data[i] = (v > upper) ? lower : v;
    }

  return carry;
} /* inc_matrix() */

void
output_matrix(const gsl_matrix *m)
{
  size_t i, j;
  size_t N = m->size1;
  size_t M = m->size2;

  for (i = 0; i < N; ++i)
  {
    for (j = 0; j < M; ++j)
      {
        printf("%10.18e%s",
        /*printf("%10.18e%s",*/
               gsl_matrix_get(m, i, j),
               (j < M - 1) ? "," : ";\n");
      }
  }
}

void
print_matrix(const gsl_matrix *m, const char *str)

{
  size_t i, j;
  size_t N = m->size1;
  size_t M = m->size2;
  size_t rows, cols;
  size_t r, c;
  char buf[100];

  /*print_octave(m, str);
  return;*/

  /*rows = GSL_MIN(15, N);*/
  rows = N;
  cols = GSL_MIN(15, N);
  /*cols = N;*/

  for (i = 0; i < N; i += rows)
  {
    for (j = 0; j < M; j += cols)
    {
      r = GSL_MIN(rows, N - i);
      c = GSL_MIN(cols, N - j);

      sprintf(buf, "%s(%u:%u,%u:%u)",
              str,
              i + 1,
              i + r,
              j + 1,
              j + c);

      printf("%s = [\n", buf);

      {
        gsl_matrix_const_view v = gsl_matrix_const_submatrix(m, i, j, r, c);
        output_matrix(&v.matrix);
      }

      printf("]\n");
    }
  }
} /* print_matrix() */

int
test_evals(gsl_vector_complex *obs, gsl_vector_complex *expected,
           gsl_matrix *A, gsl_matrix *B, const char *obsname,
           const char *expname)
{
  size_t N = expected->size;
  size_t i, k;
  double max, max_abserr, max_relerr;

  max = 0.0;
  max_abserr = 0.0;
  max_relerr = 0.0;
  k = 0;

  for (i = 0; i < N; ++i)
    {
      gsl_complex z = gsl_vector_complex_get(expected, i);
      max = GSL_MAX_DBL(max, gsl_complex_abs(z));
    }

  for (i = 0; i < N; ++i)
    {
      gsl_complex z_obs = gsl_vector_complex_get(obs, i);
      gsl_complex z_exp = gsl_vector_complex_get(expected, i);

      double x_obs = GSL_REAL(z_obs);
      double y_obs = GSL_IMAG(z_obs);
      double x_exp = GSL_REAL(z_exp);
      double y_exp = GSL_IMAG(z_exp);

      double abserr_x = fabs(x_obs - x_exp);
      double abserr_y = fabs(y_obs - y_exp);
      double noise = max * GSL_DBL_EPSILON * N * N;

      max_abserr = GSL_MAX_DBL(max_abserr, abserr_x + abserr_y);

      if (abserr_x < noise && abserr_y < noise)
        continue;

      if (abserr_x > 1.0e-6 || abserr_y > 1.0e-6)
        ++k;
    }

    if (k)
      {
        printf("==== CASE %lu ===========================\n\n", count);

        print_matrix(A, "A");
        print_matrix(B, "B");

        printf("=== eval - %s ===\n", expname);
        print_vector(expected, expname);

        printf("=== eval - %s ===\n", obsname);
        print_vector(obs, obsname);

        printf("max abserr = %g  max relerr = %g\n", max_abserr, max_relerr);

        printf("=========================================\n\n");
      }

    return k;
} /* test_evals() */

/* test if A = Q S Z^t */
void
test_schur(gsl_matrix *A, gsl_matrix *S, gsl_matrix *Q, gsl_matrix *Z)
{
  const size_t N = A->size1;
  gsl_matrix *T1, *T2;
  size_t i, j, k;
  double lhs, rhs;
  double abserr;

  T1 = gsl_matrix_alloc(N, N);
  T2 = gsl_matrix_alloc(N, N);

  /* compute T1 = S Z^t */
  gsl_blas_dgemm(CblasNoTrans,
                 CblasTrans,
                 1.0,
                 S,
                 Z,
                 0.0,
                 T1);

  /* compute T2 = Q T1 = Q S Z^t */
  gsl_blas_dgemm(CblasNoTrans,
                 CblasNoTrans,
                 1.0,
                 Q,
                 T1,
                 0.0,
                 T2);

  k = 0;
  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          lhs = gsl_matrix_get(A, i, j);
          rhs = gsl_matrix_get(T2, i, j);

          abserr = fabs(lhs - rhs);

          if (abserr > 1.0e-6)
            ++k;
        }
    }

  if (k)
    {
      printf("==== CASE %lu ===========================\n\n", count);

      print_matrix(A, "A");

      printf("=== Schur Form matrix ===\n");
      print_matrix(S, "S");

      printf("=== Left Schur matrix ===\n");
      print_matrix(Q, "Q");

      printf("=== Right Schur matrix ===\n");
      print_matrix(Z, "Z");

      printf("=== Q S Z^t ===\n");
      print_matrix(T1, "Q S Z^t");

      printf("=== A - Q S Z^t ===\n");
      gsl_matrix_sub(T2, A);
      print_matrix(T1, "A - Q S Z^t");

      printf("=========================================\n\n");
    }

  gsl_matrix_free(T1);
  gsl_matrix_free(T2);
} /* test_schur() */

void
print_vector(const gsl_vector_complex *eval, const char *str)
{
  size_t N = eval->size;
  size_t i;
  gsl_complex z;

  printf("%s = [\n", str);

  for (i = 0; i < N; ++i)
    {
      z = gsl_vector_complex_get(eval, i);
      printf("%.18e + %.18ei;\n", GSL_REAL(z), GSL_IMAG(z));
    }

  printf("]\n");
} /* print_vector() */

int
cmp(double a, double b)
{
  return ((a > b) ? 1 : ((a < b) ? -1 : 0));
} /* cmp() */

int
compare(const void *a, const void *b)
{
  const double *x = a;
  const double *y = b;
  int r1 = cmp(y[0], x[0]);
  int r2 = cmp(y[1], x[1]);

  if (!gsl_finite(x[0]))
    return 1;
  if (!gsl_finite(y[0]))
    return -1;

  if (fabs(x[0] - y[0]) < 1.0e-8)
    {
      /* real parts are very close to each other */
      return r2;
    }
  else
    {
      return r1 ? r1 : r2;
    }
} /* compare() */

void
sort_complex_vector(gsl_vector_complex *v)
{
  qsort(v->data, v->size, 2 * sizeof(double), &compare);
} /* sort_complex_vector() */

int
test_eigenvectors(const gsl_matrix *A, const gsl_matrix *B,
                  const gsl_vector_complex *alpha, const gsl_vector *beta,
                  const gsl_matrix_complex *evec)
{
  const size_t N = A->size1;
  size_t i, j;
  int k, s;
  gsl_matrix_complex *ma, *mb;
  gsl_vector_complex *x, *y;
  gsl_complex z_one, z_zero;

  ma = gsl_matrix_complex_alloc(N, N);
  mb = gsl_matrix_complex_alloc(N, N);
  y = gsl_vector_complex_alloc(N);
  x = gsl_vector_complex_alloc(N);

  /* ma <- A, mb <- B */
  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          gsl_complex z;

          GSL_SET_COMPLEX(&z, gsl_matrix_get(A, i, j), 0.0);
          gsl_matrix_complex_set(ma, i, j, z);

          GSL_SET_COMPLEX(&z, gsl_matrix_get(B, i, j), 0.0);
          gsl_matrix_complex_set(mb, i, j, z);
        }
    }

  GSL_SET_COMPLEX(&z_one, 1.0, 0.0);
  GSL_SET_COMPLEX(&z_zero, 0.0, 0.0);

  s = 0;

  /* check eigenvalues */
  for (i = 0; i < N; ++i)
    {
      gsl_vector_complex_const_view vi = gsl_matrix_complex_const_column(evec, i);
      gsl_complex ai = gsl_vector_complex_get(alpha, i);
      double bi = gsl_vector_get(beta, i);
      double norm = gsl_blas_dznrm2(&vi.vector);

      /* check that eigenvector is normalized */
      gsl_test_rel(norm, 1.0, N * GSL_DBL_EPSILON, "case %u, normalized",
                   count);

      /* compute x = alpha * B * v */
      gsl_blas_zgemv(CblasNoTrans, z_one, mb, &vi.vector, z_zero, x);
      gsl_blas_zscal(ai, x);

      /* compute y = beta * A v */
      gsl_blas_zgemv(CblasNoTrans, z_one, ma, &vi.vector, z_zero, y);
      gsl_blas_zdscal(bi, y);

      k = 0;

      /* now test if y = x */
      for (j = 0; j < N; ++j)
        {
          gsl_complex z;
          double lhs_r, lhs_i;
          double rhs_r, rhs_i;

          z = gsl_vector_complex_get(y, j);
          lhs_r = GSL_REAL(z);
          lhs_i = GSL_IMAG(z);

          z = gsl_vector_complex_get(x, j);
          rhs_r = GSL_REAL(z);
          rhs_i = GSL_IMAG(z);

          if (fabs(lhs_r - rhs_r) > 1e10 * GSL_DBL_EPSILON)
            ++k;
          if (fabs(lhs_i - rhs_i) > 1e10 * GSL_DBL_EPSILON)
            ++k;
        }

      if (k)
        {
          s++;

          printf("==== CASE %lu ===========================\n\n", count);

          print_matrix(A, "A");
          print_matrix(B, "B");

          printf("alpha = %.10e + %.10ei\n", GSL_REAL(ai), GSL_IMAG(ai));
          printf("beta = %.10e\n", bi);
          printf("alpha/beta = %.10e + %.10ei\n", GSL_REAL(ai)/bi, GSL_IMAG(ai)/bi);

          print_vector(&vi.vector, "v");

          print_vector(y, "beta*A*v");

          print_vector(x, "alpha*B*v");

          printf("=========================================\n\n");
        }
    }

  gsl_matrix_complex_free(ma);
  gsl_matrix_complex_free(mb);
  gsl_vector_complex_free(y);
  gsl_vector_complex_free(x);

  return s;
} /* test_eigenvectors() */

int
main(int argc, char *argv[])
{
  gen_workspace *gen_workspace_p;
  size_t N;
  int c;
  int lower;
  int upper;
  int incremental;
  size_t nmat;
  gsl_matrix *A, *B;
  gsl_rng *r;
  int s;
  int compute_schur;
  size_t i;

  gsl_ieee_env_setup();
  gsl_rng_env_setup();
  /*gsl_set_error_handler_off();*/

  N = 30;
  lower = -10;
  upper = 10;
  incremental = 0;
  nmat = 0;
  compute_schur = 0;

  while ((c = getopt(argc, argv, "ic:n:l:u:z")) != (-1))
    {
      switch (c)
        {
          case 'i':
            incremental = 1;
            break;

          case 'n':
            N = strtol(optarg, NULL, 0);
            break;

          case 'l':
            lower = strtol(optarg, NULL, 0);
            break;

          case 'u':
            upper = strtol(optarg, NULL, 0);
            break;

          case 'c':
            nmat = strtoul(optarg, NULL, 0);
            break;

          case 'z':
            compute_schur = 1;
            break;

          case '?':
          default:
            printf("usage: %s [-i] [-z] [-n size] [-l lower-bound] [-u upper-bound] [-c num]\n", argv[0]);
            exit(1);
            break;
        } /* switch (c) */
    }

  A = gsl_matrix_alloc(N, N);
  B = gsl_matrix_alloc(N, N);
  gen_workspace_p = gen_alloc(N, compute_schur);

  r = gsl_rng_alloc(gsl_rng_default);

  if (incremental)
    {
      make_start_matrix(A, lower);

      /* we need B to be non-singular */
      make_random_integer_matrix(B, r, lower, upper);
    }

  fprintf(stderr, "testing N = %d", N);
  if (incremental)
    fprintf(stderr, " incrementally");
  else
    fprintf(stderr, " randomly");

  fprintf(stderr, " on element range [%d, %d]", lower, upper);

  if (compute_schur)
    fprintf(stderr, ", with Schur vectors");

  fprintf(stderr, "\n");

  while (1)
    {
      if (nmat && (count >= nmat))
        break;

      ++count;

      if (!incremental)
        {
          make_random_matrix(A, r, lower, upper);
          make_random_matrix(B, r, lower, upper);
        }
      else
        {
          s = inc_matrix(A, lower, upper);
          if (s)
            break; /* all done */

          make_random_integer_matrix(B, r, lower, upper);
        }

      /*if (count != 23441872)
        continue;*/

      /* make copies of matrices */
      gsl_matrix_memcpy(gen_workspace_p->A, A);
      gsl_matrix_memcpy(gen_workspace_p->B, B);
      gsl_matrix_memcpy(gen_workspace_p->Av, A);
      gsl_matrix_memcpy(gen_workspace_p->Bv, B);

      /* compute eigenvalues with GSL */
      s = gen_proc(gen_workspace_p);

      if (s != GSL_SUCCESS)
        {
          printf("=========== CASE %lu ============\n", count);
          printf("Failed to converge: found %u eigenvalues\n",
                 gen_workspace_p->n_evals);
          print_matrix(A, "A");
          print_matrix(B, "B");
          print_matrix(gen_workspace_p->Av, "S");
          print_matrix(gen_workspace_p->Bv, "T");
          continue;
        }

      /* compute alpha / beta vectors */
      for (i = 0; i < N; ++i)
        {
          double beta = gsl_vector_get(gen_workspace_p->beta, i);
          gsl_complex alpha =
            gsl_vector_complex_get(gen_workspace_p->alpha, i);
          gsl_complex z;

          if (!gsl_finite(beta) || !gsl_finite(GSL_REAL(alpha)) ||
              !gsl_finite(GSL_IMAG(alpha)))
            {
              printf("nan/inf in element %u of (alpha,beta): alphar = %g, alphai = %g, beta = %g\n",
                     i, GSL_REAL(alpha), GSL_IMAG(alpha), beta);
            }

          if (beta == 0.0)
            GSL_SET_COMPLEX(&z, GSL_POSINF, GSL_POSINF);
          else
            z = gsl_complex_div_real(alpha, beta);

          gsl_vector_complex_set(gen_workspace_p->evals, i, z);
        }

      test_eigenvectors(A,
                        B,
                        gen_workspace_p->alphav,
                        gen_workspace_p->betav,
                        gen_workspace_p->evec);

      if (compute_schur)
        {
          test_schur(A,
                     gen_workspace_p->A,
                     gen_workspace_p->Q,
                     gen_workspace_p->Z);
          test_schur(B,
                     gen_workspace_p->B,
                     gen_workspace_p->Q,
                     gen_workspace_p->Z);
        }
    }

  gsl_matrix_free(A);
  gsl_matrix_free(B);
  gen_free(gen_workspace_p);
  gsl_rng_free(r);

  return 0;
} /* main() */

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#include <assert.h>
#include <cblas.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#include "cnn_calc.h"
#include "cnn_init.h"
#include "cnn_private.h"

inline void cnn_restrict(float* mat, int size, float limit)
{
    for (int __i = 0; __i < size; __i++)
    {
        mat[__i] = fminf(mat[__i], limit);
    }
}

void cnn_mat_update(struct CNN_MAT* matPtr, float lRate, float limit)
{
    int size = matPtr->rows * matPtr->cols;

    // Limit gradient and update weight
#ifdef CNN_WITH_CUDA
    cnn_fminf_gpu(matPtr->grad, matPtr->grad, size, limit);
    cublasSaxpy(cnnInit.blasHandle, size, &lRate, matPtr->grad, 1, matPtr->mat,
                1);
#else
    cnn_restrict(matPtr->grad, size, limit);
    cblas_saxpy(size, lRate, matPtr->grad, 1, matPtr->mat, 1);
#endif

    // Clear gradient
#ifdef CNN_WITH_CUDA
    cudaMemset
#else
    memset
#endif
        (matPtr->grad, 0, size * sizeof(float));
}

void cnn_update(cnn_t cnn, float lRate, float gradLimit)
{
    int i;
    struct CNN_CONFIG* cfgRef;
    union CNN_LAYER* layerRef;

    // Set reference
    layerRef = cnn->layerList;
    cfgRef = &cnn->cfg;

    // Update network and clear gradient
    for (i = cfgRef->layers - 1; i > 0; i--)
    {
        // Clear layer gradient
#ifdef CNN_WITH_CUDA
        cudaMemset
#else
        memset
#endif
            (layerRef[i].outMat.data.grad, 0,
             sizeof(float) * layerRef[i].outMat.data.rows *
                 layerRef[i].outMat.data.cols);

        switch (cfgRef->layerCfg[i].type)
        {
            // Fully connected
            case CNN_LAYER_FC:
                cnn_mat_update(&layerRef[i].fc.weight, lRate, gradLimit);
                cnn_mat_update(&layerRef[i].fc.bias, lRate, gradLimit);
                break;

            // Convolution
            case CNN_LAYER_CONV:
                cnn_mat_update(&layerRef[i].conv.kernel, lRate, gradLimit);

#if defined(CNN_CONV_BIAS_FILTER) || defined(CNN_CONV_BIAS_LAYER)
                cnn_mat_update(&layerRef[i].conv.bias, lRate, gradLimit);
#endif
                break;

            // Batch normalization
            case CNN_LAYER_BN:
                cnn_mat_update(&layerRef[i].bn.bnScale, lRate, gradLimit);
                cnn_mat_update(&layerRef[i].bn.bnBias, lRate, gradLimit);
                break;

            case CNN_LAYER_INPUT:
            case CNN_LAYER_ACTIV:
            case CNN_LAYER_POOL:
            case CNN_LAYER_DROP:
                break;
        }
    }
}

static inline void cnn_backward_kernel(cnn_t cnn, float* errGrad)
{
    int i;

    struct CNN_CONFIG* cfgRef;
    union CNN_LAYER* layerRef;

    // Set reference
    layerRef = cnn->layerList;
    cfgRef = &cnn->cfg;

    // Backpropagation
    for (i = cfgRef->layers - 1; i > 0; i--)
    {
        switch (cfgRef->layerCfg[i].type)
        {
            // Fully connected
            case CNN_LAYER_FC:
                cnn_backward_fc(layerRef, cfgRef, i);
                break;

            // Activation function
            case CNN_LAYER_ACTIV:
                cnn_backward_activ(layerRef, cfgRef, i);
                break;

            // Convolution
            case CNN_LAYER_CONV:
                cnn_backward_conv(layerRef, cfgRef, i);
                break;

            // Pooling
            case CNN_LAYER_POOL:
                cnn_backward_pool(layerRef, cfgRef, i);
                break;

            // Dropout
            case CNN_LAYER_DROP:
                cnn_backward_drop(layerRef, cfgRef, i);
                break;

            // Batch normalization
            case CNN_LAYER_BN:
                cnn_backward_bn(layerRef, cfgRef, i);
                break;

            default:
                assert(!"Invalid layer type");
        }
    }
}

void cnn_backward(cnn_t cnn, float* errGrad)
{
    int size;

    struct CNN_CONFIG* cfgRef;
    union CNN_LAYER* layerRef;

    // Set reference
    layerRef = cnn->layerList;
    cfgRef = &cnn->cfg;

    // Copy gradient vector
    size = sizeof(float) * layerRef[cfgRef->layers - 1].outMat.data.rows *
           layerRef[cfgRef->layers - 1].outMat.data.cols;

#ifdef CNN_WITH_CUDA
    cudaMemcpy(layerRef[cfgRef->layers - 1].outMat.data.grad, errGrad, size,
               cudaMemcpyHostToDevice);
#else
    memcpy(layerRef[cfgRef->layers - 1].outMat.data.grad, errGrad, size);
#endif

    // Backpropagation
    cnn_backward_kernel(cnn, errGrad);
}

#ifdef CNN_WITH_CUDA
void cnn_backward_gpu(cnn_t cnn, float* errGrad)
{
    int size;

    struct CNN_CONFIG* cfgRef;
    union CNN_LAYER* layerRef;

    // Set reference
    layerRef = cnn->layerList;
    cfgRef = &cnn->cfg;

    // Copy gradient vector
    size = sizeof(float) * layerRef[cfgRef->layers - 1].outMat.data.rows *
           layerRef[cfgRef->layers - 1].outMat.data.cols;
    cudaMemcpy(layerRef[cfgRef->layers - 1].outMat.data.grad, errGrad, size,
               cudaMemcpyDeviceToDevice);

    // Backpropagation
    cnn_backward_kernel(cnn, errGrad);
}
#endif

static inline void cnn_forward_kernel(cnn_t cnn, float* inputMat,
                                      float* outputMat)
{
    int i;

    struct CNN_CONFIG* cfgRef;
    union CNN_LAYER* layerRef;

    // Set reference
    layerRef = cnn->layerList;
    cfgRef = &cnn->cfg;

    // Forward computation
    for (i = 1; i < cfgRef->layers; i++)
    {
        switch (cfgRef->layerCfg[i].type)
        {
            // Fully connected
            case CNN_LAYER_FC:
                cnn_forward_fc(layerRef, cfgRef, i);
                break;

            // Activation function
            case CNN_LAYER_ACTIV:
                cnn_forward_activ(layerRef, cfgRef, i);
                break;

            // Convolution
            case CNN_LAYER_CONV:
                cnn_forward_conv(layerRef, cfgRef, i);
                break;

            // Pooling
            case CNN_LAYER_POOL:
                cnn_forward_pool(layerRef, cfgRef, i);
                break;

            // Dropout
            case CNN_LAYER_DROP:
                if (cnn->opMode == CNN_OPMODE_TRAIN)
                {
                    cnn_forward_drop(layerRef, cfgRef, i);
                }
                else
                {
                    cnn_recall_drop(layerRef, cfgRef, i);
                }

                break;

            // Batch normalization
            case CNN_LAYER_BN:
                if (cnn->opMode == CNN_OPMODE_TRAIN)
                {
                    cnn_forward_bn(layerRef, cfgRef, i);
                }
                else
                {
                    cnn_recall_bn(layerRef, cfgRef, i);
                }

                break;

            default:
                assert(!"Invalid layer type");
        }
    }
}

void cnn_forward(cnn_t cnn, float* inputMat, float* outputMat)
{
    int size;

    struct CNN_CONFIG* cfgRef;
    union CNN_LAYER* layerRef;

    // Set reference
    layerRef = cnn->layerList;
    cfgRef = &cnn->cfg;

    // Copy input
    size = sizeof(float) * layerRef[0].outMat.data.rows *
           layerRef[0].outMat.data.cols;
#ifdef CNN_WITH_CUDA
    cudaMemcpy(layerRef[0].outMat.data.mat, inputMat, size,
               cudaMemcpyHostToDevice);
#else
    memcpy(layerRef[0].outMat.data.mat, inputMat, size);
#endif

    // Forward computation
    cnn_forward_kernel(cnn, inputMat, outputMat);

    // Copy output
    if (outputMat != NULL)
    {
        size = sizeof(float) * layerRef[cfgRef->layers - 1].outMat.data.rows *
               layerRef[cfgRef->layers - 1].outMat.data.cols;
#ifdef CNN_WITH_CUDA
        cudaMemcpy(outputMat, layerRef[cfgRef->layers - 1].outMat.data.mat,
                   size, cudaMemcpyDeviceToHost);
#else
        memcpy(outputMat, layerRef[cfgRef->layers - 1].outMat.data.mat, size);
#endif
    }
#ifdef CNN_WITH_CUDA
    else
    {
        cudaDeviceSynchronize();
    }
#endif
}

#ifdef CNN_WITH_CUDA
void cnn_forward_gpu(cnn_t cnn, float* inputMat, float* outputMat)
{
    int size;

    struct CNN_CONFIG* cfgRef;
    union CNN_LAYER* layerRef;

    // Set reference
    layerRef = cnn->layerList;
    cfgRef = &cnn->cfg;

    // Copy input
    size = sizeof(float) * layerRef[0].outMat.data.rows *
           layerRef[0].outMat.data.cols;
    cudaMemcpy(layerRef[0].outMat.data.mat, inputMat, size,
               cudaMemcpyDeviceToDevice);

    // Forward computation
    cnn_forward_kernel(cnn, inputMat, outputMat);

    // Copy output
    if (outputMat != NULL)
    {
        size = sizeof(float) * layerRef[cfgRef->layers - 1].outMat.data.rows *
               layerRef[cfgRef->layers - 1].outMat.data.cols;
        cudaMemcpy(outputMat, layerRef[cfgRef->layers - 1].outMat.data.mat,
                   size, cudaMemcpyDeviceToDevice);
    }
    else
    {
        cudaDeviceSynchronize();
    }
}
#endif

max_stars_repo_name: jamesljlster/cnn

max_issues_repo_name: jamesljlster/cnn

max_forks_repo_name: jamesljlster/cnn


text:
#pragma once

#include <gsl\gsl>
#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include <map>
#include "DrawableGameComponent.h"
#include "FullScreenRenderTarget.h"
#include "FullScreenQuad.h"

namespace Library
{
	class Texture2D;
}

namespace Rendering
{
	enum class DistortionMaps
	{
		Glass,
		Text,
		NoDistortion,
		End
	};

	class DiffuseLightingDemo;

	class DistortionMappingDemo final : public Library::DrawableGameComponent
	{
	public:
		DistortionMappingDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);
		DistortionMappingDemo(const DistortionMappingDemo&) = delete;
		DistortionMappingDemo(DistortionMappingDemo&&) = default;
		DistortionMappingDemo& operator=(const DistortionMappingDemo&) = default;
		DistortionMappingDemo& operator=(DistortionMappingDemo&&) = default;
		~DistortionMappingDemo();

		std::shared_ptr<DiffuseLightingDemo> DiffuseLighting() const;

		float DisplacementScale() const;
		void SetDisplacementScale(float displacementScale);

		DistortionMaps DistortionMap() const;
		const std::string& DistortionMapString() const;
		void SetDistortionMap(DistortionMaps distortionMap);

		virtual void Initialize() override;
		virtual void Update(const Library::GameTime& gameTime) override;
		virtual void Draw(const Library::GameTime& gameTime) override;

	private:
		static const std::map<DistortionMaps, std::string> DistortionMapNames;

		struct PixelCBufferPerObject
		{
			float DisplacementScale{ 1.0f };
			DirectX::XMFLOAT3 Padding{ 0.0f, 0.0f, 0.0f };
		};

		std::shared_ptr<DiffuseLightingDemo> mDiffuseLightingDemo;
		Library::FullScreenRenderTarget mRenderTarget;
		Library::FullScreenQuad mFullScreenQuad;
		winrt::com_ptr<ID3D11Buffer> mPixelCBufferPerObject;
		PixelCBufferPerObject mPixelCBufferPerObjectData;
		std::map<DistortionMaps, std::shared_ptr<Library::Texture2D>> mDistortionMaps;
		DistortionMaps mActiveDistortionMap{ DistortionMaps::Glass };
	};
}


text:
/*--------------------------------------------------------------------
 * $Id$
 *
 * This file is part of libRadtran.
 * Copyright (c) 1997-2012 by Arve Kylling, Bernhard Mayer,
 *                            Claudia Emde, Robert Buras
 *
 * ######### Contact info: http://www.libradtran.org #########
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 *--------------------------------------------------------------------*/

/* 03.04.2002: Fixed problem with finding correct grid cell for direct beam, AKy. */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#ifndef _AIX
#include <getopt.h>
#else
#include <unistd.h>
#endif

#if HAVE_LIBGSL
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_vegas.h>
#endif

#include "ascii.h"

#define PROGRAM "angres"
#define VERSION "0.3"
#define DATE "3 Feb, 2003"

#define DIM_ANG_INT 2 /* Integration to be made over 2 dimensions */

#ifndef PI
#define PI 3.14159265358979323846264338327
#endif

typedef struct {
  int      nphi;
  int      numu;
  double*  umu;
  double*  phi;
  double** prod;
} mcinput;

/*************************************************************/
/* print usage information                                   */
/*************************************************************/

/*
  Documentation in latex for the User's Guide. Extracted by sdoc2.awk
  and input in the doc/tools.tex.
*/
/*
  <lpdoc>

\subsection{Angular response and tilted surfaces - \codeidx{angres}}
The \code{angres} tool takes a precalculated radiance field and integrates
it over a given angular area using any angular response. Typical usages of
\code{angres} are calculation of radiation on tilted surfaces and estimation
of effects of imperfect angular response functions.

The \code{angres} tool is invoked as follows:
\begin{Verbatim}[fontsize=\footnotesize]
 angres angres_file raddis_file
\end{Verbatim}
The two required input files will be read by the \code{angres} tool.
\begin{description}
\item[angres\_file] is a two column file with the first column holding
     the angle and the second column the angular response, e.g. a
     measured cosine response. To generate standard angular response
     function see the \code{make\_angres\_func} tool.

\item[raddis\_file] holds the radiance distribution as output from uvspec
     with the disort solvers for one single wavelength.
\end{description}
After reading the two input files the angular response will be tilted
and rotated if specified with the \code{-t} and \code{-r} options respectively. Finally
the product of the resulting angular response and radiance distribution
field are integrated using Monte Carlo methods to yield the effective
response. The integration is done for the diffuse radiation field only. To
include the direct contribution the -s and -z options must be set to give
the direction of the sun.

Output is 3 numbers:
\begin{enumerate}
\item The integral of the diffuse radiation field times angular response.
\item Estimated absolute error of the above integral.
\item The integral of the diffuse+direct radiation field times angular response
       (requires that \code{-s} and \code{-z} are specified, otherwise same as first number
\end{enumerate}

The angles in the \file{angres\_file} must be in radians if not
the \code{-a} option is used. The \file{raddis\_file} must contain
output from \code{uvspec} run for one single wavelength with one of
the disort solvers and with \code{phi} and \code{umu} set.
Note that the angles in the \file{angres\_file} must follow the same
conventions as for the disort algorithm. This is different from
that typically used when reporting measurements of the angular response.

The \code{angres} tool accepts the following command line options:

\begin{description}
\item[-h]  show this page.
\item[-c]  number of random points used for Monte Carlo integration.
\item[-i]  The diffuse radiation is assumed to be isotropic.
\item[-a]  angular response angle given in degrees and not cosine of angle.
\item[-r]  rotation angle in degrees.
\item[-t]  tilt angle in degrees.
\item[-s]  solar zenith angle in degrees.
\item[-z]  solar azimuth angle in degrees.
\item[-p]  pgm files are made of the angular response before
      and after tilt and rotate.
\end{description}

Sample \code{angres} input and output files are found in the \file{examples} directory.
The following
\begin{Verbatim}[fontsize=\footnotesize]
 angres examples/ANGRES_1_ANG.DAT  \
        examples/ANGRES_RADDIS_1.DAT -a -t -r 0 -s 32 -z 0
\end{Verbatim}
calculates the radiation on a horisontal surface given the angular response
in \file{examples/ANGRES_1_ANG.DAT}. The input used to calculate the radiance file is
given in the start of \file{examples/ANGRES_RADDIS_1.DAT}.

An example of the use of \code{angres} together with \code{uvspec} is given in
\citet[section 4.6]{mayer2005}.
  </lpdoc>
*/

static void print_usage (char* filename) {
  fprintf (stderr, "%s - %s   \n\n", PROGRAM, VERSION);
  fprintf (stderr, "written by  Arve Kylling,\n");
  fprintf (stderr, "            Norwegian Institute for Air Research (NILU)\n");
  fprintf (stderr, "            P.O. Box 100, 2027 Kjeller, Norway\n");
  fprintf (stderr, "Version %s finished %s\n\n", VERSION, DATE);
  fprintf (stderr, "Be aware that this program is a beta version! Please report any\n");
  fprintf (stderr, "kind of error to the author, including the error message and the\n");
  fprintf (stderr, "database being processed when the error occured. Thanks!\n\n");

  fprintf (stderr, "USAGE: %s [options] angres_file raddis_file \n\n", filename);
  fprintf (stderr, "%s will read the two required input files:\n", PROGRAM);
  fprintf (stderr, "  angres_file is a two column file with the first column holding\n");
  fprintf (stderr, "     the angle and the second column the angular response, e.g. a\n");
  fprintf (stderr, "     measured cosine response.\n");
  fprintf (stderr, "  raddis_file holds the radiance distribution as output from uvspec\n");
  fprintf (stderr, "     with the disort solvers for one single wavelength.\n");
  fprintf (stderr, "After reading the two input files the angular response will be tilted\n");
  fprintf (stderr, "and rotated if specified with the -t and -r options respectively. Finally\n");
  fprintf (stderr, "the product of the resulting angular response and radiance distribution\n");
  fprintf (stderr, "field are integrated using Monte Carlo methods to yield the effective\n");
  fprintf (stderr, "response. The integration is done for the diffuse radiation field only. To\n");
  fprintf (stderr, "include the direct contribution the -s and -z options must be set to give\n");
  fprintf (stderr, "the direction of the sun. \n");
  fprintf (stderr, "\n");
  fprintf (stderr, "Output is 3 numbers:\n");
  fprintf (stderr, "   1: The integral of the diffuse radiation field times angular response.\n");
  fprintf (stderr, "   2: Estimated absolute error of the above integral.\n");
  fprintf (stderr, "   3: The integral of the diffuse+direct radiation field times angular response\n");
  fprintf (stderr, "       (requires that -s and -z are specified, otherwise same as 1).\n");
  fprintf (stderr, "\n");
  fprintf (stderr, "The angles in the angres_file must be in radians if not\n");
  fprintf (stderr, "the -a option is used. The raddis_file must contain\n");
  fprintf (stderr, "output from uvspec run for one single wavelength with one of\n");
  fprintf (stderr, "the disort solvers and with phi and umu set.\n");
  fprintf (stderr, "Note that the angles in the angres_file must follow the same\n");
  fprintf (stderr, "conventions as for the disort algorithm. This is different from\n");
  fprintf (stderr, "that typically used when reporting measurements of the angular response.\n");
  fprintf (stderr, "\n");
  fprintf (stderr, " command line options:\n");
  fprintf (stderr, "  -h  show this page\n");
  fprintf (stderr, "  -c  number of random points used for Monte Carlo integration\n");
  fprintf (stderr, "  -i  The diffuse radiation is assumed to be isotropic\n");
  fprintf (stderr, "  -a  angular response angle given in degrees and not cosine of angle\n");
  fprintf (stderr, "  -r  rotation angle in degrees\n");
  fprintf (stderr, "  -t  tilt angle in degrees\n");
  fprintf (stderr, "  -s  solar zenith angle in degrees\n");
  fprintf (stderr, "  -z  solar azimuth angle in degrees\n");
  fprintf (stderr, "  -p  pgm files are made of the angular response before\n");
  fprintf (stderr, "      and after tilt and rotate\n");
}

static int get_options (int     argc,
                        char**  argv,
                        char*   programname,
                        char*   angresfilename,
                        char*   raddisfilename,
                        int*    pgm,
                        int*    dev,
                        int*    ang,
                        size_t* calls,
                        double* sza,
                        double* phi0,
                        double* theta_t,
                        double* phi_r,
                        int*    incdir,
                        int*    isotropic) {
  int c = 0;

  /* default settings */

  /* save name of program */
  strncpy (programname, argv[0], FILENAME_MAX - 1);

  while ((c = getopt (argc, argv, "ahipc:d:s:t:r:z:")) != EOF) {
    switch (c) {
    case 'h': /* help */
      print_usage (programname);
      return (-1);
      break;
    case 'd':
      *dev = atoi (optarg);
      break;
    case 'a':
      *ang = 1;
      break;
    case 'i':
      *isotropic = 1;
      break;
    case 'p':
      *pgm = 1;
      break;
    case 'c':
      *calls = atoi (optarg);
      break;
    case 'r':
      *phi_r = atof (optarg);
      break;
    case 's':
      *sza    = atof (optarg);
      *incdir = 1;
      break;
    case 'z':
      *phi0   = atof (optarg);
      *incdir = 1;
      break;
    case 't':
      *theta_t = atof (optarg);
      break;
    case '?':
      print_usage (programname);
      return (-1);
      break;
    default:
      print_usage (programname);
      return (-1);
    }
  }

  /* check number of remaining command line arguments */
  if (argc - optind != 2) {
    print_usage (programname);
    return -1;
  }

  /* save command line arguments */
  strncpy (angresfilename, argv[optind + 0], FILENAME_MAX - 1);
  strncpy (raddisfilename, argv[optind + 1], FILENAME_MAX - 1);

  return 0; /* if o.k. */
}

double**
tilt_and_rotate (double* aangres, int nphi, double* phi, int numu, double* umu, int pgm, double theta_t, double phi_r, int dev) {
  int status = 0, i = 0, j = 0, k = 0, it = 0, jt = 0, ia = 0, ja = 0, nt;
  /* double umu_t; */
  double   theta = 0, x = 0, y = 0, z = 0, xp = 0, yp = 0, zp = 0, r = 0;
  double   a11, a12, a13, a21, a22, a23, a31, a32, a33;
  double   rn = 0, theta_n = 0, phi_n = 0;
  double **angres = NULL, **angrestilted = NULL;
  double** angresindcnt = NULL;
  double   delphirad    = 0;
  FILE*    pgmfile      = NULL;

  /* allocate memory for the angular response function and the tilted and rotated one */
  status = ASCII_calloc_double (&angres, numu, nphi);
  if (status != 0) {
    fprintf (stderr, "tilt_and_rotate: error allocating memory for angular response function\n");
    /*    return status; */
  }
  status = ASCII_calloc_double (&angrestilted, numu, nphi);
  if (status != 0) {
    fprintf (stderr, "tilt_and_rotate: error allocating memory for angular response function\n");
    /*    return status; */
  }

  status = ASCII_calloc_double (&angresindcnt, numu, nphi);
  if (status != 0) {
    fprintf (stderr, "tilt_and_rotate: error allocating memory for angular response function\n");
    /*    return status; */
  }

  /* The read angular response is assumed to have no azimuth dependence, */
  /* hence it is the same for all phi                                    */
  for (i = 0; i < numu; i++)
    for (j = 0; j < nphi; j++)
      angres[i][j] = aangres[i];

  if (pgm) {
    if ((pgmfile = fopen ("angres_orginal.pgm", "w")) == NULL) {
      fprintf (stderr, "Error opening angres_orginal.pgm for writing\n");
      return NULL;
    }

    fprintf (pgmfile, "P2\n");
    fprintf (pgmfile, "%d %d\n", nphi, numu);
    fprintf (pgmfile, "256\n");
    for (i = 0; i < numu; i++) {
      for (j = 0; j < nphi; j++) {
        k = (int)(angres[i][j] * 256);
        fprintf (pgmfile, "%d\n", k);
      }
    }
    fclose (pgmfile);
  }

  if (dev == 1)
    fprintf (stdout, "Tilt: %f Rotate: %f\n", theta_t, phi_r);
  /* commented out by RB, because was not used but caused compiler warnings */
  /* umu_t = cos(theta_t*PI/180.0); */

  /* Rotation matrix. */
  /* Eq 4-47., p. 147, Goldstein, Classical Mechanics. psi=0 for our purposes */

  /*  a11 =  cos(phi_r*PI/180.0);                      */
  /*  a12 = -cos(theta_t*PI/180.0)*sin(phi_r*PI/180.0); */
  /*  a13 =  sin(theta_t*PI/180.0)*sin(phi_r*PI/180.0); */
  /*  a21 =  sin(phi_r*PI/180.0);                      */
  /*  a22 =  cos(theta_t*PI/180.0)*cos(phi_r*PI/180.0); */
  /*  a23 = -sin(theta_t*PI/180.0)*cos(phi_r*PI/180.0); */
  /*  a31 = 0;                                        */
  /*  a32 =  sin(theta_t*PI/180.0);                    */
  /*  a33 =  cos(theta_t*PI/180.0);                    */

  /* First only tilt */
  a11 = 1;
  a12 = 0;
  a13 = 0;
  a21 = 0;
  a22 = cos (theta_t * PI / 180.0);
  a23 = -sin (theta_t * PI / 180.0);
  a31 = 0;
  a32 = sin (theta_t * PI / 180.0);
  a33 = cos (theta_t * PI / 180.0);

  if (dev == 1) {
    printf ("a11: %7.3f, a12: %7.3f a13: %7.3f\n", a11, a12, a13);
    printf ("a21: %7.3f, a22: %7.3f a23: %7.3f\n", a21, a22, a23);
    printf ("a31: %7.3f, a32: %7.3f a33: %7.3f\n", a31, a32, a33);
  }
  if (dev == 1)
    fprintf (stdout, "Tilt: %f Rotate: %f\n", theta_t, phi_r);

  /* First tilt the angular response */
  for (i = 0; i < numu; i++) {
    theta = acos (umu[i]);
    for (j = 0; j < nphi; j++) {
      rn = 1;
      xp = rn * cos (phi[j] * PI / 180) * sin (theta);
      yp = rn * sin (phi[j] * PI / 180) * sin (theta);
      zp = rn * cos (theta);
      if (dev == 1) {
        printf ("\nxp:  %7.3f, yp:  %7.3f zp:  %7.3f theta: %7.4f phi: %5.1f r: %6.3f umu: %f\n",
                xp,
                yp,
                zp,
                theta,
                phi[j],
                r,
                umu[i]);
      }

      if (dev == 1)
        fprintf (stdout, "Tilt: %f Rotate: %f\n", theta_t, phi_r);

      x = a11 * xp + a12 * yp + a13 * zp;
      y = a21 * xp + a22 * yp + a23 * zp;
      z = a31 * xp + a32 * yp + a33 * zp;

      rn = sqrt (x * x + y * y + z * z);
      if (rn > 0)
        theta_n = acos (z / rn);
      else
        theta_n = 0;

      if (rn > 0) {
        phi_n = acos (x / (rn * sin (theta_n)));
      } else {
        phi_n = 0;
      }

      if (dev == 1) {
        printf ("x: %7.3f, y: %7.3f z: %7.3f theta: %7.4f phi: %5.1f\n", x, y, z, theta * 180 / PI, phi[j]);
        printf ("rn: %7.3f, theta_n: %7.4f phi_n: %5.1f phi_r: %f\n", rn, theta_n * 180 / PI, phi_n * 180 / PI, phi_r);
      }

      if (dev == 1)
        fprintf (stdout, "Tilt: %f Rotate: %f\n", theta_t, phi_r);

      it = 0;
      ia = 1;

      while ((theta_n * 180 / PI - acos (umu[ia++]) * 180 / PI) < 0 && ia < numu)
        ;
      it = ia - 1; /* This one is fun:-) */

      /* Make sure we get the closest grid point */
      if (fabs (acos (umu[it]) * 180 / PI - theta_n * 180 / PI) > fabs (acos (umu[it - 1]) * 180 / PI - theta_n * 180 / PI))
        it = it - 1;

      jt        = 0;
      ja        = 0;
      delphirad = fabs (phi[1] - phi[0]) * PI / 180;
      while (fabs (phi_n - phi[ja++] * PI / 180) > delphirad) {
        if (dev == 1) {
          printf ("%d %d %d %f %f %f %f\n",
                  j,
                  ja,
                  jt,
                  phi_n,
                  phi[ja - 1] * PI / 180,
                  fabs (phi_n - phi[ja - 1] * PI / 180),
                  delphirad);
        }
      }
      jt                 = ja - 1;
      angrestilted[i][j] = angres[it][jt];
    }
  }

  /* Now we rotate */
  nt = (phi_r * PI / 180) / delphirad;
  if (jt >= nphi)
    jt -= nphi;
  for (i = 0; i < numu; i++) {
    theta = acos (umu[i]);
    for (j = 0; j < nphi; j++) {
      jt = j + nt;
      if (jt >= nphi)
        jt -= nphi;
      if (dev == 1 || dev == 10)
        fprintf (stdout,
                 "%3d %3d nt: %d %3d jt: %3d %3d, nphi: %d, phi_r: %5.1f %3.1f\n",
                 i,
                 j,
                 nt,
                 it,
                 jt,
                 ja - 1,
                 nphi,
                 phi_r,
                 angres[it][jt]);

      angres[i][j] = angrestilted[i][jt];
      if (dev == 1) {
        printf ("angrestilted: %7.4f %7.4f %7.4f\n", acos (umu[it]) * 180 / PI, phi[jt], angrestilted[i][j]);
      }
    }
  }

  if (pgm) {

    if ((pgmfile = fopen ("angres_tilted.pgm", "w")) == NULL) {
      fprintf (stderr, "Error opening angres_tilted.pgm for writing\n");
      return NULL;
    }

    fprintf (pgmfile, "P2\n");
    fprintf (pgmfile, "%d %d\n", nphi, numu);
    fprintf (pgmfile, "256\n");
    for (i = 0; i < numu; i++) {
      for (j = 0; j < nphi; j++) {
        k = (int)(angres[i][j] * 256);
        fprintf (pgmfile, "%d\n", k);
      }
    }
    fclose (pgmfile);
  }

  return angres;
}

double linpol (double x1, double x2, double y1, double y2, double x) {
  /* Linearly interpolate between two points to get wanted y value for x. */
  double y;
  double a = 0, b = 0;

  if (x2 > 0 && x1 > 0 && fabs (x2 - x1) < 0.0000001) {
    y = y1;
  } else {
    a = (y2 - y1) / (x2 - x1);
    b = y1 - a * x1;
    y = a * x + b;
  }
  return y;
}

double mc_func (double* k, size_t dim, void* params) {
  int     i1 = 0, i2 = 0, j1 = 0, j2 = 0;
  double  p1, p2, p3;
  double  umu = 0, phi = 0;
  mcinput lmcinp;

  lmcinp = *(mcinput*)params;

  umu = k[0];
  phi = k[1];
  i2  = 0;
  while (umu > lmcinp.umu[i2++])
    ;
  i1 = i2 - 1;
  if (i2 == 0) {
    i2 = 1;
    i1 = 0;
  }
  if (i2 >= lmcinp.numu) {
    i2 = lmcinp.numu - 1;
    i1 = i2 - 1;
  }

  j2 = 0;
  while (phi > lmcinp.phi[j2++])
    ;
  j1 = j2 - 1;
  if (j2 == 0) {
    j2 = 1;
    j1 = 0;
  }
  if (j2 >= lmcinp.nphi) {
    j2 = lmcinp.nphi - 1;
    j1 = j2 - 1;
  }

  p1 = linpol (lmcinp.phi[j1], lmcinp.phi[j2], lmcinp.prod[i1][j1], lmcinp.prod[i1][j2], phi);
  p2 = linpol (lmcinp.phi[j1], lmcinp.phi[j2], lmcinp.prod[i2][j1], lmcinp.prod[i2][j2], phi);

  p3 = linpol (lmcinp.umu[i1], lmcinp.umu[i2], p1, p2, umu);

  return p3;
}

double mc_funciso (double* k, size_t dim, void* params) {
  int     i1 = 0, i2 = 0, j1 = 0, j2 = 0;
  double  p1, p2, p3;
  double  umu = 0, phi = 0;
  mcinput lmcinp;

  lmcinp = *(mcinput*)params;

  umu = k[0];
  phi = k[1];
  i2  = 0;
  while (umu > lmcinp.umu[i2++])
    ;
  i1 = i2 - 1;
  if (i2 == 0) {
    i2 = 1;
    i1 = 0;
  }
  if (i2 >= lmcinp.numu) {
    i2 = lmcinp.numu - 1;
    i1 = i2 - 1;
  }

  j2 = 0;
  while (phi > lmcinp.phi[j2++])
    ;
  j1 = j2 - 1;
  if (j2 == 0) {
    j2 = 1;
    j1 = 0;
  }
  if (j2 >= lmcinp.nphi) {
    j2 = lmcinp.nphi - 1;
    j1 = j2 - 1;
  }

  p1 = linpol (lmcinp.phi[j1], lmcinp.phi[j2], lmcinp.prod[i1][j1], lmcinp.prod[i1][j2], phi);
  p2 = linpol (lmcinp.phi[j1], lmcinp.phi[j2], lmcinp.prod[i2][j1], lmcinp.prod[i2][j2], phi);

  p3 = linpol (lmcinp.umu[i1], lmcinp.umu[i2], p1, p2, umu);

  return p3;
}

int main (int argc, char** argv) {
  int    i = 0, j = 0, k = 0, status = 0, pgm = 0, dev = 0, ang = 0, incdir = 0, isotropic = 0;
  int    ia = 0, it = 0, ja = 0, jt = 0;
  double delphirad = 0, dir = 0;

  char programname[FILENAME_MAX]    = "";
  char angresfilename[FILENAME_MAX] = ""; /* File holding angular response      */
  char raddisfilename[FILENAME_MAX] = ""; /* File holding radiance distribution */

  double *phi = NULL, *rphi = NULL, *aumu = NULL, *rumu = NULL, *aangres = NULL;
  double  theta_t = 0, phi_r = 0, sza = 0, phi0 = 0, max = 0;
  double  dphi = 0, dtheta = 0, tmpprod = 0;
  double  uavgso;
  /* commented out by RB, because was not used but caused compiler warnings */
  /* double rfldir, rfldn, flup, uavgso, uavgdn, uavgup; */

  int arows = 0, amin_columns = 0, amax_columns = 0; /* Angular response file      */
  int rrows = 0, rmin_columns = 0, rmax_columns = 0; /* Radiance distribution file */
  int nphi = 0, numu = 0;

  mcinput mcinp;
  /* void *pmcinp; */
  mcinput mcinpiso;
  /* void *pmcinpiso; */

  double **angres = NULL, **avalue = NULL, **rvalue = NULL, **raddis = NULL, **prod = NULL;

  FILE* pgmfile = NULL;

  /* Monte Carlo integration stuff */
  double res = 0, err = 0;

  double xl[DIM_ANG_INT];
  double xu[DIM_ANG_INT];

  /*  size_t calls = 10; */
  size_t calls = 5000000;

#if HAVE_LIBGSL
  const gsl_rng_type* T;
  const gsl_rng_type* Tiso;
  gsl_rng*            r;
  gsl_rng*            riso;

  gsl_monte_function G    = {&mc_func, DIM_ANG_INT, &mcinp};
  gsl_monte_function Giso = {&mc_funciso, DIM_ANG_INT, &mcinpiso};

  gsl_monte_plain_state* s    = gsl_monte_plain_alloc (DIM_ANG_INT);
  gsl_monte_plain_state* siso = gsl_monte_plain_alloc (DIM_ANG_INT);
#endif
  /* commented out by RB, because was not used but caused compiler warnings */
  /* pmcinp    = &mcinp; */
  /* pmcinpiso = &mcinpiso; */

  /* get command line options  */
  status = get_options (argc,
                        argv,
                        programname,
                        angresfilename,
                        raddisfilename,
                        &pgm,
                        &dev,
                        &ang,
                        &calls,
                        &sza,
                        &phi0,
                        &theta_t,
                        &phi_r,
                        &incdir,
                        &isotropic);

  if (status != 0)
    return status;

#if !HAVE_LIBGSL
  fprintf (stderr, "%s was built without the gsl-library. Thus NO SUPPORT\n", programname);
  fprintf (stderr, "for integration of product of angular response and radiance field.\n");
#endif

  /* read angular response file  */
  status = ASCII_file2double (angresfilename, &arows, &amax_columns, &amin_columns, &avalue);

  if (status != 0) {
    fprintf (stderr, "error %d reading input file %s\n", status, angresfilename);
    return status;
  }

  if (amin_columns < 2) {
    fprintf (stderr, "error, input file must have at least two columns\n");
    return -1;
  }

  if (amin_columns != amax_columns) {
    fprintf (stderr, " WARNING! %s is not a rectangular matrix; only the first\n", angresfilename);
    fprintf (stderr, " %d columns will be used for the analysis\n", amin_columns);
  }
  numu = arows;

  /* Read radiance field file  */
  status = ASCII_file2double (raddisfilename, &rrows, &rmax_columns, &rmin_columns, &rvalue);

  if (status != 0) {
    fprintf (stderr, "error %d reading input file %s\n", status, raddisfilename);
    return status;
  }

  /* Get umu and phi from radiance field  */
  rphi = ASCII_row (rvalue, rmax_columns, 1);
  j    = 0;
  for (i = 0; i < rmax_columns; i++) {
    if (rphi[i] >= 0) {
      nphi++;
    }
  }
  phi = (double*)calloc (nphi, sizeof (double));
  j   = 0;
  for (i = 0; i < nphi; i++) {
    if (rphi[i] != NAN) { /* ??? is that correct? Not even NAN == NAN ??? */
      phi[j++] = rphi[i];
    }
  }
  rumu = ASCII_column (rvalue, rrows, 0);
  for (i = 0; i < rrows - 2; i++) {
    rumu[i] = rumu[i + 2];
  } /* Shift to skip two first lines */
  rumu[rrows - 2] = NAN;
  rumu[rrows - 1] = NAN;

  /* allocate memory for the radiance field  */
  status = ASCII_calloc_double (&raddis, numu, nphi);
  if (status != 0) {
    fprintf (stderr, "angres: allocating memory for produce\n");
    return status;
  }
  for (i = 0; i < numu; i++) {
    for (j = 0; j < nphi; j++) {
      raddis[i][j] = rvalue[i + 2][j + 2];
    }
  }
  /* commented out by RB, because was not used but caused compiler warnings */
  /* rfldir = rvalue[0][1]; */
  /* rfldn   = rvalue[0][2]; */
  /* flup   = rvalue[0][3]; */
  uavgso = rvalue[0][4];
  /* uavgdn = rvalue[0][5]; */
  /* uavgup = rvalue[0][6]; */

  aumu = ASCII_column (avalue, arows, 0);
  if (ang)
    for (i = 0; i < numu; i++)
      aumu[i] = cos (aumu[i] * PI / 180);

  aangres = ASCII_column (avalue, arows, 1);

  /* compare umu with similar from angular response file, differences should be small,  */
  /* otherwise interpolate angular response to radiance field.                          */
  for (i = 0; i < numu; i++)
    if (fabs (aumu[i] - rumu[i]) > 0.0001)
      fprintf (stderr, "umu %f %f %f\n", aumu[i], rumu[i], aumu[i] - rumu[i]);

  /* Make diffuse radiation field isotropic */
  if (isotropic) {
    /* Integrate (angular response)*radiance over all angles */

    /* Upper and lower integration bounds for each dimension */
    xl[0] = aumu[0];
    xu[0] = aumu[numu - 1];
    xl[1] = phi[0] * PI / 180;
    xu[1] = phi[nphi - 1] * PI / 180;

    mcinpiso.nphi = nphi;
    mcinpiso.numu = numu;
    mcinpiso.phi  = (double*)calloc (nphi, sizeof (double));
    mcinpiso.umu  = (double*)calloc (numu, sizeof (double));
    status        = ASCII_calloc_double (&mcinpiso.prod, numu, nphi);
    for (i = 0; i < nphi; i++)
      mcinpiso.phi[i] = phi[i] * PI / 180;
    for (i = 0; i < numu; i++)
      mcinpiso.umu[i] = aumu[i];
    for (i = 0; i < numu; i++) {
      for (j = 0; j < nphi; j++) {
        mcinpiso.prod[i][j] = raddis[i][j];
      }
    }

#if HAVE_LIBGSL
    gsl_rng_env_setup();

    Tiso = gsl_rng_default;
    riso = gsl_rng_alloc (Tiso);

    gsl_monte_plain_integrate (&Giso, xl, xu, DIM_ANG_INT, calls, riso, siso, &res, &err);
    gsl_monte_plain_free (siso);
#endif

    dphi    = phi[1] - phi[0];
    dtheta  = fabs (acos (aumu[1]) * 180 / PI - acos (aumu[0]) * 180 / PI);
    tmpprod = dtheta * dphi;

    for (i = 0; i < numu; i++) {
      for (j = 0; j < nphi; j++) {
        raddis[i][j] = res / tmpprod;
      }
    }
  }

  /* Tilt and rotate angular response */
  angres = tilt_and_rotate (aangres, nphi, phi, numu, aumu, pgm, theta_t, phi_r, dev);
  if (angres == NULL) {
    fprintf (stderr, "Error calling tilt_and_rotate()\n");
    return -1;
  }

  /* Multiply angular response and radiance field */
  /* allocate memory for the product              */

  status = ASCII_calloc_double (&prod, numu, nphi);
  if (status != 0) {
    fprintf (stderr, "angres: error allocating memory for product\n");
    return status;
  }
  for (i = 0; i < numu; i++) {
    for (j = 0; j < nphi; j++) {
      prod[i][j] = angres[i][j] * raddis[i][j];
    }
  }

  if (pgm) {

    if ((pgmfile = fopen ("prod.pgm", "w")) == NULL) {
      fprintf (stderr, "Error opening prod.pgm for writing\n");
      return -1;
    }

    fprintf (pgmfile, "P2\n");
    fprintf (pgmfile, "%d %d\n", nphi, numu);
    fprintf (pgmfile, "256\n");
    max = -9999;
    for (i = 0; i < numu; i++)
      for (j = 0; j < nphi; j++)
        if (prod[i][j] > max)
          max = prod[i][j];

    for (i = 0; i < numu; i++)
      for (j = 0; j < nphi; j++) {
        k = (int)(prod[i][j] * 256 / max);
        fprintf (pgmfile, "%d\n", k);
      }

    fclose (pgmfile);

    if ((pgmfile = fopen ("raddis.pgm", "w")) == NULL) {
      fprintf (stderr, "Error opening raddis.pgm for writing\n");
      return -1;
    }

    fprintf (pgmfile, "P2\n");
    fprintf (pgmfile, "%d %d\n", nphi, numu);
    fprintf (pgmfile, "256\n");
    max = -9999;
    for (i = 0; i < numu; i++)
      for (j = 0; j < nphi; j++)
        if (raddis[i][j] > max)
          max = raddis[i][j];

    for (i = 0; i < numu; i++)
      for (j = 0; j < nphi; j++) {
        k = (int)(raddis[i][j] * 256 / max);
        fprintf (pgmfile, "%d\n", k);
      }

    fclose (pgmfile);
  }

  /* Find angular response for direct beam contribution */
  it = 0;
  ia = 0;

  while (((180 - sza) - acos (aumu[ia++]) * 180 / PI) < 0 && ia < numu)
    ;

  it = ia - 1; /* This one is fun:-) */

  /* Make sure we get the closest grid point */
  if (fabs (acos (aumu[it]) * 180 / PI - (180 - sza)) > fabs (acos (aumu[it - 1]) * 180 / PI - (180 - sza)))
    it = it - 1;

  jt        = 0;
  ja        = 0;
  delphirad = fabs (phi[1] - phi[0]) * PI / 180;
  while (fabs (phi0 - phi[ja++]) * PI / 180 > delphirad)
    ;

  jt = ja - 1;

  if (incdir) {
    dir = angres[it][jt] * uavgso * 4 * PI;
  } else {
    dir = 0;
  }

  ASCII_free_double (avalue, arows);
  ASCII_free_double (rvalue, rrows);
  ASCII_free_double (angres, numu);

  /* Integrate (angular response)*radiance over all angles */

  /* Upper and lower integration bounds for each dimension */
  xl[0] = aumu[0];
  xu[0] = aumu[numu - 1];
  xl[1] = phi[0] * PI / 180;
  xu[1] = phi[nphi - 1] * PI / 180;

  mcinp.nphi = nphi;
  mcinp.numu = numu;
  mcinp.phi  = (double*)calloc (nphi, sizeof (double));
  mcinp.umu  = (double*)calloc (numu, sizeof (double));
  status     = ASCII_calloc_double (&mcinp.prod, numu, nphi);
  for (i = 0; i < nphi; i++)
    mcinp.phi[i] = phi[i] * PI / 180;
  for (i = 0; i < numu; i++)
    mcinp.umu[i] = aumu[i];
  for (i = 0; i < numu; i++) {
    for (j = 0; j < nphi; j++) {
      mcinp.prod[i][j] = prod[i][j];
    }
  }

#if HAVE_LIBGSL
  gsl_rng_env_setup();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  gsl_monte_plain_integrate (&G, xl, xu, DIM_ANG_INT, calls, r, s, &res, &err);
  gsl_monte_plain_free (s);
#endif
  /*  fprintf(stdout,"Result: %f,  error: %f\n", res, err); */
  fprintf (stdout, "%f %f %f\n", res, err, res + dir);

  ASCII_free_double (raddis, numu);

  return 0; /* if o.k. */
}


text:
#ifndef libcv_ctf
#define libcv_ctf

#include "Array.h"
#import "Ellipse.h"
#include "MRC.h"
#include "Image.h"
#include "util.h"
#include "cvtypes.h"
#include <fftw3.h>
#include <math.h>
#include "geometry.h"
#include <gsl/gsl_math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_multimin.h>

typedef struct CTFST {

	char * imgpath;

	f64 defocus_x;
	f64 defocus_y;
	f64 astig_angle;
	f64 amplitude;
	f64 volts;
	f64 lambda;
	f64 size;
	f64 apix;
	f64 scale;

	f64 * ctf1d;
	f64 * ctfsynth1d;
	f64 * background1d;
	f64 * envelope1d;
	f64 ctf1dscale;
	f64 ctf1dapix;
	u32 ctf1dsize;

	f64 * ctf2d;
	f64 * ctfsynth2d;
	f64 * background2d;
	f64 * envelope2d;
	f64 ctf2dapix;
	u32 ctf2drows;
	u32 ctf2dcols;

} * CTF;

@interface Array ( CTF_Functions )

-(id) ellipse1DAvg:(EllipseP)ellipse;

@end

void ctf2_calcv( f64 c[], f64 ctf[], u64 size );
f64 ctf_norm( f64 fit_data[], f64 ctf_p[], f64 ctf[], f64 norm[], u32 size );
f64 calculate_score(f64 c1[], f64 c2[], u32 lcut, u32 rcut );
ArrayP ctfNormalize( ArrayP fit_data, ArrayP ctf_params );
void estimateDefocus( ArrayP fit_data, ArrayP ctf_params );
void peakNormalize( f64 values[], u32 size );
ArrayP fitNoise( ArrayP fit_data );
ArrayP fitEnvelope( ArrayP fit_data );
f64 getTEMLambda( f64 volts );
f64 getTEMVoltage( f64 lambda );
f64 ctf_calc( f64 c[], f64 x );
void fit2DCTF( ArrayP image, f64 d1, f64 d2, f64 th, f64 apix, f64 kv, f64 cs, f64 ac );
f64 positionForPeak( f64 c[], u32 peak_pos );
f64 defocusForPeak( f64 c[], f64 peak_pos, u32 peak_num );
ArrayP generate2DCTF( f64 df1, f64 df2, f64 theta, u32 rows, u32 cols, f64 apix, f64 cs, f64 kv, f64 ac );
ArrayP g2DCTF( f64 df1, f64 df2, f64 theta, u32 rows, u32 cols, f64 apix, f64 cs, f64 kv, f64 ac );

#endif