text:
/* interpolation/linear.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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.
 */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>

static int
linear_init (void * vstate,
             const double x_array[],
             const double y_array[],
             size_t size)
{
  return GSL_SUCCESS;
}

static
int
linear_eval (const void * vstate,
             const double x_array[], const double y_array[], size_t size,
             double x,
             gsl_interp_accel * a,
             double *y)
{
  double x_lo, x_hi;
  double y_lo, y_hi;
  double dx;
  size_t index;

  if (a != 0)
    {
      index = gsl_interp_accel_find (a, x_array, size, x);
    }
  else
    {
      index = gsl_interp_bsearch (x_array, x, 0, size - 1);
    }

  /* evaluate */
  x_lo = x_array[index];
  x_hi = x_array[index + 1];
  y_lo = y_array[index];
  y_hi = y_array[index + 1];
  dx = x_hi - x_lo;
  if (dx > 0.0)
    {
      *y = y_lo + (x - x_lo) / dx * (y_hi - y_lo);
      return GSL_SUCCESS;
    }
  else
    {
      *y = 0.0;
      return GSL_EINVAL;
    }
}

static
int
linear_eval_deriv (const void * vstate,
                   const double x_array[], const double y_array[], size_t size,
                   double x,
                   gsl_interp_accel * a,
                   double *dydx)
{
  double x_lo, x_hi;
  double y_lo, y_hi;
  double dx;
  double dy;
  size_t index;

  if (a != 0)
    {
      index = gsl_interp_accel_find (a, x_array, size, x);
    }
  else
    {
      index = gsl_interp_bsearch (x_array, x, 0, size - 1);
    }

  /* evaluate */
  x_lo = x_array[index];
  x_hi = x_array[index + 1];
  y_lo = y_array[index];
  y_hi = y_array[index + 1];
  dx = x_hi - x_lo;
  dy = y_hi - y_lo;
  if (dx > 0.0)
    {
      *dydx = dy / dx;;
      return GSL_SUCCESS;
    }
  else
    {
      *dydx = 0.0;
      return GSL_EINVAL;
    }
}

static
int
linear_eval_deriv2 (const void * vstate,
                    const double x_array[], const double y_array[], size_t size,
                    double x,
                    gsl_interp_accel * a,
                    double *y_pp)
{
  *y_pp = 0.0;

  return GSL_SUCCESS;
}

static
int
linear_eval_integ (const void * vstate,
                   const double x_array[], const double y_array[], size_t size,
                   gsl_interp_accel * acc,
                   double a, double b,
                   double * result)
{
  size_t i, index_a, index_b;

  if (acc != 0)
    {
      index_a = gsl_interp_accel_find (acc, x_array, size, a);
      index_b = gsl_interp_accel_find (acc, x_array, size, b);
    }
  else
    {
      index_a = gsl_interp_bsearch (x_array, a, 0, size - 1);
      index_b = gsl_interp_bsearch (x_array, b, 0, size - 1);
    }

    /* endpoints span more than one interval */

  *result = 0.0;

  /* interior intervals */
  for(i=index_a; i<=index_b; i++) {
    const double x_hi = x_array[i + 1];
    const double x_lo = x_array[i];
    const double y_lo = y_array[i];
    const double y_hi = y_array[i + 1];
    const double dx = x_hi - x_lo;

    if(dx != 0.0) {
      if (i == index_a || i == index_b)
        {
          double x1 = (i == index_a) ? a : x_lo;
          double x2 = (i == index_b) ? b : x_hi;
          const double D = (y_hi-y_lo)/dx;
          *result += (x2-x1) * (y_lo + 0.5*D*((x2-x_lo)+(x1-x_lo)));
        }
      else
        {
          *result += 0.5 * dx * (y_lo + y_hi);
        }
    }
  }

  return GSL_SUCCESS;
}

static const gsl_interp_type linear_type =
{
  "linear",
  2,
  NULL, /* alloc, not applicable */
  &linear_init,
  &linear_eval,
  &linear_eval_deriv,
  &linear_eval_deriv2,
  &linear_eval_integ,
  NULL, /* free, not applicable */
};

const gsl_interp_type * gsl_interp_linear = &linear_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

#define STRICT
#define NOMINMAX
#include <phnt/phnt_windows.h>
#include <phnt/phnt.h>
#pragma comment( lib, "ntdll.lib" )

#include <delayimp.h>
#ifndef _delayimp_h
#define _delayimp_h
#endif
#pragma comment( lib, "delayimp.lib" )

EXTERN_C const IMAGE_DOS_HEADER __ImageBase;

#include <stdlib.h>
#include <stdio.h>
#include <strsafe.h>
#include <stdint.h>

#include <string>
#include <filesystem>
namespace fs = std::filesystem;

#include <codecvt>
#include <gsl/gsl>
#include <wil/token_helpers.h>
#include <wil/resource.h>

#define DETOURS_INTERNAL
#include <detours/detours.h>

#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_DEBUG
#define SPDLOG_WCHAR_TO_UTF8_SUPPORT
#include <spdlog/spdlog.h>
#include <spdlog/sinks/msvc_sink.h>
#include <spdlog/fmt/bin_to_hex.h>

ext: h

max_stars_repo_name: LowerCode/bnspatch

max_issues_repo_name: LowerCode/bnspatch

max_forks_repo_name: LowerCode/bnspatch


text:
#ifndef TOOLS_H
#define TOOLS_H

#include <iostream>
#include <string>
#include <vector>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <math.h>

using namespace std;

//required for TRACE / DEBUG
extern int indentLevel;
extern string BLANKLINE;

//how many blanks to indent
#define INDENTATION 2

#define INDENT BLANKLINE.substr(0,(indentLevel)*INDENTATION)
#define INDENT_INC BLANKLINE.substr(0,(indentLevel++)*INDENTATION)
#define INDENT_DEC BLANKLINE.substr(0,(--indentLevel)*INDENTATION)

#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define ABS(x) ((x)>0 ? (x) : -(x))

#ifdef STDOUT
#define OUTPUT(x) cout << INDENT << (x) <<"\n";
#define OUTPUT1(x,arg1) cout << INDENT << (x) <<arg1<<"\n";
#define OUTPUT2(x,arg1,arg2) cout << INDENT << (x) <<arg1<<" / " <<arg2<<"\n";

#else
#define OUTPUT(x) ""
#define OUTPUT1(x,arg1) ""
#define OUTPUT2(x,arg1,arg2) ""

#endif //STDOUT

#ifdef TRACE
#define TRACE_IN(x) cout << INDENT_INC << "***TRACE_IN: "<< (x) <<"()\n";
#define TRACE_IN1(x,arg1) cout << INDENT_INC << "***TRACE_IN: "<< (x) <<"("<<arg1<<")\n";
#define TRACE_OUT(x,ok) cout << INDENT_DEC << "***TRACE_OUT: "<< (x) <<"() = " << ok << "\n";

#else
#define TRACE_IN(x) ""
#define TRACE_IN1(x,arg1) ""
#define TRACE_OUT(x,ok) ""

#endif //TRACE

#ifdef TRACEV
#define TRACEV_IN(x) cout << INDENT_INC << "***TRACE_IN: "<< (x) <<"()\n";
#define TRACEV_IN1(x,arg1) cout << INDENT_INC << "***TRACE_IN: "<< (x) <<"("<<arg1<<")\n";
#define TRACEV_OUT(x,ok) cout << INDENT_DEC << "***TRACE_OUT: "<< (x) <<"() = " << ok << "\n";

#else
#define TRACEV_IN(x) ""
#define TRACEV_IN1(x,arg1) ""
#define TRACEV_OUT(x,ok) ""

#endif //TRACEV

#ifdef DBUG
#define DEBUG(x) cout << INDENT << "***DEBUG: "<< (x) <<"\n";
#define DEBUG1(x,arg1) cout << INDENT << "***DEBUG: "<< (x) <<arg1<<"\n";
#define DEBUG2(x,arg1,arg2) cout << INDENT << "***DEBUG: "<< (x) <<arg1<<" / " <<arg2<<"\n";

#else
#define DEBUG(x) ""
#define DEBUG1(x,arg1) ""
#define DEBUG2(x,arg1,arg2) ""

#endif //DBUG

//verbose debug information
//mainly exact positions (quite big with dim=30)
#ifdef DBUGV
#define DEBUGV(x) cout << "***DEBUG: "<< (x) <<"\n";
#define DEBUGV1(x,arg1) cout << "***DEBUG: "<< (x) <<arg1<<"\n";
#define DEBUGV2(x,arg1,arg2) cout << "***DEBUG: "<< (x) <<arg1<<" / " <<arg2<<"\n";

#else
#define DEBUGV(x) ""
#define DEBUGV1(x,arg1) ""
#define DEBUGV2(x,arg1,arg2) ""

#endif //DBUGV

///Little helper function calculating the maximum height
/**of a regularly built pyramid of
  *@param minDegree minimum degree*/
int getMaxHeight(int minDegree, int swarmsize);

//fix the seed for the global RNG
void fixSeed(unsigned seed);

///Free the RNG ressources
void freeRng();
///Return gaussian distributed double from [0;1)
double randDoubleGaussian(double sigma);
///Return double from [0;1)
double randDouble();
///Wrapper to GSL Function call
/**Returns double from range [min,max)*/
double randDoubleRange(double min, double max);
///Wrapper to GSL Function call
/**Returns int from range [min,max]*/
int randIntRange(int min, int max);

string printVec(const vector<double>&);

double euclideanDistance(const vector<double>& v1, const vector<double>& v2);

///Generate a reproducable sequence
/**The seed is passed and thus the sequence is repeatable*/
class SequenceGenerator {
public:
	///Pass the seed
	SequenceGenerator(unsigned seed_in);
	///Destructor
	~SequenceGenerator();
	///Return exponentially distributed random variable, with mean mu
	double nextDoubleExponential(double mu);

	///Wrapper to GSL Function call
	/**Returns double from range [min,max)*/
	double nextDoubleRange(double min, double max);
	///Wrapper to GSL Function call
	/**Returns int from range [min,max]*/
	int nextIntRange(int min, int max);
private:
	///The random number generator
	gsl_rng*  randGen;

	///The random seed defining the sequence
	unsigned seed;

};

#endif //ifndef TOOLS_H

max_stars_repo_name: andrejadd/ABC-bee-opt

max_issues_repo_name: andrejadd/ABC-bee-opt

max_forks_repo_name: andrejadd/ABC-bee-opt


text:
#pragma once

#include <gsl/gsl_math.h>
#include <gsl/gsl_linalg.h>

#include "../vector.h"

namespace gsl_wrapper::bits
{
  class MatrixRow
  {
  public:
    // Consructor
    MatrixRow(gsl_vector_view view);

    // Operators
    operator ::gsl_wrapper::Vector() const;

    auto operator[](const size_t index) -> double &;
    auto operator[](const size_t index) const -> const double &;

  private:
    gsl_vector_view m_view;
  };

  inline MatrixRow::MatrixRow(gsl_vector_view view)
      : m_view{view}
  {
  }

  inline MatrixRow::operator ::gsl_wrapper::Vector() const
  {
    auto space = gsl_vector_calloc(m_view.vector.size);
    gsl_vector_memcpy(space, &m_view.vector);
    return ::gsl_wrapper::Vector(space);
  }

  inline auto MatrixRow::operator[](const size_t index) -> double &
  {
    return *gsl_vector_ptr(&m_view.vector, index);
  }

  inline auto MatrixRow::operator[](const size_t index) const -> const double &
  {
    return *gsl_vector_const_ptr(&m_view.vector, index);
  }

}


text:
#include <stdio.h>
#include <gsl/gsl_vector.h>

int
main (void)
{
  int i;
  gsl_vector * v = gsl_vector_alloc (10);

  {
     FILE * f = fopen ("test.dat", "r");
     gsl_vector_fscanf (f, v);
     fclose (f);
  }

  for (i = 0; i < 10; i++)
    {
      printf ("%g\n", gsl_vector_get(v, i));
    }

  gsl_vector_free (v);
  return 0;
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified


text:
/* rng/random.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

/* This file provides support for random() generators. There are three
   versions in widespread use today,

   - The original BSD version, e.g. on SunOS 4.1 and FreeBSD.

   - The Linux libc5 version, which is differs from the BSD version in
   its seeding procedure, possibly due to the introduction of a typo
   in the multiplier.

   - The GNU glibc2 version, which has a new (and better) seeding
   procedure.

   They all produce different numbers, due to the different seeding
   algorithms, but the algorithm for the generator is the same in each
   case.

 */

static inline long int random_get (int * i, int * j, int n, long int * x);

static inline unsigned long int random8_get (void *vstate);
static inline unsigned long int random32_get (void *vstate);
static inline unsigned long int random64_get (void *vstate);
static inline unsigned long int random128_get (void *vstate);
static inline unsigned long int random256_get (void *vstate);

static double random8_get_double (void *vstate);
static double random32_get_double (void *vstate);
static double random64_get_double (void *vstate);
static double random128_get_double (void *vstate);
static double random256_get_double (void *vstate);

static void random8_glibc2_set (void *state, unsigned long int s);
static void random32_glibc2_set (void *state, unsigned long int s);
static void random64_glibc2_set (void *state, unsigned long int s);
static void random128_glibc2_set (void *state, unsigned long int s);
static void random256_glibc2_set (void *state, unsigned long int s);

static void random8_libc5_set (void *state, unsigned long int s);
static void random32_libc5_set (void *state, unsigned long int s);
static void random64_libc5_set (void *state, unsigned long int s);
static void random128_libc5_set (void *state, unsigned long int s);
static void random256_libc5_set (void *state, unsigned long int s);

static void random8_bsd_set (void *state, unsigned long int s);
static void random32_bsd_set (void *state, unsigned long int s);
static void random64_bsd_set (void *state, unsigned long int s);
static void random128_bsd_set (void *state, unsigned long int s);
static void random256_bsd_set (void *state, unsigned long int s);

static void bsd_initialize (long int * x, int n, unsigned long int s);
static void libc5_initialize (long int * x, int n, unsigned long int s);
static void glibc2_initialize (long int * x, int n, unsigned long int s);

typedef struct
  {
    long int x;
  }
random8_state_t;

typedef struct
  {
    int i, j;
    long int x[7];
  }
random32_state_t;

typedef struct
  {
    int i, j;
    long int x[15];
  }
random64_state_t;

typedef struct
  {
    int i, j;
    long int x[31];
  }
random128_state_t;

typedef struct
  {
    int i, j;
    long int x[63];
  }
random256_state_t;

static inline unsigned long int
random8_get (void *vstate)
{
  random8_state_t *state = (random8_state_t *) vstate;

  state->x = (1103515245 * state->x + 12345) & 0x7fffffffUL;
  return state->x;
}

static inline long int
random_get (int * i, int * j, int n, long int * x)
{
  long int k ;

  x[*i] += x[*j] ;
  k = (x[*i] >> 1) & 0x7FFFFFFF ;

  (*i)++ ;
  if (*i == n)
    *i = 0 ;

  (*j)++ ;
  if (*j == n)
    *j = 0 ;

  return k ;
}

static inline unsigned long int
random32_get (void *vstate)
{
  random32_state_t *state = (random32_state_t *) vstate;
  unsigned long int k = random_get (&state->i, &state->j, 7, state->x) ;
  return k ;
}

static inline unsigned long int
random64_get (void *vstate)
{
  random64_state_t *state = (random64_state_t *) vstate;
  long int k = random_get (&state->i, &state->j, 15, state->x) ;
  return k ;
}

static inline unsigned long int
random128_get (void *vstate)
{
  random128_state_t *state = (random128_state_t *) vstate;
  unsigned long int k = random_get (&state->i, &state->j, 31, state->x) ;
  return k ;
}

static inline unsigned long int
random256_get (void *vstate)
{
  random256_state_t *state = (random256_state_t *) vstate;
  long int k = random_get (&state->i, &state->j, 63, state->x) ;
  return k ;
}

static double
random8_get_double (void *vstate)
{
  return random8_get (vstate) / 2147483648.0 ;
}

static double
random32_get_double (void *vstate)
{
  return random32_get (vstate) / 2147483648.0 ;
}

static double
random64_get_double (void *vstate)
{
  return random64_get (vstate) / 2147483648.0 ;
}

static double
random128_get_double (void *vstate)
{
  return random128_get (vstate) / 2147483648.0 ;
}

static double
random256_get_double (void *vstate)
{
  return random256_get (vstate) / 2147483648.0 ;
}

static void
random8_bsd_set (void *vstate, unsigned long int s)
{
  random8_state_t *state = (random8_state_t *) vstate;

  if (s == 0)
    s = 1;

  state->x = s;
}

static void
random32_bsd_set (void *vstate, unsigned long int s)
{
  random32_state_t *state = (random32_state_t *) vstate;
  int i;

  bsd_initialize (state->x, 7, s) ;

  state->i = 3;
  state->j = 0;

  for (i = 0 ; i < 10 * 7 ; i++)
    random32_get (state) ;
}

static void
random64_bsd_set (void *vstate, unsigned long int s)
{
  random64_state_t *state = (random64_state_t *) vstate;
  int i;

  bsd_initialize (state->x, 15, s) ;

  state->i = 1;
  state->j = 0;

  for (i = 0 ; i < 10 * 15 ; i++)
    random64_get (state) ;
}

static void
random128_bsd_set (void *vstate, unsigned long int s)
{
  random128_state_t *state = (random128_state_t *) vstate;
  int i;

  bsd_initialize (state->x, 31, s) ;

  state->i = 3;
  state->j = 0;

  for (i = 0 ; i < 10 * 31 ; i++)
    random128_get (state) ;
}

static void
random256_bsd_set (void *vstate, unsigned long int s)
{
  random256_state_t *state = (random256_state_t *) vstate;
  int i;

  bsd_initialize (state->x, 63, s) ;

  state->i = 1;
  state->j = 0;

  for (i = 0 ; i < 10 * 63 ; i++)
    random256_get (state) ;
}

static void
bsd_initialize (long int * x, int n, unsigned long int s)
{
  int i;

  if (s == 0)
    s = 1 ;

  x[0] = s;

  for (i = 1 ; i < n ; i++)
    x[i] = 1103515245 * x[i-1] + 12345 ;
}

static void
libc5_initialize (long int * x, int n, unsigned long int s)
{
  int i;

  if (s == 0)
    s = 1 ;

  x[0] = s;

  for (i = 1 ; i < n ; i++)
    x[i] = 1103515145 * x[i-1] + 12345 ;
}

static void
glibc2_initialize (long int * x, int n, unsigned long int s)
{
  int i;

  if (s == 0)
    s = 1 ;

  x[0] = s;

  for (i = 1 ; i < n ; i++)
    {
      const long int h = s / 127773;
      const long int t = 16807 * (s - h * 127773) - h * 2836;
      if (t < 0)
        {
          s = t + 2147483647 ;
        }
      else
        {
          s = t ;
        }

    x[i] = s ;
    }
}

static void
random8_glibc2_set (void *vstate, unsigned long int s)
{
  random8_state_t *state = (random8_state_t *) vstate;

  if (s == 0)
    s = 1;

  state->x = s;
}

static void
random32_glibc2_set (void *vstate, unsigned long int s)
{
  random32_state_t *state = (random32_state_t *) vstate;
  int i;

  glibc2_initialize (state->x, 7, s) ;

  state->i = 3;
  state->j = 0;

  for (i = 0 ; i < 10 * 7 ; i++)
    random32_get (state) ;
}

static void
random64_glibc2_set (void *vstate, unsigned long int s)
{
  random64_state_t *state = (random64_state_t *) vstate;
  int i;

  glibc2_initialize (state->x, 15, s) ;

  state->i = 1;
  state->j = 0;

  for (i = 0 ; i < 10 * 15 ; i++)
    random64_get (state) ;
}

static void
random128_glibc2_set (void *vstate, unsigned long int s)
{
  random128_state_t *state = (random128_state_t *) vstate;
  int i;

  glibc2_initialize (state->x, 31, s) ;

  state->i = 3;
  state->j = 0;

  for (i = 0 ; i < 10 * 31 ; i++)
    random128_get (state) ;
}

static void
random256_glibc2_set (void *vstate, unsigned long int s)
{
  random256_state_t *state = (random256_state_t *) vstate;
  int i;

  glibc2_initialize (state->x, 63, s) ;

  state->i = 1;
  state->j = 0;

  for (i = 0 ; i < 10 * 63 ; i++)
    random256_get (state) ;
}

static void
random8_libc5_set (void *vstate, unsigned long int s)
{
  random8_state_t *state = (random8_state_t *) vstate;

  if (s == 0)
    s = 1;

  state->x = s;
}

static void
random32_libc5_set (void *vstate, unsigned long int s)
{
  random32_state_t *state = (random32_state_t *) vstate;
  int i;

  libc5_initialize (state->x, 7, s) ;

  state->i = 3;
  state->j = 0;

  for (i = 0 ; i < 10 * 7 ; i++)
    random32_get (state) ;
}

static void
random64_libc5_set (void *vstate, unsigned long int s)
{
  random64_state_t *state = (random64_state_t *) vstate;
  int i;

  libc5_initialize (state->x, 15, s) ;

  state->i = 1;
  state->j = 0;

  for (i = 0 ; i < 10 * 15 ; i++)
    random64_get (state) ;
}

static void
random128_libc5_set (void *vstate, unsigned long int s)
{
  random128_state_t *state = (random128_state_t *) vstate;
  int i;

  libc5_initialize (state->x, 31, s) ;

  state->i = 3;
  state->j = 0;

  for (i = 0 ; i < 10 * 31 ; i++)
    random128_get (state) ;
}

static void
random256_libc5_set (void *vstate, unsigned long int s)
{
  random256_state_t *state = (random256_state_t *) vstate;
  int i;

  libc5_initialize (state->x, 63, s) ;

  state->i = 1;
  state->j = 0;

  for (i = 0 ; i < 10 * 63 ; i++)
    random256_get (state) ;
}

static const gsl_rng_type random_glibc2_type =
{"random-glibc2",                       /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random128_state_t),
 &random128_glibc2_set,
 &random128_get,
 &random128_get_double};

static const gsl_rng_type random8_glibc2_type =
{"random8-glibc2",                      /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random8_state_t),
 &random8_glibc2_set,
 &random8_get,
 &random8_get_double};

static const gsl_rng_type random32_glibc2_type =
{"random32-glibc2",                     /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random32_state_t),
 &random32_glibc2_set,
 &random32_get,
 &random32_get_double};

static const gsl_rng_type random64_glibc2_type =
{"random64-glibc2",                     /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random64_state_t),
 &random64_glibc2_set,
 &random64_get,
 &random64_get_double};

static const gsl_rng_type random128_glibc2_type =
{"random128-glibc2",                    /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random128_state_t),
 &random128_glibc2_set,
 &random128_get,
 &random128_get_double};

static const gsl_rng_type random256_glibc2_type =
{"random256-glibc2",                    /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random256_state_t),
 &random256_glibc2_set,
 &random256_get,
 &random256_get_double};

static const gsl_rng_type random_libc5_type =
{"random-libc5",                        /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random128_state_t),
 &random128_libc5_set,
 &random128_get,
 &random128_get_double};

static const gsl_rng_type random8_libc5_type =
{"random8-libc5",                       /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random8_state_t),
 &random8_libc5_set,
 &random8_get,
 &random8_get_double};

static const gsl_rng_type random32_libc5_type =
{"random32-libc5",                      /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random32_state_t),
 &random32_libc5_set,
 &random32_get,
 &random32_get_double};

static const gsl_rng_type random64_libc5_type =
{"random64-libc5",                      /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random64_state_t),
 &random64_libc5_set,
 &random64_get,
 &random64_get_double};

static const gsl_rng_type random128_libc5_type =
{"random128-libc5",                     /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random128_state_t),
 &random128_libc5_set,
 &random128_get,
 &random128_get_double};

static const gsl_rng_type random256_libc5_type =
{"random256-libc5",                     /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random256_state_t),
 &random256_libc5_set,
 &random256_get,
 &random256_get_double};

static const gsl_rng_type random_bsd_type =
{"random-bsd",                  /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random128_state_t),
 &random128_bsd_set,
 &random128_get,
 &random128_get_double};

static const gsl_rng_type random8_bsd_type =
{"random8-bsd",                 /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random8_state_t),
 &random8_bsd_set,
 &random8_get,
 &random8_get_double};

static const gsl_rng_type random32_bsd_type =
{"random32-bsd",                        /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random32_state_t),
 &random32_bsd_set,
 &random32_get,
 &random32_get_double};

static const gsl_rng_type random64_bsd_type =
{"random64-bsd",                        /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random64_state_t),
 &random64_bsd_set,
 &random64_get,
 &random64_get_double};

static const gsl_rng_type random128_bsd_type =
{"random128-bsd",               /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random128_state_t),
 &random128_bsd_set,
 &random128_get,
 &random128_get_double};

static const gsl_rng_type random256_bsd_type =
{"random256-bsd",               /* name */
 0x7fffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (random256_state_t),
 &random256_bsd_set,
 &random256_get,
 &random256_get_double};

const gsl_rng_type *gsl_rng_random_libc5    = &random_libc5_type;
const gsl_rng_type *gsl_rng_random8_libc5   = &random8_libc5_type;
const gsl_rng_type *gsl_rng_random32_libc5  = &random32_libc5_type;
const gsl_rng_type *gsl_rng_random64_libc5  = &random64_libc5_type;
const gsl_rng_type *gsl_rng_random128_libc5 = &random128_libc5_type;
const gsl_rng_type *gsl_rng_random256_libc5 = &random256_libc5_type;

const gsl_rng_type *gsl_rng_random_glibc2    = &random_glibc2_type;
const gsl_rng_type *gsl_rng_random8_glibc2   = &random8_glibc2_type;
const gsl_rng_type *gsl_rng_random32_glibc2  = &random32_glibc2_type;
const gsl_rng_type *gsl_rng_random64_glibc2  = &random64_glibc2_type;
const gsl_rng_type *gsl_rng_random128_glibc2 = &random128_glibc2_type;
const gsl_rng_type *gsl_rng_random256_glibc2 = &random256_glibc2_type;

const gsl_rng_type *gsl_rng_random_bsd    = &random_bsd_type;
const gsl_rng_type *gsl_rng_random8_bsd   = &random8_bsd_type;
const gsl_rng_type *gsl_rng_random32_bsd  = &random32_bsd_type;
const gsl_rng_type *gsl_rng_random64_bsd  = &random64_bsd_type;
const gsl_rng_type *gsl_rng_random128_bsd = &random128_bsd_type;
const gsl_rng_type *gsl_rng_random256_bsd = &random256_bsd_type;

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
//-*-C++-*-
/***************************************************************************
 *
 *   Copyright (C) 2010 by Paul Demorest
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/

#ifndef __Pulsar_SplineFit_h
#define __Pulsar_SplineFit_h

#include "Estimate.h"
#include "Reference.h"

#include <gsl/gsl_bspline.h>
#include <gsl/gsl_matrix.h>

//! Spline fitting for smoothing and/or interpolation, using
//! GSL's bspline routines.
class SplineFit : public Reference::Able {

public:

  //! Default constructor
  SplineFit ();

  //! Destructor
  virtual ~SplineFit ();

  //! Clear all current data, results
  void reset();

  //! Set the degree of the fit
  void set_order(int n) { order=n; calculated=false; }

  //! Get the current degree
  int get_order() { return order; }

  //! Set uniform breakpoints to span the data
  void set_uniform_breaks(int nint);

  //! Add a data point
  void add_data(double x, Estimate<double> y);

  //! Compute the fit using current data
  void compute();

  //! Evaluate the fit solution at the given x
  double evaluate(double x);

  //! Evaluate the fit solution's derivative at the given x
  double evaluate_deriv(double x);

  //! Get the reduced chi2 of the fit
  double get_rchi2();

protected:

  //! The x values for the fit
  std::vector<double> x;

  //! The y values/errors for the fit
  std::vector< Estimate<double> > y;

  //! The spline breakpoints
  std::vector<double> bp;

  //! Spline order (0=const, 3=cubic, etc)
  int order;

  //! The fit chi2
  double chi2;

  //! Fit NDOF
  int ndof;

  //! Has the fit been calculated?
  bool calculated;

  //! The fitted coeffs
  gsl_vector *coeffs;

  //! The fit cov matrix
  gsl_matrix *cov;

  //! Check if a requested x val is in the fit range
  bool check_range(double x);

  //! Check if spline intervals and data make sense
  void interval_check(bool fix=false);

  //! Free spline workspaces
  void free_workspaces();

  //! bspline temp space
  gsl_bspline_workspace *bwork;

private:

};

#endif


text:
/* -*- Mode: C; c-basic-offset:4 ; indent-tabs-mode:nil ; -*- */
/** \file
 * jitconv "doRef" calculation - a simplified gemm impl.
 *
 * This should be correct, and may even be somewhat fast "in general".
 * Do not expect "fastest" performance.
 *
 * libvednn has a more sophisticated GEMM convolution.
 */
#include "vednn_helper.h"
#include "convolution_gemm.h"
#include <string.h>
#include <stdlib.h>
#include <cblas.h>

#define LOCAL_FTRACE 1
#if LOCAL_FTRACE
#include "conv_test_param.h" // just for FTRACE macros
#define LFTRACE_BEGIN(...) FTRACE_BEGIN(__VA_ARGS__)
#define LFTRACE_END(...) FTRACE_END(__VA_ARGS__)
#define LFTRACE_IF(...) FTRACE_IF(__VA_ARGS__)
#else
#define LFTRACE_BEGIN(...) do{}while(0)
#define LFTRACE_END(...) do{}while(0)
#define LFTRACE_IF(...) do{}while(0)
#endif // LOCAL_FTRACE

// Scratchpad: none -- hardwired to malloc/free im2col buffers
//                     slower, but safer.

#define GEMM_PARA_THRESH 32768

// Around Jan. 2022, ncc-3.4.20, SGEMM began to segfault for M==1,
// so a workaround writes out the trivial matrix multiply.
/// 0 : SGEMM works great
/// 1 : safe macros
/// 2: dev code (generic)
#define SGEMM_M1_SEGFAULTS 1
/// 0 : old code, bias via sgemm (w/ problems)
/// 1 : just use cblas_saxpy alternate (simpler and avoids sgemm bug!)
#define BIAS_SAXPY 0

#if SGEMM_M1_SEGFAULTS==0 // issues ncc-3.4.20 and M==1 ???
#define SGEMM sgemm_
#define SGEMM_A1B0 sgemm_
#define SGEMM_A1B1K1 sgemm_
#define SGEMM_A1B0t sgemm_
#define SGEMM_A1tB1 sgemm_
#else
#define SGEMM sgemm_ // dangerous now - require bug workaround for M=1
#define SGEMM_A1B0   SGEMM_SAFE_A1B0
#define SGEMM_A1B1K1 SGEMM_SAFE_A1B1K1
#define SGEMM_A1B0t SGEMM_SAFE_A1B0t
#define SGEMM_A1tB1 SGEMM_SAFE_A1tB1
#endif

/// A workaround for SGEMM M==1 segfaults. (circa ncc-3.4.20, Jan 2022)
/// For alpha=1, beta=0 (main gemm calculation)
#define SGEMM_SAFE_A1B0(TRANSA,TRANSB, N,M,K, ALPHA,A,LDA, B,LDB, BETA,C,LDC) \
  do { \
    if(*(M) > 1){ \
      sgemm_(TRANSA,TRANSB, N,M,K, ALPHA,A,LDA, B,LDB, BETA,C,LDC); \
    }else{ \
      int const NN = *(N); \
      /*int const MM = *(M);*/ \
      int const KK = *(K); \
      if(*(M) == 1 && *(K) > 1){ /* using just M==1 */ \
        _Pragma("omp parallel if(NN * KK > 32768)") /* C99 */ \
        for (int n=0; n < (NN); ++n) { \
          float acc = 0.0f; \
          for (int k=0; k < (KK); ++k) { \
            acc += (A)[k * (NN) + n] * (B)[k]; \
          } \
          (C)[n] = acc; /* M==1 && beta==0.0 : no accumulation into C */ \
        } \
      }else{ /* M=1, K=1 */ \
        _Pragma("omp parallel if((NN) > 32768)") \
        for (int n=0; n < *(N); ++n) { \
          (C)[n] = (A)[n] * (B)[0]; \
        } \
      } \
    } \
  }while(0)
/// Backward Data also has alpha=1, beta=0, but B[] is transposed
// XXX test with jitconv -T BackwardData
#define SGEMM_SAFE_A1B0t(TRANSA,TRANSB, N,M,K, ALPHA,A,LDA, B,LDB, BETA,C,LDC) \
  do { \
    if(*(M) > 1){ \
      sgemm_(TRANSA,TRANSB, N,M,K, ALPHA,A,LDA, B,LDB, BETA,C,LDC); \
    }else{ \
      /* for M=1, B is K x 1, so vector ignoring the transpose is OK */ \
      int const NN = *(N); \
      /*int const MM = *(M);*/ \
      int const KK = *(K); \
      if(*(M) == 1 && *(K) > 1){ /* using just M==1 */ \
        _Pragma("omp parallel if(NN * KK > 32768)") /* C99 */ \
        for (int n=0; n < (NN); ++n) { \
          float acc = 0.0f; \
          for (int k=0; k < (KK); ++k) { \
            acc += (A)[k * (NN) + n] * (B)[k]; \
          } \
          (C)[n] = acc; /* M==1 && beta==0.0 : no accumulation into C */ \
        } \
      }else{ /* M=1, K=1 */ \
        _Pragma("omp parallel if((NN) > 32768)") \
        for (int n=0; n < *(N); ++n) { \
          (C)[n] = (A)[n] * (B)[0]; \
        } \
      } \
    } \
  }while(0)
/// for BackwardFilter, alpha=1, beta=1 and A is transposed
// XXX test with jitconv -T BackwardFilter
// try only a single omp || ?
// needs testing of ALL impls (update jitconv testBackwardFilter!!!)
#define SGEMM_SAFE_A1tB1_0(TRANSA,TRANSB, N,M,K, ALPHA,A,LDA, B,LDB, BETA,C,LDC) \
  do { \
    /*printf(" A1tB1 N=%d M=%d K=%d\n",*(N),*(M),*(K)); fflush(stdout);*/ \
    if(*(M) > 1){ \
      sgemm_(TRANSA,TRANSB, N,M,K, ALPHA,A,LDA, B,LDB, BETA,C,LDC); \
    }else{ \
      /* for M=1, A is N x K, so need A transpose wrt. Forward impl */ \
      int const NN = *(N); \
      int const KK = *(K); \
      int const NNKK = NN*KK; \
      /* M==1, any K */ \
      _Pragma("omp parallel if(NNKK > 32768)") /* C99 */ \
      for (int n=0; n < (NN); ++n) { \
        float acc = 0.0f; \
        for (int k=0; k < (KK); ++k) { \
          acc += (A)[n * (KK) + k] * (B)[k]; \
        } \
        (C)[n] += acc; /* beta=1 accumulation into C */ \
      } \
    } \
  }while(0)
#define SGEMM_SAFE_A1tB1(TRANSA,TRANSB, N,M,K, ALPHA,A,LDA, B,LDB, BETA,C,LDC) \
  do { \
    /*printf(" A1tB1 N=%d M=%d K=%d\n",*(N),*(M),*(K)); fflush(stdout);*/ \
    if(*(M) > 1){ \
      sgemm_(TRANSA,TRANSB, N,M,K, ALPHA,A,LDA, B,LDB, BETA,C,LDC); \
    }else{ \
      /* for M=1, A is N x K, so need A transpose wrt. Forward impl */ \
      int const NN = *(N); \
      int const KK = *(K); \
      if(*(M) == 1 && *(K) > 1){ /* using just M==1 */ \
        /*_Pragma("omp parallel if(NN * KK > 32768)")*/ /* C99 */ \
        for (int n=0; n < (NN); ++n) { \
          float acc = 0.0f; \
          for (int k=0; k < (KK); ++k) { \
            acc += (A)[n * (KK) + k] * (B)[k]; \
          } \
          (C)[n] += acc; /* beta=1 accumulation into C */ \
        } \
      }else{ /* M=1, K=1 */ \
        /*_Pragma("omp parallel if((NN) > 32768)")*/ \
        for (int n=0; n < (NN); ++n) { \
          (C)[n] += (A)[n] * (B)[0]; /* beta=1 accum */ \
        } \
      } \
    } \
  }while(0)
// using M==1 and K==1
// here A[N] is 1.0
/// workaround for M=1 bias segfault.
/// Here K=1, alpha=1, beta=1, \b and A[] is all-1.0  (bias accumulation)
#define SGEMM_SAFE_A1B1K1(TRANSA,TRANSB, N,M,K, ALPHA,A,LDA, B,LDB, BETA,C,LDC) \
  do { \
    if(1 && *(M) > 1) /* always elide? */ \
    { \
      sgemm_(TRANSA,TRANSB, N,M,K, ALPHA,A,LDA, B,LDB, BETA,C,LDC); \
    }else{ \
      int const MN = *(M) * *(N); \
      float const B_0 = *(B);/* (B)[0] */ \
      /* wrong output if try to parallelize? */ \
      /* _Pragma("omp parallel if(MN > 32768)") */ \
      for (int mn=0; mn < MN; ++mn) { \
        (C)[mn] += B_0; \
      } \
    } \
  }while(0)

#if 0
#define DBG(...) do{printf(__VA_ARGS__);fflush(stdout);}while(0)
#else
#define DBG(...)
#endif

#if 1
void sgemm_(char *TRANSA, char *TRANSB, int *M, int *N, int *K,
    float *ALPHA, float *A,  int *LDA, float *B, int *LDB,
    float *BETA, float *C, int *LDC ) ;
//void cblas_saxpy(const int N, const float alpha, const float *X,
//                 const int incX, float *Y, const int incY);
#endif

static char  TRANS   = 'T';
static char  NOTRANS = 'N';
static float FONE    = 1.0f;
static float FZERO   = 0.0f;
static int   IONE    = 1;

/* ----------------------------------------------------------------------- */
static inline int is_a_ge_zero_and_a_lt_b(int a, int b) {
  //return (unsigned)a < (unsigned)b;
  return a>=0  && a<b; // for ncc auto vectorization, this is better
}

  static void
#if 0
im2col_cpu(const float * restrict data_im, const int channels,
    const int height, const int width, const int kernel_h, const int kernel_w,
    const int output_h, const int output_w,
    const int pad_h, const int pad_w,
    const int stride_h, const int stride_w,
    const int dilation_h, const int dilation_w,
    float * restrict data_col)
#else
im2col_cpu(const float * data_im, const int channels,
    const int height, const int width, const int kernel_h, const int kernel_w,
    const int output_h, const int output_w,
    const int pad_h, const int pad_w,
    const int stride_h, const int stride_w,
    const int dilation_h, const int dilation_w,
    float * data_col)
#endif
{
  LFTRACE_BEGIN("im2col_cpu");
#if 0
  const int output_h = (height + 2 * pad_h - (dilation_h * (kernel_h - 1) + 1)) / stride_h + 1;
  const int output_w = (width + 2 * pad_w -  (dilation_w * (kernel_w - 1) + 1)) / stride_w + 1;
#endif
  const int channel_size = height * width;

  int channel;
  if (0){
    // Note: alloc is for ic * ih * kw, but probably only need (ic/g) * kh * kw here?
    DBG("im2col_cpu _alloc needs %llu floats\n",(long long)channels*channel_size);
    for(size_t i=0; i<(size_t)channels*(size_t)channel_size; ++i){
      data_col[i] = 0.0f;
    }
    DBG("data_col / pColBuf accessible");
  }

//#pragma omp parallel for if(channels>=3)
  for (channel = 0 ; channel < channels; channel++) {       // inChannel
    //printf(" i2c c%d/%d",channel,channels); fflush(stdout);
    int kernel_row, kernel_col, output_rows, output_cols, output_col;

    int inOffset = channel * channel_size;
    int outOffset = channel * output_h * output_w * kernel_h * kernel_w;

    for (kernel_row = 0; kernel_row < kernel_h; kernel_row++) {   // kernHeight
      for (kernel_col = 0; kernel_col < kernel_w; kernel_col++) {   // kernWidth
        int input_row = -pad_h + kernel_row * dilation_h;
        for (output_rows = output_h; output_rows; output_rows--) {  // outHeight
          if (!is_a_ge_zero_and_a_lt_b(input_row, height)) {
            for (output_cols = output_w; output_cols; output_cols--) {
              // *(data_col++) = 0;
              data_col[outOffset++] = 0.f;
            }
          } else {
            int input_col = -pad_w + kernel_col * dilation_w;
            // following still bombed
            //#pragma _NEC novector
            for (output_col = output_w; output_col; output_col--) { // outWidth
#if 1 // newer
              data_col[outOffset++] //*(data_col++)
                = (is_a_ge_zero_and_a_lt_b(input_col, width)
                    ? data_im[inOffset + input_row * width + input_col]
                    : 0.f);
#else // older
              if (outOffset < 0 || outOffset >= channels*kernel_h*kernel_w*output_h*output_w){
                printf("ERROR: outOffset"); fflush(stdout); exit(-1);
              }
              if (is_a_ge_zero_and_a_lt_b(input_col, width)) {
                // *(data_col++) = data_im[input_row * width + input_col];
                data_col[outOffset] = data_im[inOffset + input_row * width + input_col];
              } else {
                // *(data_col++) = 0;
                data_col[outOffset] = 0;
              }
              ++outOffset;
#endif
              input_col += stride_w;
            }
          }
          input_row += stride_h;
        }
      }
    }
  }
  LFTRACE_END("im2col_cpu");
}

static void
col2im_cpu(
    const float* data_col, const int channels,
    const int height, const int width, const int kernel_h, const int kernel_w,
    const int output_h, const int output_w,
    const int pad_h, const int pad_w,
    const int stride_h, const int stride_w,
    const int dilation_h, const int dilation_w,
    float* data_im) {

  LFTRACE_BEGIN("col2im_cpu");
  memset(data_im, 0, sizeof(float)*height*width*channels) ;

#if 0
  const int output_h = (height + 2 * pad_h - (dilation_h * (kernel_h - 1) + 1)) / stride_h + 1;
  const int output_w = (width + 2 * pad_w -  (dilation_w * (kernel_w - 1) + 1)) / stride_w + 1;
#endif

  const int channel_size = height * width;

  int channel;

#pragma omp parallel for if(channels>=3)
  for (channel = 0 ; channel < channels; channel++) {       // inChannel
    int kernel_row, kernel_col, output_rows, output_cols, output_col;

    int inOffset = channel * channel_size;
    int outOffset = channel * output_h * output_w * kernel_h * kernel_w;

    for (kernel_row = 0; kernel_row < kernel_h; kernel_row++) {   // kernHeight
      for (kernel_col = 0; kernel_col < kernel_w; kernel_col++) {   // kernWidth
        int input_row = -pad_h + kernel_row * dilation_h;
        for (output_rows = output_h; output_rows; output_rows--) {  // outHeight
          if (!is_a_ge_zero_and_a_lt_b(input_row, height)) {
            for (output_cols = output_w; output_cols; output_cols--) {
              data_col[outOffset++] ;
            }
          } else {
            int input_col = -pad_w + kernel_col * dilation_w;
            for (output_col = output_w; output_col; output_col--) { // outWidth
              if (is_a_ge_zero_and_a_lt_b(input_col, width)) {
                data_im[inOffset + input_row * width + input_col] += data_col[outOffset++] ;
              } else {
                outOffset++ ;
              }
              input_col += stride_w;
            }
          }
          input_row += stride_h;
        }
      }
    }
  }
  LFTRACE_END("col2im_cpu");
}

// pOne is oh*ow of 1.0f
// pColBuff is scratch of ic*kw*kh * ow*oh * iw*ih (huge, in this version) see conv_test_param.c
  vednnError_t
convolution_forward_gemm(
    const vednnTensorParam_t * restrict pParamIn, const void * restrict pDataIn,
    const vednnFilterParam_t * restrict pParamKernel, const void * restrict pDataKernel,
    const vednnBiasParam_t * restrict pParamBias, const void * restrict pDataBias,
    const vednnTensorParam_t * restrict pParamOut, void * restrict pDataOut,
    //const float * restrict pOne,  float * restrict pColBuff,
    const float * pOne,  float * pColBuff,
    const vednnConvolutionParam_t * restrict pParamConv )
{
  LFTRACE_BEGIN("convolution_forward_gemm");

  int batch   = pParamIn->batch;
  int inChannel = pParamIn->channel;
  int inWidth   = pParamIn->width;
  int inHeight  = pParamIn->height;
  int outChannel  = pParamOut->channel;
  int outWidth  = pParamOut->width;
  int outHeight = pParamOut->height;
  int kernWidth = pParamKernel->width;
  int kernHeight  = pParamKernel->height;

  int group   = pParamConv->group;
  int strideWidth = pParamConv->strideWidth;;
  int strideHeight  = pParamConv->strideHeight;
  int padWidth  = pParamConv->padWidth;
  int padHeight = pParamConv->padHeight;
  int dilationWidth = pParamConv->dilationWidth;
  int dilationHeight  = pParamConv->dilationHeight;

  int inChannelGroup  = inChannel  / group; // pParamKernel->inChannel
  int outChannelGroup = outChannel / group; // pParamKernel->outChannel

  int no_im2col = (kernWidth == 1 && kernHeight == 1 && strideWidth == 1 && strideHeight == 1 && padWidth == 0 && padHeight == 0);
  if (0){
    printf("mb=%d g=%d ic=%d ic/g=%d noi2c=%d", batch,group,inChannel,inChannelGroup, no_im2col);
    fflush(stdout);
    float const* pk = (float const*)(pDataKernel); // assume float (debug) XXX
    int const ksz = group * outChannelGroup * inChannelGroup * kernHeight * kernWidth;
    size_t const ksz2 = getKernelSize(pParamKernel) * pParamConv->group;
    if (ksz != ksz2) {
      printf("ksz %d != ksz2 %d\n");
    }
    for(int i=0; i<ksz; ++i){
      if (isnan(pk[i])) {
        printf("generateRandomData --> nans!\n");
        printf("ksz=%d i=%d\n");
        exit(-1);
      }
      if (pk[i] < -5.0 || pk[i] > +5.0){
        printf("generateRandomData --> outside [-5.0,5.0]\n");
        printf("ksz=%d i=%d\n");
        exit(-1);
      }
    }
    printf("input pDataKernel[0..%d - 1] looks good\n",ksz);
    printf("input pDataKernel[0..%lu - 1] looks good\n",(long unsigned)ksz2);
  }

  float * transformed_filter = NULL ;
  if( pParamKernel->layout == VEDNN_FILTER_LAYOUT_HWCN ) { // only support group=1
    if (group!=1){
      printf("Unsupported ref calc: HWCN wants group==1");
      exit(-1);
      //return VEDNN_ERROR_INVALID_PARAM;
    }

    const int N = outChannel ;
    const int C = inChannel ;
    const int H = kernHeight ;
    const int W = kernWidth ;

    float * filter = (float *) pDataKernel ;
    transformed_filter = (float *) malloc(sizeof(float)*N*C*H*W) ;
#pragma omp parallel for
    for(int n=0; n<N ; n++) {
      for(int c=0; c<C ; c++) {
        for(int hw=0; hw<H*W ; hw++) {
          transformed_filter[((n*C+c)*H)*W+hw] = filter[((hw)*C+c)*N+n] ;
        }
      }
    }
  }

  const float * restrict pIn     = pDataIn;
  const float * restrict pBias   = pDataBias;
  const float * restrict pKernel = transformed_filter == NULL ? pDataKernel : transformed_filter ;
  float * restrict pOut    = pDataOut;

  for (int n = 0; n < batch; n++) { // this->num_
    int inBatchOffset  = n * inChannel  * inWidth  * inHeight;
    int outBatchOffset = n * outChannel * outWidth * outHeight;

    for (int g = 0; g < group; g++) {
      int inGroupOffset   = g * inChannelGroup                   * inHeight   * inWidth;
      int outGroupOffset  = g * outChannelGroup                  * outHeight  * outWidth;
      int kernGroupOffset = g * outChannelGroup * inChannelGroup * kernHeight * kernWidth;
      int biasGroupOffset = g * outChannelGroup;

      int inOffset  = inBatchOffset  + inGroupOffset;
      int outOffset = outBatchOffset + outGroupOffset;

      if (no_im2col) {
        int M = outChannelGroup;
        int N = outWidth * outHeight;
        int K = inChannelGroup;
        int LDA = inWidth * inHeight;
        //printf(" M=%d N=%d K=%d ",M,N,K);

#if SGEMM_M1_SEGFAULTS==0 // issues ncc-3.4.20 and M==1 ???
        SGEMM(&NOTRANS, &NOTRANS, &N, &M, &K,
            &FONE, (float *) &pIn[inOffset], &LDA,
            (float *) &pKernel[kernGroupOffset], &K,
            &FZERO, &pOut[outOffset], &N);
#elif SGEMM_M1_SEGFAULTS==1 // issues ncc-3.4.20 and M==1 ???
        SGEMM_A1B0(&NOTRANS, &NOTRANS, &N, &M, &K,
            &FONE, (float *) &pIn[inOffset], &LDA,
            (float *) &pKernel[kernGroupOffset], &K,
            &FZERO, &pOut[outOffset], &N);
#else // ncc-3.4.20 fixup...  Here I show general equivalent for the SGEMM...
        if (M>1) { // || K>1) {
          SGEMM(&NOTRANS, &NOTRANS, &N, &M, &K,
              &FONE, (float *) &pIn[inOffset], &LDA, // LDA=N for no_im2col
              (float *) &pKernel[kernGroupOffset], &K, // LDB=K for no_im2col
              &FZERO, &pOut[outOffset], &N);
        }else{
          // M==1 has some BLAS segv !!! just write it out for now...
          // At -O3 and -O4, ncc should have -fassociative-math -fmatrix-multiply
          // and should emit matrix-multiply code for these
          // consulting sgemm docs
          //  3rd dim is the summation index
          float const* A = (float const*)&pIn[inOffset];            // size N x K
          float const* B = (float const*)&pKernel[kernGroupOffset]; // size K x M
          float * C = &pOut[outOffset];                             // size N x M
#if 0 // generic, long-hand matrix multiply
          if (M>1) { // actually this is NOT quite right yet :(
            for (int n=0; n<N; n++) {
              for (int m=0; m<M; m++) {
                C[m*N + n] = 0.0f;
              }
            }
            for (int n=0; n<N; n++) {
              for (int m=0; m<M; m++) {
                float acc = 0.0f;
                for (int k=0; k<K; k++) {
                  acc += A[k*N + n] * B[m*K + k];
                }
                C[m*N + n] += acc; // beta=0
              }
            }
          }else
#endif
#if 1
            if(M==1 && K>1){ // using just M==1
#pragma omp parallel if(N*K>GEMM_PARA_THRESH)
              for (int n=0; n<N; n++) {
                float acc = 0.0f;
                for (int k=0; k<K; k++) {
                  acc += A[k*N+n] * B[k];
                }
                C[n] = acc; // beta=0, no accumulation into C
              }
            }else
#endif
          { // using M==1 and K==1
            float b0 = B[0];
#pragma omp parallel if(N*K>GEMM_PARA_THRESH)
            for (int n=0; n<N; n++) {
              C[n] = A[n] * b0;
            }
          }
        }
#endif

        if (pBias) {
          //printf("noi2c pBias M=%d N=%d K=%d ", M,N,K);
#if SGEMM_M1_SEGFAULTS==0 // issues ncc-3.4.20 and M==1 ???
          SGEMM(&NOTRANS, &NOTRANS, &N, &M, &IONE,
              &FONE, (float *) pOne, &N,
              (float *) &pBias[biasGroupOffset], &IONE,
              &FONE, &pOut[outOffset], &N);
#elif SGEMM_M1_SEGFAULTS==1 // issues ncc-3.4.20 and M==1 ???
#if BIAS_SAXPY==0
          SGEMM_A1B1K1(&NOTRANS, &NOTRANS, &N, &M, &IONE,
              &FONE, (float *) pOne, &N,                // N x 1
              (float *) &pBias[biasGroupOffset], &IONE, // 1 x M
              &FONE, &pOut[outOffset], &N);             // N x M
#else
          // note that this might be formulated as a
          // SAXPY( N=MN,
          //        ALPHA=1.0,
          //        X=&pBias[biasGroupOffset],
          //        INCX=0, /* <-- "add constant" */
          //        Y=&pOut[outOffset], /* add to  this */
          //        INCY=1
          //      )
          // Unfortunately, this only handles for M=1
          if (M>1) {
            SGEMM(&NOTRANS, &NOTRANS, &N, &M, &IONE,
                &FONE, (float *) pOne, &N,
                (float *) &pBias[biasGroupOffset], &IONE,
                &FONE, &pOut[outOffset], &N);
          }else{
            //printf(" saxpy M=%d N=%d\n", M,N);
            cblas_saxpy( N, 1.0, &pBias[biasGroupOffset], 0,
                &pOut[outOffset], 1);
          }
#endif
#elif 0 // debug
          {
            // M==1 has some BLAS segv !!! just write it out for now...
            // K=1 summation index is a huge simplification
            //   We might want to never fully call the SGEMM!
            //float const* A = (float const*)pOne;            // size N x K
            float const* B = (float const*)&pBias[biasGroupOffset]; // size K x M
            float * C = &pOut[outOffset];                   // size N x M
            if (1) { // further simplificcation only elides 1 scalar multiply
              //for (int n=0; n<N; n++) for (int m=0; m<M; m++) C[m*N + n] = 0.0f;
#if 0 // version 0, all loops before simplification
              printf("x7");
              int const K = 1;
              for (int n=0; n<N; n++) {
                for (int m=0; m<M; m++) {
                  float acc = 0.0f;
                  for (int k=0; k<K; k++) {
                    acc += 1.0/* A[k*N + n] */ * B[m*K + k];
                  }
                  C[m*N + n] += acc; // beta=0
                }
              }
#elif 0 // K=1 is a drastic simplification
              printf("x8");
              for (int m=0; m<M; ++m) {
                for (int n=0; n<N; ++n) {
                  C[m*N + n] += B[0]; // beta=0
                }
              }
#else // not working with omp parllel
              DBG("x9");
              int const MN = M * N;
//#pragma omp parallel if(MN > GEMM_PARA_THRESH) /* this cause wrong output */
              for (int mn=0; mn < MN; ++mn) {
                C[mn] += B[0]; // beta=0
              }
            }
          }
#endif

#else // dev code, summarized
          // maybe it's faster to always elide the SGEMM?
          if (M>1) {
            SGEMM(&NOTRANS, &NOTRANS, &N, &M, &IONE,
                &FONE, (float *) pOne, &N,
                (float *) &pBias[biasGroupOffset], &IONE,
                &FONE, &pOut[outOffset], &N);
          }else{
            // workaround for M=1 segfault circa ncc 3.4.20
            // using M==1 and K==1
            // here A[N] is 1.0
            float b0 = pBias[biasGroupOffset+0]; // B[0]
            int const MN = M * N;
            /* _Pragma("omp parallel if(MN > 32768)") //wrong output? */
            for (int mn=0; mn < MN; ++mn) {
              pOut[outOffset+mn] += b0; // C[n]
            }
          }
#endif
        }// if pBias

      } else {

        int M = outChannelGroup;
        int N = outWidth * outHeight;
        int K = inChannelGroup * kernWidth * kernHeight;

        im2col_cpu(&pIn[inOffset],
            inChannelGroup, inHeight, inWidth, kernHeight, kernWidth, outHeight, outWidth,
            padHeight, padWidth, strideHeight, strideWidth, dilationHeight, dilationWidth,
            pColBuff);

#if SGEMM_M1_SEGFAULTS==0 // issues ncc-3.4.20 and M==1 ???
        SGEMM(&NOTRANS, &NOTRANS, &N, &M, &K,
            &FONE, pColBuff, &N,
            (float *)&pKernel[kernGroupOffset], &K,
            &FZERO, &pOut[outOffset], &N);
        // segfault if M==1, at least w/ ncc 3.4.20 etc.
#elif SGEMM_M1_SEGFAULTS==1 // issues ncc-3.4.20 and M==1 ???
        SGEMM_A1B0(&NOTRANS, &NOTRANS, &N, &M, &K,
            &FONE, pColBuff, &N,
            (float *)&pKernel[kernGroupOffset], &K,
            &FZERO, &pOut[outOffset], &N);
#else
        if (M>1) {
          SGEMM(&NOTRANS, &NOTRANS, &N, &M, &K,
              &FONE, pColBuff, &N,                    // N x K
              (float *)&pKernel[kernGroupOffset], &K, // K x M
              &FZERO, &pOut[outOffset], &N);          // N x M
        }else{
          // M==1 has some BLAS segv !!! just write it out for now...
          // At -O3 and -O4, ncc should have -fassociative-math -fmatrix-multiply
          // and should emit matrix-multiply code for these
          float const* A = pColBuff;                                  // osz x icg*ksz (?)
          float const* B = (float const*)&pKernel[kernGroupOffset];   // icg*ksz x ocg
          float      * C = &pOut[outOffset];                          // osz x ocg (?)
#if 0
          if (M>1) {
            for (int n=0; n<N; n++) {
              for (int m=0; m<M; m++) {
                C[m*N + n] = 0.0f;
              }
            }
            for (int n=0; n<N; n++) {
              for (int m=0; m<M; m++) {
                float acc = 0.0f;
                for (int k=0; k<K; k++) {
                  acc += A[k*N + n] * B[m*K + k];
                }
                C[m*N + n] += acc; // beta=0, no accumulation into C
              }
            }
          }else
#endif
          { // M==1
            for (int n=0; n<N; n++) {
              float acc = 0.0f;
              for (int k=0; k<K; k++) {
                acc += A[k*N+n] * B[k];
              }
              C[n] = acc;
            }
          }
        }
#endif
        //printf(" back from SGEMM..."); fflush(stdout);

        if (pBias) {
          //printf("i2c bias...\n"); fflush(stdout);
#if SGEMM_M1_SEGFAULTS==0 // issues ncc-3.4.20 and M==1 ???
          SGEMM(&NOTRANS, &NOTRANS, &N, &M, &IONE,
              &FONE, (float *)pOne, &N,
              (float *) &pBias[biasGroupOffset], &IONE,
              &FONE, &pOut[outOffset], &N);
#elif SGEMM_M1_SEGFAULTS==1 // issues ncc-3.4.20 and M==1 ???
#if 1 //BIAS_SAXPY==0
          SGEMM_A1B1K1(&NOTRANS, &NOTRANS, &N, &M, &IONE,
              &FONE, (float *)pOne, &N,                 // size N x 1
              (float *) &pBias[biasGroupOffset], &IONE, // size 1 x M
              &FONE, &pOut[outOffset], &N);             // size N x M
#else
          cblas_saxpy( M*N, 1.0, &pBias[biasGroupOffset], 0,
              &pOut[outOffset], 1);
#endif
#else
          if (M>1) {
            SGEMM(&NOTRANS, &NOTRANS, &N, &M, &IONE,
                &FONE, (float *) pOne, &N,
                (float *) &pBias[biasGroupOffset], &IONE,
                &FONE, &pOut[outOffset], &N);
          }else{
            // workaround for M=1 segfault circa ncc 3.4.20
            // using M==1 and K==1
            // here A[N] is 1.0
            float b0 = pBias[biasGroupOffset+0]; // B[0]
            for (int n=0; n<N; n++) {
              pOut[outOffset+n] += b0; // C[n]
            }
          }
#endif
        }
      } // no_im2col?
    } // group
  } // batch

  if( transformed_filter != NULL ) free(transformed_filter) ;

  LFTRACE_END("convolution_forward_gemm");
  return VEDNN_SUCCESS;
}

  vednnError_t
convolution_backward_data_gemm(
    const vednnTensorParam_t * restrict pParamGradOut, const void * restrict pDataGradOut,
    const vednnFilterParam_t * restrict pParamKernel, const void * restrict pDataKernel,
    const vednnTensorParam_t * restrict pParamGradIn, void * restrict pDataGradIn,
    float * restrict pColBuff,
    const vednnConvolutionParam_t * restrict pParamConv )
{
  LFTRACE_BEGIN("convolution_backward_data_gemm");
  int n, g;

  int batch   = pParamGradOut->batch;
  int gOutChannel = pParamGradOut->channel;
  int gOutWidth = pParamGradOut->width;
  int gOutHeight  = pParamGradOut->height;
  int gInChannel  = pParamGradIn->channel;
  int gInWidth  = pParamGradIn->width;
  int gInHeight = pParamGradIn->height;
  int kernWidth = pParamKernel->width;
  int kernHeight  = pParamKernel->height;

  int group   = pParamConv->group;
  int strideWidth = pParamConv->strideWidth;;
  int strideHeight  = pParamConv->strideHeight;
  int padWidth  = pParamConv->padWidth;
  int padHeight = pParamConv->padHeight;
  int dilationWidth = pParamConv->dilationWidth;
  int dilationHeight  = pParamConv->dilationHeight;

  int gOutChannelGroup = gOutChannel  / group;
  int gInChannelGroup  = gInChannel / group;

  int no_im2col = (kernWidth == 1 && kernHeight == 1 && strideWidth == 1 && strideHeight == 1 && padWidth == 0 && padHeight == 0);

  float * transformed_filter = NULL ;
  if( pParamKernel->layout == VEDNN_FILTER_LAYOUT_HWCN ) { // only support group=1

    const int N = gOutChannel ;
    const int C = gInChannel ;
    const int H = kernHeight ;
    const int W = kernWidth ;

    float * filter = (float *) pDataKernel ;
    transformed_filter = (float *) malloc(sizeof(float)*N*C*H*W) ;
#pragma omp parallel for
    for(int n=0; n<N ; n++) {
      for(int c=0; c<C ; c++) {
        for(int hw=0; hw<H*W ; hw++) {
          transformed_filter[((n*C+c)*H)*W+hw] = filter[((hw)*C+c)*N+n] ;
        }
      }
    }
  }

  const float * restrict pGradOut = pDataGradOut;
  const float * restrict pKernel  = transformed_filter == NULL ? pDataKernel : transformed_filter ;
  float * restrict pGradIn  = pDataGradIn;

  for (n = 0; n < batch; n++) { // this->num_
    int gOutBatchOffset  = n * gOutChannel  * gOutWidth  * gOutHeight;
    int gInBatchOffset = n * gInChannel * gInWidth * gInHeight;

    for (g = 0; g < group; g++) {
      int gOutGroupOffset = g *                   gOutChannelGroup * gOutHeight * gOutWidth;
      int gInGroupOffset  = g * gInChannelGroup                    * gInHeight  * gInWidth;
      int kernGroupOffset = g * gInChannelGroup * gOutChannelGroup * kernHeight * kernWidth;

      int gOutOffset = gOutBatchOffset + gOutGroupOffset;
      int gInOffset  = gInBatchOffset  + gInGroupOffset;

      int M = gInChannelGroup * kernWidth * kernHeight;
      int N = gOutWidth * gOutHeight;
      int K = gOutChannelGroup;

      if( no_im2col ) {
        SGEMM_A1B0t(&NOTRANS, &TRANS, &N, &M, &K,
            &FONE, (float *) &pGradOut[gOutOffset], &N, // N x K
            (float *) &pKernel[kernGroupOffset], &M,    // M x K (trans!)
            &FZERO, &pGradIn[gInOffset], &N);           // N x M
      }
      else {
        SGEMM_A1B0t(&NOTRANS, &TRANS, &N, &M, &K,
            &FONE, (float *) &pGradOut[gOutOffset], &N,
            (float *) &pKernel[kernGroupOffset], &M,
            &FZERO, pColBuff, &N);

        col2im_cpu(pColBuff,
            gInChannelGroup, gInHeight, gInWidth, kernHeight, kernWidth, gOutHeight, gOutWidth,
            padHeight, padWidth, strideHeight, strideWidth, dilationHeight, dilationWidth,
            &pGradIn[gInOffset]);
      }
    } // group
  } // batch

  if( transformed_filter != NULL ) free(transformed_filter) ;

  LFTRACE_END("convolution_backward_data_gemm");
  return VEDNN_SUCCESS;
}

  vednnError_t
convolution_backward_filter_gemm(
    const vednnTensorParam_t * restrict pParamIn, const void * restrict pDataIn,
    const vednnTensorParam_t * restrict pParamGradOut, const void * restrict pDataGradOut,
    const vednnFilterParam_t * restrict pParamGradKernel, void * restrict pDataGradKernel,
    float * restrict pColBuff,
    const vednnConvolutionParam_t * restrict pParamConv )
{
  LFTRACE_BEGIN("convolution_backward_filter_gemm");
  int n, g;

  int batch   = pParamIn->batch;
  int inChannel = pParamIn->channel;
  int inWidth   = pParamIn->width;
  int inHeight  = pParamIn->height;
  int outChannel  = pParamGradOut->channel;
  int outWidth  = pParamGradOut->width;
  int outHeight = pParamGradOut->height;
  int kernWidth = pParamGradKernel->width;
  int kernHeight  = pParamGradKernel->height;

  int group   = pParamConv->group;
  int strideWidth = pParamConv->strideWidth;;
  int strideHeight  = pParamConv->strideHeight;
  int padWidth  = pParamConv->padWidth;
  int padHeight = pParamConv->padHeight;
  int dilationWidth = pParamConv->dilationWidth;
  int dilationHeight  = pParamConv->dilationHeight;

  int inChannelGroup  = inChannel  / group; // pParamKernel->inChannel
  int outChannelGroup = outChannel / group; // pParamKernel->outChannel

  int no_im2col = (kernWidth == 1 && kernHeight == 1 && strideWidth == 1 && strideHeight == 1 && padWidth == 0 && padHeight == 0);

  float * transformed_filter = NULL ;
  if( pParamGradKernel->layout == VEDNN_FILTER_LAYOUT_HWCN ) { // only support group=1
    const int N = outChannel ;
    const int C = inChannel ;
    const int H = kernHeight ;
    const int W = kernWidth ;

    transformed_filter = (float *) malloc(sizeof(float)*N*C*H*W) ;
#pragma omp parallel for
    for(int i=0; i<N*C*H*W; i++) transformed_filter[i] = 0.f ;
  }

  const float * restrict pIn     = pDataIn;
  const float * restrict pOut    = pDataGradOut;
  float * restrict pKernel = transformed_filter == NULL ? pDataGradKernel : transformed_filter ;

  for (n = 0; n < batch; n++) { // this->num_
    int inBatchOffset  = n * inChannel  * inWidth  * inHeight;
    int outBatchOffset = n * outChannel * outWidth * outHeight;

    for (g = 0; g < group; g++) {
      int inGroupOffset   = g * inChannelGroup                   * inHeight   * inWidth;
      int outGroupOffset  = g * outChannelGroup                  * outHeight  * outWidth;
      int kernGroupOffset = g * outChannelGroup * inChannelGroup * kernHeight * kernWidth;

      int inOffset  = inBatchOffset  + inGroupOffset;
      int outOffset = outBatchOffset + outGroupOffset;

      if( no_im2col ) {
        int M = outChannelGroup;
        int N = inChannelGroup * kernWidth * kernHeight;
        int K = outWidth * outHeight;

        SGEMM_A1tB1(&TRANS, &NOTRANS, &N, &M, &K,
            &FONE,  (float*)&pIn[inOffset], &K,
            (float*)&pOut[outOffset], &K,
            &FONE, &pKernel[kernGroupOffset], &N);
      }
      else {
        im2col_cpu(&pIn[inOffset],
            inChannelGroup, inHeight, inWidth, kernHeight, kernWidth, outHeight, outWidth,
            padHeight, padWidth, strideHeight, strideWidth, dilationHeight, dilationWidth,
            pColBuff);

        int M = outChannelGroup;
        int N = inChannelGroup * kernWidth * kernHeight;
        int K = outWidth * outHeight;

        SGEMM_A1tB1(&TRANS, &NOTRANS, &N, &M, &K,
            &FONE,  pColBuff, &K,
            (float*)&pOut[outOffset], &K,
            &FONE, &pKernel[kernGroupOffset], &N);
      }
    } // group
  } // batch

  if( transformed_filter != NULL ) {
    const int N = outChannel ;
    const int C = inChannel ;
    const int H = kernHeight ;
    const int W = kernWidth ;

    float * filter = (float *) pDataGradKernel ;
#pragma omp parallel for
    for(int n=0; n<N ; n++) {
      for(int c=0; c<C ; c++) {
        for(int hw=0; hw<H*W ; hw++) {
          filter[((hw)*C+c)*N+n] += transformed_filter[((n*C+c)*H)*W+hw] ;
        }
      }
    }

    free(transformed_filter) ;
  }

  LFTRACE_END("convolution_backward_filter_gemm");
  return VEDNN_SUCCESS;
}
// vim: et ts=2 sw=2 cindent cino=^0,=0,l0,\:0,N-s syntax=cpp.doxygen

max_stars_repo_name: sx-aurora-dev/vednn

max_issues_repo_name: sx-aurora-dev/vednn

max_forks_repo_name: sx-aurora-dev/vednn


text:
#include "indexes/utils/Utils.h"

#include <doctest/doctest.h>
#include <gsl/span>

#include <inttypes.h>
#include <limits>
#include <map>
#include <random>
#include <string>
#include <thread>
#include <vector>

enum class ConcurrentMapTestWorkload {
  WL_CONTENTED,
  WL_CONTENTED_SWAP,
  WL_RANDOM,
};

enum class LookupType {
  LT_DEFAULT,
  LT_CUSTOM,
};

std::vector<int64_t> generateUniqueValues(int num_threads, int perthread_count,
                                          ConcurrentMapTestWorkload workload);

template <typename MapType> void MixedMapTest() {
  indexes::utils::ThreadRegistry::RegisterThread();

  MapType map;

  int num_operations = 1024 * 1024;
  int cardinality = num_operations * 0.1;

  constexpr auto INSERT_OP = 1;
  constexpr auto LOOKUP_OP = 2;
  constexpr auto DELETE_OP = 3;

  std::random_device r;
  std::seed_seq seed{r(), r(), r(), r(), r(), r(), r(), r()};
  std::mt19937 rnd(seed);

  std::uniform_int_distribution<int> key_dist{1, cardinality};
  std::uniform_int_distribution<int> val_dist;
  std::uniform_int_distribution<> op_dist{INSERT_OP, DELETE_OP};

  std::map<int, int> key_values;

  for (int i = 0; i < num_operations; i++) {
    const auto key = key_dist(rnd);
    const auto op = op_dist(rnd);
    int val;

    switch (op) {
    case INSERT_OP: {
      val = val_dist(rnd);

      if (key_values.count(key)) {
        REQUIRE(*map.Search(key) == key_values[key]);
        REQUIRE(*map.Update(key, val) == key_values[key]);
      } else {
        REQUIRE(map.Insert(key, val) == true);
        REQUIRE(*map.Search(key) == val);
      }

      key_values[key] = val;
      break;
    }

    case LOOKUP_OP: {
      if (key_values.count(key))
        REQUIRE(*map.Search(key) == key_values[key]);
      else
        REQUIRE(map.Search(key).has_value() == false);

      break;
    }

    case DELETE_OP: {
      if (key_values.count(key)) {
        REQUIRE(*map.Search(key) == key_values[key]);
        REQUIRE(*map.Delete(key) == key_values[key]);
        key_values.erase(key);
      } else {
        REQUIRE(map.Delete(key).has_value() == false);
      }
      break;
    }

    default:
      continue;
    }
  }

  REQUIRE(map.size() == key_values.size());

  for (const auto &kv : key_values) {
    REQUIRE(*map.Delete(kv.first) == kv.second);
  }

  REQUIRE(map.size() == 0);

  indexes::utils::ThreadRegistry::UnregisterThread();
}

template <typename MapType, LookupType LkType, typename LookupOp>
static void lookup_worker(MapType &map, gsl::span<const int64_t> vals,
                          int64_t min_val, int64_t max_val,
                          const LookupOp &op) {
  indexes::utils::ThreadRegistry::RegisterThread();

  if constexpr (LkType == LookupType::LT_DEFAULT) {
    (void)op;
    for (auto val : vals) {
      map.Search(val);
    }
  } else {
    op(map, min_val, max_val, vals.size());
  }
  indexes::utils::ThreadRegistry::UnregisterThread();
}

template <typename MapType>
static void insert_worker(MapType &map, gsl::span<const int64_t> vals) {
  indexes::utils::ThreadRegistry::RegisterThread();

  for (auto val : vals) {
    REQUIRE(map.Insert(val, val) == true);
    REQUIRE(*map.Search(val) == val);
  }

  indexes::utils::ThreadRegistry::UnregisterThread();
}

enum class OpType { INSERT, DELETE, DELETE_AND_INSERT };

template <typename MapType>
static void delete_worker(MapType &map, gsl::span<const int64_t> vals,
                          OpType op) {
  indexes::utils::ThreadRegistry::RegisterThread();

  for (auto val : vals) {
    switch (op) {
    case OpType::DELETE:
      REQUIRE(*map.Delete(val) == val);
      REQUIRE(map.Search(val).has_value() == false);
      break;

    case OpType::DELETE_AND_INSERT:
      REQUIRE(*map.Delete(val) == val);
      REQUIRE(map.Insert(val, val) == true);
      break;

    case OpType::INSERT:
      break;
    }
  }

  indexes::utils::ThreadRegistry::UnregisterThread();
}

template <typename MapType, LookupType LkType, typename LookupOp>
void ConcurrentMapTest(ConcurrentMapTestWorkload workload, LookupOp &&lop) {
  MapType map;
  constexpr int PER_THREAD_OP_COUNT = 256 * 1024;
  const int NUM_THREADS = std::thread::hardware_concurrency();
  const std::vector<int64_t> vals =
      generateUniqueValues(NUM_THREADS, PER_THREAD_OP_COUNT, workload);
  auto [min_it, max_it] = std::minmax_element(vals.begin(), vals.begin());
  int64_t min_val = *min_it, max_val = *max_it;

  indexes::utils::ThreadRegistry::RegisterThread();

  auto run_test = [&](OpType op) {
    std::vector<std::thread> workers;
    int startval = 0;
    int quantum = vals.size() / NUM_THREADS;

    for (int i = 0; i < NUM_THREADS; i++) {
      if (op == OpType::INSERT) {
        workers.emplace_back(
            insert_worker<MapType>, std::ref(map),
            gsl::span<const int64_t>{vals.data() + startval, quantum});
      } else {
        workers.emplace_back(
            delete_worker<MapType>, std::ref(map),
            gsl::span<const int64_t>{vals.data() + startval, quantum}, op);
      }

      startval += quantum;
    }

    workers.emplace_back(lookup_worker<MapType, LkType, LookupOp>,
                         std::ref(map), gsl::span<const int64_t>{vals}, min_val,
                         max_val, std::cref(lop));

    for (auto &worker : workers) {
      worker.join();
    }
  };

  // Insert check
  {
    run_test(OpType::INSERT);

    REQUIRE(map.size() == vals.size());

    for (auto val : vals) {
      REQUIRE(*map.Search(val) == val);
    }
  }

  // Delete And Insert check
  {
    run_test(OpType::DELETE_AND_INSERT);

    for (auto val : vals) {
      REQUIRE(*map.Search(val) == val);
    }

    REQUIRE(map.size() == vals.size());
  }

  // Delete check
  {
    run_test(OpType::DELETE);

    REQUIRE(map.size() == 0);
  }

  indexes::utils::ThreadRegistry::UnregisterThread();
}


text:
#ifndef QDM_KNOTS_H
#define QDM_KNOTS_H 1

#include <gsl/gsl_rng.h>
#include <gsl/gsl_vector.h>

gsl_vector *
qdm_knots_vector(
  size_t spline_df,
  gsl_vector *knots_inter
);

int
qdm_knots_optimize(
  gsl_vector *result,

  gsl_rng *rng,

  gsl_vector *sorted_data,
  gsl_vector *middle,
  gsl_vector *possible_knots,

  size_t iterate_n,
  size_t spline_df
);

#endif /* QDM_KNOTS_H */

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
/*
Multisteps: a program to get optime Runge-Kutta and multi-steps methods.

Copyright 2011-2019, Javier Burguete Tolosa.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

	1. Redistributions of source code must retain the above copyright notice,
		this list of conditions and the following disclaimer.

	2. Redistributions in binary form must reproduce the above copyright notice,
		this list of conditions and the following disclaimer in the
		documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``AS IS'' AND ANY EXPRESS OR
IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL Javier Burguete Tolosa OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/**
 * \file rk_4_4.c
 * \brief Source file to optimize Runge-Kutta 5 steps 4th order methods.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2011-2019.
 */
#define _GNU_SOURCE
#include <string.h>
#include <math.h>
#include <libxml/parser.h>
#include <glib.h>
#include <libintl.h>
#include <gsl/gsl_rng.h>
#include "config.h"
#include "utils.h"
#include "optimize.h"
#include "rk.h"
#include "rk_4_4.h"

#define DEBUG_RK_4_4 0          ///< macro to debug.

/**
 * Function to obtain the coefficients of a 4 steps 4th order Runge-Kutta
 * method.
 */
int
rk_tb_4_4 (Optimize * optimize) ///< Optimize struct.
{
  long double *tb, *r;
#if DEBUG_RK_4_4
  fprintf (stderr, "rk_tb_4_4: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t4 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  t3 (tb) = 1.L;
  b43 (tb) = (0.25L - 1.L / 3.L * t1 (tb)
              - (1.L / 3.L - 0.5L * t1 (tb)) * t2 (tb))
    / (t3 (tb) * (t3 (tb) - t2 (tb)) * (t3 (tb) - t1 (tb)));
  if (isnan (b43 (tb)))
    return 0;
  b42 (tb) = (1.L / 3.L - 0.5L * t1 (tb)
              - b43 (tb) * t3 (tb) * (t3 (tb) - t1 (tb)))
    / (t2 (tb) * (t2 (tb) - t1 (tb)));
  if (isnan (b42 (tb)))
    return 0;
  b41 (tb) = (0.5L - b42 (tb) * t2 (tb) - b43 (tb) * t3 (tb)) / t1 (tb);
  if (isnan (b41 (tb)))
    return 0;
  b32 (tb) = (1.L / 12.L - 1.L / 6.L * t1 (tb))
    / (b43 (tb) * t2 (tb) * (t2 (tb) - t1 (tb)));
  if (isnan (b32 (tb)))
    return 0;
  b31 (tb) = ((0.125L - 1.L / 6.L * t2 (tb)) / (b43 (tb) * (t3 (tb) - t2 (tb)))
              - b32 (tb) * t2 (tb)) / t1 (tb);
  if (isnan (b31 (tb)))
    return 0;
  b21 (tb) = 1.L / 24.L / (t1 (tb) * b43 (tb) * b32 (tb));
  if (isnan (b21 (tb)))
    return 0;
  rk_b_4 (tb);
#if DEBUG_RK_4_4
  fprintf (stderr, "rk_tb_4_4: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 4 steps 4th order, 5th order in
 * equations depending only in time, Runge-Kutta method.
 */
int
rk_tb_4_4t (Optimize * optimize)        ///< Optimize struct.
{
  long double *tb, *r;
#if DEBUG_RK_4_4
  fprintf (stderr, "rk_tb_4_4t: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t4 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = 0.5L * (t1 (tb) - 0.6L) / (t1 (tb) - 0.5L);
  t3 (tb) = 1.L;
  b43 (tb) = (0.25L - 1.L / 3.L * t1 (tb)
              - (1.L / 3.L - 0.5L * t1 (tb)) * t2 (tb))
    / (t3 (tb) * (t3 (tb) - t2 (tb)) * (t3 (tb) - t1 (tb)));
  if (isnan (b43 (tb)))
    return 0;
  b42 (tb) = (1.L / 3.L - 0.5L * t1 (tb)
              - b43 (tb) * t3 (tb) * (t3 (tb) - t1 (tb)))
    / (t2 (tb) * (t2 (tb) - t1 (tb)));
  if (isnan (b42 (tb)))
    return 0;
  b41 (tb) = (0.5L - b42 (tb) * t2 (tb) - b43 (tb) * t3 (tb)) / t1 (tb);
  if (isnan (b41 (tb)))
    return 0;
  b32 (tb) = (1.L / 12.L - 1.L / 6.L * t1 (tb))
    / (b43 (tb) * t2 (tb) * (t2 (tb) - t1 (tb)));
  if (isnan (b32 (tb)))
    return 0;
  b31 (tb) = ((0.125L - 1.L / 6.L * t2 (tb)) / (b43 (tb) * (t3 (tb) - t2 (tb)))
              - b32 (tb) * t2 (tb)) / t1 (tb);
  if (isnan (b31 (tb)))
    return 0;
  b21 (tb) = 1.L / 24.L / (t1 (tb) * b43 (tb) * b32 (tb));
  if (isnan (b21 (tb)))
    return 0;
  rk_b_4 (tb);
#if DEBUG_RK_4_4
  fprintf (stderr, "rk_tb_4_4t: end\n");
#endif
  return 1;
}

/**
 * Function to calculate the objective function of a 4 steps 4th order
 * Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_4_4 (RK * rk)   ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_4_4
  fprintf (stderr, "rk_objective_tb_4_4: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b21 (tb) < 0.L)
    o += b21 (tb);
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b31 (tb) < 0.L)
    o += b31 (tb);
  if (b32 (tb) < 0.L)
    o += b32 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b41 (tb) < 0.L)
    o += b41 (tb);
  if (b42 (tb) < 0.L)
    o += b42 (tb);
  if (b43 (tb) < 0.L)
    o += b43 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L + fmaxl (1.L, fmaxl (t1 (tb), fmaxl (t2 (tb), t3 (tb))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_4_4
  fprintf (stderr, "rk_objective_tb_4_4: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_4_4: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 4 steps 4th order, 5th
 * order in equations depending only in time, Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_4_4t (RK * rk)  ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_4_4
  fprintf (stderr, "rk_objective_tb_4_4: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b21 (tb) < 0.L)
    o += b21 (tb);
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b31 (tb) < 0.L)
    o += b31 (tb);
  if (b32 (tb) < 0.L)
    o += b32 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b41 (tb) < 0.L)
    o += b41 (tb);
  if (b42 (tb) < 0.L)
    o += b42 (tb);
  if (b43 (tb) < 0.L)
    o += b43 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L + fmaxl (1.L, fmaxl (t1 (tb), fmaxl (t2 (tb), t3 (tb))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_4_4
  fprintf (stderr, "rk_objective_tb_4_4: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_4_4: end\n");
#endif
  return o;
}

max_stars_repo_name: jburguete/ode

max_issues_repo_name: jburguete/ode

max_forks_repo_name: jburguete/ode


text:
#include <string.h>
#include <mpi.h>
#include <math.h>
#include <alloca.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>

#include <fastpm/libfastpm.h>
#include <fastpm/logging.h>
#include <fastpm/string.h>

#include "pmpfft.h"

void
fastpm_powerspectrum_init(FastPMPowerSpectrum * ps, const size_t size)
{
    fastpm_funck_init(&ps->base, size);
    ps->pm = NULL;
    ps->edges = malloc(sizeof(ps->edges[0]) * (ps->base.size + 1));
    ps->Nmodes = malloc(sizeof(ps->Nmodes[0]) * ps->base.size);
}

void
fastpm_powerspectrum_init_from(FastPMPowerSpectrum * ps, const FastPMPowerSpectrum * other)
{
    fastpm_powerspectrum_init(ps, other->base.size);
    memcpy(ps->base.k, other->base.k, sizeof(ps->base.k[0]) * ps->base.size);
    memcpy(ps->base.f, other->base.f, sizeof(ps->base.f[0]) * ps->base.size);
    memcpy(ps->edges, other->edges, sizeof(ps->edges[0]) * (ps->base.size + 1));
    memcpy(ps->Nmodes, other->Nmodes, sizeof(ps->Nmodes[0]) * ps->base.size);
}

void
fastpm_powerspectrum_init_from_delta(FastPMPowerSpectrum * ps, PM * pm, const FastPMFloat * delta1_k, const FastPMFloat * delta2_k)
{
    /* This function measures powerspectrum from two overdensity or 1+overdensity fields */
    /* normalize them with fastpm_apply_normalize_transfer if needed before using this function.*/
    /* N is used to store metadata -- the shot-noise level. */
    fastpm_powerspectrum_init(ps, pm_nmesh(pm)[0] / 2);

    ps->pm = pm;
    double Volume = 1.0;
    int d;
    double k0 = 2 * M_PI / pm_boxsize(ps->pm)[0];
    for (d = 0; d < 3; d ++) {
        Volume *= pm_boxsize(pm)[d];
    }
    ps->Volume = Volume;
    ps->k0 = k0;

    memset(ps->base.f, 0, sizeof(ps->base.f[0]) * ps->base.size);
    memset(ps->base.k, 0, sizeof(ps->base.k[0]) * ps->base.size);
    memset(ps->edges, 0, sizeof(ps->edges[0]) * (ps->base.size + 1));
    memset(ps->Nmodes, 0, sizeof(ps->Nmodes[0]) * ps->base.size);

    int i;
    for(i = 0; i < ps->base.size + 1; i ++) {
        ps->edges[i] = i * k0;
    }

#pragma omp parallel
    {
        PMKIter kiter;
        for(pm_kiter_init(ps->pm, &kiter);
            !pm_kiter_stop(&kiter);
            pm_kiter_next(&kiter)) {
            int d;
            ptrdiff_t kk = 0.;
            for(d = 0; d < 3; d++) {
                double ik = kiter.iabs[d];
                if(ik > pm->Nmesh[d] / 2) ik -= pm->Nmesh[d];
                kk += ik * ik;
            }

            ptrdiff_t ind = kiter.ind;

            ptrdiff_t bin = ((ptrdiff_t)floor(sqrt(kk))) - 2;
            if(bin < 0) bin = 0;
            while((bin + 1) * (bin + 1) <= kk) {
                bin ++;
            }

            double k = sqrt(kk) * k0;

            if(bin >= 0 && bin < ps->base.size) {
                double real1 = delta1_k[ind + 0];
                double imag1 = delta1_k[ind + 1];
                double real2 = delta2_k[ind + 0];
                double imag2 = delta2_k[ind + 1];
                double value = real1 * real2 + imag1 * imag2;
                int w = 2;

                if(kiter.i[2] == 0 || kiter.i[2] == pm->Nmesh[2] / 2) w = 1;

                if(kiter.iabs[0] == 0 &&
                   kiter.iabs[1] == 0 &&
                   kiter.iabs[2] == 0) {

                } else {
                    #pragma omp atomic
                    ps->Nmodes[bin] += w;
                    #pragma omp atomic
                    ps->base.f[bin] += w * value; /// cic;
                    #pragma omp atomic
                    ps->base.k[bin] += w * k;
                }
            }
        }
    }

    MPI_Allreduce(MPI_IN_PLACE, ps->base.f, ps->base.size, MPI_DOUBLE, MPI_SUM, ps->pm->Comm2D);
    MPI_Allreduce(MPI_IN_PLACE, ps->Nmodes, ps->base.size, MPI_DOUBLE, MPI_SUM, ps->pm->Comm2D);
    MPI_Allreduce(MPI_IN_PLACE, ps->base.k, ps->base.size, MPI_DOUBLE, MPI_SUM, ps->pm->Comm2D);

    ptrdiff_t ind;
    for(ind = 0; ind < ps->base.size; ind++) {
        if(ps->Nmodes[ind] == 0) continue;
        ps->base.k[ind] /= ps->Nmodes[ind];
        ps->base.f[ind] /= ps->Nmodes[ind];
        ps->base.f[ind] *= ps->Volume;
    }
}

void
fastpm_transferfunction_init(FastPMPowerSpectrum * ps, PM * pm, FastPMFloat * src_k, FastPMFloat * dest_k)
{
    FastPMPowerSpectrum * ps2 = alloca(sizeof(*ps2));

    fastpm_powerspectrum_init_from_delta(ps, pm, src_k, src_k);
    fastpm_powerspectrum_init_from_delta(ps2, pm, dest_k, dest_k);

    ptrdiff_t i;
    for(i = 0; i < ps->base.size; i ++) {
        ps->base.f[i] = sqrt(ps2->base.f[i] / ps->base.f[i]);
    }

    fastpm_powerspectrum_destroy(ps2);
}

void
fastpm_powerspectrum_destroy(FastPMPowerSpectrum * ps) {
    free(ps->edges);
    free(ps->Nmodes);
    fastpm_funck_destroy(&ps->base);
}

void
fastpm_powerspectrum_write(FastPMPowerSpectrum * ps, char * filename, double N)
{
    FILE * fp = fopen(filename, "w");
    int i;
    fprintf(fp, "# k p N \n");
    for(i = 0; i < ps->base.size; i ++) {
        fprintf(fp, "%g %g %g\n", ps->base.k[i], ps->base.f[i], ps->Nmodes[i]);
    }
    double * BoxSize = pm_boxsize(ps->pm);
    fprintf(fp, "# metadata 7\n");
    fprintf(fp, "# volume %g float64\n", ps->Volume);
    fprintf(fp, "# shotnoise %g float64\n", ps->Volume / N);
    fprintf(fp, "# N1 %g int\n", N);
    fprintf(fp, "# N2 %g int\n", N);
    fprintf(fp, "# Lz %g float64\n", BoxSize[2]);
    fprintf(fp, "# Lx %g float64\n", BoxSize[0]);
    fprintf(fp, "# Ly %g float64\n", BoxSize[1]);
    fclose(fp);
}

double
fastpm_powerspectrum_large_scale(FastPMPowerSpectrum * ps, int Nmax)
{
    double kmax = Nmax * ps->k0;
    ptrdiff_t i;
    double Plin = 0;
    double Nmodes = 0;
    /* ignore zero mode ! */
    for(i = 0; (i == 0) || (i < ps->base.size && ps->base.k[i] <= kmax); i ++) {
        Plin += ps->base.f[i] * ps->Nmodes[i];
        Nmodes += ps->Nmodes[i];
    }
    Plin /= Nmodes;
    return Plin;
}

/* inverted calling signature for callbacks */
double
fastpm_powerspectrum_eval2(double k, FastPMPowerSpectrum * ps)
{
    return fastpm_powerspectrum_eval(ps, k);
}

double
fastpm_powerspectrum_eval(FastPMPowerSpectrum * ps, double k)
{
    return fastpm_funck_eval(&ps->base, k);
}

double
fastpm_powerspectrum_get(FastPMPowerSpectrum * ps, double k)
{
    if(k == 0) return 1;

    /* ignore the 0 mode */

    int l = 0;
    int r = ps->base.size;

    while(r - l > 1) {
        int m = (r + l) / 2;
        /* if we are on the exact k, return value from there.*/
        if(k <= ps->edges[m])
            r = m;
        else
            l = m;
    }

    return ps->base.f[l];
}

double
fastpm_powerspectrum_get2(double k, FastPMPowerSpectrum * ps)
{
    return fastpm_powerspectrum_get(ps, k);
}

struct sigma2_int {
    FastPMPowerSpectrum * ps;
    double R;
};
static
double sigma2_int(double k, struct sigma2_int * param)
{
    double kr, kr3, kr2, w, x;
    double r_tophat = param->R;

    kr = r_tophat * k;
    kr2 = kr * kr;
    kr3 = kr2 * kr;

    if(kr < 1e-8)
        return 0;

    w = 3 * (sin(kr) / kr3 - cos(kr) / kr2);
    x = 4 * M_PI * k * k * w * w * fastpm_powerspectrum_eval(param->ps, k);

    return x / pow(2 * M_PI, 3);
}

double
fastpm_powerspectrum_sigma(FastPMPowerSpectrum * ps, double R)
{
    const int WORKSIZE = 81920;

    double result, abserr;
    gsl_integration_workspace *workspace;
    gsl_function F;

    struct sigma2_int param;
    param.ps = ps;
    param.R = R;

    workspace = gsl_integration_workspace_alloc(WORKSIZE);

    F.function = (void*) &sigma2_int;
    F.params = &param;

    void * handler = gsl_set_error_handler_off ();
    //
    // note: 500/R is here chosen as (effectively) infinity integration boundary
    gsl_integration_qag(&F, 0, 500.0 * 1 / R,
            0, 1.0e-4, WORKSIZE, GSL_INTEG_GAUSS41, workspace, &result, &abserr);

    // high precision caused error
    gsl_integration_workspace_free(workspace);
    gsl_set_error_handler (handler);

    return sqrt(result);
}

void
fastpm_powerspectrum_scale(FastPMPowerSpectrum * ps, double factor)
{
    ptrdiff_t i;
    /* neglect the zero mode */
    for(i = 1; i < ps->base.size; i ++) {
        ps->base.f[i] *= factor;
    }
}

void
fastpm_powerspectrum_rebin(FastPMPowerSpectrum * ps, size_t newsize)
{
    /* this doesn't work */
    double * k1 = malloc(newsize * sizeof(k1[0]));
    double * p1 = malloc(newsize * sizeof(p1[0]));
    double * Nmodes1 = malloc(newsize * sizeof(Nmodes1[0]));
    double * edges1 = malloc((newsize + 1) * sizeof(edges1[0]));

    ptrdiff_t i;

    for(i = 0; i < newsize; i ++) {
        ptrdiff_t j1 = (i    ) * ps->base.size / newsize;
        ptrdiff_t j2 = (i + 1) * ps->base.size / newsize;
        ptrdiff_t j;
        k1[i] = 0;
        p1[i] = 0;
        Nmodes1[i] = 0;
        edges1[i] = ps->edges[j1];
        edges1[i + 1] = ps->edges[j2];
        for(j = j1; j < j2; j ++) {
            k1[i] += ps->base.k[j] * ps->Nmodes[j];
            p1[i] += ps->base.f[j] * ps->Nmodes[j];
            Nmodes1[i] += ps->Nmodes[j];
        }
        if(Nmodes1[i] > 0) {
            k1[i] /= Nmodes1[i];
            p1[i] /= Nmodes1[i];
        }
    }
    free(ps->base.k);
    free(ps->base.f);
    free(ps->Nmodes);
    free(ps->edges);
    ps->base.k = k1;
    ps->base.f = p1;
    ps->Nmodes = Nmodes1;
    ps->edges = edges1;
    ps->base.size = newsize;
}

int
fastpm_powerspectrum_init_from_string(FastPMPowerSpectrum * ps, const char * string)
{
    int r = fastpm_funck_init_from_string(&ps->base, string);
    ps->edges = malloc(sizeof(ps->edges[0]) * (ps->base.size + 1));
    ps->Nmodes = malloc(sizeof(ps->Nmodes[0]) * ps->base.size);
    return r;
}

void
fastpm_funck_init(FastPMFuncK * fk, const size_t size)
{
    fk->size = size;
    fk->k = malloc(sizeof(fk->k[0]) * fk->size);
    fk->f = malloc(sizeof(fk->f[0]) * fk->size);
}

int
fastpm_funck_init_from_string(FastPMFuncK * fk, const char * string)
{
    char ** list = fastpm_strsplit(string, "\n");
    char ** line;
    ptrdiff_t i;
    int pass = 0;
    /* two pass parsing, first pass for counting */
    /* second pass for assignment */
    while(pass < 2) {
        i = 0;
        for (line = list; *line; line++) {
            double k, f;
            if(2 == sscanf(*line, "%lg\t%lg", &k, &f)) {   // note tab. Could make more general?
                if(pass == 1) {
                    fk->k[i] = k;
                    fk->f[i] = f;
                }
                i ++;
            }
        }
        if(pass == 0) {
            fastpm_funck_init(fk, i);
        }
        pass ++;
    }

    free(list);
    if(fk->size == 0) {
        /* Nothing is savagable in the file.*/
        return 0;
    }
    return 0;
}

/* inverted calling signature for callbacks */
double
fastpm_funck_eval2(double k, FastPMFuncK * fk)
{
    return fastpm_funck_eval(fk, k);
}

double
fastpm_funck_eval(FastPMFuncK * fk, double k)
{

    /* ignore the 0 mode */
    if(k == 0) return 1;

    int l = 0;
    int r = fk->size - 1;

    /* determine the right and left values around k */
    while(r - l > 1) {
        int m = (r + l) / 2;
        if(k < fk->k[m])
            r = m;
        else
            l = m;
    }
    double k2 = fk->k[r],
           k1 = fk->k[l];
    double f2 = fk->f[r],
           f1 = fk->f[l];

    if(l == r) {
        return fk->f[l];
    }

    if(f1 <= 0 || f2 <= 0 || k1 == 0 || k2 == 0) {
        /* if any of the f is zero, or negative, use linear interpolation */
        double f = (k - k1) * f2 + (k2 - k) * f1;
        f /= (k2 - k1);
        return f;
    } else {
        k = log(k);
        f1 = log(f1);
        f2 = log(f2);
        k1 = log(k1);
        k2 = log(k2);
        double f = (k - k1) * f2 + (k2 - k) * f1;
        f /= (k2 - k1);
        return exp(f);
    }
}

void
fastpm_funck_destroy(FastPMFuncK * fk)
{
    free(fk->f);
    free(fk->k);
}

max_stars_repo_name: sbird/FastPMRunner

max_issues_repo_name: sbird/FastPMRunner

max_forks_repo_name: sbird/FastPMRunner


text:
#ifndef INCLUDED_TELEPATH_BLAS_3_GEMM_H
#define INCLUDED_TELEPATH_BLAS_3_GEMM_H

#include <cblas.h>

namespace telepath{

    template< typename T >
    void gemm( CBLAS_ORDER order, CBLAS_TRANSPOSE transA,
        CBLAS_TRANSPOSE transB, std::size_t m std::size_t n, std::size_t k,
        T alpha, const T* A, std::size_t lda, const T* B, std::size_t ldb,
        T beta, T* C, std::size_t ldc );
    template<>
    void gemm<float>( CBLAS_ORDER order, CBLAS_TRANSPOSE transA,
        CBLAS_TRANSPOSE transB, std::size_t m, std::size_t n, std::size_t k,
        float alpha, const float* A, std::size_t lda, const float* B,
        std::size_t ldb, float beta, float* C, std::size_t ldc ){
        return cblas_sgemm( order, transA, transB, m, n, k, alpha, A, lda, B,
            ldb, beta, C, ldc );
    }
    template<>
    void gemm<double>( CBLAS_ORDER order, CBLAS_TRANSPOSE transA,
        CBLAS_TRANSPOSE transB, std::size_t m, std::size_t n, std::size_t k,
        double alpha, const double* A, std::size_t lda, const double* B,
        std::size_t ldb, double beta, double* C, std::size_t ldc ){
        return cblas_dgemm( order, transA, transB, m, n, k, alpha, A, lda, B,
            ldb, beta, C, ldc );
    }
    template<>
    void gemm<std::complex<float>>( CBLAS_ORDER order, CBLAS_TRANSPOSE transA,
        CBLAS_TRANSPOSE transB, std::size_t m, std::size_t n, std::size_t k,
        const std::complex<float>* alpha, const std::complex<float>* A,
        std::size_t lda, std::const complex<float>* B, std::size_t ldb,
        const std::complex<float>* beta, std::complex<float>* C,
        std::size_t ldc ){
        return cblas_cgemm( order, transA, transB, m, n, k, (float*) alpha,
            (float*) A, lda, (float*) B, ldb, (float*) beta, (float*) C, ldc );
    }
    template<>
    void gemm<std::complex<double>>( CBLAS_ORDER order, CBLAS_TRANSPOSE transA,
        CBLAS_TRANSPOSE transB, std::size_t m, std::size_t n, std::size_t k,
        const std::complex<double>* alpha, const std::complex<double>* A,
        std::size_t lda, const std::complex<double>* B, std::size_t ldb,
        const std::complex<double>* beta, std::complex<double>* C,
        std::size_t ldc ){
        return cblas_zgemm( order, transA, transB, m, n, k, (double*) alpha,
            (double*) A, lda, (double*) B, ldb, (double*) beta, (double*) C,
            ldc );
    }

} //namespace telepath

#endif

max_stars_repo_name: tbepler/telepath

max_issues_repo_name: tbepler/telepath

max_forks_repo_name: tbepler/telepath


text:
/**
 *@file CosanBO.h
 *@author Xinyu Zhang
 *@brief CosanBO
 *@date 2021-04-10
 */

#ifndef __COSANBO_H_INCLUDED__
#define __COSANBO_H_INCLUDED__
#include <string>
#ifndef FMT_HEADER_ONLY
#define FMT_HEADER_ONLY
#include <fmt/format.h>
#endif
#include <gsl/gsl>
#include <type_traits>
#include <cosan/utils/utils.h>
#include <cosan/utils/Exceptions.h>
#include <random>

template<typename NumericType>
concept Numeric = std::is_arithmetic<NumericType>::value ;

template <class T, class U>
concept Derived = std::is_base_of<U, T>::value;

namespace Cosan
{

//    template<typename NumericType,
//            typename = typename std::enable_if<std::is_arithmetic<NumericType>::value,NumericType>::type
//    >

    /** Description: data structure for matrix. Will be used throughout the project. */
    template<Numeric NumericType>
    using CosanMatrix = Eigen::Matrix<NumericType, Eigen::Dynamic, Eigen::Dynamic> ;

//    template<typename NumericType,
//            typename = typename std::enable_if<std::is_arithmetic<NumericType>::value,NumericType>::type
//    >
    /** Description: data structure for column vector. Will be used throughout the project. */
    template<Numeric NumericType>
    using CosanColVector =  Eigen::Matrix<NumericType,  Eigen::Dynamic, 1> ;

//    template<typename NumericType,
//            typename = typename std::enable_if<std::is_arithmetic<NumericType>::value,NumericType>::type
//    >

/** Description: data structure for row vector. Will be used throughout the project. */
    template<Numeric NumericType>
    using CosanRowVector =  Eigen::Matrix<NumericType,  1, Eigen::Dynamic> ;

//    typedef std::variant<CosanMatrix<unsigned int>, CosanMatrix<unsigned long>,
//                         CosanMatrix<unsigned long long>, >

    /**
     * @brief Cosan Base Object
     * @details: Description: Base object CosanBO which all other Cosan objects inherited from.
     * */
    class CosanBO
    {
    public:
        /**
         *  @brief Default constructor.
         **/
        CosanBO(){
        }
//        virtual ~CosanBO();
//        virtual CosanBO *Shallow_copy() const;
//        virtual CosanBO *Deep_copy() const;

        /**
         *  @brief Get the name of the objects.
         *  @details Description: Get the name of the object. It should return "Abstract Object"
         **/
        virtual const std::string  GetName() const {
            return "Abstract Object";}

//        virtual bool SaveFile(const std::string &path ,const std::string & prefix = "");
//        virtual bool LoadFile(const std::string &path);
//        void PrintModel();
//        virtual bool Equals(CosanBO* other, float accuracy = 0.0);
//        virtual CosanBO* Clone();

    protected:
    };

}

#endif

max_stars_repo_name: zhxinyu/cosan

max_issues_repo_name: zhxinyu/cosan

max_forks_repo_name: zhxinyu/cosan


text:
#pragma once
#include <wrl/client.h>
#include <d3d11.h>
#include <gsl/span>

namespace vrutil
{
template <typename T> using ComPtr = Microsoft::WRL::ComPtr<T>;
using SharedTexture = std::pair<ComPtr<ID3D11Texture2D>, HANDLE>;

class D3DManager
{
    class D3D11ManagerImpl *m_impl;

public:
    D3DManager();
    ~D3DManager();
    bool Initialize();
    ComPtr<ID3D11Texture2D> CreateStaticTexture(int w, int h, const gsl::span<uint8_t> &bytes);
    SharedTexture CreateSharedStaticTexture(int w, int h, const gsl::span<uint8_t> &bytes);
    void Copy(const ComPtr<ID3D11Texture2D> &src, const ComPtr<ID3D11Texture2D> &dst);
};

} // namespace vrutil

max_stars_repo_name: ousttrue/VROverlaySample

max_issues_repo_name: ousttrue/VROverlaySample

max_forks_repo_name: ousttrue/VROverlaySample


text:
/* Author:  G. Jungman */

#ifndef __GSL_SF_H__
#define __GSL_SF_H__

#include <gsl/gsl_sf_result.h>

#include <gsl/gsl_sf_airy.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_sf_clausen.h>
#include <gsl/gsl_sf_coupling.h>
#include <gsl/gsl_sf_coulomb.h>
#include <gsl/gsl_sf_dawson.h>
#include <gsl/gsl_sf_debye.h>
#include <gsl/gsl_sf_dilog.h>
#include <gsl/gsl_sf_elementary.h>
#include <gsl/gsl_sf_ellint.h>
#include <gsl/gsl_sf_elljac.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_expint.h>
#include <gsl/gsl_sf_fermi_dirac.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_gegenbauer.h>
#include <gsl/gsl_sf_hyperg.h>
#include <gsl/gsl_sf_laguerre.h>
#include <gsl/gsl_sf_lambert.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf_psi.h>
#include <gsl/gsl_sf_synchrotron.h>
#include <gsl/gsl_sf_transport.h>
#include <gsl/gsl_sf_trig.h>
#include <gsl/gsl_sf_zeta.h>

#endif /* __GSL_SF_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws