text:
#pragma once

#include <chrono>
#include <fmt/printf.h>
#include <gsl/gsl-lite.hpp>

namespace angonoka::cli {
/**
    Pretty-print values.

    This class is mostly used for specializing user-defined
    formatters for the fmt library. We don't want to specialize
    the formatter for the original type because we might not
    always want a human-readable output.

    Instead of specializing fmt::formatter<std::chrono::seconds>
    we specialize fmt::formatter<humanize<std::chrono::seconds>>
    thus not interfering with the original formatter.

    @var value Value to be pretty-printed
*/
template <typename T> struct humanize {
    T value;
};
template <typename T> humanize(T) -> humanize<T>;

namespace detail {
    /**
        Enum that tells if a human-readable duration
        printing should continue parsing durations of a
        smaller scale or stop.
    */
    enum DurationParsing : bool { continue_, stop };

    /**
        Helper function to print durations in human-readable form.

        @param article  "a" or "an"
        @param name     Duration (seconds, minutes, etc)

        @return Formatter function
    */
    template <typename T>
    constexpr auto
    humanize_duration(gsl::czstring article, gsl::czstring name)
    {
        return [=](auto total, auto& ctx) {
            const auto dur = std::chrono::round<T>(total);
            if (dur == dur.zero()) return DurationParsing::continue_;
            const auto ticks = dur.count();
            if (ticks == 1) {
                fmt::format_to(
                    ctx.out(),
                    "about {} {}",
                    article,
                    name);
            } else {
                fmt::format_to(ctx.out(), "{} {}s", ticks, name);
            }
            return DurationParsing::stop;
        };
    }
} // namespace detail
} // namespace angonoka::cli

namespace fmt {
using angonoka::cli::humanize;
/**
    User-defined formatter for std::chrono durations.
*/
template <typename... Ts>
struct fmt::formatter<humanize<std::chrono::duration<Ts...>>> {
    using value_type = humanize<std::chrono::duration<Ts...>>;
    static constexpr auto parse(format_parse_context& ctx)
    {
        return ctx.end();
    }

    template <typename FormatContext>
    constexpr auto format(const value_type& obj, FormatContext& ctx)
    {
        using angonoka::cli::detail::humanize_duration;
        using namespace std::literals::chrono_literals;
        using namespace std::chrono;
        constexpr auto min_threshold = 5s;
        if (obj.value <= min_threshold)
            return format_to(ctx.out(), "a few seconds");
        [&](auto&&... fns) {
            (fns(obj.value, ctx) || ...);
        }(humanize_duration<months>("a", "month"),
          humanize_duration<days>("a", "day"),
          humanize_duration<hours>("an", "hour"),
          humanize_duration<minutes>("a", "minute"),
          humanize_duration<seconds>("a", "second"));
        return ctx.out();
    }
};
} // namespace fmt

max_stars_repo_name: coffee-lord/angonoka

max_issues_repo_name: coffee-lord/angonoka

max_forks_repo_name: coffee-lord/angonoka


text:
#pragma once

#include "RTTI.h"
#include <stack>
#include <gsl\gsl>
#include <GL/gl3w.h>
#include "Framebuffer.h"

namespace Library
{
	class RenderTarget : public RTTI
	{
		RTTI_DECLARATIONS(RenderTarget, RTTI)

	public:
		RenderTarget() = default;
		RenderTarget(const RenderTarget&) = delete;
		RenderTarget(RenderTarget&&) = default;
		RenderTarget& operator=(const RenderTarget&) = delete;
		RenderTarget& operator=(RenderTarget&&) = default;
		virtual ~RenderTarget() = default;

		void GenTexture(GLuint width, GLuint height);

		GLuint GetTextureID();

		void Begin();
		void End();
		void RebindCurrentRenderTargets();

	protected:
		std::shared_ptr<Framebuffer> mDrawTarget;

	private:
		static std::stack<std::shared_ptr<Framebuffer>> sRenderTargetStack;
	};
}

ext: h

max_stars_repo_name: DakkyDaWolf/OpenGL

max_issues_repo_name: DakkyDaWolf/OpenGL

max_forks_repo_name: DakkyDaWolf/OpenGL


text:
/*
Copyright 2016, D. E. Shaw Research.
All rights reserved.

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

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

* 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.

* Neither the name of D. E. Shaw Research nor the names of its
  contributors may be used to endorse or promote products derived from
  this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"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 THE COPYRIGHT
OWNER 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.
*/

#ifdef __GNUC__
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wfloat-equal"
#pragma GCC diagnostic ignored "-Wconversion"
#endif

#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wexpansion-to-defined"
#pragma clang diagnostic ignored "-Wreserved-id-macro"
#endif

#ifdef _MSC_FULL_VER
// 4311 - 'type cast': pointer truncation from 'char[54]' to 'long'
// 4202 - nonstandard extension used: non-constant aggregate initializer
// 4127 - conditional expression is constant
// 4302 - 'type cast': truncation from 'const char *' to 'long'
#pragma warning(push)
#pragma warning(disable : 4311 4204 4127 4302)
#endif

#include "rng/config.h"

#include "Random123/conventional/gsl_cbrng.h"
#include "ut_gsl.h"
#include <assert.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>

/* Exercise the GSL_CBRNG macro */

GSL_CBRNG(cbrng, threefry4x64); /* creates gsl_rng_cbrng */

int main(int argc, char **argv) {
  int i;
  gsl_rng *r;
  gsl_rng *rcopy;
  unsigned long x;
  rngRemember(unsigned long save);
  rngRemember(unsigned long saved[5]);
  double sum = 0.;
  (void)argc;
  (void)argv; /* unused */

  /* Silence an unused-parameter warning. */
  (void)argc;

  r = gsl_rng_alloc(gsl_rng_cbrng);
  assert(gsl_rng_min(r) == 0);
  assert(gsl_rng_max(r) == 0xffffffffUL); /* Not necessarily ~0UL */
  assert(gsl_rng_size(r) > 0);

  printf("%s\nulongs from %s in initial state\n", argv[0], gsl_rng_name(r));
  for (i = 0; i < 5; i++) {
    x = gsl_rng_get(r);
    rngRemember(saved[i] = x);
    printf("%d: 0x%lx\n", i, x);
    assert(x != 0);
  }
  printf("uniforms from %s\n", gsl_rng_name(r));
  for (i = 0; i < 5; i++) {
    double z = gsl_rng_uniform(r);
    sum += z;
    printf("%d: %.4g\n", i, z);
  }
  assert(sum < 0.9 * 5 && sum > 0.1 * 5 &&
         (long)"sum must be reasonably close  to 0.5*number of trials");
  rngRemember(save =) gsl_rng_get(r);

  gsl_rng_set(r, 0xdeadbeef); /* set a non-zero seed */
  printf("ulongs from %s after seed\n", gsl_rng_name(r));
  for (i = 0; i < 5; i++) {
    x = gsl_rng_get(r);
    printf("%d: 0x%lx\n", i, x);
    assert(x != 0);
  }
  /* make a copy of the total state */
  rcopy = gsl_rng_alloc(gsl_rng_cbrng);
  gsl_rng_memcpy(rcopy, r);
  printf("uniforms from %s\n", gsl_rng_name(r));
  sum = 0.;
  for (i = 0; i < 5; i++) {
    double x2 = gsl_rng_uniform(r);
    rngRemember(double y = gsl_rng_uniform(rcopy));
    printf("%d: %.4g\n", i, x2);
    sum += x2;
    assert(x2 == y);
  }
  assert(gsl_rng_get(r) != save);
  assert(sum < 0.9 * 5 && sum > 0.1 * 5 &&
         (long)"sum must be reasonably close  to 0.5*number of trials");

  /* gsl_rng_set(*, 0) is supposed to recover the default seed */
  gsl_rng_set(r, 0);
  printf("ulongs from %s after restore to initial\n", gsl_rng_name(r));
  for (i = 0; i < 5; i++) {
    x = gsl_rng_get(r);
    assert(x == saved[i]);
    printf("%d: 0x%lx\n", i, x);
    assert(x != 0);
  }
  printf("uniforms from %s\n", gsl_rng_name(r));
  for (i = 0; i < 5; i++) {
    printf("%d: %.4g\n", i, gsl_rng_uniform(r));
  }

  gsl_rng_free(rcopy);
  gsl_rng_free(r);

  printf("ut_gsl: OK\n");
  return 0;
}

#ifdef _MSC_FULL_VER
#pragma warning(pop)
#endif

#ifdef __clang__
// Restore clang diagnostics to previous state.
#pragma clang diagnostic pop
#endif

#ifdef __GNUC__
// Restore GCC diagnostics to previous state.
#pragma GCC diagnostic pop
#endif

ext: c

max_stars_repo_name: jhchang-lanl/Draco

max_issues_repo_name: jhchang-lanl/Draco

max_forks_repo_name: jhchang-lanl/Draco


text:
#include <linux/types.h>
#include "sweeny_bfs.h"
#include <stdio.h>
#include "timeseries.h"
#include "extract_args.h"
#include <math.h>
#include <gsl/gsl_rng.h>

#ifdef SEQUENTIAL
char impl_title[] = "Sequential BFS implementation";
#define TS_FILE_D "../time_series/sbfs/simulation_l%u_q%.4f_b%.4f_c%.4f_s%u.hdf5"
#define INFO_STRING "*Sequential BFS implementation of Sweeny's algorithm*\n"
#define INIT init_sweeny_sbfs
#define SIMULATE  simulate_sweeny_sbfs
#define DESTROY destroy_sweeny_sbfs
#else
char impl_title[] = "Interleaved BFS implementation";
#define TS_FILE_D "../time_series/ibfs/simulation_l%u_q%.4f_b%.4f_c%.4f_s%u.hdf5"
#define INFO_STRING "*Interleaved BFS implementation of Sweeny's algorithm*\n"
#define INIT init_sweeny_ibfs
#define SIMULATE  simulate_sweeny_ibfs
#define DESTROY destroy_sweeny_ibfs
#endif
#define BUF_LEN 512
char fn[BUF_LEN];

char verbose=0;
__u32 DX;
__u32 N;
__u32 seed;
double q;
double coupling;
double beta;
double v;
double K;
__u32 steps;
__s8 *bonds;
gsl_rng * r;
char *filename;
__u32 cutoff;
__u32 *num_bonds , *num_cluster, *size_giant;
__u64 *sec_cs_moment,*four_cs_moment;
int main(int argc, char *argv[])
{

    if(!extractArgs(argc, argv,impl_title))
      return EXIT_FAILURE;
    snprintf(fn,BUF_LEN,TS_FILE_D,DX,q,beta,coupling,seed);

    if(!init_observables())
        return EXIT_FAILURE;
    if(!INIT(q,DX,beta,coupling,cutoff,steps,seed,num_bonds,num_cluster,size_giant,sec_cs_moment,four_cs_moment))
        return EXIT_FAILURE;
    if(!SIMULATE())
       return EXIT_FAILURE;
    DESTROY();
    save_timeseries();
    destroy_observables();
    return EXIT_SUCCESS;
}
/******************************************************************************
 *****************************************************************************/

max_stars_repo_name: ernmeel/sweeny

max_issues_repo_name: ernmeel/sweeny


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sort.h>

#include "cosmocalc.h"

static double linear_powspec_norm;
static gsl_spline *cosmocalc_sigmaR_spline = NULL;
static gsl_interp_accel *cosmocalc_sigmaR_acc = NULL;
static gsl_spline *cosmocalc_Rsigma_spline = NULL;
static gsl_interp_accel *cosmocalc_Rsigma_acc = NULL;

static void init_cosmocalc_peakheight_table(void);

static void init_cosmocalc_peakheight_table(void)
{
  static int initFlag = 1;
  static int currCosmoNum;

  double sigmar_table[COSMOCALC_PEAKHEIGHT_TABLE_LENGTH];
  double topHatRad_table[COSMOCALC_PEAKHEIGHT_TABLE_LENGTH];
  size_t sindex[COSMOCALC_PEAKHEIGHT_TABLE_LENGTH];
  double tmpDouble[COSMOCALC_PEAKHEIGHT_TABLE_LENGTH];
  long i;
  double rad;

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      initFlag = 0;
      currCosmoNum = cosmoData.cosmoNum;

      linear_powspec_norm = cosmoData.Sigma8*cosmoData.Sigma8/tophatradnorm_linear_powspec_exact_nonorm(8.0);

      for(i=0;i<COSMOCALC_PEAKHEIGHT_TABLE_LENGTH;++i)
        {
          rad = log(PH_R_MAX/PH_R_MIN)/(COSMOCALC_PEAKHEIGHT_TABLE_LENGTH-1.0)*((double) i) + log(PH_R_MIN);
          topHatRad_table[i] = rad;
          sigmar_table[i] = log(sqrt(linear_powspec_norm*tophatradnorm_linear_powspec_exact_nonorm(exp(rad))));
        }

      //init the spline and accelerators
      if(cosmocalc_sigmaR_spline != NULL)
        gsl_spline_free(cosmocalc_sigmaR_spline);
      cosmocalc_sigmaR_spline = gsl_spline_alloc(GSL_SPLINE_TYPE,(size_t) (COSMOCALC_PEAKHEIGHT_TABLE_LENGTH));
      gsl_spline_init(cosmocalc_sigmaR_spline,topHatRad_table,sigmar_table,(size_t) (COSMOCALC_PEAKHEIGHT_TABLE_LENGTH));
      if(cosmocalc_sigmaR_acc != NULL)
        gsl_interp_accel_reset(cosmocalc_sigmaR_acc);
      else
        cosmocalc_sigmaR_acc = gsl_interp_accel_alloc();

      //reverse the table
      gsl_sort_index(sindex,sigmar_table,(size_t) 1,(size_t) (COSMOCALC_PEAKHEIGHT_TABLE_LENGTH));
      for(i=0;i<COSMOCALC_PEAKHEIGHT_TABLE_LENGTH;++i)
        tmpDouble[i] = sigmar_table[sindex[i]];
      for(i=0;i<COSMOCALC_PEAKHEIGHT_TABLE_LENGTH;++i)
        sigmar_table[i] = tmpDouble[i];
      for(i=0;i<COSMOCALC_PEAKHEIGHT_TABLE_LENGTH;++i)
        tmpDouble[i] = topHatRad_table[sindex[i]];
      for(i=0;i<COSMOCALC_PEAKHEIGHT_TABLE_LENGTH;++i)
        topHatRad_table[i] = tmpDouble[i];

      //init the spline and accelerators
      if(cosmocalc_Rsigma_spline != NULL)
        gsl_spline_free(cosmocalc_Rsigma_spline);
      cosmocalc_Rsigma_spline = gsl_spline_alloc(GSL_SPLINE_TYPE,(size_t) (COSMOCALC_PEAKHEIGHT_TABLE_LENGTH));
      gsl_spline_init(cosmocalc_Rsigma_spline,sigmar_table,topHatRad_table,(size_t) (COSMOCALC_PEAKHEIGHT_TABLE_LENGTH));
      if(cosmocalc_Rsigma_acc != NULL)
        gsl_interp_accel_reset(cosmocalc_Rsigma_acc);
      else
        cosmocalc_Rsigma_acc = gsl_interp_accel_alloc();
    }
}

double sigmaRtophat_exact(double topHatRad, double a)
{
  static int initFlag = 1;
  static int currCosmoNum;
  static double linear_powspec_norm;
  double gf;

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      initFlag = 0;
      currCosmoNum = cosmoData.cosmoNum;

      linear_powspec_norm = cosmoData.Sigma8*cosmoData.Sigma8/tophatradnorm_linear_powspec_exact_nonorm(8.0);
    }

  gf = growth_function(a);

  return sqrt(linear_powspec_norm*tophatradnorm_linear_powspec_exact_nonorm(topHatRad))*gf;
}

double sigmaRtophat(double topHatRad, double a)
{
  static int initFlag = 1;
  static int currCosmoNum;

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      initFlag = 0;
      currCosmoNum = cosmoData.cosmoNum;
      init_cosmocalc_peakheight_table();
    }

  return exp(gsl_spline_eval(cosmocalc_sigmaR_spline,log(topHatRad),cosmocalc_sigmaR_acc))*growth_function(a);
}

double sigmaMtophat(double m, double a)
{
  return sigmaRtophat(pow(m/(4.0/3.0*M_PI*RHO_CRIT*cosmoData.OmegaM),1.0/3.0),a);
}

double sigmaMtophat_exact(double m, double a)
{
  return sigmaRtophat_exact(pow(m/(4.0/3.0*M_PI*RHO_CRIT*cosmoData.OmegaM),1.0/3.0),a);
}

double inverse_sigmaRtophat(double sigmaR, double a)
{
  static int initFlag = 1;
  static int currCosmoNum;

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      initFlag = 0;
      currCosmoNum = cosmoData.cosmoNum;
      init_cosmocalc_peakheight_table();
    }

  return exp(gsl_spline_eval(cosmocalc_Rsigma_spline,log(sigmaR/growth_function(a)),cosmocalc_Rsigma_acc));
}

double inverse_sigmaMtophat(double sigmaR, double a)
{
  return (4.0/3.0*M_PI*RHO_CRIT*cosmoData.OmegaM)*pow(inverse_sigmaRtophat(sigmaR,a),3.0);
}

double inverse_nuMtophat(double nu, double a)
{
  return (4.0/3.0*M_PI*RHO_CRIT*cosmoData.OmegaM)*pow(inverse_sigmaRtophat(DELTAC/nu,a),3.0);
}

double inverse_nuRtophat(double nu, double a)
{
  return inverse_sigmaRtophat(DELTAC/nu,a);
}

max_stars_repo_name: beckermr/cosmocalc

max_issues_repo_name: beckermr/cosmocalc

max_forks_repo_name: beckermr/cosmocalc


text:
#ifdef _MSC_VER
#if _MSC_VER >= 1910
#include <CppCoreCheck/warnings.h>
#pragma warning(disable: ALL_CPPCORECHECK_WARNINGS)
#endif // _MSC_VER >= 1910
#endif // _MSC_VER

#include <gsl/gsl>

#ifdef _MSC_VER
#if _MSC_VER >= 1910
#pragma warning(default: ALL_CPPCORECHECK_WARNINGS)
#endif // _MSC_VER >= 1910
#endif // _MSC_VER

max_stars_repo_name: TomoyukiAota/CoreCheckerAndGsl

max_issues_repo_name: TomoyukiAota/CoreCheckerAndGsl

max_forks_repo_name: TomoyukiAota/CoreCheckerAndGsl


text:
#include "odeint.h"

#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <assert.h>
#include <float.h>
#include <stdio.h>
#if HAVE_MKL_CBLAS
#include <mkl_cblas.h>
#elif HAVE_OPENBLAS
#include <openblas/cblas.h>
#else
#include <cblas.h>
#endif

/* TODO:
 * - implicit methods
 * - off-GPU
 */

/**
 * The maximum number of steps.
 */
#define MAX_STEP 8

/**
 * @name: Name of the method
 * @steps: Number of steps
 * @last_step_is_update: Indicate if the last step correspond to the update, in
 * which case no final update is computed and the last transient 'y' is used.
 * @c: Whole step weight
 * @a: Previous steps weights
 * @b: Final step weights
 * @order: Order of the final update used to adjust the step size
 * @e: Error estimate
 * @error_estimate_order: Order of the error estimate (unused)
 */
typedef struct _OdeIntMeta
{
    const char *name;
    int         steps;
    int         last_step_is_update;
    int         error_estimate;
    double      c[MAX_STEP];
    double      a[(MAX_STEP*(MAX_STEP-1))/2];
    double      b[MAX_STEP];
    int         order;
    double      e[MAX_STEP];
    int         error_estimate_order;
} OdeIntMeta;

/*
 * Butcher table
 *
 * For method that provide no error estimate, we reuse the last step, which
 * will yield no error.
 */
const OdeIntMeta INTEGRATOR_META[] =
{
    {"euler", 1, 0, 0,
        {0.0},
        {0},
        {1.0}, 1,
        {0},   0},
    {"heuns", 2, 0, 0,
        {0.0,     1.0},
        {1.0},
        {1.0/2.0, 1.0/2.0}, 2,
        {0},                0},
    {"heuns-euler", 2, 0, 1,
        {0.0,     1.0},
        {1.0},
        {1.0/2.0, 1.0/2.0}, 2,
        {1.0,     0.0},     1},
    {"bogacki-shampine", 4, 1, 1,
        {0.0,      1.0/2.0, 3.0/4.0, 1.0},
        {1.0/2.0,
         0.0,      3.0/4.0,
         2.0/9.0,  1.0/3.0, 4.0/9.0},
        {2.0/9.0,  1.0/3.0, 4.0/9.0, 0.0},     3,
        {7.0/24.0, 1.0/4.0, 1.0/3.0, 1.0/8.0}, 2},
    {"runge-kutta", 4, 0, 0,
        {0.0,     1.0/2.0, 1.0/2.0, 1.0},
        {1.0/2.0,
         0.0,     1.0/2.0,
         0.0,     0.0,     1.0},
        {1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0}, 4,
        {0},                                  0},
    {"runge-kutta-fehlberg", 6, 0, 1,
        {0.0,           1.0/4.0,        3.0/8.0,        12.0/13.0,       1.0,       1.0/2.0},
        {1.0/4.0,
         3.0/32.0,      9.0/32.0,
         1932.0/2197.0, -7200.0/2197.0, 7296.0/2197.0,
         439.0/216.0,   -8.0,           3680.0/513.0,   -845/4104.0,
         -8.0/27.0,     2.0,            -3544.0/2565.0, 1859.0/4104.0,   -11.0/40.0},
        {16.0/135.0,    0.0,            6656.0/12825.0, 28561.0/56430.0, -9.0/50.0, 2.0/55.0}, 5,
        {25.0/216.0,    0.0,            1408.0/2565.0,  2197.0/4104.0,   -1.0/5.0,  0.0},      4},
    {"cash-karp", 6, 0, 1,
        {0.0,            1.0/5.0,     3.0/10.0,        3.0/5.0,          1.0,           7.0/8.0},
        {1.0/5.0,
         3.0/40.0,       9.0/40.0,
         3.0/10.0,       -9.0/10.0,   6.0/5.0,
         -11.0/54.0,     5.0/2.0,     -70.0/27.0,      35.0/27.0,
         1631.0/55296.0, 175.0/512.0, 575.0/13824.0,   44275.0/110592.0, 253.0/4096.0},
        {37.0/378.0,     0.0,         250.0/621.0,     125.0/594.0,      0.0,           512.0/1771.0}, 5,
        {2825.0/27648.0, 0.0,         18575.0/48384.0, 13525.0/55296.0,  277.0/14336.0, 1.0/4.0},      4},
    {"dormand-prince", 7, 1, 1,
        {0.0,            1.0/5.0,         3.0/10.0,       4.0/5.0,     8.0/9.0,           1.0,          1.0},
        {1.0/5.0,
         3.0/40.0,       9.0/40.0,
         44.0/45.0,      -56.0/15.0,      32.0/9.0,
         19372.0/6561.0, -25360.0/2187.0, 64448.0/6561.0, -212.0/729.0,
         9017.0/3168.0,  -355.0/33.0,     46732.0/5247.0, 49.0/176.0,  -5103.0/18656.0,
         35.0/384.0,     0.0,             500.0/1113.0,   125.0/192.0, -2187.0/6784.0,    11.0/84.0},
        {35.0/384.0,     0.0,             500.0/1113.0,   125.0/192.0, -2187.0/6784.0,    11.0/84.0,    0.0},      5,
        {5179.0/57600.0, 0.0,             7571.0/16695.0, 393.0/640.0, -92097.0/339200.0, 187.0/2100.0, 1.0/40.0}, 4}
};

struct _OdeIntIntegrator
{
    const OdeIntMeta *integrator_meta;
    double *t;
    double *y;
    size_t  n;
    double  rtol;
    double  atol;
    double *transient_y;
    double *transient_ye;
    double *F; // dy/dt
    double  h;
};

/**
 * odeint_integrator_new:
 * @t0: Initial timestep
 * @y0: Initial value for the system which will be updated throughout the
 * integration
 * @n: Number of equations and variables (each variable has a differential equation w.r.t. time 't')
 * @rtol: Relative tolerance (i.e. number of significant digits)
 * @atol: Absolute tolerance
 */
OdeIntIntegrator *
odeint_integrator_new (OdeIntMethod  method,
                       double       *t0,
                       double       *y0,
                       size_t        n,
                       double        rtol,
                       double        atol)
{
    OdeIntIntegrator *ret = malloc (sizeof (OdeIntIntegrator));

    ret->integrator_meta = &INTEGRATOR_META[method];

    ret->t = t0;
    ret->y = y0;
    ret->n = n;
    /* we need at least 2 steps for the initial step size estimate  */
    ret->F = calloc ((ret->integrator_meta->steps + 1) * n, sizeof (double));

    /* transient states for multi-step methods */
    ret->transient_y = malloc ( n * sizeof (double) );
    ret->transient_ye = malloc ( n * sizeof (double) );

    ret->rtol = rtol;
    ret->atol = atol;

    ret->h = 0;

    assert (ret->F != NULL);

    return ret;
}

/**
 * odeint_integrator_integrate:
 * @func: The function to integrate
 * @user_data:
 * @tw: The upper integration bound, which may also be smaller than the current
 * time.
 *
 * Integrate dy/dt = @func(t, y) over [t, @tw] or [@tw, t] if @tw < t.
 */
void
odeint_integrator_integrate (OdeIntIntegrator *self,
                             OdeIntFunc        func,
                             void             *user_data,
                             double            tw)
{
    double t0 = *self->t;
    double t;

    /* transient state for the multi-step method */
    double *y  = self->transient_y;
    double *ye = self->transient_ye;

    // Estimate initial step size using order, local curvature and tolerance
    // information
    //
    // reference:
    // H. A. Watts, “Starting Step Size for an ODE Solver,” Journal of
    // Computational and Applied Mathematics 9, no. 2 (June 1, 1983): 177–91,
    // https://doi.org/10.1016/0377-0427(83)90040-7.
    if (self->h == 0)
    {
        double epsilon = sqrt (DBL_EPSILON);

        func (t0, self->y, self->F, user_data);

        /* simple Euler step */
        cblas_dcopy (self->n, self->y, 1, y, 1);
        cblas_daxpy (self->n,
                     epsilon,
                     self->F,
                     1,
                     y,
                     1);

        func (t0 + epsilon, y, self->F + self->n, user_data);

        cblas_daxpy (self->n,
                     -1,
                     self->F,
                     1,
                     self->F + self->n,
                     1);

        double y_norm = cblas_dnrm2 (self->n, self->F + self->n, 1) / epsilon;

        assert (isfinite (y_norm));

        // if the has no curvature whatsoever, the step size will be infinite,
        self->h = sqrt (2.0) * pow (self->atol, 1.0 / (self->integrator_meta->order + 1)) / sqrt (y_norm);
    }

    size_t iteration = 0;

    while (fabs (tw - *self->t) >= self->rtol * fabs (tw) + self->atol)
    {
        iteration++;

        t = *self->t;

        // ensure we always land exactly on the upper integration bound
        // if we went too far, the step size will be negative and still point
        // toward the integration bound
        #define fsign(x)((x) < 0. ? -1. : 1.)
        self->h = fsign (tw - t) * fmin (fabs (self->h), fabs (tw - t));

        int step;
        for (step = 0; step < self->integrator_meta->steps; step++)
        {
            /* restore state at the beginning of the step */
            cblas_dcopy (self->n, self->y, 1, y, 1);

            int prev_step;
            for (prev_step = 0; prev_step < step; prev_step++)
            {
                cblas_daxpy (self->n,
                             self->h * self->integrator_meta->a[(step * step - step)/2 + prev_step],
                             self->F + prev_step * self->n,
                             1,
                             y,
                             1);
            }

            /* update from previous steps */
            func (t + self->h * self->integrator_meta->c[step], y, self->F + (step * self->n), user_data);
        }

        // final update
        if (!self->integrator_meta->last_step_is_update)
        {
            /* restore state at the beginning of the step */
            cblas_dcopy (self->n, self->y, 1, y, 1);

            int step;
            for (step = 0; step < self->integrator_meta->steps; step++)
            {
                cblas_daxpy (self->n,
                             self->h * self->integrator_meta->b[step],
                             self->F + step * self->n,
                             1,
                             y,
                             1);
            }
        }

        // error estimate
        if (self->integrator_meta->error_estimate)
        {
            /* restore state at the beginning of the step */
            cblas_dcopy (self->n, self->y, 1, ye, 1);

            int step;
            for (step = 0; step < self->integrator_meta->steps; step++)
            {
                cblas_daxpy (self->n,
                             self->h * self->integrator_meta->e[step],
                             self->F + step * self->n,
                             1,
                             ye,
                             1);
            }

            // component-wise error
            double error_ratio = 0;
            size_t i;
            for (i = 0; i < self->n; i++)
            {
                double tol = self->rtol * fmax (fabs (ye[i]), fabs(y[i])) + self->atol;
                double error = fabs (ye[i] - y[i]);
                error_ratio = fmax (error_ratio, error / tol);
            }

            assert (isfinite (error_ratio));

            // check if the step is accepted
            if (error_ratio <= 1)
            {
                *self->t = t + self->h;
                cblas_dcopy (self->n, y, 1, self->y, 1);
            }

            // compute optimal step size
            self->h = self->h * pow (1.0 / error_ratio, 1.0 / self->integrator_meta->order);
        }
        else
        {
            *self->t = t + self->h;
            cblas_dcopy (self->n, y, 1, self->y, 1);
        }
    }
}

void odeint_integrator_free (OdeIntIntegrator *self)
{
    free (self->transient_y);
    free (self->transient_ye);
    free (self->F);
    free (self);
}

max_stars_repo_name: major-lab/mirbooking

max_issues_repo_name: major-lab/mirbooking

max_forks_repo_name: major-lab/mirbooking


text:
#include <config.h>
#include <stddef.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <stdio.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>

#include <gsl/gsl_fft_complex.h>

#define BASE_DOUBLE
#include "templates_on.h"
#include "compare_source.c"
#include "templates_off.h"
#undef  BASE_DOUBLE

#define BASE_FLOAT
#include "templates_on.h"
#include "compare_source.c"
#include "templates_off.h"
#undef  BASE_FLOAT

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/**
 *
 * @file testing_zgetmi.c
 *
 *  PLASMA testings module
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 *  This program tests the implementation of the inplace format
 *  conversion based on the GKK algorithm by Gustavson, Karlsson,
 *  Kagstrom.
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 *
 * @precisions normal z -> c d s
 *
 * Purpose :
 *    Test the in place matrix transposition with the 6 different
 *    formats of input.
 **/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <lapacke.h>
#include <plasma.h>
#include "testing_zmain.h"

static int check_solution(int m, int n, int mb, int nb,
                   PLASMA_Complex64_t *A, PLASMA_Complex64_t *B,
                   int (*mapA)(int, int, int, int, int, int)) {
    int i, j;

    for( j=0; j<n; j++) {
        for (i=0; i<m; i++) {
            if (A[ mapA(m, n, mb, nb, i, j) ] != B[ mapA(n, m, mb, nb, j, i) ] ) {
                return -1;
            }
        }
    }
    return 0;
}

int testing_zgetmi(int argc, char **argv){

    PLASMA_Complex64_t *A, *B;
    int m, n, mb, nb;
    int i, ret, size;

    /* Check for number of arguments*/
    if (argc != 4){
        USAGE("GETMI", "M N MB NB ntdbypb with \n",
              "   - M       : the number of rows of the matrix    \n"
              "   - N       : the number of columns of the matrix \n"
              "   - MB      : the number of rows of each block    \n"
              "   - NB      : the number of columns of each block \n");
        return -1;
    }

    m      = atoi(argv[0]);
    n      = atoi(argv[1]);
    mb     = atoi(argv[2]);
    nb     = atoi(argv[3]);

    size = m*n*sizeof(PLASMA_Complex64_t);
    A = (PLASMA_Complex64_t *)malloc(size);
    B = (PLASMA_Complex64_t *)malloc(size);
    LAPACKE_zlarnv_work(1, ISEED, m*n, A);

    for(i=0; i<6; i++) {
        memcpy(B, A, size);

        printf(" - TESTING ZGETMI (%4s) ...", formatstr[i]);
        ret = PLASMA_zgetmi( m, n, A, format[i], mb, nb );

        if (ret != PLASMA_SUCCESS) {
            printf("Failed\n");
            continue;
        }

        if ( check_solution(m, n, mb, nb, B, A,
                            (int (*)(int, int, int, int, int, int))formatmap[i]) == 0 )
            printf("............ PASSED !\n");
        else
            printf("... FAILED !\n");
    }

    free( A ); free( B );

    return 0;
}


text:
#include<iostream>
#include<string.h>
#include<vector>
#include<set>
#include<algorithm>
#include<fstream>
#include<cstdlib>
#include<time.h>
#include<iterator>
#include<sys/time.h>
#include<stdio.h>
#include<numeric>
#include <stdexcept>
#include<math.h>
#include<limits>
#include<random>
#include<ctime>
#include<map>
#include<omp.h>
#include <gsl/gsl_rng.h>
#include<pthread.h>
#include<chrono>

extern int GenNum;
extern std::vector<int> ArgPop;
//double PopulationFitnessTable = new double[numberofSNPs*diploidSize];
extern int popSize;
extern double RecombRate;
extern double MutationRate;
extern int m;
extern int epiflag;
extern int ChromLen;

max_stars_repo_name: SimRA-ARG/SimRA

max_issues_repo_name: SimRA-ARG/SimRA

max_forks_repo_name: SimRA-ARG/SimRA