text:
/*
 * Command-line numerical polynomial equation solver using the GNU Scientific Library (GSL).
 * GSL is available from: "http://www.gnu.org/software/gsl/".
 *
 * Copyright (C) 2008-2011 George Gesslein II

This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.

This library 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
Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA

The chief copyright holder can be contacted at gesslein@mathomatic.org, or
George Gesslein II, P.O. Box 224, Lansing, NY  14882-0224  USA.

 */

/*
 * To compile:

./compile.roots

 * or

gcc -O3 -Wall -o roots roots.c -lgsl -lgslcblas

 * This program nicely solves any degree polynomial
 * with real coefficients by calling the
 * GSL function gsl_poly_complex_solve().
 * This file is also useful for testing
 * and as an example of using the GSL from C.
 * The results are double precision floating point values
 * that are sometimes accurate to 14 digits.
 * Complex numbers are output if successful.
 * Here is an example of it solving a cubic polynomial:

$ roots 1 1 1 1
The 3 approximate floating point solutions of:
+x^3 +x^2 +x +1 = 0
are:

x = +0 +1*i
x = +0 -1*i
x = -1
$

Try this for a large amount of error:
roots -1 -1 1 1

Proof the GSL isn't perfect.

 */

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <float.h>
#include <math.h>
#include <errno.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_errno.h>

#define EPSILON	0.00000000000005	/* a good value for doubles */
#define HELP	1			/* display helpful text */

void	display_root(int i);

char	*prog_name = "roots";		/* name of this program */
double	*a, *z;				/* input and output arrays, respectively */
int	precision = DBL_DIG - 1;	/* display precision, it is not useful to set this higher than 15 */

void
usage(void)
{
  printf("\n%s version 1.0 - numerical polynomial equation solver\n", prog_name);
  printf("Uses the GNU Scientific Library (GSL).\n");
  printf("\nSolves polynomial = 0 when given all real coefficients of the polynomial.\n");
  printf("Double precision floating point math is used, accurate to about 14 digits.\n");
  printf("\nUsage: %s highest-power-coefficient ... constant-term\n", prog_name);
  printf("\nThe coefficients must be decimal, floating point, real numbers.\n");
  printf("For example, if 4 real numbers are given, there will be 3 complex number\n");
  printf("results or \"roots\" that are all valid solutions to polynomial = 0.\n");
}

int
main(int argc, char **argv)
{
  int i, highest_power;
  char *cp, *arg;
  gsl_poly_complex_workspace *w;

  if (argc <= 1) {
    fprintf(stderr, "%s: The polynomial coefficients must be specified on the command line.\n", prog_name);
    usage();
    exit(2);
  }

  highest_power = argc - 2;
  a = calloc(highest_power + 1, sizeof(double)); /* allocate real double input array */
  z = calloc(2 * highest_power, sizeof(double)); /* allocate complex double output array */

/* parse the command line into the coefficient array a[] */
  for (i = 0; i < argc - 1; i++) {
    arg = argv[argc-i-1];
    errno = 0;
    a[i] = strtod(arg, &cp);
    if (arg == cp || *cp) {
      fprintf(stderr, "%s: Argument \"%s\" is not a floating point number.\n", prog_name, arg);
      usage();
      exit(2);
    }
    if (errno) {
      fprintf(stderr, "%s: Argument \"%s\" is out of range.\n", prog_name, arg);
      exit(2);
    }
  }

#if	HELP
/* nicely display the actual polynomial equation we are solving */
  printf("The %d approximate floating point solutions of:\n", highest_power);
  for (i = highest_power; i >= 0; i--) {
    if (a[i]) {
      if (i && a[i] == 1.0) {
	printf("+x");
      } else {
        printf("%+.*g", precision, a[i]);
        if (i) {
	  printf("*x");
        }
      }
      if (i > 1) {
        printf("^%d", i);
      }
      printf(" ");
    }
  }
  printf("= 0\nare:\n\n");
#endif

/* solve the polynomial equation */
  w = gsl_poly_complex_workspace_alloc(highest_power + 1);
  if (gsl_poly_complex_solve(a, highest_power + 1, w, z) != GSL_SUCCESS) {
    fprintf(stderr, "%s: GSL approximation failed.\n", prog_name);
    exit(1);
  }
  gsl_poly_complex_workspace_free(w);

/* display all solutions */
  for (i = 0; i < highest_power; i++) {
#ifdef	EPSILON /* zero out relatively very small values (which are floating point error) */
    if (fabs(z[2*i] * EPSILON) > fabs(z[2*i+1]))
      z[2*i+1] = 0.0;
    else if (fabs(z[2*i+1] * EPSILON) > fabs(z[2*i]))
      z[2*i] = 0.0;
#endif
#if	HELP
    printf("x = ");
#endif
    display_root(i);
    printf("\n");
  }

  exit(0);
}

void
display_root(int i)
{
  printf("%+.*g", precision, z[2*i]);		/* output real part */
  if (z[2*i+1])
    printf(" %+.*g*i", precision, z[2*i+1]);	/* output imaginary part */
}


text:
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlin.h>

struct data {
  size_t n;
  double * y;
  double * sigma;
};

int
expb_f (const gsl_vector * x, void *params,
        gsl_vector * f)
{
  size_t n = ((struct data *)params)->n;
  double *y = ((struct data *)params)->y;
  double *sigma = ((struct data *) params)->sigma;

  double A = gsl_vector_get (x, 0);
  double lambda = gsl_vector_get (x, 1);
  double b = gsl_vector_get (x, 2);

  size_t i;

  for (i = 0; i < n; i++)
    {
      /* Model Yi = A * exp(-lambda * i) + b */
      double t = i;
      double Yi = A * exp (-lambda * t) + b;
      gsl_vector_set (f, i, (Yi - y[i])/sigma[i]);
    }

  return GSL_SUCCESS;
}

int
expb_df (const gsl_vector * x, void *params,
         gsl_matrix * J)
{
  size_t n = ((struct data *)params)->n;
  double *sigma = ((struct data *) params)->sigma;

  double A = gsl_vector_get (x, 0);
  double lambda = gsl_vector_get (x, 1);

  size_t i;

  for (i = 0; i < n; i++)
    {
      /* Jacobian matrix J(i,j) = dfi / dxj, */
      /* where fi = (Yi - yi)/sigma[i],      */
      /*       Yi = A * exp(-lambda * i) + b  */
      /* and the xj are the parameters (A,lambda,b) */
      double t = i;
      double s = sigma[i];
      double e = exp(-lambda * t);
      gsl_matrix_set (J, i, 0, e/s);
      gsl_matrix_set (J, i, 1, -t * A * e/s);
      gsl_matrix_set (J, i, 2, 1/s);

    }
  return GSL_SUCCESS;
}

int
expb_fdf (const gsl_vector * x, void *params,
          gsl_vector * f, gsl_matrix * J)
{
  expb_f (x, params, f);
  expb_df (x, params, J);

  return GSL_SUCCESS;
}

#define N 40

int
main (void)
{
  const gsl_multifit_fdfsolver_type *T;
  gsl_multifit_fdfsolver *s;

  int status;
  size_t i, iter = 0;

  const size_t n = N;
  const size_t p = 3;

  gsl_matrix *covar = gsl_matrix_alloc (p, p);

  double y[N], sigma[N];

  struct data d = { n, y, sigma};

  gsl_multifit_function_fdf f;

  double x_init[3] = { 1.0, 0.0, 0.0 };

  gsl_vector_view x = gsl_vector_view_array (x_init, p);

  const gsl_rng_type * type;
  gsl_rng * r;

  gsl_rng_env_setup();

  type = gsl_rng_default;
  r = gsl_rng_alloc (type);

  f.f = &expb_f;
  f.df = &expb_df;
  f.fdf = &expb_fdf;
  f.n = n;
  f.p = p;
  f.params = &d;

  /* This is the data to be fitted */

  for (i = 0; i < n; i++)
    {
      double t = i;
      y[i] = 1.0 + 5 * exp (-0.1 * t)
                 + gsl_ran_gaussian (r, 0.1);
      sigma[i] = 0.1;
      printf ("data: %d %g %g\n", i, y[i], sigma[i]);
    };

  T = gsl_multifit_fdfsolver_lmsder;
  s = gsl_multifit_fdfsolver_alloc (T, n, p);
  gsl_multifit_fdfsolver_set (s, &f, &x.vector);

  print_state (iter, s);

  do
    {
      iter++;
      status = gsl_multifit_fdfsolver_iterate (s);

      printf ("status = %s\n", gsl_strerror (status));

      print_state (iter, s);

      if (status)
        break;

      status = gsl_multifit_test_delta (s->dx, s->x,
                                        1e-4, 1e-4);
    }
  while (status == GSL_CONTINUE && iter < 500);

  gsl_multifit_covar (s->J, 0.0, covar);

  gsl_matrix_fprintf (stdout, covar, "%g");

#define FIT(i) gsl_vector_get(s->x, i)
#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))

  printf ("A      = %.5f +/- %.5f\n", FIT(0), ERR(0));
  printf ("lambda = %.5f +/- %.5f\n", FIT(1), ERR(1));
  printf ("b      = %.5f +/- %.5f\n", FIT(2), ERR(2));

  {
    double chi = gsl_blas_dnrm2(s->f);
    printf("chisq/dof = %g\n",  pow(chi, 2.0)/ (n - p));
  }

  printf ("status = %s\n", gsl_strerror (status));

  gsl_multifit_fdfsolver_free (s);
  return 0;
}

int
print_state (size_t iter, gsl_multifit_fdfsolver * s)
{
  printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f "
          "|f(x)| = %g\n",
          iter,
          gsl_vector_get (s->x, 0),
          gsl_vector_get (s->x, 1),
          gsl_vector_get (s->x, 2),
          gsl_blas_dnrm2 (s->f));
}

ext: c

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 <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <sys/time.h>
#include <gsl/gsl_cblas.h>

#define N 5000;

void mm_bruteforce_ijk(double *a, double *b, double *c, int I, int K, int J) {
    for(int i = 0; i < I; i++) {
        for(int j = 0; j < J; j++) {
            for(int k = 0; k < K; k++) {
                c[i * J + j] += a[i * K + k] * b[k * J + j];
            }
        }
    }
}

void mm_bruteforce_ikj(double *a, double *b, double *c, int I, int K, int J) {
    for(int i = 0; i < I; i++) {
        for(int k = 0; k < K; k++) {
            double dv = a[i * K + k];
            for(int j = 0; j < J; j++) {
                c[i * J + j] += dv * b[k * J + j];
            }
        }
    }
}

void mm_omp(double *a, double *b, double *c, int I, int K, int J) {
    #pragma omp parallel for
    for(int i = 0; i < I; i++) {
        for(int k = 0; k < K; k++) {
            register double dv = a[i * K + k];
            for(int j = 0; j < J; j++) {
                c[i * J + j] += dv * b[k * J + j];
            }
        }
    }
}

void mm_cblas_dgemm(double *a, double *b, double *c, int p, int q, int r) {
    int l = p;
    int m = q;
    int n = r;

    int lda = m;
    int ldb = n;
    int ldc = n;

    double alpha = 1.0;
    double beta = 0.0;

    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
            l, n, m, alpha, a, lda, b, ldb, beta, c, ldc);
}

void init_arange(double *mat, int a, int b) {
    for(int i = 0; i < (a*b); i++) {
        mat[i] = i + 1;
    }
}

void init_value(double *mat, int a, int b, double value) {
    for(int i = 0; i < (a*b); i++) {
        mat[i] = value;
    }
}

void init_zero(double *mat, int a, int b) {
    init_value(mat, a, b, 0);
}

void init_one(double *mat, int a, int b) {
    init_value(mat, a, b, 1);
}

void timer_start(struct timeval *pstv) {
    gettimeofday(pstv, NULL);
}

void timer_end(struct timeval *petv) {
    gettimeofday(petv, NULL);
}

void timer_print(struct timeval *pstv, struct timeval *petv) {
    time_t sec;
    suseconds_t usec;
    sec = petv->tv_sec - pstv->tv_sec;
    usec = petv->tv_usec - pstv->tv_usec;
    if(usec < 0) {
        sec--;
        usec += 1000000;
    }
    printf("elapsed time : %ld.%ld\n", sec, usec);
}

typedef void (*fptr_mm)(double *a, double *b, double *c, int I, int K, int J);

void check_etime_mm(double *a, double *b, double *c, int I, int K, int J, fptr_mm mm) {
    struct timeval stv;
    struct timeval etv;
    init_zero(c, I, J);
    timer_start(&stv);
    mm(a, b, c, I, K, J);
    timer_end(&etv);
    timer_print(&stv, &etv);
}

int main(int argc, char *argv[]) {
    int I, K, J;
    I = K = J = N;

    double *a = (double*)malloc(sizeof(double) * I * K);
    double *b = (double*)malloc(sizeof(double) * K * J);
    double *c = (double*)malloc(sizeof(double) * I * J);

    init_arange(a, I, K);
    init_arange(b, K, J);

    check_etime_mm(a, b, c, I, K, J, mm_bruteforce_ijk);
    check_etime_mm(a, b, c, I, K, J, mm_bruteforce_ikj);
    check_etime_mm(a, b, c, I, K, J, mm_omp);
    check_etime_mm(a, b, c, I, K, J, mm_cblas_dgemm);

    return 0;
}


text:
#ifndef _ALLVARS_H
#define _ALLVARS_H

#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <math.h>

#define PI ((M_PI))

#define max(a, b) ((a>=b?a:b))
#define min(a, b) ((a>=b?b:a))

// file names
extern char fname_param[100];
extern char fname_con[100], fname_line[100], fname_results[100];

extern FILE *fp_results;

// transfer function parameters
extern double tau_lim_low, tau_lim_up;
extern int nc, nc_lim_low, nc_lim_up, nc_best;
extern double *grid_tau;
extern const int ntau;
extern double *TF_tau, *TF;
// mcmc
extern const int ntheta_max;
extern int ntheta;
extern int nmcmc, nbuilt;
extern double *workspace;
extern int flag_detrend, flag_sim, flag_mcmc;
extern double *Cmat, *ICmat, *Smat, *Nmat, *USmat, *ASmat;
extern double *ICvmat, *Tmat1, *Tmat2, *Tmat3, *Tmat4;

extern double *cov_matrix;
extern double **theta_range, *theta_best, *theta_best_var;
extern double *theta_input, *sigma_input;
extern int *theta_fixed;

extern double *theta_best_con, *theta_best_var_con;

extern const gsl_rng_type * gsl_T;
extern const gsl_rng * gsl_r;

// light curve data
extern const int ndata_max;
extern int ncon_data, nline_data, nall_data;
extern double *Tcon_data, *Tline_data, *Fcon_data, *Fcerrs_data, *Fline_data, *Flerrs_data;
extern double *Fall_data;
extern double scale_con, scale_line;
extern double len_con, len_line, cad_con, cad_line;

extern int ncon, nline, nall;
extern double *Tcon, *Fcon, *Fcerrs, *Tline, *Fline, *Flerrs;

// error exit
extern char str_error_exit[200];

// mathematic functions
extern int *workspace_ipiv;
#endif

ext: h

max_stars_repo_name: LiyrAstroph/MICA

max_issues_repo_name: LiyrAstroph/MICA

max_forks_repo_name: LiyrAstroph/MICA


text:
#pragma once

#include <CesiumAsync/IAssetResponse.h>

#include <gsl/span>

#include <string>
#include <vector>

namespace CesiumCpp {

namespace CesiumNativeImpl {
/**
 * @brief Implementation of an IAssetResponse for data from a file.
 */
class FileAssetResponse : public CesiumAsync::IAssetResponse {
public:
  /**
   * @brief Creates a new instance.
   *
   * @param fileName The file name.
   * @param statusCode The HTTP status code
   * @param contentType The content type.
   * @param fileData The file data.
   */
  FileAssetResponse(const std::string &fileName, uint16_t statusCode,
                    const std::string &contentType,
                    const CesiumAsync::HttpHeaders &headers,
                    const std::vector<std::byte> &fileData);

  uint16_t statusCode() const noexcept override;

  virtual const CesiumAsync::HttpHeaders &headers() const override {
    return this->_headers;
  }

  std::string contentType() const noexcept override;

  gsl::span<const std::byte> data() const noexcept override;

private:
  std::string _fileName;
  uint16_t _statusCode;
  std::string _contentType;
  CesiumAsync::HttpHeaders _headers;

  std::vector<std::byte> _fileData;
};
} // namespace CesiumNativeImpl
} // namespace CesiumCpp

max_stars_repo_name: javagl/cesium-cpp

max_issues_repo_name: javagl/cesium-cpp

max_forks_repo_name: javagl/cesium-cpp


text:
#pragma once
#include <utility>
#include <algorithm>
#include <gsl/span>
#include <array>
#include <vector>
#include <memory>
#include <chrono>

struct patternbyte
{
    uint8_t value;
    uint8_t mask;

    patternbyte()
    {
        this->clear();
    }

    patternbyte(const uint8_t value, const uint8_t mask)
    {
        this->value = value;
        this->mask = mask;
    }

    bool opaque() const
    {
        return this->mask == 0xffui8;
    }

    void clear()
    {
        this->value = 0ui8;
        this->mask = 0xffui8;
    }

    friend inline bool operator==(
        const patternbyte& left,
        const uint8_t& right)
    {
        return (right & left.mask) == (left.value & left.mask); // value should already be masked
    }

    friend inline bool operator==(
        const uint8_t& left,
        const patternbyte& right)
    {
        return (left & right.mask) == (right.value & right.mask);
    }
};

inline int hexchtoint(char ch)
{
    if (ch >= '0' && ch <= '9')
        return ch - '0';
    else if (ch >= 'A' && ch <= 'F')
        return ch - 'A' + 10;
    else if (ch >= 'a' && ch <= 'f')
        return ch - 'a' + 10;
    return -1;
}

static auto compile_pattern(
    const char* pattern)
{
    auto compiled_pattern = std::make_unique<std::vector< patternbyte>>();
    patternbyte b;

    char c;
    auto shift = 0x4UL;
    while (c = *pattern++) {
        if (c == '?' || c == '.') {
            b.mask &= ~(0xf << shift);
        }
        else if (auto num = hexchtoint(c); num != -1) {
            b.value |= (num & 0xf) << shift;
        }
        else {
            continue;
        }

        if (!shift) {
            compiled_pattern->push_back(b);
            b.clear();
        }
        shift ^= 0x4UL;
    }
    if (!shift) {
        b.mask &= 0xf;
        compiled_pattern->push_back(b);
    }
    return compiled_pattern;
}

class pattern_searcher
{
    std::unique_ptr<std::vector< patternbyte>> _cp;

public:
    pattern_searcher(const char* pattern)
    {
        _cp = compile_pattern(pattern);
    }

    template <class Iterator>
    std::pair<Iterator, Iterator> operator()(Iterator first, Iterator last) const
    {
        if (first == last)
            return std::make_pair(last, last);
        if (_cp->empty())
            return std::make_pair(first, first);

        auto it = std::search(first, last, std::begin(*_cp), std::end(*_cp));

        if (it == last)
            return std::make_pair(last, last);
        return std::make_pair(it, it + std::size(*_cp));
    }
};

max_stars_repo_name: bnsmodpolice/loginhelper

max_issues_repo_name: bnsmodpolice/loginhelper

max_forks_repo_name: bnsmodpolice/loginhelper


text:
#ifndef TYPEDEFS
#define TYPEDEFS

#define MH_REPS 100
#define TRUE 1
#define FALSE 0

#include <gsl/gsl_vector.h>
#include "utils.h"

typedef struct topic
{
    double w_tot;                // total # of words assigned to the topic
    gsl_vector* w_cnt;           // vector of word counts
    gsl_vector* log_prob_w;      // vector of log probabilities
    gsl_vector* lgam_w_plus_eta; // to prevent many repeated computations

    double doc_tot;         // total number of doc-level instances
    double log_doc_tot;     // precomputed log of doc_tot
    int id;                 // a unique ID for this topic

    int level;              // level in the tree
    int nchild;             // number of children
    double scaling;         // scaling factor for my DP
    struct topic** child;   // array of pointers to child topics
    struct topic* parent;   // pointer to the parent topic
    struct tree* tr;        // pointer to my tree

    double prob;            // probability (used to sample a path)
} topic;

typedef struct tree
{
    int depth;              // depth of the tree
    gsl_vector* eta;        // topic dirichlet parameter
    gsl_vector* gam;        // scaling parameter (!!! not used; see G prior)
    double scaling_shape;   // shape parameter for the G prior on scaling
    double scaling_scale;   // scale parameter for the G prior on scaling
    topic* root;            // root topic of the tree
    int next_id;            // the next id for a new topic
} tree;

typedef struct doc
{
    int_vector* word;            // each word
    int_vector* levels;          // level assigned to each word

    int id;
    topic** path;                // path of topics
    gsl_vector* tot_levels;      // level counts (convenience, for sampling)
    gsl_vector* log_p_level;     // log p(level) [ unnormalized ]
    double* gem_mean;
    double* gem_scale;
    double score;
} doc;

typedef struct corpus
{
    double gem_mean;
    double gem_scale;
    int ndoc;
    int nterms;
    doc** doc;
} corpus;

typedef struct gibbs_state
{
    // data and hidden variables
    corpus* corp;
    tree* tr;

    // current scores and iteration
    double score;
    double gem_score;
    double eta_score;
    double gamma_score;
    double max_score;
    int iter;

    // log files
    char* run_dir;
    FILE* score_log;
    FILE* tree_structure_log;

    // sampling parameters
    int shuffle_lag;
    int hyper_lag;
    int level_lag;
    int output_lag;
    int sample_eta;
    int sample_gem;
    int sample_gam;
} gibbs_state;

#endif

max_stars_repo_name: thatindiandude/PGM-Project

max_issues_repo_name: thatindiandude/PGM-Project

max_forks_repo_name: thatindiandude/PGM-Project


text:
/* ode-initval/rkf45.c
 *
 * Copyright (C) 2001 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 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Runge-Kutta-Fehlberg 4(5)*/

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include "gsl_odeiv.h"

#include "odeiv_util.h"

/* Runge-Kutta-Fehlberg constants */
static const double ah[] = { 1.0/4.0, 3.0/8.0, 12.0/13.0, 1.0, 1.0/2.0 };

static const double b3[] = { 3.0/32.0, 9.0/32.0 };

static const double b4[] = { 1932.0/2197.0, -7200.0/2197.0, 7296.0/2197.0};

static const double b5[] = { 8341.0/4104.0, -32832.0/4104.0, 29440.0/4104.0, -845.0/4104.0};

static const double b6[] = { -6080.0/20520.0, 41040.0/20520.0, -28352.0/20520.0, 9295.0/20520.0, -5643.0/20520.0};

static const double c1 = 902880.0/7618050.0;
static const double c3 = 3953664.0/7618050.0;
static const double c4 = 3855735.0/7618050.0;
static const double c5 = -1371249.0/7618050.0;
static const double c6 = 277020.0/7618050.0;

static const double ec[] = { 0.0,
  1.0 / 360.0,
  0.0,
  -128.0 / 4275.0,
  -2197.0 / 75240.0,
  1.0 / 50.0,
  2.0 / 55.0
};

typedef struct
{
  double *k1;
  double *k2;
  double *k3;
  double *k4;
  double *k5;
  double *k6;
  double *y0;
  double *ytmp;
}
rkf45_state_t;

static void *
rkf45_alloc (size_t dim)
{
  rkf45_state_t *state = (rkf45_state_t *) malloc (sizeof (rkf45_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for rkf45_state", GSL_ENOMEM);
    }

  state->k1 = (double *) malloc (dim * sizeof (double));

  if (state->k1 == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k1", GSL_ENOMEM);
    }

  state->k2 = (double *) malloc (dim * sizeof (double));

  if (state->k2 == 0)
    {
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k2", GSL_ENOMEM);
    }

  state->k3 = (double *) malloc (dim * sizeof (double));

  if (state->k3 == 0)
    {
      free (state->k2);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k3", GSL_ENOMEM);
    }

  state->k4 = (double *) malloc (dim * sizeof (double));

  if (state->k4 == 0)
    {
      free (state->k3);
      free (state->k2);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k4", GSL_ENOMEM);
    }

  state->k5 = (double *) malloc (dim * sizeof (double));

  if (state->k5 == 0)
    {
      free (state->k4);
      free (state->k3);
      free (state->k2);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k5", GSL_ENOMEM);
    }

  state->k6 = (double *) malloc (dim * sizeof (double));

  if (state->k6 == 0)
    {
      free (state->k5);
      free (state->k4);
      free (state->k3);
      free (state->k2);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k6", GSL_ENOMEM);
    }

  state->y0 = (double *) malloc (dim * sizeof (double));

  if (state->y0 == 0)
    {
      free (state->k6);
      free (state->k5);
      free (state->k4);
      free (state->k3);
      free (state->k2);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
    }

  state->ytmp = (double *) malloc (dim * sizeof (double));

  if (state->ytmp == 0)
    {
      free (state->y0);
      free (state->k6);
      free (state->k5);
      free (state->k4);
      free (state->k3);
      free (state->k2);
      free (state->k1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
    }

  return state;
}

static int
rkf45_apply (void *vstate,
	    size_t dim,
	    double t,
	    double h,
	    double y[],
	    double yerr[],
	    const double dydt_in[],
	    double dydt_out[], const gsl_odeiv_system * sys)
{
  rkf45_state_t *state = (rkf45_state_t *) vstate;

  size_t i;
  int status = 0;

  double *const k1 = state->k1;
  double *const k2 = state->k2;
  double *const k3 = state->k3;
  double *const k4 = state->k4;
  double *const k5 = state->k5;
  double *const k6 = state->k6;
  double *const ytmp = state->ytmp;

  /* k1 step */
  if (dydt_in != NULL)
    {
      DBL_MEMCPY (k1, dydt_in, dim);
    }
  else
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t, y, k1);
      GSL_STATUS_UPDATE (&status, s);
    }

  for (i = 0; i < dim; i++)
    ytmp[i] = y[i] +  ah[0] * h * k1[i];

  /* k2 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[0] * h, ytmp, k2);
    GSL_STATUS_UPDATE (&status, s);
  }
  for (i = 0; i < dim; i++)
    ytmp[i] = y[i] + h * (b3[0] * k1[i] + b3[1] * k2[i]);

  /* k3 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[1] * h, ytmp, k3);
    GSL_STATUS_UPDATE (&status, s);
  }
  for (i = 0; i < dim; i++)
    ytmp[i] = y[i] + h * (b4[0] * k1[i] + b4[1] * k2[i] + b4[2] * k3[i]);

  /* k4 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[2] * h, ytmp, k4);
    GSL_STATUS_UPDATE (&status, s);
  }
  for (i = 0; i < dim; i++)
    ytmp[i] =
      y[i] + h * (b5[0] * k1[i] + b5[1] * k2[i] + b5[2] * k3[i] +
		  b5[3] * k4[i]);

  /* k5 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[3] * h, ytmp, k5);
    GSL_STATUS_UPDATE (&status, s);
  }
  for (i = 0; i < dim; i++)
    ytmp[i] =
      y[i] + h * (b6[0] * k1[i] + b6[1] * k2[i] + b6[2] * k3[i] +
		  b6[3] * k4[i] + b6[4] * k5[i]);

  /* k6 step and final sum */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[4] * h, ytmp, k6);
    GSL_STATUS_UPDATE (&status, s);
  }

  for (i = 0; i < dim; i++)
    {
      const double d_i = c1 * k1[i] + c3 * k3[i] + c4 * k4[i] + c5 * k5[i] + c6 * k6[i];
      y[i] += h * d_i;
      if (dydt_out != NULL)
	dydt_out[i] = d_i;
    }

  /* difference between 4th and 5th order */
  for (i = 0; i < dim; i++)
    yerr[i] =
      h * (ec[1] * k1[i] + ec[3] * k3[i] + ec[4] * k4[i] + ec[5] * k5[i] +
	   ec[6] * k6[i]);

  return status;
}

static int
rkf45_reset (void *vstate, size_t dim)
{
  rkf45_state_t *state = (rkf45_state_t *) vstate;

  DBL_ZERO_MEMSET (state->k1, dim);
  DBL_ZERO_MEMSET (state->k2, dim);
  DBL_ZERO_MEMSET (state->k3, dim);
  DBL_ZERO_MEMSET (state->k4, dim);
  DBL_ZERO_MEMSET (state->k5, dim);
  DBL_ZERO_MEMSET (state->k6, dim);
  DBL_ZERO_MEMSET (state->ytmp, dim);

  return GSL_SUCCESS;
}

static unsigned int
rkf45_order (void *vstate)
{
  rkf45_state_t *state = (rkf45_state_t *) vstate;
  state = 0; /* prevent warnings about unused parameters */
  return 5;
}

static void
rkf45_free (void *vstate)
{
  rkf45_state_t *state = (rkf45_state_t *) vstate;

  free (state->ytmp);
  free (state->y0);
  free (state->k6);
  free (state->k5);
  free (state->k4);
  free (state->k3);
  free (state->k2);
  free (state->k1);
  free (state);
}

static const gsl_odeiv_step_type rkf45_type = { "rkf45",	/* name */
  1,				/* can use dydt_in */
  0,				/* gives exact dydt_out */
  &rkf45_alloc,
  &rkf45_apply,
  &rkf45_reset,
  &rkf45_order,
  &rkf45_free
};

const gsl_odeiv_step_type *gsl_odeiv_step_rkf45 = &rkf45_type;


text:
#include <gsl/gsl_sf_hyperg.h>

double potential(double x) { return pow(x, 4) - pow(x, 2) + 1; }
double ground_state_energy() { return 0.0; }
double exact(double x) { return 0.0; }

// Works with 70 70 70 70 & 1000 points of int


text:
//! @file
//! Various functions used during experiences (dgemmsy)
//!
//! @author
//!    Copyright (C) 2009-2011 BigDFT group
//!    This file is distributed under the terms of the
//!    GNU General Public License, see ~/COPYING file
//!    or http://www.gnu.org/copyleft/gpl.txt .
//!    For the list of contributors, see ~/AUTHORS
//! Eric Bainville, Mar 2010

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <malloc.h>
#include <math.h>
#include <unistd.h>
#include <xmmintrin.h>

// Set to 1 to use BLAS axpy call to accumulate partial results.
#ifndef USE_AXPY
#define USE_AXPY 1
/*#define ACML_BLAS 1*/
#endif

#if USE_AXPY
#if defined(MKL_BLAS) || defined(WIN32)
#include <mkl_cblas.h>
#else
#if defined(ACML_BLAS)
#include <acml.h>
#define cblas_daxpy daxpy
#else
#include <cblas.h>
#endif
#endif
#endif

#include "utils.h"

// System-specific functions

#include <unistd.h>
void * allocPages(size_t sz)
{
  size_t page_size = sysconf(_SC_PAGESIZE);
  size_t extra = sz % page_size;
  if (extra != 0) sz = (sz - extra) + page_size; // round to next page
  void * address = 0;
  int status = posix_memalign(&address,page_size,sz);
  if (status != 0) return 0; // Failed
  return address;
}
void freePages(void * address)
{
  free(address);
}
double * allocDoubles(int n)
{
  void * address = 0;
  int status = posix_memalign(&address,16,n*sizeof(double));
  if (status != 0) return 0; // Failed
  return address;
}
float * allocFloats(int n)
{
  void * address = 0;
  int status = posix_memalign(&address,16,n*sizeof(float));
  if (status != 0) return 0; // Failed
  return address;
}
void freeDoubles(double * x)
{
  free(x);
}
void freeFloats(float * x)
{
  free(x);
}

// Reference implementations

void gemm_block_2x2_ref(const double * a,const double * b,long long int n,double * y,long long int ldy)
{
  int k;
  for (k=0;k<n;k+=2)
  {
    y[0] += a[2*k] * b[4*k];
    y[0] += a[2*k+1] * b[4*k+1];
    y[1] += a[2*k] * b[4*k+2];
    y[1] += a[2*k+1] * b[4*k+3];
    y[ldy] += a[2*k+2] * b[4*k];
    y[ldy] += a[2*k+3] * b[4*k+1];
    y[ldy+1] += a[2*k+2] * b[4*k+2];
    y[ldy+1] += a[2*k+3] * b[4*k+3];
  }
}

void gemm_block_2x4_ref(const double * a,const double * b,long long int n,double * y,long long int ldy)
{
  int k;
  double * work;
  for (k=0;k<n;k+=2)
  {
    work = y;
    work[0] += a[2*k]   * b[4*k];
    work[0] += a[2*k+1] * b[4*k+1];
    work += ldy;
    work[0] += a[2*k]   * b[4*k+2];
    work[0] += a[2*k+1] * b[4*k+3];
    work += ldy;
    work[0] += a[2*k]   * b[4*k+4];
    work[0] += a[2*k+1] * b[4*k+5];
    work += ldy;
    work[0] += a[2*k]   * b[4*k+6];
    work[0] += a[2*k+1] * b[4*k+7];

    work = y + 1;
    work[0] += a[2*k+2] * b[4*k];
    work[0] += a[2*k+3] * b[4*k+1];
    work += ldy;
    work[0] += a[2*k+2] * b[4*k+2];
    work[0] += a[2*k+3] * b[4*k+3];
    work += ldy;
    work[0] += a[2*k+2] * b[4*k+4];
    work[0] += a[2*k+3] * b[4*k+5];
    work += ldy;
    work[0] += a[2*k+2] * b[4*k+6];
    work[0] += a[2*k+3] * b[4*k+7];
  }
}

int dgemmsy_base(dgemmsyBaseArgs * args)
{
  int status = 0; // return value
  int row,col;
  size_t slice_size,ywork_size;
  int p2; // P rounded to next multiple of 4
  double *a_slice,*b_slice,*y_work;
  ComputeData cdata;
  TransposeData tdata;
  int slice_n,n,p;
  const double *a,*b;
  double *y;
  int lda,ldb,ldy;
  int transa,transb;
  BlockPattern_2x4_Proc pattern;
  void * pattern_arg;
  double alpha;

  if (args == 0) return -2;

  pattern = args->params.pattern;
  pattern_arg = args->params.pattern_arg;
  slice_n = args->params.slice_n;

  transa = args->transa;
  transb = args->transb;
  n = args->n;
  p = args->p;
  a = args->a;
  lda = args->lda;
  b = args->b;
  ldb = args->ldb;
  y = args->y;
  ldy = args->ldy;
  alpha = args->alpha;

  // Check dimensions
  if (slice_n < 8) return -1;
  if ( n <= 0 || p <= 0 ) return -1;
  p2 = (p+3) & 0x7FFFFFFC;

  // printf("N=%d P=%d LDA=%d LDB=%d LDY=%d  P2=%d\n",n,p,lda,ldb,ldy,p2);

  // Check other arguments
  if (a == 0 || b == 0 || y == 0 || pattern == 0) return -5;

  // Allocate memory
  slice_size = slice_n * p2 * sizeof(double); // Slice size in bytes
  ywork_size = p2 * p2 * sizeof(double); // Result size in bytes
  a_slice = (double *)allocPages(slice_size);
  b_slice = (double *)allocPages(slice_size);
  y_work = (double *)allocPages(ywork_size);
  if (a_slice == 0 || b_slice == 0 || y_work == 0) { status = -3; goto END; }
  memset(y_work,0,ywork_size);

  // Loop on all slices and accumulate products
  initComputeData(&cdata,p2,slice_n,a_slice,b_slice,y_work);
  for (row=0;row<n;row+=slice_n)
  {
    int s = slice_n;
    if (row+s > n) s = n-row; // Limit size of last slice if needed

    // 2-pack A slice
    if (transa) tpack_2(s,p,a+lda*row,lda, slice_n,p2,a_slice);
    else npack_2(s,p,a+row,lda, slice_n,p2,a_slice);

    // 4-pack B slice
    if (transb) tpack_4(s,p,b+ldb*row,ldb, slice_n,p2,b_slice);
    else npack_4(s,p,b+row,ldb, slice_n,p2,b_slice);

    pattern(pattern_arg,p2,Compute_visitor,&cdata);
  }
  cleanupComputeData(&cdata);

  // Complete result by symmetry
  initTransposeData(&tdata,p2,y_work);
  pattern(pattern_arg,p2,Transpose_visitor,&tdata);
  cleanupTransposeData(&tdata);

  // Combine and store (untransposed) result. If we are multithreading,
  // we must protect the update with a mutex, since all threads
  // will update the same Y.
  if (args->yMutex != 0)
    {
      int locked = pthread_mutex_lock(args->yMutex);
      if (locked != 0) status = -4;
    }
  for (col=0;col<p;col++)
    {
      double * yy = y+ldy*col;
      double * yy_work = y_work+p2*col;
#if USE_AXPY
      cblas_daxpy(p,alpha,yy_work,1,yy,1);
#else
      if (alpha == 1)
	{
	  for (row=0;row<p;row++) yy[row] += yy_work[row];
	}
      else
	{
	  for (row=0;row<p;row++) yy[row] += alpha * yy_work[row];
	}
#endif
    }
  if (args->yMutex != 0)
    {
      int unlocked = pthread_mutex_unlock(args->yMutex);
      if (unlocked != 0) status = -4;
    }

END:
  // Cleanup
  freePages(y_work);
  freePages(a_slice);
  freePages(b_slice);
  return status;
}

void npack_2(int n,int p,const double * a,int lda,int s,int q,double * slice)
{
  int row,col,icol;
  double * x = 0;
  int oddn = 0;
  if (n&1) { n--; oddn=1; }
  for (col=0;col<q;col++)
    {
      icol = col&1;
      x = slice + s*(col-icol) + (icol<<1);
      row = 0;
      if (col < p)
	{
	  for ( ;row<n;row+=2)
	    {
	      x[2*row] = a[row];
	      x[2*row+1] = a[row+1];
	    }
	  if (oddn) { x[2*row] = a[row]; x[2*row+1] = 0; row+=2; }
	}
      for ( ;row<s;row+=2)
	{
	  x[2*row] = 0;
	  x[2*row+1] = 0;
	}
      a += lda;
    }
}

void npack_4(int n,int p,const double * a,int lda,int s,int q,double * slice)
{
  int row,col,icol;
  double * x = 0;
  int oddn = 0;
  if (n&1) { n--; oddn=1; }
  for (col=0;col<q;col++)
    {
      icol = col&3;
      x = slice + s*(col-icol) + (icol<<1);
      row = 0;
      if (col < p)
	{
	  for ( ;row<n;row+=2)
	    {
	      x[4*row] = a[row];
	      x[4*row+1] = a[row+1];
	    }
	  if (oddn) { x[4*row] = a[row]; x[4*row+1] = 0; row+=2; }
	}
      for ( ;row<s;row+=2)
	{
	  x[4*row] = 0;
	  x[4*row+1] = 0;
	}
      a += lda;
    }
}

void tpack_2(int n,int p,const double * a,int lda,int s,int q,double * slice)
{
  int row,col,icol;
  double * x = 0;
  const double * aa;
  int da = lda<<1;
  int oddn = 0;
  if (n&1) { n--; oddn=1; }
  for (col=0;col<q;col++)
    {
      icol = col&1;
      x = slice + s*(col-icol) + (icol<<1);
      row = 0;
      if (col < p)
	{
	  aa = a;
	  for ( ;row<n;row+=2)
	    {
	      x[2*row] = aa[0];
	      x[2*row+1] = aa[lda];
	      aa += da;
	    }
	  if (oddn) { x[2*row] = aa[0]; x[2*row+1] = 0; row += 2; }
	}
      for ( ;row<s;row+=2)
	{
	  x[2*row] = 0;
	  x[2*row+1] = 0;
	}
      a++;
    }
}

void tpack_4(int n,int p,const double * a,int lda,int s,int q,double * slice)
{
  int row,col,icol;
  double * x = 0;
  const double * aa;
  int da = lda<<1;
  int oddn = 0;
  if (n&1) { n--; oddn=1; }
  for (col=0;col<q;col++)
    {
      icol = col&3;
      x = slice + s*(col-icol) + (icol<<1);
      row = 0;
      if (col < p)
	{
	  aa = a;
	  for ( ;row<n;row+=2)
	    {
	      x[4*row] = aa[0];
	      x[4*row+1] = aa[lda];
	      aa += da;
	    }
	  if (oddn) { x[4*row] = aa[0]; x[4*row+1] = 0; row += 2; }
	}
      for ( ;row<s;row+=2)
	{
	  x[4*row] = 0;
	  x[4*row+1] = 0;
	}
      a++;
    }
}

max_stars_repo_name: gimunu/octopus-metric

max_issues_repo_name: gimunu/octopus-metric

max_forks_repo_name: gimunu/octopus-metric


text:
//This does CELL (~soma) stage of minimal GRU (gated recurrent unit) model.
//This requires each neuron to have 2 input time series, X and Xf,
//where X is the usual input and Xf the input for the forget-gate.
//Both X and Xf are the output of a linear IN stage (weights and baises).

//For dim=0: F[:,t] = sig(Xf[:,t] + Uf*Y[:,t-1])
//           H[:,t] = F[:,t].*Y[:,t-1]
//           Y[:,t] = H[:,t] + (1-f[:,t]).*tanh(X[:,t] + U*H[:,t])
//
//For dim=1: F[t,:] = sig(Xf[t,:] + Y[t-1,:]*Uf)
//           H[t,:] = F[t,:].*Y[t-1,:]
//           Y[t,:] = H[t,:] + (1-F[t,:]).*tanh(X[t,:] + H[t,:]*U)
//
//where sig is the logistic (sigmoid) nonlinearity = 1/(1+exp(-x)),
//F is the forget gate signal, H is an intermediate vector,
//U and Uf are NxN update matrices, and Y is the output.

//Note that, the neurons of a layer are independent only if U and Uf are diagonal matrices.
//This is only really a CELL (~soma) stage in that case.

//Unlike v2, this version takes only one X input: X = [X; Xf] for dim=0
//                                                X = [X Xf]  for dim=1
//This is essentially to allow treatment as a layer with Linux piping conventions.

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

#ifdef __cplusplus
namespace codee {
extern "C" {
#endif

int gru_min_s (float *Y, const float *X, const float *U, const float *Uf, const size_t N, const size_t T, const char iscolmajor, const size_t dim);
int gru_min_d (double *Y, const double *X, const double *U, const double *Uf, const size_t N, const size_t T, const char iscolmajor, const size_t dim);

int gru_min_inplace_s (float *X, const float *U, const float *Uf, const size_t N, const size_t T, const char iscolmajor, const size_t dim);
int gru_min_inplace_d (double *X, const double *U, const double *Uf, const size_t N, const size_t T, const char iscolmajor, const size_t dim);

int gru_min_s (float *Y, const float *X, const float *U, const float *Uf, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const float o = 1.0f;
    const size_t NT = N*T;
    size_t nT, tN, tN2;

    float *F, *H;
    if (!(F=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in gru_min_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in gru_min_s: problem with malloc. "); perror("malloc"); return 1; }

    if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                F[n] = 1.0f / (1.0f+expf(-X[N+n]));
                Y[n] = (1.0f-F[n]) * tanhf(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN2 = 2*tN;
                cblas_scopy((int)N,&X[tN2+N],1,F,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&Y[tN-N],1,o,F,1);
                for (size_t n=0u; n<N; ++n)
                {
                    F[n] = 1.0f / (1.0f+expf(-F[n]));
                    H[n] = F[n] * Y[tN-N+n];
                }
                cblas_scopy((int)N,&X[tN2],1,&Y[tN],1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[tN],1);
                for (size_t n=0u; n<N; ++n)
                {
                    Y[tN+n] = H[n] + (1.0f-F[n])*tanhf(Y[tN+n]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                F[n] = 1.0f / (1.0f+expf(-X[NT+nT]));
                Y[nT] = (1.0f-F[n]) * tanhf(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&X[NT+t],(int)T,F,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&Y[t-1],(int)T,o,F,1);
                for (size_t n=0u; n<N; ++n)
                {
                    F[n] = 1.0f / (1.0f+expf(-F[n]));
                    H[n] = F[n] * Y[t-1+n*T];
                }
                cblas_scopy((int)N,&X[t],(int)T,&Y[t],(int)T);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    Y[t+nT] = H[n] + (1.0f-F[n])*tanhf(Y[t+nT]);
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                F[n] = 1.0f / (1.0f+expf(-X[NT+nT]));
                Y[nT] = (1.0f-F[n]) * tanhf(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&X[NT+t],(int)T,F,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&Y[t-1],(int)T,o,F,1);
                for (size_t n=0u; n<N; ++n)
                {
                    F[n] = 1.0f / (1.0f+expf(-F[n]));
                    H[n] = F[n] * Y[t-1+n*T];
                }
                cblas_scopy((int)N,&X[t],(int)T,&Y[t],(int)T);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    Y[t+nT] = H[n] + (1.0f-F[n])*tanhf(Y[t+nT]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                F[n] = 1.0f / (1.0f+expf(-X[N+n]));
                Y[n] = (1.0f-F[n]) * tanhf(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN2 = 2*tN;
                cblas_scopy((int)N,&X[tN2+N],1,F,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&Y[tN-N],1,o,F,1);
                for (size_t n=0u; n<N; ++n)
                {
                    F[n] = 1.0f / (1.0f+expf(-F[n]));
                    H[n] = F[n] * Y[tN-N+n];
                }
                cblas_scopy((int)N,&X[tN2],1,&Y[tN],1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[tN],1);
                for (size_t n=0u; n<N; ++n)
                {
                    Y[tN+n] = H[n] + (1.0f-F[n])*tanhf(Y[tN+n]);
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in gru_min_s: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int gru_min_d (double *Y, const double *X, const double *U, const double *Uf, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const double o = 1.0;
    const size_t NT = N*T;
    size_t nT, tN, tN2;

    double *F, *H;
    if (!(F=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in gru_min_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in gru_min_d: problem with malloc. "); perror("malloc"); return 1; }

    if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                F[n] = 1.0 / (1.0+exp(-X[N+n]));
                Y[n] = (1.0-F[n]) * tanh(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN2 = 2*tN;
                cblas_dcopy((int)N,&X[tN2+N],1,F,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&Y[tN-N],1,o,F,1);
                for (size_t n=0u; n<N; ++n)
                {
                    F[n] = 1.0 / (1.0+exp(-F[n]));
                    H[n] = F[n] * Y[tN-N+n];
                }
                cblas_dcopy((int)N,&X[tN2],1,&Y[tN],1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[tN],1);
                for (size_t n=0u; n<N; ++n)
                {
                    Y[tN+n] = H[n] + (1.0-F[n])*tanh(Y[tN+n]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                F[n] = 1.0 / (1.0+exp(-X[NT+nT]));
                Y[nT] = (1.0-F[n]) * tanh(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&X[NT+t],(int)T,F,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&Y[t-1],(int)T,o,F,1);
                for (size_t n=0u; n<N; ++n)
                {
                    F[n] = 1.0 / (1.0+exp(-F[n]));
                    H[n] = F[n] * Y[t-1+n*T];
                }
                cblas_dcopy((int)N,&X[t],(int)T,&Y[t],(int)T);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    Y[t+nT] = H[n] + (1.0-F[n])*tanh(Y[t+nT]);
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                F[n] = 1.0 / (1.0+exp(-X[NT+nT]));
                Y[nT] = (1.0-F[n]) * tanh(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&X[NT+t],(int)T,F,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&Y[t-1],(int)T,o,F,1);
                for (size_t n=0u; n<N; ++n)
                {
                    F[n] = 1.0 / (1.0+exp(-F[n]));
                    H[n] = F[n] * Y[t-1+n*T];
                }
                cblas_dcopy((int)N,&X[t],(int)T,&Y[t],(int)T);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    Y[t+nT] = H[n] + (1.0-F[n])*tanh(Y[t+nT]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                F[n] = 1.0 / (1.0+exp(-X[N+n]));
                Y[n] = (1.0-F[n]) * tanh(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN2 = 2*tN;
                cblas_dcopy((int)N,&X[tN2+N],1,F,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&Y[tN-N],1,o,F,1);
                for (size_t n=0u; n<N; ++n)
                {
                    F[n] = 1.0 / (1.0+exp(-F[n]));
                    H[n] = F[n] * Y[tN-N+n];
                }
                cblas_dcopy((int)N,&X[tN2],1,&Y[tN],1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&Y[tN],1);
                for (size_t n=0u; n<N; ++n)
                {
                    Y[tN+n] = H[n] + (1.0-F[n])*tanh(Y[tN+n]);
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in gru_min_d: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int gru_min_inplace_s (float *X, const float *U, const float *Uf, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const float o = 1.0f;
    const size_t NT = N*T, N2 = 2*N;
    size_t nT, tN, tN2;

    float *F, *H;
    if (!(F=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in gru_inplace_min2_s: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(float *)malloc(N*sizeof(float)))) { fprintf(stderr,"error in gru_inplace_min2_s: problem with malloc. "); perror("malloc"); return 1; }

    if (dim==0u)
    {
        if (iscolmajor)
        {
            fprintf(stderr,"d0, colmajor \n");
            for (size_t n=0u; n<N; ++n)
            {
                F[n] = 1.0f / (1.0f+expf(-X[N+n]));
                X[n] = (1.0f-F[n]) * tanhf(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN2 = 2*tN;
                cblas_scopy((int)N,&X[tN2+N],1,F,1);
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&X[tN2-N2],1,o,F,1);
                for (size_t n=0u; n<N; ++n)
                {
                    F[n] = 1.0f / (1.0f+expf(-F[n]));
                    H[n] = F[n]*X[tN2-N2+n];
                }
                cblas_sgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[tN2],1);
                for (size_t n=0u; n<N; ++n)
                {
                    X[tN2+n] = F[n]*X[tN2-N2+n] + (1.0f-F[n])*tanhf(X[tN2+n]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                F[n] = 1.0f / (1.0f+expf(-X[NT+nT]));
                X[nT] = (1.0f-F[n]) * tanhf(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&X[NT+t],(int)T,F,1);
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&X[t-1],(int)T,o,F,1);
                for (size_t n=0u; n<N; ++n)
                {
                    F[n] = 1.0f / (1.0f+expf(-F[n]));
                    H[n] = F[n]*X[t-1+n*T];
                }
                cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    X[t+nT] = F[n]*X[t-1+nT] + (1.0f-F[n])*tanhf(X[t+nT]);
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                F[n] = 1.0f / (1.0f+expf(-X[NT+nT]));
                X[nT] = (1.0f-F[n]) * tanhf(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_scopy((int)N,&X[NT+t],(int)T,F,1);
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&X[t-1],(int)T,o,F,1);
                for (size_t n=0u; n<N; ++n)
                {
                    F[n] = 1.0f / (1.0f+expf(-F[n]));
                    H[n] = F[n]*X[t-1+n*T];
                }
                cblas_sgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    X[t+nT] = F[n]*X[t-1+nT] + (1.0f-F[n])*tanhf(X[t+nT]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                F[n] = 1.0f / (1.0f+expf(-X[N+n]));
                X[n] = (1.0f-F[n]) * tanhf(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN2 = 2*tN;
                cblas_scopy((int)N,&X[tN2+N],1,F,1);
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&X[tN2-N2],1,o,F,1);
                for (size_t n=0u; n<N; ++n)
                {
                    F[n] = 1.0f / (1.0f+expf(-F[n]));
                    H[n] = F[n]*X[tN2-N2+n];
                }
                cblas_sgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[tN2],1);
                for (size_t n=0u; n<N; ++n)
                {
                    X[tN2+n] = F[n]*X[tN2-N2+n] + (1.0f-F[n])*tanhf(X[tN2+n]);
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in gru_inplace_min2_s: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int gru_min_inplace_d (double *X, const double *U, const double *Uf, const size_t N, const size_t T, const char iscolmajor, const size_t dim)
{
    const double o = 1.0;
    const size_t NT = N*T, N2 = 2*N;
    size_t nT, tN, tN2;

    double *F, *H;
    if (!(F=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in gru_inplace_min2_d: problem with malloc. "); perror("malloc"); return 1; }
    if (!(H=(double *)malloc(N*sizeof(double)))) { fprintf(stderr,"error in gru_inplace_min2_d: problem with malloc. "); perror("malloc"); return 1; }

    if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                F[n] = 1.0 / (1.0+exp(-X[N+n]));
                X[n] = (1.0-F[n]) * tanh(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN2 = 2*tN;
                cblas_dcopy((int)N,&X[tN2+N],1,F,1);
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&X[tN2-N2],1,o,F,1);
                for (size_t n=0u; n<N; ++n)
                {
                    F[n] = 1.0 / (1.0+exp(-F[n]));
                    H[n] = F[n]*X[tN2-N2+n];
                }
                cblas_dgemv(CblasColMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[tN2],1);
                for (size_t n=0u; n<N; ++n)
                {
                    X[tN2+n] = F[n]*X[tN2-N2+n] + (1.0-F[n])*tanh(X[tN2+n]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                F[n] = 1.0 / (1.0+exp(-X[NT+nT]));
                X[nT] = (1.0-F[n]) * tanh(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&X[NT+t],(int)T,F,1);
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,Uf,(int)N,&X[t-1],(int)T,o,F,1);
                for (size_t n=0u; n<N; ++n)
                {
                    F[n] = 1.0 / (1.0+exp(-F[n]));
                    H[n] = F[n]*X[t-1+n*T];
                }
                cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    X[t+nT] = F[n]*X[t-1+nT] + (1.0-F[n])*tanh(X[t+nT]);
                }
            }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=0u; n<N; ++n)
            {
                nT = n*T;
                F[n] = 1.0 / (1.0+exp(-X[NT+nT]));
                X[nT] = (1.0-F[n]) * tanh(X[nT]);
            }
            for (size_t t=1; t<T; ++t)
            {
                cblas_dcopy((int)N,&X[NT+t],(int)T,F,1);
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&X[t-1],(int)T,o,F,1);
                for (size_t n=0u; n<N; ++n)
                {
                    F[n] = 1.0 / (1.0+exp(-F[n]));
                    H[n] = F[n]*X[t-1+n*T];
                }
                cblas_dgemv(CblasColMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[t],(int)T);
                for (size_t n=0u; n<N; ++n)
                {
                    nT = n*T;
                    X[t+nT] = F[n]*X[t-1+nT] + (1.0-F[n])*tanh(X[t+nT]);
                }
            }
        }
        else
        {
            for (size_t n=0u; n<N; ++n)
            {
                F[n] = 1.0 / (1.0+exp(-X[N+n]));
                X[n] = (1.0-F[n]) * tanh(X[n]);
            }
            for (size_t t=1; t<T; ++t)
            {
                tN = t*N; tN2 = 2*tN;
                cblas_dcopy((int)N,&X[tN2+N],1,F,1);
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,Uf,(int)N,&X[tN2-N2],1,o,F,1);
                for (size_t n=0u; n<N; ++n)
                {
                    F[n] = 1.0 / (1.0+exp(-F[n]));
                    H[n] = F[n]*X[tN2-N2+n];
                }
                cblas_dgemv(CblasRowMajor,CblasTrans,(int)N,(int)N,o,U,(int)N,H,1,o,&X[tN2],1);
                for (size_t n=0u; n<N; ++n)
                {
                    X[tN2+n] = F[n]*X[tN2-N2+n] + (1.0-F[n])*tanh(X[tN2+n]);
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in gru_inplace_min2_d: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif


text:
/**
 * @file utilities.h
 * @author Matthew Nicely (mnicely@nvidia.com)
 * @date 2020-01-06
 * @version 1.0
 * @brief Contains additional functions used throughout
 * the application.
 *
 * @copyright Copyright (c) 2020
 *
 * @license This project is released under the GNU Public License
 *
 * @note * Target Processor: Intel x86
 * @n * Target Compiler: GCC 7.4.0
 * @n * NVCC Compiler: CUDA Toolkit 10.0 or later
 */
#ifndef UTILITIES_H_
#define UTILITIES_H_

#include <algorithm>  // std::remove_if, std::nth_element
#include <cassert>    // assert
#include <cblas.h>    //cblas_sgemm
#include <cmath>      // std::sqrt,
#include <fstream>    // std::ofstream
#include <iterator>   // std::vector<>::const_iterator
#include <lapacke.h>  // LAPACKE_sgetrf, LAPACKE_sgetri
#include <map>        // std::map
#include <numeric>    // std::accumulate, std::inner_product
#include <random>     // std::default_random_engine, std::random_device, std::uniform_real_distribution
#include <stdexcept>  // std::runtime_error
#include <string>     // std::string
#include <vector>     // std::vector

// *************** FOR ERROR CHECKING *******************
#ifndef CUDA_RT_CALL
#define CUDA_RT_CALL( call )                                                                                \
    {                                                                                \
        auto status = static_cast<cudaError_t>( call );                                                                \
        if ( status != cudaSuccess )                                                                                \
            fprintf( stderr,                                                                                \
                     "ERROR: CUDA RT call \"%s\" in line %d of file %s failed "                                        \
                     "with "                                                                                \
                     "%s (%d).\n",                                                                                \
                     #call,                                                                                \
                     __LINE__,                                                                                \
                     __FILE__,                                                                                \
                     cudaGetErrorString( status ),                                                                     \
                     status );                                                                                \
    }
#endif  // CUDA_RT_CALL
// *************** FOR ERROR CHECKING *******************

/**
 * @brief Contains additional functions used throughout application
 */
namespace utility {

constexpr int kSysDim { 4 };  /**< Number of states in system model */
constexpr int kMeasDim { 2 }; /**< Number of states in measurement model */

/**
 * @enum Processor
 * @brief Enumeration of processor types
 */
enum class Processor : bool { kCpu, kGpu };

using ProcessorMap = std::map<Processor, std::string>;

static const ProcessorMap processor_map = { { Processor::kCpu, "CPU" }, { Processor::kGpu, "GPU" } };

/**
 * @enum Timing
 * @brief Enumeration of timing calculations
 */
enum class Timing : int { kMedian, kMean, kStdDev, kCount };

using TimingMap = std::map<Timing, std::string>;

static const TimingMap timing_map = { { Timing::kMedian, "Median" },
                                      { Timing::kMean, "Mean" },
                                      { Timing::kStdDev, "StdDev" } };

/**
 * @enum Filter
 * @brief Enumeration of possible particle filters
 */
enum class Filter : int { kBootstrap, kCount };

using FilterMap = std::map<Filter, std::string>;

static const FilterMap filters_map = { { Filter::kBootstrap, "Bootstrap" } };

/**
 * @enum Method
 * @brief Enumeration of possible resampling method
 */
enum class Method : int { kSystematic, kStratified, kMetropolisC2, kCount };

using MethodMap = std::map<Method, std::string>;

static const MethodMap method_map = { { Method::kSystematic, "Systematic" },
                                      { Method::kStratified, "Stratified" },
                                      { Method::kMetropolisC2, "MetropolisC2" } };

/**
 * @struct FilterInfo
 * @brief Structure to store user-specified data and send to filter classes
 */
typedef struct _filterInfo {
    int samples;     /**< Quantity of samples specified from input file */
    int particles;   /**< Quantity of particles specified from input file */
    int num_mcs;     /**< Quantity of Monte Carlos specified from input file */
    int resampling;  /**< Systmatic = 0, Stratified = 1, Metropolis = 2 */
    int filter_type; /**< Type of filter: Bootstrap = 0 */
} FilterInfo;

/**
 * @struct Matrix
 * @brief Structure to store matrix and information
 *
 * @param[in] row Stores quantity of rows
 * @param[in] col Stores quantity of columns
 * @param[in] val Store matrix data in a vector
 */
typedef struct _matrix {

    int                row;
    int                col;
    std::vector<float> val;

    _matrix( const _matrix &a ) : row( a.row ), col( a.col ), val( a.val ) {};

    _matrix( int const &row, int const &col, std::initializer_list<float> const &list ) :
        row( row ), col( col ), val( list ) {}

    _matrix( int const &row, int const &col ) : row( row ), col( col ) {
        val.resize( row * col, 0.0f );
    }

    _matrix( int const &row, int const &col, std::vector<float> const val ) : row( row ), col( col ), val( val ) {}
} Matrix;

/**
 * @brief Template function to compute square root of a matrix using floats
 *
 * @see http://www.netlib.org/lapack/lapacke.html
 * @see
 * http://www.netlib.org/lapack/explore-html/d7/d75/lapacke__ssyevr_8c_source.html
 * @see http://www.openblas.net/
 * @see http://www.seehuhn.de/pages/matrixfn
 *
 * @param[in,out] matrix_a Matrix to process
 */
inline void MatrixSqrt( utility::Matrix &matrix_a ) {

    int    n { matrix_a.row };
    float *a { matrix_a.val.data( ) };

    /* allocate space for the output parameters and workspace arrays */
    int   m {};
    float abstol { -1.0f };

    std::vector<float> w( n, 0.0f );
    std::vector<float> z( n * n, 0.0f );
    std::vector<int>   suppz( 2 * n, 0 );

    /* get the eigenvalues and eigenvectors */
    int info = LAPACKE_ssyevr(
        LAPACK_ROW_MAJOR, 'V', 'A', 'L', n, a, n, 0, 0, 0, 0, abstol, &m, w.data( ), z.data( ), n, suppz.data( ) );

    if ( info != 0 ) {
        throw std::runtime_error( "Issue with LAPACKE_ssyevr.\n" );
    }

    /* allocate and Initialize a new matrix B=Z*D */
    std::vector<float> b( n * n, 0.0f );
    for ( int j = 0; j < n; ++j ) {
        float lambda = std::sqrt( w[j] );
        for ( int i = 0; i < n; ++i ) {
            b[i * n + j] = z[i * n + j] * lambda;
        }
    }

    float alpha { 1.0f };
    float beta { 0.0f };

    /* calculate the square root A=B*Z^T */
    cblas_sgemm( CblasRowMajor, CblasNoTrans, CblasTrans, n, n, n, alpha, b.data( ), n, z.data( ), n, beta, a, n );
}

/**
 * @brief Matrix-Matrix multiplication function using floats
 *
 * @see http://www.openblas.net/
 *
 * @param[out] matrix_c Structure to store A * B
 * @param[in] matrix_a Structure reference to matrix A
 * @param[in] matrix_b Structure reference to matrix B
 */

inline void MatrixMult( utility::Matrix &matrix_c, utility::Matrix const &matrix_a, utility::Matrix const &matrix_b ) {

    assert( matrix_a.col == matrix_b.row );

    float const *a { matrix_a.val.data( ) };
    float const *b { matrix_b.val.data( ) };
    float *      c { matrix_c.val.data( ) };

    float alpha { 1.0f };
    float beta { 0.0f };

    cblas_sgemm( CblasRowMajor,
                 CblasNoTrans,
                 CblasNoTrans,
                 matrix_a.row,
                 matrix_b.col,
                 matrix_b.row,
                 alpha,
                 a,
                 matrix_b.row,
                 b,
                 matrix_b.col,
                 beta,
                 c,
                 matrix_b.col );
}

/**
 * @brief Matrix-Matrix-Matrix multiplication function using floats
 *
 * @see http://www.openblas.net/
 *
 * @param[out] matrix_d Structure to store A * B * C
 * @param[in] matrix_a Structure reference to matrix A
 * @param[in] matrix_b Structure reference to matrix B
 * @param[in] matrix_c Structure reference to matrix C
 */
inline void MatrixMult( utility::Matrix &      matrix_d,
                        utility::Matrix const &matrix_a,
                        utility::Matrix const &matrix_b,
                        utility::Matrix const &matrix_c ) {

    assert( matrix_a.col == matrix_b.row );
    assert( matrix_b.col == matrix_c.row );

    float const *a { matrix_a.val.data( ) };
    float const *b { matrix_b.val.data( ) };
    float const *c { matrix_c.val.data( ) };
    float *      d { matrix_d.val.data( ) };

    std::vector<float> temp_storage( matrix_a.row * matrix_b.col, 0.0f );
    float *            ab { temp_storage.data( ) };  // Temp store for A*B

    float alpha { 1.0f };
    float beta { 0.0f };

    cblas_sgemm( CblasRowMajor,
                 CblasNoTrans,
                 CblasNoTrans,
                 matrix_a.row,
                 matrix_b.col,
                 matrix_b.row,
                 alpha,
                 a,
                 matrix_b.row,
                 b,
                 matrix_b.col,
                 beta,
                 ab,
                 matrix_b.col );

    cblas_sgemm( CblasRowMajor,
                 CblasNoTrans,
                 CblasNoTrans,
                 matrix_a.row,
                 matrix_c.col,
                 matrix_c.row,
                 alpha,
                 ab,
                 matrix_c.row,
                 c,
                 matrix_c.col,
                 beta,
                 d,
                 matrix_c.col );
}

/**
 * @brief Matrix-Vector multiplication function using floats
 *
 * @see http://www.openblas.net/
 *
 * @param[out] vector_c Vector to store A * B
 * @param[in] matrix_a Structure reference to matrix A
 * @param[in] vector_b Vector reference to vector B
 */
inline void
MatrixMult( std::vector<float> &vector_c, utility::Matrix const &matrix_a, std::vector<float> const &vector_b ) {

    assert( matrix_a.col == vector_b.size( ) );

    float const *a { matrix_a.val.data( ) };
    float const *b { vector_b.data( ) };
    float *      c { vector_c.data( ) };

    float alpha { 1.0f };
    float beta { 0.0f };

    cblas_sgemv( CblasRowMajor, CblasNoTrans, matrix_a.row, matrix_a.col, alpha, a, matrix_a.col, b, 1, beta, c, 1 );
}

/**
 * @brief Function to compute the inversion of a matrix using floats
 *
 * @see http://www.netlib.org/lapack/lapacke.html
 * @see
 * http://www.netlib.org/lapack/explore-html/d2/d96/lapacke__dgetrf_8c_source.html
 * @see
 * http://www.netlib.org/lapack/explore-html/da/d0e/lapacke__dgetri_8c_source.html
 * @see
 * https://stackoverflow.com/questions/3519959/computing-the-inverse-of-a-matrix-using-lapack-in-c
 *
 * @param[in,out] matrix_a Matrix to be inverted
 */
inline void ComputeInverse( utility::Matrix &matrix_a ) {

    // A must be a square matrix???
    int              n { matrix_a.row };
    float *          a = matrix_a.val.data( );
    std::vector<int> piv( n + 1, 0 );

    int info {};
    info = LAPACKE_sgetrf( LAPACK_ROW_MAJOR, n, n, a, n, piv.data( ) );

    if ( info != 0 ) {
        throw std::runtime_error( "Issue with LAPACKE_sgetrf.\n" );
    }

    info = LAPACKE_sgetri( LAPACK_ROW_MAJOR, n, a, n, piv.data( ) );

    if ( info != 0 ) {
        throw std::runtime_error( "Issue with LAPACKE_sgetri.\n" );
    }
}

/**
 * @brief Normal distribution random number generator
 *
 * @param[in] random_numbers Stores random numbers
 * @param[in] lBound Lower bound of distribution
 * @param[in] uBound Upper bound of distribution
 */
inline void GenerateRandomNum( std::vector<float> &random_numbers, float const &lBound, float const &uBound ) {

    // Random device class instance, source of 'true' randomness for
    // initializing random seed
    std::default_random_engine gen( std::random_device {}( ) );

    // Instance of class std::normal_distribution with specific mean and stddev
    // Lower and Upper limits are declared in base class
    std::normal_distribution<float> distr( lBound, uBound );

    // Compute random numbers for all particles
    for_each( random_numbers.begin( ), random_numbers.end( ), [&distr, &gen]( float &a ) { a = distr( gen ); } );
}

/**
 * @brief Compute median for Monte Carlos times
 *
 * @param[in/out] time_vector Vector holding times for all Monte Carlos
 * @param[out] median Median time for all Monte Carlos
 */
inline void ComputeMedianTime( std::vector<float> &time_vector, float &median ) {
    // Compute median
    std::nth_element( time_vector.begin( ), time_vector.begin( ) + time_vector.size( ) / 2, time_vector.end( ) );
    median = time_vector[time_vector.size( ) / 2];

    // Remove outliers
    time_vector.erase( std::remove_if( time_vector.begin( ),
                                       time_vector.end( ),
                                       [&median]( float const &n ) { return ( n > ( median * 1.05f ) ); } ),
                       time_vector.end( ) );
}

/**
 * @brief Compute mean for Monte Carlos times
 *
 * @param[in] time_vector Vector holding times for all Monte Carlos
 * @param[out] mean Mean time for all Monte Carlos
 */
inline void ComputeMeanTime( std::vector<float> const &time_vector, float &mean ) {
    // Compute mean
    mean = std::accumulate( time_vector.begin( ), time_vector.end( ), 0.0f ) / time_vector.size( );
}

/**
 * @brief Compute standard deviation for Monte Carlos times
 *
 * @param[in] time_vector Vector holding times for all Monte Carlos
 * @param[in] mean Mean time for all Monte Carlos
 * @param[out] stdDev
 */
inline void ComputeStdDevTime( std::vector<float> const &time_vector, float const &mean, float &stdDev ) {

    // Compute standard deviation
    std::vector<float> diff( time_vector.size( ) );
    std::transform(
        time_vector.begin( ), time_vector.end( ), diff.begin( ), [&mean]( float const &x ) { return ( x - mean ); } );

    // Compute difference
    float sqrtSum { std::inner_product( diff.begin( ), diff.end( ), diff.begin( ), 0.0f ) };  // Compute product
    stdDev = std::sqrt( sqrtSum / time_vector.size( ) );
}

/**
 * @brief Compute median, mean, standard deviation for Monte Carlo times
 *
 * @param[in/out] time_vector
 * @param[in/out] median
 * @param[in/out] mean
 * @param[in/out] stdDev
 */
inline void ComputeOverallTiming( std::vector<float> &timeVec, float &median, float &mean, float &stdDev ) {
    ComputeMedianTime( timeVec, median );
    ComputeMeanTime( timeVec, mean );
    ComputeStdDevTime( timeVec, mean, stdDev );
}

/**
 * @brief Compute overall timing across all Monte Carlos
 *
 * @param[in] timing_results 2D vector containing timing results for each Monte
 * Carlo
 */
inline void PrintTimingResults( std::vector<std::vector<float>> const &timing_results ) {

    for ( int i = 0; i < static_cast<int>( utility::Timing::kCount ); i++ ) {
        std::printf( "%s\t ", utility::timing_map.find( utility::Timing( i ) )->second.c_str( ) );
    }
    std::printf( "\n" );

    for ( int i = 0; i < static_cast<int>( utility::Timing::kCount ); i++ ) {
        float sum { std::accumulate( timing_results[i].begin( ), timing_results[i].end( ), 0.0f ) };
        float mean { sum / timing_results[i].size( ) };
        std::printf( "%0.0f\t ", mean );
    }
    std::printf( "\n\n" );
}

/**
 * @brief Writes data to output file for generated system
 *        and measurement truth data
 *
 * @param[in] filename Output filename
 * @param[in] dim Number of system dimensions
 * @param[in] samples Number of samples
 * @param[in] it Iterator to beginning of vector to be printed
 */
inline void WriteToFile( std::string const &                          filename,
                         int const &                                  dim,
                         int const &                                  samples,
                         typename std::vector<float>::const_iterator &it ) {

    // We jump to the end to add each Monte Carlo
    std::ofstream fout( filename, std::ios_base::app );
    fout.precision( 15 );
    fout.setf( std::ios::fixed, std::ios::floatfield );
    // fout << dim << " " << samples << "\n";
    for ( int i = 0; i < samples; i++ ) {
        int idxS { i * dim };
        int idxE { idxS + dim };
        std::copy( it + idxS, it + idxE, std::ostream_iterator<float>( fout, " " ) );
        fout << "\n";
    }
    fout.close( );
}

/**
 * @brief Writes data to output file for filtered estimates
 *
 * @param[in] filename Output filename
 * @param[in] type Filter type
 * @param[in] samples Number of samples
 * @param[in] median Timing median across all Monte Carlos
 * @param[in] mean Timing mean across all Monte Carlos
 * @param[in] stdDev Timing standard deviation across all Monte Carlos
 * @param[in] it Iterator to beginning of vector to be printed
 */
inline void WriteToFile( std::string const &                          filename,
                         std::string const &                          type,
                         int const &                                  samples,
                         float const &                                median,
                         float const &                                mean,
                         float const &                                stdDev,
                         typename std::vector<float>::const_iterator &it ) {

    // We jump to the end to add each Monte Carlo
    std::ofstream fout( filename, std::ios_base::app );
    fout.setf( std::ios::fixed, std::ios::floatfield );
    fout << type << " " << median << " " << mean << " " << stdDev << "\n";
    fout.precision( 15 );
    for ( int i = 0; i < samples; i++ ) {
        int idxS { i * utility::kSysDim };
        int idxE { idxS + utility::kSysDim };
        std::copy( it + idxS, it + idxE, std::ostream_iterator<float>( fout, " " ) );
        fout << "\n";
    }
    fout.close( );
}

/**
 * @brief Write header to output file for generated system
 *        and measurement truth data
 *
 * @param[in] filename Output filename
 * @param[in] info Struct containing filter information
 */
inline void WriteTruthHeader( std::string const &filename, int const &num_mcs, int const &dim, int const &samples ) {
    std::ofstream fout( filename, std::ios_base::trunc );
    fout.setf( std::ios::fixed, std::ios::floatfield );
    fout << num_mcs << " " << dim << " " << samples << "\n";
    fout.close( );
}

/**
 * @brief Write header to output file.
 *
 * @param[in] filename Output filename
 * @param[in] info Struct containing filter information
 */
inline void WriteDataHeader( std::string const &filename, utility::FilterInfo const &filter_info ) {
    std::ofstream fout( filename, std::ios_base::trunc );
    fout.setf( std::ios::fixed, std::ios::floatfield );
    fout << filter_info.num_mcs << " " << kSysDim << " " << filter_info.samples << " " << filter_info.particles << "\n";
    fout.close( );
}

#ifdef USE_NVTX
#include "nvtx3/nvToolsExt.h"

uint32_t const colors[] { 0x0000ff00, 0x0f0000ff, 0x000000ff, 0x00ffff00, 0x00ff00ff,
                          0x0000ffff, 0x00ff0000, 0x00ffffff, 0x0f0f0f0f };
int const      num_colors { sizeof( colors ) / sizeof( uint32_t ) };

/**
 * class Tracer
 */
class Tracer {
  public:
    /**
     * @brief Construct a new Tracer object
     *
     * @see
     * https://docs.nvidia.com/gameworks/content/gameworkslibrary/nvtx/nvidia_tools_extension_library_nvtx.htm
     * @see
     * https://devblogs.nvidia.com/cuda-pro-tip-generate-custom-application-profile-timelines-nvtx/
     *
     * @param[in] name Function name
     * @param[in] cid color id
     */

    Tracer( char const *name, int const &cid ) {
        int color_id                      = cid;
        color_id                          = color_id % num_colors;
        nvtxEventAttributes_t eventAttrib = {};
        eventAttrib.version               = NVTX_VERSION;
        eventAttrib.size                  = NVTX_EVENT_ATTRIB_STRUCT_SIZE;
        eventAttrib.colorType             = NVTX_COLOR_ARGB;
        eventAttrib.color                 = colors[color_id];
        eventAttrib.messageType           = NVTX_MESSAGE_TYPE_ASCII;
        eventAttrib.message.ascii         = name;
        nvtxRangePushEx( &eventAttrib );
    }

    /**
     * @brief Destroy the Tracer object
     *
     */
    ~Tracer( ) {
        nvtxRangePop( );
    }
};

#define RANGE( name, cid ) utility::Tracer uniq_name_using_macros( name, cid );
#else
#define RANGE( name, cid )
#endif

} /* namespace utility */

#endif /* UTILITIES_H_ */


text:
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_sf_legendre.h>
#include "sphere.h"
#include "sphSolv.h"
#include <petsc.h>

#undef __FUNCT__
#define __FUNCT__ "CalcEnm"
PetscErrorCode CalcEnm(Vec chargeSph, Vec chargeMag, PetscInt n, PetscInt m, PetscReal *rval, PetscReal *ival)
{
  PetscErrorCode ierr;
  PetscInt numCharges;
  PetscScalar *chargeMagPtr, *chargeSphPtr;
  PetscFunctionBegin;

  ierr = VecGetSize(chargeMag, &numCharges);CHKERRQ(ierr);

  PetscInt absm = abs(m);
  PetscReal coef = 1;
  for(PetscInt i=n-absm+1; i <= n + absm; ++i)
    coef *= i;
  coef = 1./coef;

  *rval = 0;
  *ival = 0;

  PetscReal qk, rk, rkton, thetak, phik, lterm, rterm, iterm;

  ierr = VecGetArrayRead(chargeMag, &chargeMagPtr);CHKERRQ(ierr);
  ierr = VecGetArrayRead(chargeSph, &chargeSphPtr);CHKERRQ(ierr);
  for(PetscInt k=0; k<numCharges; ++k) {
    qk = chargeMagPtr[k];
    rk = chargeSphPtr[3*k+0];
    thetak = chargeSphPtr[3*k+1];
    phik = chargeSphPtr[3*k+2];

    rkton = 1;
    for(PetscInt i=0; i<n; ++i)
      rkton = rkton*rk;

    //seems like we don't have to worry about (-1)^m? if(m%2 ==1) lterm*=-1; not necessary
    lterm = gsl_sf_legendre_Plm(n, absm, cos(thetak));
    rterm = PetscCosReal(m*phik);
    iterm = -PetscSinReal(m*phik);

    //update val
    *rval += qk*coef*rkton*lterm*rterm;
    *ival += qk*coef*rkton*lterm*iterm;
  }
  ierr = VecRestoreArrayRead(chargeMag, &chargeMagPtr);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(chargeSph, &chargeSphPtr);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

void calcEnm(SPoint *positions, double *charges, int nCharges, int n, int m, double *rval, double *ival)
{
  int print = 0;
  if(n < 2) {
    printf("\n\nn = %d\nm = %d\n", n, m);
    print = 0;
  }

  int absm = abs(m);
  double coef = 1;
  //compute factorial coefficient
  for(int i = n - absm + 1; i <= n + absm; ++i)
    coef *= i;
  coef = 1./coef;

  *rval = 0;
  *ival = 0;

  double qk, rk, rkton, thetak, phik, lterm, rterm, iterm;
  //double *legvals = (double*) malloc(sizeof(double)*(2*n+1));
  //sum over charges
  for(int k=0; k<nCharges; ++k) {

    qk = charges[k];
    rk = positions[k].x1;
    thetak = positions[k].x2;
    phik = positions[k].x3;

    rkton = 1;
    for(int i=0; i<n; ++i)
      rkton = rkton*rk;

    //seems like we don't have to worry about (-1)^m? if(m%2 ==1) lterm*=-1; not necessary
    lterm = gsl_sf_legendre_Plm(n, absm, cos(thetak));//boost::math::legendre_p(n, absm, cos(thetak));

    rterm = cos(m*phik);
    iterm = -sin(m*phik);

    if(print) {
      printf("QK: %15.15f\n", qk);
      printf("rkton: %15.15f\n", rkton);
      printf("ltern: %15.15f\n", lterm);
    }

    //update val
    *rval += qk * coef * rkton * lterm * rterm;
    *ival += qk * coef * rkton * lterm * iterm;
  }

}

void calcBnmFromEnm(int n, int m, double Enmrval, double Enmival, SProblem *problem, double *Bnmrval, double *Bnmival) {
  double e1 = problem->e1;
  double e2 = problem->e2;
  double b  = problem->b;

  double val = 1.0;
  for(int i=0; i<2*n+1; ++i) {
    val /= b;
  }

  val /= e1;

  val *= (e1-e2)*(n+1);
  val /= (e1*n + e2*(n+1));
  //printf("VAL: %15.15f\n", val);
  *Bnmrval = val*Enmrval;
  *Bnmival = val*Enmival;

}

void calcCnmFromEnm(int n, int m, double Enmrval, double Enmival, SProblem *problem, double *Cnmrval, double *Cnmival) {
  double e1 = problem->e1;
  double e2 = problem->e2;

  double val;

  val = n*(e1-e2)/(n*e1+e2*(n+1));
  val = (1./e2)*(1-val);

  *Cnmrval = val * Enmrval;
  *Cnmival = val * Enmival;

}

#undef __FUNCT__
#define __FUNCT__ "CalcSphericalCoulombPotential"
PetscErrorCode CalcSphericalCoulombPotential(PetscReal b, PetscReal eps1, PetscReal eps2, PetscInt nCharges, Vec chargeXYZ, Vec chargeMag, PetscInt nSol, Vec solXYZ, PetscInt Nmax, Vec targetSol)
{
  PetscErrorCode ierr;
  PetscInt numSol;
  Vec chargeSph;
  Vec solSph;
  Vec coulCoefs;
  PetscScalar *vecPtr, *solSphPtr;
  PetscReal solRad, solTheta, solPhi;
  PetscReal lterm, rterm, iterm;
  PetscReal Enmrval, Enmival;
  PetscReal radpow = 1.0;
  PetscFunctionBegin;

  //get number of solution points
  ierr = VecGetSize(targetSol, &numSol);CHKERRQ(ierr);

  // create charge spherical vec and convert from xyz
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*nCharges, &chargeSph);CHKERRQ(ierr);
  ierr = CartesianToSpherical(chargeXYZ, chargeSph);CHKERRQ(ierr);
  //create solution spherical vec and convert from xyz
  ierr = VecCreateSeq(PETSC_COMM_SELF, 3*nSol, &solSph);CHKERRQ(ierr);
  ierr = CartesianToSpherical(solXYZ, solSph);CHKERRQ(ierr);

  //zero solutions vec
  ierr = VecSet(targetSol, 0.0);CHKERRQ(ierr);

  ierr = VecGetArrayRead(solSph, &solSphPtr);CHKERRQ(ierr);
  ierr = VecGetArray(targetSol, &vecPtr);CHKERRQ(ierr);

  for(PetscInt n=0; n<=Nmax; ++n) {
    for(PetscInt m=-n; m <= n; ++m) {

      ierr = CalcEnm(chargeSph, chargeMag, n, m, &Enmrval, &Enmival);CHKERRQ(ierr);
      for(PetscInt k=0; k<numSol; ++k) {
	solRad = solSphPtr[3*k+0];
	solTheta = solSphPtr[3*k+1];
	solPhi = solSphPtr[3*k+2];
	lterm = gsl_sf_legendre_Plm(n, abs(m), cos(solTheta));//boost::math::legendre_p(n, abs(m), cos(theta));//legvals[n+abs(m)];//boost::math::legendre_p(n, m, cos(theta));
	rterm = cos(m*solPhi);
	iterm = sin(m*solPhi);
	vecPtr[k] += ((Enmrval*rterm - Enmival*iterm)*lterm)/(eps1*pow(solRad,n+1));
      }
    }
  }

  ierr = VecRestoreArray(targetSol, &vecPtr);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(solSph, &solSphPtr);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

double calcCoulombSpherical(SProblem *problem, int maxDeg, SPoint *r) {
  double val = 0;
  double e1 = problem->e1;
  double rad = r->x1;
  double theta = r->x2;
  double phi = r->x3;
  double radpow = 1.0;
  double lterm;
  double Enmrterm, Enmiterm;
  double rterm, iterm;

  for(int n=0; n<=maxDeg; ++n) {
    radpow *= rad;
    for(int m=-n; m<=n; ++m) {
      calcEnm(problem->positions, problem->charges, problem->nCharges, n, m, &Enmrterm, &Enmiterm);
      //printf("Enm %d %d %15.15f\n", n, m, Enm);

      //printf("cos(theta): %15.15f\n", theta);
      lterm = gsl_sf_legendre_Plm(n, abs(m), cos(theta));//boost::math::legendre_p(n, abs(m), cos(theta));//legvals[n+abs(m)];//boost::math::legendre_p(n, m, cos(theta));

      rterm = cos(m*phi);
      iterm = sin(m*phi);

      val += ((Enmrterm*rterm - Enmiterm*iterm)*lterm)/(e1*radpow);
      //printf("lterm: %15.15f\n", lterm);
      //printf("cos(mphi) = %15.15f\n\n", cos(m*phi));
    }
    printf("n = %d\nval = %15.15f\n", n, val);
  }

  return 0.0;
}

double calcCoulombSphericalGrid(SProblem *problem, int maxDeg, SPoint *r, int nPoints, double *solution) {
  double e1 = problem->e1;

  double rad, theta, phi;
  double lterm;
  double Enmrval, Enmival;
  double Bnmrval, Bnmival;
  double Cnmrval, Cnmival;
  double rval, ival;

  for(int i=0; i<nPoints; ++i)
    solution[i] = 0.0;
  for(int n=0; n<=maxDeg; ++n) {
    for(int m=-n; m<=n; ++m) {

      calcEnm(problem->positions, problem->charges, problem->nCharges, n, m, &Enmrval, &Enmival);
      calcBnmFromEnm(n, m, Enmrval, Enmival, problem, &Bnmrval, &Bnmival);
      calcCnmFromEnm(n, m, Enmrval, Enmival, problem, &Cnmrval, &Cnmival);

      for(int k=0; k<nPoints; ++k) {

	rad   = r[k].x1;
	theta = r[k].x2;
	phi   = r[k].x3;

	lterm = gsl_sf_legendre_Plm(n, abs(m), cos(theta));//boost::math::legendre_p(n, abs(m), cos(theta));

	rval = cos(m*phi);
	ival = sin(m*phi);

	if(rad < problem->b) {

	  solution[k] += (Bnmrval*rval - Bnmival*ival)*lterm*pow(rad,n);

	}
	else if(rad >  problem->b) {

	  solution[k] += ((Cnmrval*rval - Cnmival*ival)*lterm)/pow(rad,n+1);

	}
	else {

	  printf("CALCULATING ON BOUNDARY\n");
	  solution[k] += ((Enmrval*rval - Enmival*ival)*lterm)/(e1*pow(rad,n+1));
	}
      }
    }
    printf("n = %d\nval = %15.15f\n", n, solution[0]);
  }

  return 0.0;
}

max_stars_repo_name: tom-klotz/ellipsoid-solvation

max_issues_repo_name: tom-klotz/ellipsoid-solvation

max_forks_repo_name: tom-klotz/ellipsoid-solvation


text:
/** \file filters.h
 * A collection of filters and FFT related routines.
 */
#ifndef __FILTERS_H__
#define __FILTERS_H__
#include <fftw3.h>
#include <gsl/gsl_wavelet.h>
#include "common.h"

#define FFTW_NTHREADS_DEFAULT 4
#define FFTW_FLAGS_DEFAULT FFTW_ESTIMATE

typedef struct filters_handle
{
    /* public for i/o */
    ANALYSIS_WAVEFORM_BASE_TYPE *inWav;
    ANALYSIS_WAVEFORM_BASE_TYPE *outWav;
    ANALYSIS_WAVEFORM_BASE_TYPE *respWav;
    WAVELET_BASE_TYPE *waveletWav;

    /* private for internal work */
    size_t wavLen, respLen;
    int malloced;

    /* fft */
    int fftUsed;
    size_t fftLen;
    size_t fftwNThreads;
    unsigned fftwFlags;
    FFTW(plan) fftwPlan, fftwPlan1, fftwPlan2;
    FFT_BASE_TYPE *fftwWork, *fftwWork1;
    /* window function */
    FFT_BASE_TYPE *fftwWin;
    double fftwS1, fftwS2, dt;

    /* wavelet */
    gsl_wavelet *gslDWT;
    gsl_wavelet_workspace *gslDWTWork;
} filters_t;

/** Initialize the filter.
 * @param[in] inWav waveform input.
 * If inWav == NULL, a space is malloced and freed upon closure.
 * If inWav is given, it is not freed.
 * @param[in] n waveform length.
 * The output is also of length n.  n can be odd in principle,
 * but an even number is preferred.
 */
filters_t *filters_init(ANALYSIS_WAVEFORM_BASE_TYPE *inWav, size_t n);
/** Initialize for convolution and fft
 * @param[in] np length of response function.  Must be an odd number.
 * np=0 is reserved for fft spectrum calculations.
 */
filters_t *filters_init_for_convolution(ANALYSIS_WAVEFORM_BASE_TYPE *inWav, size_t n, size_t np);
/** Destroy the filter */
int filters_close(filters_t *fHdl);
/** Savitzky-Golay filter.
 * @param[in] m order of polynomial.
 * @param[in] np number of points (response length)
 * @param[in] ld degree of derivative
 */
int filters_SavitzkyGolay(filters_t *fHdl, int m, int ld);
/** Raised cosine as convolution kernel.
 * Outermost bin, i=(fHdl->respLen-1)/2, is set to 0.0
 **/
int filters_raisedCosine(filters_t *fHdl);
int filters_convolute(filters_t *fHdl);
int filters_hanning_window(filters_t *fHdl);
/** Compute the spectrum in Fourier space.  Requires init_for_convolution with np = 0.
 * Computed power spectrum is stored in fftwWork(s), normalized.
 * This function handles both odd and even number of input points.
 * The resulting power spectrum has the length (int)(n/2)+1 for both even and odd n.
 * Spectra density (linearized) is stored in fftwWork,
 * Spectrum (linearized) is stored in fftwWork1.
 */
int filters_fft_spectrum(filters_t *fHdl);
/** Discrete wavelet transform. */
int filters_DWT(filters_t *fHdl);
/** Median filter with moving window size n. */
int filters_median(filters_t *fHdl, size_t n);
/** Trapezoidal filter as in Knoll NIMA 345(1994) 337-345.
 * @param[in] k rise time.
 * @param[in] l delay of peak.
 * l-k is the flat-top duration.
 * @param[in] M decay time constant (in number of samples) of the input
 * pulse.  Set M=-1.0 to deal with a step-like input function.
 */
int filters_trapezoidal(filters_t *fHdl, size_t k, size_t l, double M);

#endif // __FILTERS_H__