text:
/* block/gsl_block_ushort.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#ifndef __GSL_BLOCK_USHORT_H__
#define __GSL_BLOCK_USHORT_H__

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

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

struct gsl_block_ushort_struct
{
  size_t size;
  unsigned short *data;
};

typedef struct gsl_block_ushort_struct gsl_block_ushort;

gsl_block_ushort *gsl_block_ushort_alloc (const size_t n);
gsl_block_ushort *gsl_block_ushort_calloc (const size_t n);
void gsl_block_ushort_free (gsl_block_ushort * b);

int gsl_block_ushort_fread (FILE * stream, gsl_block_ushort * b);
int gsl_block_ushort_fwrite (FILE * stream, const gsl_block_ushort * b);
int gsl_block_ushort_fscanf (FILE * stream, gsl_block_ushort * b);
int gsl_block_ushort_fprintf (FILE * stream, const gsl_block_ushort * b, const char *format);

int gsl_block_ushort_raw_fread (FILE * stream, unsigned short * b, const size_t n, const size_t stride);
int gsl_block_ushort_raw_fwrite (FILE * stream, const unsigned short * b, const size_t n, const size_t stride);
int gsl_block_ushort_raw_fscanf (FILE * stream, unsigned short * b, const size_t n, const size_t stride);
int gsl_block_ushort_raw_fprintf (FILE * stream, const unsigned short * b, const size_t n, const size_t stride, const char *format);

size_t gsl_block_ushort_size (const gsl_block_ushort * b);
unsigned short * gsl_block_ushort_data (const gsl_block_ushort * b);

__END_DECLS

#endif /* __GSL_BLOCK_USHORT_H__ */

ext: h


text:
static char help[] =
"Solve the linear biharmonic equation in 2D.  Equation is\n"
"  Lap^2 u = f\n"
"where Lap = - grad^2 is the positive Laplacian, equivalently\n"
"  u_xxxx + 2 u_xxyy + u_yyyy = f(x,y)\n"
"Domain is unit square  S = (0,1)^2.   Boundary conditions are homogeneous\n"
"simply-supported:  u = 0,  Lap u = 0.  The equation is rewritten as a\n"
"2x2 block system with SPD Laplacian blocks on the diagonal:\n"
"   | Lap |  0  |  | v |   | f | \n"
"   |-----|-----|  |---| = |---| \n"
"   | -I  | Lap |  | u |   | 0 | \n"
"Includes manufactured, polynomial exact solution.  The discretization is\n"
"structured-grid (DMDA) finite differences.  Includes analytical Jacobian.\n"
"Recommended preconditioning combines fieldsplit:\n"
"   -pc_type fieldsplit -pc_fieldsplit_type multiplicative|additive \n"
"with multigrid as the preconditioner for the diagonal blocks:\n"
"   -fieldsplit_v_pc_type mg|gamg -fieldsplit_u_pc_type mg|gamg\n"
"(GMG requires setting levels and Galerkin coarsening.)  One can also do\n"
"monolithic multigrid (-pc_type mg|gamg).\n\n";

#include <petsc.h>

typedef struct {
    PetscReal  v, u;
} Field;

typedef struct {
    PetscReal  (*f)(PetscReal x, PetscReal y);  // right-hand side
} BiharmCtx;

static PetscReal c(PetscReal x) {
    return x*x*x * (1.0-x)*(1.0-x)*(1.0-x);
}

static PetscReal ddc(PetscReal x) {
    return 6.0 * x * (1.0-x) * (1.0 - 5.0 * x + 5.0 * x*x);
}

static PetscReal d4c(PetscReal x) {
    return - 72.0 * (1.0 - 5.0 * x + 5.0 * x*x);
}

static PetscReal u_exact_fcn(PetscReal x, PetscReal y) {
    return c(x) * c(y);
}

static PetscReal lap_u_exact_fcn(PetscReal x, PetscReal y) {
    return - ddc(x) * c(y) - c(x) * ddc(y);  // Lap u = - grad^2 u
}

static PetscReal f_fcn(PetscReal x, PetscReal y) {
    return d4c(x) * c(y) + 2.0 * ddc(x) * ddc(y) + c(x) * d4c(y);  // Lap^2 u = grad^4 u
}

extern PetscErrorCode FormExactWLocal(DMDALocalInfo*, Field**, BiharmCtx*);
extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*, Field**, Field **FF, BiharmCtx*);
extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*, Field**, Mat, Mat, BiharmCtx*);

int main(int argc,char **argv) {
    PetscErrorCode ierr;
    DM             da;
    SNES           snes;
    Vec            w, w_initial, w_exact;
    BiharmCtx      user;
    Field          **aW;
    PetscReal      normv, normu, errv, erru;
    DMDALocalInfo  info;

    ierr = PetscInitialize(&argc,&argv,NULL,help); if (ierr) return ierr;

    user.f = &f_fcn;
    ierr = DMDACreate2d(PETSC_COMM_WORLD,
                        DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR,
                        3,3,PETSC_DECIDE,PETSC_DECIDE,
                        2,1,              // degrees of freedom, stencil width
                        NULL,NULL,&da); CHKERRQ(ierr);
    ierr = DMSetApplicationContext(da,&user); CHKERRQ(ierr);
    ierr = DMSetFromOptions(da); CHKERRQ(ierr);
    ierr = DMSetUp(da); CHKERRQ(ierr);  // this must be called BEFORE SetUniformCoordinates
    ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,-1.0,-1.0); CHKERRQ(ierr);
    ierr = DMDASetFieldName(da,0,"v"); CHKERRQ(ierr);
    ierr = DMDASetFieldName(da,1,"u"); CHKERRQ(ierr);

    ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
    ierr = SNESSetDM(snes,da); CHKERRQ(ierr);
    ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,
               (DMDASNESFunction)FormFunctionLocal,&user); CHKERRQ(ierr);
    ierr = DMDASNESSetJacobianLocal(da,
               (DMDASNESJacobian)FormJacobianLocal,&user); CHKERRQ(ierr);
    ierr = SNESSetType(snes,SNESKSPONLY); CHKERRQ(ierr);
    ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);

    ierr = DMGetGlobalVector(da,&w_initial); CHKERRQ(ierr);
    ierr = VecSet(w_initial,0.0); CHKERRQ(ierr);
    ierr = SNESSolve(snes,NULL,w_initial); CHKERRQ(ierr);
    ierr = DMRestoreGlobalVector(da,&w_initial); CHKERRQ(ierr);
    ierr = DMDestroy(&da); CHKERRQ(ierr);

    ierr = SNESGetSolution(snes,&w); CHKERRQ(ierr);
    ierr = SNESGetDM(snes,&da); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);

    ierr = DMCreateGlobalVector(da,&w_exact); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(da,w_exact,&aW); CHKERRQ(ierr);
    ierr = FormExactWLocal(&info,aW,&user); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(da,w_exact,&aW); CHKERRQ(ierr);
    ierr = VecStrideNorm(w_exact,0,NORM_INFINITY,&normv); CHKERRQ(ierr);
    ierr = VecStrideNorm(w_exact,1,NORM_INFINITY,&normu); CHKERRQ(ierr);
    ierr = VecAXPY(w,-1.0,w_exact); CHKERRQ(ierr);
    ierr = VecStrideNorm(w,0,NORM_INFINITY,&errv); CHKERRQ(ierr);
    ierr = VecStrideNorm(w,1,NORM_INFINITY,&erru); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,
        "done on %d x %d grid ...\n"
        "  errors |v-vex|_inf/|vex|_inf = %.5e, |u-uex|_inf/|uex|_inf = %.5e\n",
        info.mx,info.my,errv/normv,erru/normu); CHKERRQ(ierr);

    ierr = VecDestroy(&w_exact); CHKERRQ(ierr);
    ierr = SNESDestroy(&snes); CHKERRQ(ierr);
    return PetscFinalize();
}

PetscErrorCode FormExactWLocal(DMDALocalInfo *info, Field **aW, BiharmCtx *user) {
    PetscErrorCode ierr;
    PetscInt   i, j;
    PetscReal  xymin[2], xymax[2], hx, hy, x, y;
    ierr = DMGetBoundingBox(info->da,xymin,xymax); CHKERRQ(ierr);
    hx = (xymax[0] - xymin[0]) / (info->mx - 1);
    hy = (xymax[1] - xymin[1]) / (info->my - 1);
    for (j = info->ys; j < info->ys + info->ym; j++) {
        y = j * hy;
        for (i = info->xs; i < info->xs + info->xm; i++) {
            x = i * hx;
            aW[j][i].u = u_exact_fcn(x,y);
            aW[j][i].v = lap_u_exact_fcn(x,y);
        }
    }
    return 0;
}

PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field **aW,
                                 Field **FF, BiharmCtx *user) {
    PetscErrorCode ierr;
    PetscInt   i, j;
    PetscReal  xymin[2], xymax[2], hx, hy, darea, scx, scy, scdiag, x, y,
               ve, vw, vn, vs, ue, uw, un, us;
    ierr = DMGetBoundingBox(info->da,xymin,xymax); CHKERRQ(ierr);
    hx = (xymax[0] - xymin[0]) / (info->mx - 1);
    hy = (xymax[1] - xymin[1]) / (info->my - 1);
    darea = hx * hy;               // multiply FD equations by this
    scx = hy / hx;
    scy = hx / hy;
    scdiag = 2.0 * (scx + scy);    // diagonal scaling
    for (j = info->ys; j < info->ys + info->ym; j++) {
        y = xymin[1] + j * hy;
        for (i = info->xs; i < info->xs + info->xm; i++) {
            x = xymin[0] + i * hx;
            if (i==0 || i==info->mx-1 || j==0 || j==info->my-1) {
                FF[j][i].v = scdiag * aW[j][i].v;
                FF[j][i].u = scdiag * aW[j][i].u;
            } else {
                ve = (i+1 == info->mx-1) ? 0.0 : aW[j][i+1].v;
                vw = (i-1 == 0)          ? 0.0 : aW[j][i-1].v;
                vn = (j+1 == info->my-1) ? 0.0 : aW[j+1][i].v;
                vs = (j-1 == 0)          ? 0.0 : aW[j-1][i].v;
                FF[j][i].v = scdiag * aW[j][i].v - scx * (vw + ve) - scy * (vs + vn)
                             - darea * (*(user->f))(x,y);
                ue = (i+1 == info->mx-1) ? 0.0 : aW[j][i+1].u;
                uw = (i-1 == 0)          ? 0.0 : aW[j][i-1].u;
                un = (j+1 == info->my-1) ? 0.0 : aW[j+1][i].u;
                us = (j-1 == 0)          ? 0.0 : aW[j-1][i].u;
                FF[j][i].u = - darea * aW[j][i].v
                             + scdiag * aW[j][i].u - scx * (uw + ue) - scy * (us + un);
            }
        }
    }
    ierr = PetscLogFlops(18.0*info->xm*info->ym);CHKERRQ(ierr);
    return 0;
}

PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, Field **aW,
                                 Mat J, Mat Jpre, BiharmCtx *user) {
    PetscErrorCode ierr;
    PetscInt     i, j, c, ncol;
    PetscReal    xymin[2], xymax[2], hx, hy, darea, scx, scy, scdiag, val[6];
    MatStencil   col[6], row;

    ierr = DMGetBoundingBox(info->da,xymin,xymax); CHKERRQ(ierr);
    hx = (xymax[0] - xymin[0]) / (info->mx - 1);
    hy = (xymax[1] - xymin[1]) / (info->my - 1);
    darea = hx * hy;               // multiply FD equations by this
    scx = hy / hx;
    scy = hx / hy;
    scdiag = 2.0 * (scx + scy);    // diagonal scaling
    for (j = info->ys; j < info->ys + info->ym; j++) {
        row.j = j;
        for (i = info->xs; i < info->xs + info->xm; i++) {
            row.i = i;
            for (c = 0; c < 2; c++) { // v,u equations are c=0,1
                row.c = c;
                col[0].c = c;     col[0].i = i;       col[0].j = j;
                val[0] = scdiag;
                if (i==0 || i==info->mx-1 || j==0 || j==info->my-1) {
                    ierr = MatSetValuesStencil(Jpre,1,&row,1,col,val,INSERT_VALUES);
                        CHKERRQ(ierr);
                } else {
                    ncol = 1;
                    if (i+1 < info->mx-1) {
                        col[ncol].c = c;  col[ncol].i = i+1;  col[ncol].j = j;
                        val[ncol++] = -scx;
                    }
                    if (i-1 > 0) {
                        col[ncol].c = c;  col[ncol].i = i-1;  col[ncol].j = j;
                        val[ncol++] = -scx;
                    }
                    if (j+1 < info->my-1) {
                        col[ncol].c = c;  col[ncol].i = i;    col[ncol].j = j+1;
                        val[ncol++] = -scy;
                    }
                    if (j-1 > 0) {
                        col[ncol].c = c;  col[ncol].i = i;    col[ncol].j = j-1;
                        val[ncol++] = -scy;
                    }
                    if (c == 1) {  // u equation;  has off-diagonal block entry
                        col[ncol].c = 0;
                        col[ncol].i = i;  col[ncol].j = j;
                        val[ncol++] = - darea;
                    }
                    ierr = MatSetValuesStencil(Jpre,1,&row,ncol,col,val,INSERT_VALUES);
                        CHKERRQ(ierr);
                }
            }
        }
    }

    ierr = MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    if (J != Jpre) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}

ext: c


text:
#ifndef SIMPROP_UTILS_GSL_H
#define SIMPROP_UTILS_GSL_H

#include <gsl/gsl_deriv.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_roots.h>

#include <cassert>
#include <cmath>
#include <functional>
#include <iostream>
#include <vector>

namespace simprop {
namespace utils {

// Axis
std::vector<double> LinAxis(const double &min, const double &max, const size_t &size);

std::vector<double> LogAxis(const double &min, const double &max, const size_t &size);

inline bool isInside(double x, const std::vector<double> &X) {
  return (x >= X.front() && x <= X.back());
};

double interpolate(double x, const std::vector<double> &X, const std::vector<double> &Y);

double cspline(double x, const std::vector<double> &X, const std::vector<double> &Y);

double interpolateEquidistant(double x, double lo, double hi, const std::vector<double> &Y);

double interpolate2d(double x, double y, const std::vector<double> &X, const std::vector<double> &Y,
                     const std::vector<double> &Z);

template <typename T>
T deriv(std::function<T(T)> f, T x, double rel_error = 1e-4) {
  double result;
  double abs_error = 0.0;  // disabled
  gsl_function F;

  F.function = [](double x, void *vf) -> double {
    auto &func = *static_cast<std::function<double(double)> *>(vf);
    return func(x);
  };
  F.params = &f;

  gsl_deriv_central(&F, static_cast<double>(x), rel_error, &result, &abs_error);

  return T(result);
}

template <typename T>
T deriv5pt(std::function<T(T)> f, T x, T h) {
  auto result = -f(x + 2 * h) + 8 * f(x + h) - 8 * f(x - h) + f(x - 2 * h);
  return T(result) / 12 / h;
}

template <typename T>
T QAGIntegration(std::function<T(T)> f, T start, T stop, int LIMIT, double rel_error = 1e-4) {
  double a = static_cast<double>(start);
  double b = static_cast<double>(stop);
  double abs_error = 0.0;  // disabled
  int key = GSL_INTEG_GAUSS31;
  double result;
  double error;

  gsl_function F;
  F.function = [](double x, void *vf) -> double {
    auto &func = *static_cast<std::function<double(double)> *>(vf);
    return func(x);
  };
  F.params = &f;

  gsl_integration_workspace *workspace_ptr = gsl_integration_workspace_alloc(LIMIT);
  gsl_integration_qag(&F, a, b, abs_error, rel_error, LIMIT, key, workspace_ptr, &result, &error);
  gsl_integration_workspace_free(workspace_ptr);

  return T(result);
}

template <typename T>
T QAGSIntegration(std::function<T(T)> f, T start, T stop, int LIMIT, double rel_error = 1e-4) {
  double a = static_cast<double>(start);
  double b = static_cast<double>(stop);
  double abs_error = 0.0;  // disabled
  double result;
  double error;

  gsl_function F;
  F.function = [](double x, void *vf) -> double {
    auto &func = *static_cast<std::function<double(double)> *>(vf);
    return func(x);
  };
  F.params = &f;

  gsl_integration_workspace *workspace_ptr = gsl_integration_workspace_alloc(LIMIT);
  gsl_integration_qags(&F, a, b, abs_error, rel_error, LIMIT, workspace_ptr, &result, &error);
  gsl_integration_workspace_free(workspace_ptr);

  return T(result);
}

template <typename T>
T simpsonIntegration(std::function<T(T)> f, T start, T stop, int N = 100) {
  const T a = start;
  const T b = stop;

  const T h = (b - a) / N;
  const T XI0 = f(a) + f(b);

  T XI1 = 0, XI2 = 0;

  for (int i = 1; i < N; ++i) {
    const T X = a + i * h;
    if (i % 2 == 0)
      XI2 = XI2 + f(X);
    else
      XI1 = XI1 + f(X);
  }

  return h * (XI0 + 2 * XI2 + 4 * XI1) / 3.0;
}

template <typename T>
T rootFinder(std::function<T(T)> f, T xLower, T xUpper, int maxIter, double relError = 1e-4) {
  int status;
  int iter = 0;
  const gsl_root_fsolver_type *solverType;
  gsl_root_fsolver *solver;

  T r = 0;

  gsl_function F;
  F.function = [](double x, void *vf) -> double {
    auto &func = *static_cast<std::function<double(double)> *>(vf);
    return func(x);
  };
  F.params = &f;

  solverType = gsl_root_fsolver_brent;
  solver = gsl_root_fsolver_alloc(solverType);
  gsl_root_fsolver_set(solver, &F, xLower, xUpper);

  do {
    iter++;
    status = gsl_root_fsolver_iterate(solver);
    r = (T)gsl_root_fsolver_root(solver);
    xLower = gsl_root_fsolver_x_lower(solver);
    xUpper = gsl_root_fsolver_x_upper(solver);
    status = gsl_root_test_interval(xLower, xUpper, 0, relError);
  } while (status == GSL_CONTINUE && iter < maxIter);

  gsl_root_fsolver_free(solver);

  return r;
}

template <typename T>
T rk4fixed(std::function<T(T, T)> dydx, T yStart, T xStart, T xEnd, T h) {
  size_t steps = (size_t)((xEnd - xStart) / h);
  T x = xStart;
  T y = yStart;
  while (steps--) {
    auto k1 = h * dydx(x, y);
    auto k2 = h * dydx(x + 0.5 * h, y + 0.5 * k1);
    auto k3 = h * dydx(x + 0.5 * h, y + 0.5 * k2);
    auto k4 = h * dydx(x + h, y + k3);
    y += (k1 + 2. * k2 + 2. * k3 + k4) / 6.;
    x += h;
  }
  return y;
}

template <typename T>
T odeiv(std::function<T(T, T)> dydx, T yStart, T xStart, T xEnd, T rel_error = 1e-4) {
  const size_t NEQS = 1;
  double abs_error = 0.0;  // disabled

  const gsl_odeiv2_step_type *stepType = gsl_odeiv2_step_rkf45;
  gsl_odeiv2_step *s = gsl_odeiv2_step_alloc(stepType, NEQS);
  gsl_odeiv2_control *c = gsl_odeiv2_control_y_new(abs_error, static_cast<double>(rel_error));
  gsl_odeiv2_evolve *e = gsl_odeiv2_evolve_alloc(NEQS);

  double y[1] = {yStart};

  auto func = [](double t, const double *y, double *f, void *params) -> int {
    auto dydx = *static_cast<std::function<double(double, double)> *>(params);
    f[0] = dydx(t, y[0]);
    return GSL_SUCCESS;
  };

  gsl_odeiv2_system sys = {func, nullptr, NEQS, &dydx};

  double t = xStart, t1 = xEnd;
  double h = 1e-5 * (xEnd - xStart);

  while (t < t1) {
    int status = gsl_odeiv2_evolve_apply(e, c, s, &sys, &t, t1, &h, y);
    if (status != GSL_SUCCESS) break;
  }

  gsl_odeiv2_evolve_free(e);
  gsl_odeiv2_control_free(c);
  gsl_odeiv2_step_free(s);
  return y[0];
}

}  // namespace utils
}  // namespace simprop

#endif

max_stars_repo_name: carmeloevoli/SimProp-beta

max_issues_repo_name: carmeloevoli/SimProp-beta

max_forks_repo_name: carmeloevoli/SimProp-beta


text:
/**
 *
 * @file qwrapper_chegst.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:44:59 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

#undef REAL
#define COMPLEX

/***************************************************************************//**
 *
 **/
void QUARK_CORE_chegst(Quark *quark, Quark_Task_Flags *task_flags,
                       int itype, PLASMA_enum uplo, int n,
                       PLASMA_Complex32_t *A, int lda,
                       PLASMA_Complex32_t *B, int ldb,
                       PLASMA_sequence *sequence, PLASMA_request *request,
                       int iinfo)
{
    QUARK_Insert_Task(quark, CORE_chegst_quark, task_flags,
        sizeof(int),                        &itype,      VALUE,
        sizeof(PLASMA_enum),                &uplo,      VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(PLASMA_Complex32_t)*lda*n,    A,             INOUT,
        sizeof(int),                        &lda,       VALUE,
#ifdef COMPLEX
        sizeof(PLASMA_Complex32_t)*ldb*n,    B,             INOUT,
#else
        sizeof(PLASMA_Complex32_t)*ldb*n,    B,             INPUT,
#endif
        sizeof(int),                        &ldb,       VALUE,
        sizeof(PLASMA_sequence*),           &sequence,  VALUE,
        sizeof(PLASMA_request*),            &request,   VALUE,
        sizeof(int),                        &iinfo,     VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_chegst_quark = PCORE_chegst_quark
#define CORE_chegst_quark PCORE_chegst_quark
#endif
void CORE_chegst_quark(Quark *quark)
{
    int itype;
    PLASMA_enum uplo;
    int n;
    PLASMA_Complex32_t *A;
    int lda;
    PLASMA_Complex32_t *B;
    int ldb;
    PLASMA_sequence *sequence;
    PLASMA_request *request;
    int iinfo;

    int info;

    quark_unpack_args_10(quark, itype, uplo, n, A, lda, B, ldb, sequence, request, iinfo);
    info = LAPACKE_chegst_work(
        LAPACK_COL_MAJOR,
        itype,
        lapack_const(uplo),
        n, A, lda, B, ldb);
    if (sequence->status == PLASMA_SUCCESS && info != 0)
      plasma_sequence_flush(quark, sequence, request, iinfo+info);
}


text:
#pragma once

#include <cstdint>
#include <functional>
#include <memory>
#include <sstream>
#include <string>
#include <utility>
#include <vector>

#include <gsl/gsl>
#include <nonstd/optional.hpp>

#include "chainerx/array_body.h"
#include "chainerx/array_fwd.h"
#include "chainerx/array_index.h"
#include "chainerx/array_node.h"
#include "chainerx/array_repr.h"
#include "chainerx/axes.h"
#include "chainerx/constant.h"
#include "chainerx/context.h"
#include "chainerx/device.h"
#include "chainerx/dtype.h"
#include "chainerx/enum.h"
#include "chainerx/error.h"
#include "chainerx/graph.h"
#include "chainerx/scalar.h"
#include "chainerx/shape.h"
#include "chainerx/strides.h"

namespace chainerx {
namespace internal {

BackpropId GetArrayBackpropId(const Array& array, const nonstd::optional<BackpropId>& backprop_id);

Array MakeArray(const Shape& shape, const Strides& strides, Dtype dtype, Device& device, std::shared_ptr<void> data, int64_t offset = 0);

inline const std::shared_ptr<ArrayBody>& GetArrayBody(const Array& array);

inline std::shared_ptr<ArrayBody>&& MoveArrayBody(Array&& array);

}  // namespace internal

// The user interface of multi-dimensional arrays.
//
// This wraps an ArrayBody, providing accessors, an interface for graph operations and differentiable operations.
class Array {
public:
    Array() = default;

    // TODO(hvy): Consider making this contructor private and prohibit body from being null (assert that given body is not null).
    explicit Array(std::shared_ptr<internal::ArrayBody> body) : body_{std::move(body)} {
        if (body_ == nullptr) {
            throw ChainerxError{"Cannot create an array from null."};
        }
    }

    // Copy constructor that copies the pointer to the body instead of the body itself.
    //
    // Use MakeView if you want to clone the body.
    Array(const Array& other) = default;
    Array(Array&& other) = default;

    // Assign operators just replace the body. They do not copy data between arrays.
    Array& operator=(const Array&) = default;
    Array& operator=(Array&& other) = default;

    Array operator-() const;

    Array operator==(const Array& rhs) const;
    Array operator!=(const Array& rhs) const;
    Array operator>(const Array& rhs) const;
    Array operator>=(const Array& rhs) const;
    Array operator<(const Array& rhs) const;
    Array operator<=(const Array& rhs) const;

    Array& operator+=(const Array& rhs);
    Array& operator+=(Scalar rhs);
    Array& operator-=(const Array& rhs);
    Array& operator-=(Scalar rhs);
    Array& operator*=(const Array& rhs);
    Array& operator*=(Scalar rhs);
    Array& operator/=(const Array& rhs);
    Array& operator/=(Scalar rhs);

    const Array& operator+=(const Array& rhs) const;
    const Array& operator+=(Scalar rhs) const;
    const Array& operator-=(const Array& rhs) const;
    const Array& operator-=(Scalar rhs) const;
    const Array& operator*=(const Array& rhs) const;
    const Array& operator*=(Scalar rhs) const;
    const Array& operator/=(const Array& rhs) const;
    const Array& operator/=(Scalar rhs) const;

    Array operator+(const Array& rhs) const;
    Array operator+(Scalar rhs) const;
    Array operator-(const Array& rhs) const;
    Array operator-(Scalar rhs) const;
    Array operator*(const Array& rhs) const;
    Array operator*(Scalar rhs) const;
    Array operator/(const Array& rhs) const;
    Array operator/(Scalar rhs) const;

    // Returns a view selected with the indices.
    Array At(const std::vector<ArrayIndex>& indices) const;

    // Returns a transposed view of the array.
    Array Transpose(const OptionalAxes& axes = nonstd::nullopt) const;

    // Returns a reshaped array.
    // TODO(niboshi): Support shape with dimension -1.
    Array Reshape(const Shape& newshape) const;

    // Returns a squeezed array with unit-length axes removed.
    //
    // If no axes are specified, all axes of unit-lengths are removed.
    // If no axes can be removed, an array with aliased data is returned.
    Array Squeeze(const OptionalAxes& axis = nonstd::nullopt) const;

    // Broadcasts the array to the specified shape.
    // Returned array is always a view to this array.
    Array BroadcastTo(const Shape& shape) const;

    // Returns the indices of the maximum values along the given axis.
    Array ArgMax(const OptionalAxes& axis = nonstd::nullopt) const;

    // Returns a sum of the array.
    // If `axis` is set, it will be summed over the specified axes.
    // Otherwise, it will be summed over all the existing axes.
    // Note: When implementing chainerx::Sum(), be careful of the semantics of the default value of `keepdims`. See NumPy documentation.
    Array Sum(const OptionalAxes& axis = nonstd::nullopt, bool keepdims = false) const;

    // Returns the maximum value of the array.
    // If `axis` is set, the maximum value is chosen along the specified axes.
    // Otherwise, all the elements are searched at once.
    Array Max(const OptionalAxes& axis = nonstd::nullopt, bool keepdims = false) const;

    Array Mean(const OptionalAxes& axis = nonstd::nullopt, bool keepdims = false) const;

    Array Var(const OptionalAxes& axis = nonstd::nullopt, bool keepdims = false) const;

    // Returns a dot product of the array with another one.
    Array Dot(const Array& b) const;

    // Takes elements specified by indices from the array.
    //
    // TODO(niboshi): Support Scalar and StackVector as indices.
    // TODO(niboshi): Support axis=None behavior in NumPy.
    // TODO(niboshi): Support indices dtype other than int64.
    Array Take(const Array& indices, int8_t axis) const;

    // Creates a copy.
    // It will be connected to all the graphs.
    // It will be always C-contiguous.
    Array Copy() const;

    // Creates a view.
    // It creates a new array node and connects graphs.
    Array MakeView() const;

    // Transfers the array to another device. It will be connected to all the graphs.
    //
    // If the destination is the same device, an array with aliased data is returned.
    // Otherwise, a C-contiguous Array will be created on the target device.
    // TODO(niboshi): Currently control over whether to make an alias is not supported.
    Array ToDevice(Device& dst_device) const;

    // Transfer the array to the native device. It will be connected to all the graphs.
    //
    // This is a wrapper function which calls Array::ToDevice with the native:0 device.
    // See also: Array::ToDevice();
    Array ToNative() const;

    // Creates a copy or a view. It will be disconnected from all the graphs.
    // If `kind` is `CopyKind::kCopy`, the returned array will be always C-contiguous.
    Array AsGradStopped(CopyKind kind = CopyKind::kView) const;

    // Creates a copy or a view. It will be disconnected from the specified graphs.
    // If `kind` is `CopyKind::kCopy`, the returned array will be always C-contiguous.
    Array AsGradStopped(gsl::span<const BackpropId> backprop_ids, CopyKind kind = CopyKind::kView) const;
    Array AsGradStopped(std::initializer_list<const BackpropId> backprop_ids, CopyKind kind = CopyKind::kView) const {
        return AsGradStopped(gsl::span<const BackpropId>{backprop_ids.begin(), backprop_ids.end()}, kind);
    }

    // Casts to a specified type.
    // By default, always returns a newly allocated array. If `copy` is false,
    // and the dtype requirement is satisfied, the input array is returned instead of a copy.
    Array AsType(Dtype dtype, bool copy = true) const;

    void Fill(Scalar value) const;

    // Returns the gradient of the array.
    //
    // ChainerxError is thrown if the array is constant with respect to the computation for the specified backprop ID.
    // ChainerxError is thrown if the array is not flagged as requiring gradient.
    // This function ignores no/force-backprop mode.
    const nonstd::optional<Array>& GetGrad(const nonstd::optional<BackpropId>& backprop_id = nonstd::nullopt) const;

    // Sets the gradient of the array.
    // This function also flags the array as requiring gradient, so that preceding GetGrad() can return the gradient.
    //
    // ChainerxError is thrown if the array is constant with respect to the computation for the specified backprop ID.
    // This function ignores no/force-backprop mode.
    void SetGrad(Array grad, const nonstd::optional<BackpropId>& backprop_id = nonstd::nullopt) const;

    // Clears the gradient of the array if set.
    // This function does not change the state of the array other than that. For example, if the array is flagged as requiring gradient,
    // that will not change.
    //
    // ChainerxError is thrown if the array is constant with respect to the computation for the specified backprop ID.
    // This function ignores no/force-backprop mode.
    void ClearGrad(const nonstd::optional<BackpropId>& backprop_id = nonstd::nullopt) const;

    // Returns whether the array needs to backprop.
    //
    // If no-backprop mode is set with respect to the specified backprop ID, this function returns false.
    bool IsBackpropRequired(const nonstd::optional<BackpropId>& backprop_id = nonstd::nullopt) const;
    bool IsBackpropRequired(AnyGraph any_graph) const;

    // Returns whether the array is flagged to compute the gradient during backprop.
    //
    // This function ignores no/force-backprop mode.
    bool IsGradRequired(const nonstd::optional<BackpropId>& backprop_id = nonstd::nullopt) const;

    // Flags the array to compute the gradient during backprop.
    // If the array is constant with respect to the computation of the backprop ID, this function makes the array non-constant.
    //
    // This function ignores no/force-backprop mode.
    const Array& RequireGrad(const nonstd::optional<BackpropId>& backprop_id = nonstd::nullopt) const {
        return RequireGradImpl(*this, backprop_id);
    }

    Array& RequireGrad(const nonstd::optional<BackpropId>& backprop_id = nonstd::nullopt) { return RequireGradImpl(*this, backprop_id); }

    int64_t GetTotalSize() const { return body_->GetTotalSize(); }

    int64_t GetNBytes() const { return body_->GetNBytes(); }

    int64_t GetItemSize() const { return body_->GetItemSize(); }

    bool IsContiguous() const { return body_->IsContiguous(); }

    std::string ToString() const;

    Context& context() const { return body_->device().context(); }

    Dtype dtype() const { return body_->dtype(); }

    Device& device() const { return body_->device(); }

    int8_t ndim() const { return body_->ndim(); }

    const Shape& shape() const { return body_->shape(); }

    const Strides& strides() const { return body_->strides(); }

    const std::shared_ptr<void>& data() const { return body_->data(); }

    void* raw_data() const { return body_->data().get(); }

    int64_t offset() const { return body_->offset(); }

private:
    friend Array internal::MakeArray(
            const Shape& shape, const Strides& strides, Dtype dtype, Device& device, std::shared_ptr<void> data, int64_t offset);
    friend const std::shared_ptr<internal::ArrayBody>& internal::GetArrayBody(const Array& array);
    friend std::shared_ptr<internal::ArrayBody>&& internal::MoveArrayBody(Array&& array);

    Array(const Shape& shape, const Strides& strides, Dtype dtype, Device& device, std::shared_ptr<void> data, int64_t offset = 0);

    template <typename T>
    static T& RequireGradImpl(T& array, const nonstd::optional<BackpropId>& backprop_id);

    std::shared_ptr<internal::ArrayBody> body_;
};

Array operator+(Scalar lhs, const Array& rhs);
Array operator-(Scalar lhs, const Array& rhs);
Array operator*(Scalar lhs, const Array& rhs);
// TODO(hvy): Implement Scalar / Array using e.g. multiplication with reciprocal.

namespace internal {

inline const std::shared_ptr<ArrayBody>& GetArrayBody(const Array& array) { return array.body_; }

inline std::shared_ptr<ArrayBody>&& MoveArrayBody(Array&& array) { return std::move(array.body_); }

std::vector<std::shared_ptr<ArrayBody>> MoveArrayBodies(std::vector<Array>&& arrays);

std::vector<std::shared_ptr<ArrayBody>> MoveArrayBodies(std::vector<nonstd::optional<Array>>&& arrays);

}  // namespace internal

void DebugDumpComputationalGraph(
        std::ostream& os,
        const Array& array,
        const nonstd::optional<BackpropId>& backprop_id,
        int indent = 0,
        const std::vector<std::pair<ConstArrayRef, std::string>>& array_name_map = {});

}  // namespace chainerx

max_stars_repo_name: hitsgub/chainer

max_issues_repo_name: hitsgub/chainer

max_forks_repo_name: hitsgub/chainer


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

#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_expint.h>
#include <gsl/gsl_roots.h>

#include "ccl.h"

// Analytic FT of NFW profile, from Cooray & Sheth (2002; Section 3 of https://arxiv.org/abs/astro-ph/0206508)
// Normalised such that U(k=0)=1
static double u_nfw_c(ccl_cosmology *cosmo, double rv, double c, double k, int *status) {

  double rs, ks;
  double f1, f2, f3, fc;

  // Special case to prevent numerical problems if k=0,
  // the result should be unity here because of the normalisation
  if (k==0.) {
    return 1.;
  }

  // The general k case
  else{

    // Scale radius for NFW (rs=rv/c)
    rs = rv/c;

    // Dimensionless wave-number variable
    ks = k*rs;

    // Various bits for summing together to get final result
    f1 = sin(ks)*(gsl_sf_Si(ks*(1.+c))-gsl_sf_Si(ks));
    f2 = cos(ks)*(gsl_sf_Ci(ks*(1.+c))-gsl_sf_Ci(ks));
    f3 = sin(c*ks)/(ks*(1.+c));
    fc = log(1.+c)-c/(1.+c);

    return (f1+f2-f3)/fc;

  }
}

/*----- ROUTINE: ccl_halo_concentration -----
INPUT: cosmology, a halo mass [Msun], scale factor, halo definition, concentration model label
TASK: Computes halo concentration; the ratio of virial raidus to scale radius for an NFW halo.
*/
double ccl_halo_concentration(ccl_cosmology *cosmo, double halomass,
                              double a, double odelta, int *status) {

  double gz, g0, nu, delta_c, a_form;
  double Mpiv, A, B, C;

  switch(cosmo->config.halo_concentration_method){

    // Bhattacharya et al. (2011; 1005.2239; Delta = 200rho_m; Table 2)
  case ccl_bhattacharya2011:

    if (odelta != 200.) {
      *status = CCL_ERROR_CONC_DV;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halomod.c: halo_concentration(): Bhattacharya (2011) concentration "
        "relation only valid for Delta_v = 200");
      return NAN;
    }

    gz = ccl_growth_factor(cosmo,a,status);
    g0 = ccl_growth_factor(cosmo,1.0,status);
    delta_c = 1.686;
    nu = delta_c/ccl_sigmaM(cosmo, log10(halomass), a, status);
    return 9.*pow(nu,-0.29)*pow(gz/g0,1.15);

    // Duffy et al. (2008; 0804.2486; Table 1)
  case ccl_duffy2008:

    Mpiv = 2e12/cosmo->params.h; // Pivot mass in Msun (note in the paper units are Msun/h)

    if (odelta == Dv_BryanNorman(cosmo, a, status)) {

      // Duffy et al. (2008) for virial density haloes (second section in Table 1)
      A = 7.85;
      B = -0.081;
      C = -0.71;
      return A*pow(halomass/Mpiv,B)*pow(a,-C);

    } else if (odelta == 200.) {

      // Duffy et al. (2008) for x200 mean-matter-density haloes (third section in Table 1)
      A = 10.14;
      B = -0.081;
      C = -1.01;
      return A*pow(halomass/Mpiv,B)*pow(a,-C);

    } else {

      *status = CCL_ERROR_CONC_DV;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_halomod.c: halo_concentration(): Duffy (2008) virial "
        "concentration only valid for virial Delta_v or 200\n");
      return NAN;

    }

    // Constant concentration (good for tests)
  case ccl_constant_concentration:

    return 4.;

    // Something went wrong
  default:

    *status = CCL_ERROR_HALOCONC;
    ccl_raise_gsl_warning(*status, "ccl_halomod.c: concentration-mass relation specified incorrectly");
    return NAN;

  }
}

// Fourier Transforms of halo profiles
static double window_function(ccl_cosmology *cosmo, double m, double k,
                              double a, double odelta, ccl_win_label label,
                              int *status) {
  double rho_matter, c, rv;

  switch(label){

  case ccl_nfw:

    // The mean background matter density in Msun/Mpc^3
    rho_matter = ccl_rho_x(cosmo, 1., ccl_species_m_label, 1, status);

    // The halo virial radius
    rv = r_delta(cosmo, m, a, odelta, status);

    // The halo concentration for this halo mass and at this scale factor
    c = ccl_halo_concentration(cosmo, m, a, odelta, status);

    // The function U is normalised to 1 for k<<1 so multiplying by M/rho turns units to overdensity
    return m*u_nfw_c(cosmo, rv, c, k, status)/rho_matter;

    // Something went wrong
  default:

    *status = CCL_ERROR_HALOWIN;
    ccl_raise_warning(*status, "ccl_halomod.c: Window function specified incorrectly");
    return NAN;

  }

}

// Parameters structure for the one-halo integrand
typedef struct{
  ccl_cosmology *cosmo;
  double k, a;
  int *status;
} Int_one_halo_Par;

// Integrand for the one-halo integral
static double one_halo_integrand(double log10mass, void *params){

  Int_one_halo_Par *p = (Int_one_halo_Par *)params;;
  double halomass = pow(10,log10mass);
  double odelta = Dv_BryanNorman(p->cosmo, p->a, p->status); // Virial density for haloes

  // The normalised Fourier Transform of a halo density profile
  double wk = window_function(p->cosmo,halomass, p->k, p->a, odelta, ccl_nfw, p->status);

  // Fairly sure that there should be no ln(10) factor should be here since the integration is being specified in log10 range
  double dn_dlogM = ccl_massfunc(p->cosmo, halomass, p->a, odelta, p->status);

  return dn_dlogM*pow(wk,2);
}

// The one-halo term integral using gsl
static double one_halo_integral(ccl_cosmology *cosmo, double k, double a, int *status){

  int one_halo_integral_status = 0, qagstatus;
  double result = 0, eresult;
  double log10mmin = log10(cosmo->gsl_params.HM_MMIN);
  double log10mmax = log10(cosmo->gsl_params.HM_MMAX);
  Int_one_halo_Par ipar;
  gsl_function F;
  gsl_integration_workspace *w = NULL;

  w = gsl_integration_workspace_alloc(cosmo->gsl_params.HM_LIMIT);
  if (w == NULL) {
    *status = CCL_ERROR_MEMORY;
  }

  if (*status == 0) {
    // Structure required for the gsl integration
    ipar.cosmo = cosmo;
    ipar.k = k;
    ipar.a = a;
    ipar.status = &one_halo_integral_status;
    F.function = &one_halo_integrand;
    F.params = &ipar;

    // Actually does the integration
    qagstatus = gsl_integration_qag(
      &F, log10mmin, log10mmax,
      cosmo->gsl_params.HM_EPSABS,
      cosmo->gsl_params.HM_EPSREL,
      cosmo->gsl_params.HM_LIMIT,
      cosmo->gsl_params.HM_INT_METHOD, w,
      &result, &eresult);

    // Check for errors
    if (qagstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(qagstatus, "ccl_halomod.c: one_halo_integral():");
      *status = CCL_ERROR_ONE_HALO_INT;
      ccl_cosmology_set_status_message(cosmo, "ccl_halomod.c: one_halo_integral(): Integration failure\n");
      result = NAN;
    }
  }

  // Clean up
  gsl_integration_workspace_free(w);

  return result;
}

// Parameters structure for the two-halo integrand
typedef struct{
  ccl_cosmology *cosmo;
  double k, a;
  int *status;
} Int_two_halo_Par;

// Integrand for the two-halo integral
static double two_halo_integrand(double log10mass, void *params){

  Int_two_halo_Par *p = (Int_two_halo_Par *)params;
  double halomass = pow(10,log10mass);
  double odelta = Dv_BryanNorman(p->cosmo, p->a, p->status); // Virial density for haloes

  // The normalised Fourier Transform of a halo density profile
  double wk = window_function(p->cosmo, halomass, p->k, p->a, odelta, ccl_nfw, p->status);

  // Fairly sure that there should be no ln(10) factor should be here since the integration is being specified in log10 range
  double dn_dlogM = ccl_massfunc(p->cosmo, halomass, p->a, odelta, p->status);

  // Halo bias
  double b = ccl_halo_bias(p->cosmo, halomass, p->a, odelta, p->status);

  return b*dn_dlogM*wk;
}

// The two-halo term integral using gsl
static double two_halo_integral(ccl_cosmology *cosmo, double k, double a, int *status){

  int two_halo_integral_status = 0, qagstatus;
  double result = 0, eresult;
  double log10mmin = log10(cosmo->gsl_params.HM_MMIN);
  double log10mmax = log10(cosmo->gsl_params.HM_MMAX);
  Int_two_halo_Par ipar;
  gsl_function F;
  gsl_integration_workspace *w = NULL;

  w = gsl_integration_workspace_alloc(cosmo->gsl_params.HM_LIMIT);
  if (w == NULL) {
    *status = CCL_ERROR_MEMORY;
  }

  if (*status == 0) {
    // Structure required for the gsl integration
    ipar.cosmo = cosmo;
    ipar.k = k;
    ipar.a = a;
    ipar.status = &two_halo_integral_status;
    F.function = &two_halo_integrand;
    F.params = &ipar;

    // Actually does the integration
    qagstatus = gsl_integration_qag(
      &F, log10mmin, log10mmax,
      cosmo->gsl_params.HM_EPSABS,
      cosmo->gsl_params.HM_EPSREL,
      cosmo->gsl_params.HM_LIMIT,
      cosmo->gsl_params.HM_INT_METHOD, w,
      &result, &eresult);

    // Check for errors
    if (qagstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(qagstatus, "ccl_halomod.c: two_halo_integral():");
      *status = CCL_ERROR_TWO_HALO_INT;
      ccl_cosmology_set_status_message(cosmo, "ccl_halomod.c: two_halo_integral(): Integration failure\n");
      result = NAN;
    }
  }

  // Clean up
  gsl_integration_workspace_free(w);

  return result;
}

/*----- ROUTINE: ccl_twohalo_matter_power -----
INPUT: cosmology, wavenumber [Mpc^-1], scale factor
TASK: Computes the two-halo power spectrum term in the halo model assuming NFW haloes
*/
double ccl_twohalo_matter_power(ccl_cosmology *cosmo, double k, double a, int *status){

  // Get the integral
  double I2h = two_halo_integral(cosmo, k, a, status);

  // The addative correction is the missing part of the integral below the lower-mass limit
  double A = 1.-two_halo_integral(cosmo, 0., a, status);

  // Virial overdensity for haloes
  double odelta = Dv_BryanNorman(cosmo, a, status);

  // ...multiplied by the ratio of window functions
  double W1 = window_function(cosmo, cosmo->gsl_params.HM_MMIN, k,  a, odelta, ccl_nfw, status);
  double W2 = window_function(cosmo, cosmo->gsl_params.HM_MMIN, 0., a, odelta, ccl_nfw, status);
  A = A*W1/W2;

  // Add the additive correction to the calculated integral
  I2h = I2h+A;

  return ccl_linear_matter_power(cosmo, k, a, status)*I2h*I2h;

}

/*----- ROUTINE: ccl_onehalo_matter_power -----
INPUT: cosmology, wavenumber [Mpc^-1], scale factor
TASK: Computes the one-halo power spectrum term in the halo model assuming NFW haloes
*/
double ccl_onehalo_matter_power(ccl_cosmology *cosmo, double k, double a, int *status){

  return one_halo_integral(cosmo, k, a, status);

}

/*----- ROUTINE: ccl_onehalo_matter_power -----
INPUT: cosmology, wavenumber [Mpc^-1], scale factor
TASK: Computes the halo model power spectrum by summing the two- and one-halo terms
*/
double ccl_halomodel_matter_power(ccl_cosmology *cosmo, double k, double a, int *status){

  // This matter power spectrum method doesn't work if mu / Sigma parameterisation of modified gravity
  // Is turned on:
  if (fabs(cosmo->params.mu_0)>1e-14 || fabs(cosmo->params.sigma_0)>1e-14){
	  *status = CCL_ERROR_NOT_IMPLEMENTED;
	  strcpy(cosmo->status_message,
           "ccl_halofmod.c: ccl_halomodel_matter_power(): The halo model power spectrum "
           "is not implemented the mu / Sigma modified gravity parameterisation.\n");
	  return NAN;
  }

  // Standard sum of two- and one-halo terms
  return ccl_twohalo_matter_power(cosmo, k, a, status)+ccl_onehalo_matter_power(cosmo, k, a, status);

}

max_stars_repo_name: benediktdiemer/CCL

max_issues_repo_name: benediktdiemer/CCL

max_forks_repo_name: benediktdiemer/CCL


text:
/* specfunc/poch.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 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "gsl_sf_exp.h"
#include "gsl_sf_log.h"
#include "gsl_sf_pow_int.h"
#include "gsl_sf_psi.h"
#include "gsl_sf_gamma.h"

#include "error.h"

static const double bern[21] = {
   0.0   /* no element 0 */,
  +0.833333333333333333333333333333333e-01,
  -0.138888888888888888888888888888888e-02,
  +0.330687830687830687830687830687830e-04,
  -0.826719576719576719576719576719576e-06,
  +0.208767569878680989792100903212014e-07,
  -0.528419013868749318484768220217955e-09,
  +0.133825365306846788328269809751291e-10,
  -0.338968029632258286683019539124944e-12,
  +0.858606205627784456413590545042562e-14,
  -0.217486869855806187304151642386591e-15,
  +0.550900282836022951520265260890225e-17,
  -0.139544646858125233407076862640635e-18,
  +0.353470703962946747169322997780379e-20,
  -0.895351742703754685040261131811274e-22,
  +0.226795245233768306031095073886816e-23,
  -0.574472439520264523834847971943400e-24,
  +0.145517247561486490186626486727132e-26,
  -0.368599494066531017818178247990866e-28,
  +0.933673425709504467203255515278562e-30,
  -0.236502241570062993455963519636983e-31
};

/* ((a)_x - 1)/x in the "small x" region where
 * cancellation must be controlled.
 *
 * Based on SLATEC DPOCH1().
 */
/*
C When ABS(X) is so small that substantial cancellation will occur if
C the straightforward formula is used, we use an expansion due
C to Fields and discussed by Y. L. Luke, The Special Functions and Their
C Approximations, Vol. 1, Academic Press, 1969, page 34.
C
C The ratio POCH(A,X) = GAMMA(A+X)/GAMMA(A) is written by Luke as
C        (A+(X-1)/2)**X * polynomial in (A+(X-1)/2)**(-2) .
C In order to maintain significance in POCH1, we write for positive a
C        (A+(X-1)/2)**X = EXP(X*LOG(A+(X-1)/2)) = EXP(Q)
C                       = 1.0 + Q*EXPREL(Q) .
C Likewise the polynomial is written
C        POLY = 1.0 + X*POLY1(A,X) .
C Thus,
C        POCH1(A,X) = (POCH(A,X) - 1) / X
C                   = EXPREL(Q)*(Q/X + Q*POLY1(A,X)) + POLY1(A,X)
C
*/
static
int
pochrel_smallx(const double a, const double x, gsl_sf_result * result)
{
  /*
   SQTBIG = 1.0D0/SQRT(24.0D0*D1MACH(1))
   ALNEPS = LOG(D1MACH(3))
   */
  const double SQTBIG = 1.0/(2.0*M_SQRT2*M_SQRT3*GSL_SQRT_DBL_MIN);
  const double ALNEPS = GSL_LOG_DBL_EPSILON - M_LN2;

  if(x == 0.0) {
    return gsl_sf_psi_e(a, result);
  }
  else {
    const double bp   = (  (a < -0.5) ? 1.0-a-x : a );
    const int	 incr = ( (bp < 10.0) ? 11.0-bp : 0 );
    const double b    = bp + incr;
    double dpoch1;
    gsl_sf_result dexprl;
    int stat_dexprl;
    int i;

    double var    = b + 0.5*(x-1.0);
    double alnvar = log(var);
    double q = x*alnvar;

    double poly1 = 0.0;

    if(var < SQTBIG) {
      const int nterms = (int)(-0.5*ALNEPS/alnvar + 1.0);
      const double var2 = (1.0/var)/var;
      const double rho  = 0.5 * (x + 1.0);
      double term = var2;
      double gbern[24];
      int k, j;

      gbern[1] = 1.0;
      gbern[2] = -rho/12.0;
      poly1 = gbern[2] * term;

      if(nterms > 20) {
        /* NTERMS IS TOO BIG, MAYBE D1MACH(3) IS BAD */
        /* nterms = 20; */
	result->val = 0.0;
	result->err = 0.0;
	GSL_ERROR ("error", GSL_ESANITY);
      }

      for(k=2; k<=nterms; k++) {
        double gbk = 0.0;
        for(j=1; j<=k; j++) {
          gbk += bern[k-j+1]*gbern[j];
        }
        gbern[k+1] = -rho*gbk/k;

        term  *= (2*k-2-x)*(2*k-1-x)*var2;
        poly1 += gbern[k+1]*term;
      }
    }

    stat_dexprl = gsl_sf_expm1_e(q, &dexprl);
    if(stat_dexprl != GSL_SUCCESS) {
      result->val = 0.0;
      result->err = 0.0;
      return stat_dexprl;
    }
    dexprl.val = dexprl.val/q;
    poly1 *= (x - 1.0);
    dpoch1 = dexprl.val * (alnvar + q * poly1) + poly1;

    for(i=incr-1; i >= 0; i--) {
      /*
       C WE HAVE DPOCH1(B,X), BUT BP IS SMALL, SO WE USE BACKWARDS RECURSION
       C TO OBTAIN DPOCH1(BP,X).
       */
      double binv = 1.0/(bp+i);
      dpoch1 = (dpoch1 - binv) / (1.0 + x*binv);
    }

    if(bp == a) {
      result->val = dpoch1;
      result->err = 2.0 * GSL_DBL_EPSILON * (fabs(incr) + 1.0) * fabs(result->val);
      return GSL_SUCCESS;
    }
    else {
      /*
       C WE HAVE DPOCH1(BP,X), BUT A IS LT -0.5.  WE THEREFORE USE A
       C REFLECTION FORMULA TO OBTAIN DPOCH1(A,X).
       */
      double sinpxx = sin(M_PI*x)/x;
      double sinpx2 = sin(0.5*M_PI*x);
      double t1 = sinpxx/tan(M_PI*b);
      double t2 = 2.0*sinpx2*(sinpx2/x);
      double trig  = t1 - t2;
      result->val  = dpoch1 * (1.0 + x*trig) + trig;
      result->err  = (fabs(dpoch1*x) + 1.0) * GSL_DBL_EPSILON * (fabs(t1) + fabs(t2));
      result->err += 2.0 * GSL_DBL_EPSILON * (fabs(incr) + 1.0) * fabs(result->val);
      return GSL_SUCCESS;
    }
  }
}

/* Assumes a>0 and a+x>0.
 */
static
int
lnpoch_pos(const double a, const double x, gsl_sf_result * result)
{
  double absx = fabs(x);

  if(absx > 0.1*a || absx*log(GSL_MAX_DBL(a,2.0)) > 0.1) {
    if(a < GSL_SF_GAMMA_XMAX && a+x < GSL_SF_GAMMA_XMAX) {
      /* If we can do it by calculating the gamma functions
       * directly, then that will be more accurate than
       * doing the subtraction of the logs.
       */
      gsl_sf_result g1;
      gsl_sf_result g2;
      gsl_sf_gammainv_e(a,   &g1);
      gsl_sf_gammainv_e(a+x, &g2);
      result->val  = -log(g2.val/g1.val);
      result->err  = g1.err/fabs(g1.val) + g2.err/fabs(g2.val);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      return GSL_SUCCESS;
    }
    else {
      /* Otherwise we must do the subtraction.
       */
      gsl_sf_result lg1;
      gsl_sf_result lg2;
      int stat_1 = gsl_sf_lngamma_e(a,   &lg1);
      int stat_2 = gsl_sf_lngamma_e(a+x, &lg2);
      result->val  = lg2.val - lg1.val;
      result->err  = lg2.err + lg1.err;
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      return GSL_ERROR_SELECT_2(stat_1, stat_2);
    }
  }
  else if(absx < 0.1*a && a > 15.0) {
    /* Be careful about the implied subtraction.
     * Note that both a+x and and a must be
     * large here since a is not small
     * and x is not relatively large.
     * So we calculate using Stirling for Log[Gamma(z)].
     *
     *   Log[Gamma(a+x)/Gamma(a)] = x(Log[a]-1) + (x+a-1/2)Log[1+x/a]
     *                              + (1/(1+eps)   - 1) / (12 a)
     *                              - (1/(1+eps)^3 - 1) / (360 a^3)
     *                              + (1/(1+eps)^5 - 1) / (1260 a^5)
     *                              - (1/(1+eps)^7 - 1) / (1680 a^7)
     *                              + ...
     */
    const double eps = x/a;
    const double den = 1.0 + eps;
    const double d3 = den*den*den;
    const double d5 = d3*den*den;
    const double d7 = d5*den*den;
    const double c1 = -eps/den;
    const double c3 = -eps*(3.0+eps*(3.0+eps))/d3;
    const double c5 = -eps*(5.0+eps*(10.0+eps*(10.0+eps*(5.0+eps))))/d5;
    const double c7 = -eps*(7.0+eps*(21.0+eps*(35.0+eps*(35.0+eps*(21.0+eps*(7.0+eps))))))/d7;
    const double p8 = gsl_sf_pow_int(1.0+eps,8);
    const double c8 = 1.0/p8             - 1.0;  /* these need not   */
    const double c9 = 1.0/(p8*(1.0+eps)) - 1.0;  /* be very accurate */
    const double a4 = a*a*a*a;
    const double a6 = a4*a*a;
    const double ser_1 = c1 + c3/(30.0*a*a) + c5/(105.0*a4) + c7/(140.0*a6);
    const double ser_2 = c8/(99.0*a6*a*a) - 691.0/360360.0 * c9/(a6*a4);
    const double ser = (ser_1 + ser_2)/ (12.0*a);

    double term1 = x * log(a/M_E);
    double term2;
    gsl_sf_result ln_1peps;
    gsl_sf_log_1plusx_e(eps, &ln_1peps);  /* log(1 + x/a) */
    term2 = (x + a - 0.5) * ln_1peps.val;

    result->val  = term1 + term2 + ser;
    result->err  = GSL_DBL_EPSILON*fabs(term1);
    result->err += fabs((x + a - 0.5)*ln_1peps.err);
    result->err += fabs(ln_1peps.val) * GSL_DBL_EPSILON * (fabs(x) + fabs(a) + 0.5);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result poch_rel;
    int stat_p = pochrel_smallx(a, x, &poch_rel);
    double eps = x*poch_rel.val;
    int stat_e = gsl_sf_log_1plusx_e(eps, result);
    result->err  = 2.0 * fabs(x * poch_rel.err / (1.0 + eps));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_ERROR_SELECT_2(stat_e, stat_p);
  }
}

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int
gsl_sf_lnpoch_e(const double a, const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(a <= 0.0 || a+x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    return lnpoch_pos(a, x, result);
  }
}

int
gsl_sf_lnpoch_sgn_e(const double a, const double x,
                       gsl_sf_result * result, double * sgn)
{
  if(a == 0.0 || a+x == 0.0) {
    *sgn = 0.0;
    DOMAIN_ERROR(result);
  }
  else if(x == 0.0) {
    *sgn = 1.0;
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(a > 0.0 && a+x > 0.0) {
    *sgn = 1.0;
    return lnpoch_pos(a, x, result);
  }
  else if(a < 0.0 && a+x < 0.0) {
    /* Reduce to positive case using reflection.
     */
    double sin_1 = sin(M_PI * (1.0 - a));
    double sin_2 = sin(M_PI * (1.0 - a - x));
    if(sin_1 == 0.0 || sin_2 == 0.0) {
      *sgn = 0.0;
      DOMAIN_ERROR(result);
    }
    else {
      gsl_sf_result lnp_pos;
      int stat_pp   = lnpoch_pos(1.0-a, -x, &lnp_pos);
      double lnterm = log(fabs(sin_1/sin_2));
      result->val  = lnterm - lnp_pos.val;
      result->err  = lnp_pos.err;
      result->err += 2.0 * GSL_DBL_EPSILON * (fabs(1.0-a) + fabs(1.0-a-x)) * fabs(lnterm);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      *sgn = GSL_SIGN(sin_1*sin_2);
      return stat_pp;
    }
  }
  else {
    /* Evaluate gamma ratio directly.
     */
    gsl_sf_result lg_apn;
    gsl_sf_result lg_a;
    double s_apn, s_a;
    int stat_apn = gsl_sf_lngamma_sgn_e(a+x, &lg_apn, &s_apn);
    int stat_a   = gsl_sf_lngamma_sgn_e(a,   &lg_a,   &s_a);
    if(stat_apn == GSL_SUCCESS && stat_a == GSL_SUCCESS) {
      result->val  = lg_apn.val - lg_a.val;
      result->err  = lg_apn.err + lg_a.err;
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      *sgn = s_a * s_apn;
      return GSL_SUCCESS;
    }
    else if(stat_apn == GSL_EDOM || stat_a == GSL_EDOM){
      *sgn = 0.0;
      DOMAIN_ERROR(result);
    }
    else {
      result->val = 0.0;
      result->err = 0.0;
      *sgn = 0.0;
      return GSL_FAILURE;
    }
  }
}

int
gsl_sf_poch_e(const double a, const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result lnpoch;
    double sgn;
    int stat_lnpoch = gsl_sf_lnpoch_sgn_e(a, x, &lnpoch, &sgn);
    int stat_exp    = gsl_sf_exp_err_e(lnpoch.val, lnpoch.err, result);
    result->val *= sgn;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_ERROR_SELECT_2(stat_exp, stat_lnpoch);
  }
}

int
gsl_sf_pochrel_e(const double a, const double x, gsl_sf_result * result)
{
  const double absx = fabs(x);
  const double absa = fabs(a);

  /* CHECK_POINTER(result) */

  if(absx > 0.1*absa || absx*log(GSL_MAX(absa,2.0)) > 0.1) {
    gsl_sf_result lnpoch;
    double sgn;
    int stat_poch = gsl_sf_lnpoch_sgn_e(a, x, &lnpoch, &sgn);
    if(lnpoch.val > GSL_LOG_DBL_MAX) {
      OVERFLOW_ERROR(result);
    }
    else {
      const double el = exp(lnpoch.val);
      result->val  = (sgn*el - 1.0)/x;
      result->err  = fabs(result->val) * (lnpoch.err + 2.0 * GSL_DBL_EPSILON);
      result->err += 2.0 * GSL_DBL_EPSILON * (fabs(sgn*el) + 1.0) / fabs(x);
      return stat_poch;
    }
  }
  else {
    return pochrel_smallx(a, x, result);
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_lnpoch(const double a, const double x)
{
  EVAL_RESULT(gsl_sf_lnpoch_e(a, x, &result));
}

double gsl_sf_poch(const double a, const double x)
{
  EVAL_RESULT(gsl_sf_poch_e(a, x, &result));
}

double gsl_sf_pochrel(const double a, const double x)
{
  EVAL_RESULT(gsl_sf_pochrel_e(a, x, &result));
}


text:
#pragma once
#include <vector>
#include <array>
#include <gsl/gsl>

#pragma warning(push)
#include <CppCoreCheck/Warnings.h>
#pragma warning(disable: ALL_CPPCORECHECK_WARNINGS)
#include <glm/glm.hpp>
#define GLM_ENABLE_EXPERIMENTAL
#include <glm/gtx/hash.hpp>
#pragma warning(pop)

#include <vulkan/vulkan.h>

#include "../utils/Utils.h"
#include "./RenderUtils.h"
#include "../Configuration.h"

#define VMA_USE_ALLOCATOR
#ifdef VMA_USE_ALLOCATOR
#pragma warning(push)
#include <CppCoreCheck/Warnings.h>
#pragma warning(disable: ALL_CPPCORECHECK_WARNINGS)
#include "vk_mem_alloc.h"
#pragma warning(pop)

/**
Our main unsigned in type to accomodate
vulkan's needs
*/
using uint = uint32_t;

/**
Wrapps a Vulkan Buffer with allocation information tied to it
*/
struct AllocatedBuffer {
	VkBuffer buffer{};
	VmaAllocation allocation{};
	VmaAllocationInfo allocation_info{};
};

/**
Wrapps a Vulkan Image with allocation information tied to it
*/
struct AllocatedImage {
	VkImage image{};
	VmaAllocation allocation{};
	VmaAllocationInfo allocation_info{};
};

#else

#pragma message ( "This program is meant to currently use VMA allocation" )
#error Current code needs to use VMA allocation

struct AllocatedBuffer {
	VkBuffer buffer{};
	VkDeviceMemory memory{};
};

struct AllocatedImage {
	VkImage image{};
	VkDeviceMemory memory{};
};

#endif

/**
Wraps a Vulkan Command Buffer with relevant
type and state information tied to it
*/
struct WrappedCommandBuffer {
	VkCommandBuffer buffer{};
	CommandType type{};
	bool recording{ false };
};

/**
Wraps a Vulkan Render Target with all the relevant information to
render to it like the image, memory, view, width, heigth and if it has
been initializated.

We are not using an "AllocatedImage" because allocation of a render target
is done with a custom method because of the special requirements of it being
a render target.
*/
struct WrappedRenderTarget {
	VkImage image{};
	VkDeviceMemory memory{};
	VkImageView view{};
	uint width{};
	uint heigth{};
	bool init{ false };
};

/**
Struct that holds dynamic configuration parameters of the renderer
*/
struct RenderConfiguration {
	short multisampling_samples{ config::initial_multisampling_samples };
};

struct Vertex {
	glm::vec3 pos{};
	glm::vec3 color{};
	glm::vec2 tex_coord{};

	auto static getBindingDescription() noexcept ->VkVertexInputBindingDescription {

		auto binding_description = VkVertexInputBindingDescription{};

		binding_description.binding = 0;
		binding_description.stride = sizeof(Vertex);
		binding_description.inputRate = VK_VERTEX_INPUT_RATE_VERTEX;

		return binding_description;
	}

	auto static getAttributeDescriptions() noexcept->std::array<VkVertexInputAttributeDescription, 3> {

		auto attribute_descriptions = std::array<VkVertexInputAttributeDescription, 3>{};

		attribute_descriptions[0].binding = 0;
		attribute_descriptions[0].location = 0;
		attribute_descriptions[0].format = VK_FORMAT_R32G32B32_SFLOAT;
		attribute_descriptions[0].offset = offsetof(Vertex, pos);

		attribute_descriptions[1].binding = 0;
		attribute_descriptions[1].location = 1;
		attribute_descriptions[1].format = VK_FORMAT_R32G32B32_SFLOAT;
		attribute_descriptions[1].offset = offsetof(Vertex, color);

		attribute_descriptions[2].binding = 0;
		attribute_descriptions[2].location = 2;
		attribute_descriptions[2].format = VK_FORMAT_R32G32_SFLOAT;
		attribute_descriptions[2].offset = offsetof(Vertex, tex_coord);

		return attribute_descriptions;
	}

	auto operator==(const Vertex& other) const ->bool {
		return	pos == other.pos &&
				color == other.color &&
				tex_coord == other.tex_coord;
	}
};

namespace std {
	template<> struct hash<Vertex> {
		size_t operator()(Vertex const& vertex) const {
			return(
				(hash<glm::vec3>()(vertex.pos) ^
				(hash<glm::vec3>()(vertex.color) << 1)) >> 1) ^
				(hash<glm::vec2>()(vertex.tex_coord) << 1);
		}
	};
}

/*
MVP model.

http://www.opengl-tutorial.org/es/beginners-tutorials/tutorial-3-matrices/
https://solarianprogrammer.com/2013/05/22/opengl-101-matrices-projection-view-model/
*/
struct UniformBufferObject {
	glm::mat4 model{};
	glm::mat4 view{};
	glm::mat4 proj{};
};

struct SimpleObjScene {
	std::vector<Vertex> vertices;
	std::vector<uint32_t> indices;
	AllocatedImage m_texture_image{};
	VkImageView m_texture_image_view{};
};

#if 0
#pragma warning(push)
#include <CppCoreCheck/Warnings.h>
#pragma warning(disable: 26426)
const auto vertices = std::vector<Vertex>{
	{ { -0.5f, -0.5f, 0.0f },{ 1.0f, 0.0f, 0.0f }, { 1.0f, 0.0f} },
	{ {  0.5f, -0.5f, 0.0f },{ 0.0f, 1.0f, 0.0f }, { 0.0f, 0.0f } },
	{ {  0.5f,  0.5f, 0.0f },{ 0.0f, 0.0f, 1.0f }, { 0.0f, 1.0f } },
	{ { -0.5f,  0.5f, 0.0f },{ 1.0f, 1.0f, 1.0f }, { 1.0f, 1.0f } },

	{ { -0.5f, -0.5f, -0.5f },{ 1.0f, 0.0f, 0.0f },{ 1.0f, 0.0f } },
	{ {  0.5f, -0.5f, -0.5f },{ 0.0f, 1.0f, 0.0f },{ 0.0f, 0.0f } },
	{ {  0.5f,  0.5f, -0.5f },{ 0.0f, 0.0f, 1.0f },{ 0.0f, 1.0f } },
	{ { -0.5f,  0.5f, -0.5f },{ 1.0f, 1.0f, 1.0f },{ 1.0f, 1.0f } },
};

const auto indices = std::vector<uint16_t>{
	0, 1, 2, 2, 3, 0,
	4, 5, 6, 6, 7, 4
};
#pragma warning(pop)
#endif


text:
/* monte/gsl_monte_miser.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth
 *
 * 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.
 */

/* Author: MJB */

#ifndef __GSL_MONTE_MISER_H__
#define __GSL_MONTE_MISER_H__

#include <gsl/gsl_rng.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct {
  size_t min_calls;
  size_t min_calls_per_bisection;
  double dither;
  double estimate_frac;
  double alpha;
  size_t dim;
  int estimate_style;
  int depth;
  int verbose;
  double * x;
  double * xmid;
  double * sigma_l;
  double * sigma_r;
  double * fmax_l;
  double * fmax_r;
  double * fmin_l;
  double * fmin_r;
  double * fsum_l;
  double * fsum_r;
  double * fsum2_l;
  double * fsum2_r;
  size_t * hits_l;
  size_t * hits_r;
} gsl_monte_miser_state;

int gsl_monte_miser_integrate(gsl_monte_function * f,
                              const double xl[], const double xh[],
                              size_t dim, size_t calls,
                              gsl_rng *r,
                              gsl_monte_miser_state* state,
                              double *result, double *abserr);

gsl_monte_miser_state* gsl_monte_miser_alloc(size_t dim);

int gsl_monte_miser_init(gsl_monte_miser_state* state);

void gsl_monte_miser_free(gsl_monte_miser_state* state);

__END_DECLS

#endif /* __GSL_MONTE_MISER_H__ */

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/**
 * @file bblas_cutil.c
 *
 * @brief BBLAS testing utilities for float _Complex routines.
 *
 * BBLAS is a software package provided by Univ. of Manchester,
 * Univ. of Tennessee.
 *
 * @version 1.0.0
 * @author  Samuel  D. Relton
 * @author  Pedro   V. Lara
 * @author  Mawussi Zounon
 * @date    2016-02-20
 *
 * Contains routines used in the testing to modify randomly generated matrices
 * and compute the average error over an entire batch etc.
 *
 **/

#ifndef DOXYGEN_SHOULD_SKIP_THIS
/**
 * Code generation
 * @generated from bblas_zutil.c normal z -> c, Mon Jun  6 09:44:14 2016
 **/
#endif

#include "bblas_common.h"
#if defined(BBLAS_WITH_MKL)
#include <mkl_lapacke.h>
#else
#include <lapacke.h>
#endif

/** Include complex functions since using float complex precision **/
#define COMPLEX

/** Quick access to the matrix elements **/
#define A(i,j)  A[i + j*lda]

/**
 * Make a matrix symmetric/Hermitian.  Makes diagonal real.
 * Sets Aji = conj( Aij ) for j < i, that is, copy & conjugate
 * lower triangle to upper triangle.
 **/

void bblas_cmake_hermitian(int lda, int N, BBLAS_Complex32_t* A)
{
    int i, j;
    for( i=0; i < N; ++i ) {
        A(i,i) = creal( A(i,i) );
        for( j=0; j < i; ++j ) {
            A(j,i) = conj( A(i,j) );
        }
    }
}

#ifdef COMPLEX

/**
 * Make a matrix complex-symmetric
 * Does NOT make diagonal real.
 * Sets Aji = Aij for j < i, that is,
 * copy lower triangle to upper triangle.
 **/

void bblas_cmake_symmetric(int lda, int N, BBLAS_Complex32_t* A)
{
    int i, j;
    for( i=0; i < N; ++i ) {
        for( j=0; j < i; ++j ) {
            A(j,i) =  A(i,j);
        }
    }
}
#endif

/**
 * irandRange generates a random value (int)
 * in the range min_n and max_n
 **/

int irandRange(int min_n, int max_n)
{
    return rand() % (max_n - min_n + 1) + min_n;
}

/**
 * bblas_crandRange generates a random value (BBLAS_Complex32_t)
 * in the range [0,max_n]: TODO replace by generic
 */

BBLAS_Complex32_t bblas_crandRange( int max_n )
{
    return ( BBLAS_Complex32_t )rand()/( BBLAS_Complex32_t )( RAND_MAX/max_n );
}

/**
 * Computes statistics of the relative errors to summarise them for the user.
 **/

void bblas_cstatistic(bblas_ctest_t *test)
{

    enum BBLAS_ROUTINE routine  = test->routine;

    /*Compute avg(M), avg(N), avg(K) */
    if(test->batch_opts == BBLAS_VARIABLE)
    {
	if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	    (routine == BBLAS_HEMM) || (routine == BBLAS_TRMM) ||
	    (routine == BBLAS_TRSM))
	{
	    test->avgM = bblas_avgarrayI(test->M, test->batch_count);
	}

	if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYRK)  ||
	    (routine == BBLAS_HERK) || (routine == BBLAS_SYR2K) ||
	    (routine == BBLAS_HER2K))
	{
	    test->avgK = bblas_avgarrayI(test->K, test->batch_count);
	}

	test->avgN = bblas_avgarrayI(test->N, test->batch_count);

    } else if (test->batch_opts == BBLAS_FIXED)
    {

	if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	    (routine == BBLAS_HEMM) || (routine == BBLAS_TRMM) ||
	    (routine == BBLAS_TRSM))
	{
	    test->avgM = test->M[0];
	}

	if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYRK)  ||
	    (routine == BBLAS_HERK) || (routine == BBLAS_SYR2K) ||
	    (routine == BBLAS_HER2K))
	{
	    test->avgK = test->K[0];
	}
	test->avgN = test->N[0];

    } else
    {
	bblas_error("testing_cgemm_batch.c", "wrong batch_opts value");
    }

    /*Statistics on the error */
    switch(test->target)
    {
    case BBLAS_MKL:
      test->mkl_min_error = bblas_cminarrayD(test->mkl_error, test->batch_count);
      test->mkl_avg_error = bblas_cavgarrayD(test->mkl_error, test->batch_count);
      test->mkl_max_error = bblas_cmaxarrayD(test->mkl_error, test->batch_count);
      test->mkl_std_error = bblas_cstdarrayD(test->mkl_error, test->batch_count);
      break;

    case BBLAS_CUBLAS:
    case BBLAS_MAGMA:
      test->device_min_error = bblas_cminarrayD(test->device_error, test->batch_count);
      test->device_avg_error = bblas_cavgarrayD(test->device_error, test->batch_count);
      test->device_max_error = bblas_cmaxarrayD(test->device_error, test->batch_count);
      test->device_std_error = bblas_cstdarrayD(test->device_error, test->batch_count);
      break;

    case BBLAS_OTHER:
      test->other_min_error = bblas_cminarrayD(test->other_error, test->batch_count);
      test->other_avg_error = bblas_cavgarrayD(test->other_error, test->batch_count);
      test->other_max_error = bblas_cmaxarrayD(test->other_error, test->batch_count);
      test->other_std_error = bblas_cstdarrayD(test->other_error, test->batch_count);
      break;

    default:
      printf("In bblas_cstatistic(): Target no defined\n");
      exit(EXIT_FAILURE);
    }
}

/**
 * Print a matrix.
 **/
void bblas_cprintmatrix(BBLAS_Complex32_t *matrix, int row, int col)
{
    /*Local variables */
    int i, j;

    for (i=0; i < row; i++)
    {
	printf("\n\n");
	for (j=0; j < col; j++)
	{
#ifdef COMPLEX
	    printf("%1.2f + %1.2f\t", creal(matrix[i*col+j]), cimag(matrix[i*col+j]));
#else
	    printf("%1.2f",matrix[i*col+j]);
#endif
	}
    }
    printf("\n");
}

/**
 * Decide whether a batch is fixed or variable.
 **/

char* bblas_getoption(enum BBLAS_OPTS opts)
{
   /*Local variable */
    char funcname[] = "bblas_getoption";

    switch(opts)
    {
	case BBLAS_VARIABLE:
	    return "BATCH OPTION: VARIABLE";
	    break;

	case BBLAS_FIXED:
	    return "BATCH OPTION: FIXED";
	    break;

	default:
	    printf("ERROR in %s, undefined bblas routine name\n",funcname);
	    exit(EXIT_FAILURE);
    }

}

/**
 * Get the name of the current test routine.
 **/

char* bblas_getroutine(enum BBLAS_ROUTINE routine)
{
    /*Local variable */
    char funcname[] = "bblas_getroutine";

    switch(routine)
    {
	case BBLAS_GEMM:
	    return "CGEMM";
	    break;

	case BBLAS_HEMM:
	    return "CHEMM";
	    break;

	case BBLAS_HER2K:
	    return "CHER2K";
	    break;

	case BBLAS_HERK:
	    return "CHERK";
	    break;

	case BBLAS_SYMM:
	    return "CSYMM";
	    break;

	case BBLAS_SYR2K:
	    return "CSYR2K";
	    break;

	case BBLAS_SYRK:
	    return "CSYRK";
	    break;

	case BBLAS_TRMM:
	    return "CTRMM";
	    break;

	case BBLAS_TRSM:
	    return "CTRSM";
	    break;

	default:
	    printf("ERROR in %s, undefined bblas routine name\n",funcname);
	    exit(EXIT_FAILURE);
    }
}

/**
 * Computes the maximum value of an array of real floats.
 **/

float bblas_cmaxarrayD(float *myArray, int size)
{
    int iter;
    float maxValue = myArray[0];

    for (iter = 0; iter < size; ++iter)
    {
        if ( myArray[iter] > maxValue )
	{
            maxValue = myArray[iter];
        }
    }
    return maxValue;
}

/**
 * Computes the minimum value of an array of real floats.
 **/

float bblas_cminarrayD(float *myArray, int size)
{
    int iter;
    float minValue = myArray[0];

    for (iter = 0; iter < size; ++iter)
    {
        if ( myArray[iter] < minValue )
	{
            minValue = myArray[iter];
        }
    }
    return minValue;
}

/**
 * Computes the mean value of an array of real floats.
 **/

float bblas_cavgarrayD(float *myArray, int size)
{
    int iter;
    float avg = 0.;

    for (iter = 0; iter < size; ++iter)
    {
	avg += myArray[iter];
    }
    return avg/size;
}

/**
 * Computes the standard deviation of an array of real floats.
 **/

float bblas_cstdarrayD(float *myArray, int size)
{
    int iter;
    float avg, sd=0.;

    avg =  bblas_cavgarrayD(myArray, size);

    for (iter = 0; iter < size; ++iter)
    {
	sd += (myArray[iter] -avg)*(myArray[iter] -avg);
    }
    return sd/size;
}

/**
 * Computes the minimum value of an array of integers.
 **/

int bblas_minarrayI(int *myArray, int size)
{
    int iter;
    int minValue = myArray[0];

    for (iter = 0; iter < size; ++iter)
    {
        if ( myArray[iter] < minValue )
	{
            minValue = myArray[iter];
        }
    }
    return minValue;
}

/**
 * Computes the mean of an array of integers.
 **/

int bblas_avgarrayI(int *myArray, int size)
{
    int iter;
    int avg = 0;

    for (iter = 0; iter < size; ++iter)
    {
	avg += myArray[iter];
    }
    return avg/size;
}

/**
 * Transform BBLAS enum values for <tt>trans</tt>, <tt>uplo</tt> etc. to human-readable strings.
 **/

char* bblas_op2char(unsigned int op)
{
    char *opname = (char*)malloc(30*sizeof(char));
    switch(op)
    {
        case BblasNoTrans:
            strcpy(opname,"CblasNoTrans");
            break;

        case BblasTrans:
	    strcpy(opname,"CblasTrans");
            break;

        case BblasConjTrans:
            strcpy(opname,"CblasConjTrans");
            break;

        case BblasLower:
            strcpy(opname,"CblasLower");
            break;

        case BblasUpper:
	    strcpy(opname,"CblasUpper");
            break;

        case BblasNonUnit:
	    strcpy(opname,"CblasNonUnit");
            break;

        case BblasUnit:
	   strcpy(opname,"CblasUnit");
            break;

        case BblasLeft:
	   strcpy(opname,"CblasLeft");
            break;

        case BblasRight:
	    strcpy(opname,"CblasRight");
            break;

        default:
            return 0;
            exit(EXIT_FAILURE);
    }
    return opname;
}

/**
 * Get the amount of data needed.
 **/

int bblas_cnbdata(bblas_ctest_t *test)
{
    enum BBLAS_OPTS batch_opts        = test->batch_opts;
    int nb_data;
    char function_name[NAME_LENGTH]   ="bblas_getdatacount";

    switch(batch_opts)
    {
	case BBLAS_VARIABLE:
	    nb_data =  test->batch_count;
	    break;

	case BBLAS_FIXED:
	    nb_data = 1;
	    break;

	default:
	    bblas_fatal_error(function_name, "wrong batch_opts value");
    }

    return nb_data;
}

/**
 * Inject an error into the computation if this is set in the input file.
 **/

void bblas_cset_error(bblas_ctest_t *test)
{
    int   nb_data =  bblas_cnbdata(test);
    int   routine = test->routine;
    int   error_index = irandRange(0, nb_data);

    if (test->global_error){
	test->batch_count = -1;
	return;
    }

    if (test->batch_opts == BBLAS_FIXED)
    {
	error_index = 0;
    }

    if (test->set_error)
    {
	if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	    (routine == BBLAS_HEMM) || (routine == BBLAS_TRMM) ||
	    (routine == BBLAS_TRSM))
	{
	    test->M[error_index] = -1;
	}else if ((routine == BBLAS_SYRK) || (routine == BBLAS_HERK) ||
		  (routine == BBLAS_SYR2K)|| (routine == BBLAS_HER2K))
	{
	    test->K[error_index] = -1;
	}
    }
}

/**
 * Check whether a computation has passed or failed our accuracy test.
 * When new_accuracy=1 in the input file this uses an appropriate
 * forward/backward error bound,
 * otherwise this looks at the relative error.
 **/

void bblas_cpassed_failed(bblas_ctest_t  *test, float error, char *result, int info)
{
    float eps = LAPACKE_slamch_work('e');

	/* Use our new accuracy test based on forward/backward error analysis*/
    if (test->new_accuracy)
    {
	if( (error > 1) || (info))
	{
	    strcpy(result, "FAILED");

	}else
	{
	    strcpy(result, "PASSED");
	}
    }else
    {
		/* Use old accuracy test based on the relative error */
	if((error > eps*test->tolerance) || (info))
	{
	    strcpy(result, "FAILED");

	}else
	{
	    strcpy(result, "PASSED");
	}
    }
}

/**
 * Set the batch_count.
 **/

void bblas_cset_batch_count(bblas_ctest_t *test)
{

  test->batch_count = test->minbatch_count* (test->current_iter+1);
}

#undef COMPLEX

max_stars_repo_name: NLAFET/BBLAS-ref

max_issues_repo_name: NLAFET/BBLAS-ref

max_forks_repo_name: NLAFET/BBLAS-ref


text:
/* rng/g05faf.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 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 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

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

/* This is the NAG G05FAF generator. The generator returns the
   upper 32 bits from each term of the sequence,

   x_{n+1} = (a x_n) mod m

   using 59-bit unsigned arithmetic, with a = 13^{13} and m =
   2^59. The seed specifies the upper 32 bits of the initial value,
   x_1, with the lower 16 bits set to 0x330E.

   The theoretical value of x_{10001} is 244131582646046.

   The period of this generator is ? FIXME (probably around 2^48). */

static inline void g05faf_advance (void *vstate);
static unsigned long int g05faf_get (void *vstate);
static double g05faf_get_double (void *vstate);
static void g05faf_set (void *state, unsigned long int s);

static const unsigned short int a0 = 0xFD ;
static const unsigned short int a1 = 0xC5 ;
static const unsigned short int a2 = 0x23 ;
static const unsigned short int a3 = 0x9B ;
static const unsigned short int a4 = 0x76 ;
static const unsigned short int a5 = 0x13 ;
static const unsigned short int a6 = 0x01 ;
static const unsigned short int a7 = 0x00 ;

typedef struct
  {
    unsigned short int x0, x1, x2, x3, x4, x5, x6, x7 ;
  }
g05faf_state_t;

static inline void
g05faf_advance (void *vstate)
{
  g05faf_state_t *state = (g05faf_state_t *) vstate;

  const unsigned short int x0 = state->x0 ;
  const unsigned short int x1 = state->x1 ;
  const unsigned short int x2 = state->x2 ;
  const unsigned short int x3 = state->x3 ;
  const unsigned short int x4 = state->x4 ;
  const unsigned short int x5 = state->x5 ;
  const unsigned short int x6 = state->x6 ;
  const unsigned short int x7 = state->x7 ;

  unsigned long int a ;

  /* This looks like it will be pretty slow. Maybe someone can produce
     a more efficient implementation (it has to be portable though) */

  a = a0 * x0 ;
  state->x0 = (a & 0x000000FFUL) ;

  a >>= 8 ;
  a += a0 * x1 + a1 * x0 ;
  state->x1 = (a & 0x000000FFUL) ;

  a >>= 8 ;
  a += a0 * x2 + a1 * x1 + a2 * x0 ;
  state->x2 = (a & 0x000000FFUL) ;

  a >>= 8 ;
  a += a0 * x3 + a1 * x2 + a2 * x1 + a3 * x0 ;
  state->x3 = (a & 0x000000FFUL) ;

  a >>= 8 ;
  a += a0 * x4 + a1 * x3 + a2 * x2 + a3 * x1 + a4 * x0 ;
  state->x4 = (a & 0x000000FFUL) ;

  a >>= 8 ;
  a += a0 * x5 + a1 * x4 + a2 * x3 + a3 * x2 + a4 * x1 + a5 * x0 ;
  state->x5 = (a & 0x000000FFUL) ;

  a >>= 8 ;
  a += a0 * x6 + a1 * x5 + a2 * x4 + a3 * x3 + a4 * x2 + a5 * x1 + a6 * x0 ;
  state->x6 = (a & 0x000000FFUL) ;

  a >>= 8 ;
  a += (a0 * x7 + a1 * x6 + a2 * x5 + a3 * x4 + a4 * x3 + a5 * x2 + a6 * x1
        + a7 * x0) ;
  state->x7 = (a & 0x00000007UL) ;

}

static unsigned long int
g05faf_get (void *vstate)
{
  g05faf_state_t *state = (g05faf_state_t *) vstate;

  g05faf_advance (state) ;

  return ((state->x7 << 29) + (state->x6 << 21) +
          (state->x5 << 13) + (state->x4 << 5) + (state->x3 >> 3));
}

static double
g05faf_get_double (void * vstate)
{
  g05faf_state_t *state = (g05faf_state_t *) vstate;

  g05faf_advance (state) ;

  return (ldexp((double) state->x7, -3)
          + ldexp((double) state->x6, -11)
          + ldexp((double) state->x5, -19)
          + ldexp((double) state->x4, -27)
          + ldexp((double) state->x3, -35)
          + ldexp((double) state->x2, -43)
          + ldexp((double) state->x1, -51)
          + ldexp((double) state->x0, -59)) ;
}

static void
g05faf_set (void *vstate, unsigned long int s)
{
  g05faf_state_t *state = (g05faf_state_t *) vstate;

  if (s == 0)  /* default seed */
    s = 1 ;

  /* FIXME: I have no idea what the nag seeding procedure is */

  return;
}

static const gsl_rng_type g05faf_type =
{"g05faf",                      /* name */
 0xffffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (g05faf_state_t),
 &g05faf_set,
 &g05faf_get,
 &g05faf_get_double
};

const gsl_rng_type *gsl_rng_g05faf = &g05faf_type;

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:
#ifndef __GSL_STATISTICS_H__
#define __GSL_STATISTICS_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <gsl/gsl_statistics_long_double.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_statistics_float.h>

#include <gsl/gsl_statistics_ulong.h>
#include <gsl/gsl_statistics_long.h>

#include <gsl/gsl_statistics_uint.h>
#include <gsl/gsl_statistics_int.h>

#include <gsl/gsl_statistics_ushort.h>
#include <gsl/gsl_statistics_short.h>

#include <gsl/gsl_statistics_uchar.h>
#include <gsl/gsl_statistics_char.h>

#endif /* __GSL_STATISTICS_H__ */


text:
#ifndef VM_CODE_CODE_VIEW_H
#define VM_CODE_CODE_VIEW_H

#include <cstddef>
#include <cstring>
#include <cstdint>
#include <type_traits>
#include <array>
#include <bitset>
#include <optional>
#include <variant>
#include <algorithm>
#include <iosfwd>
#include <ostream>
#include <iomanip>
#include <gsl/span>
#include <gsl/gsl>
#include "wasm_base.h"

namespace wasm::opc {

namespace detail {

template <class It>
[[gnu::pure]]
std::tuple<wasm_uint32_t, wasm_uint32_t, It> read_memory_immediate(It first, It last)
{
	alignas(wasm_uint32_t) char buff1[sizeof(wasm_uint32_t)];
	alignas(wasm_uint32_t) char buff2[sizeof(wasm_uint32_t)];
	wasm_uint32_t flags;
	wasm_uint32_t offset;
	for(auto& chr: buff1)
	{
		assert(first != last);
		chr = *first++;
	}
	for(auto& chr: buff2)
	{
		assert(first != last);
		chr = *first++;
	}
	std::memcpy(&flags, buff1, sizeof(flags));
	std::memcpy(&offset, buff1, sizeof(offset));
	return std::make_tuple(flags, offset, first);
}

template <class T, class It>
[[gnu::pure]]
std::pair<T, It> read_serialized_immediate(It first, It last)
{
	static_assert(std::is_trivially_copyable_v<T>);
	T value;
	alignas(T) char buff[sizeof(T)];
	for(auto& chr: buff)
	{
		assert(first != last);
		chr = *first++;
	}
	std::memcpy(&value, buff, sizeof(value));
	return std::make_pair(value, first);
}

} /* namespace detail */

struct BadOpcodeError:
	public std::logic_error
{
	BadOpcodeError(OpCode op, const char* msg):
		std::logic_error(msg),
		opcode(op)
	{

	}

	const OpCode opcode;
};

template <class It, class Visitor>
decltype(auto) visit_opcode(Visitor visitor, It first, It last)
{
	assert(first != last);
	OpCode op;
	auto pos = first;
	using value_type = typename std::iterator_traits<It>::value_type;
	using underlying_type = std::underlying_type_t<OpCode>;
	if constexpr(not std::is_same_v<value_type, OpCode>)
	{
		value_type opcode = *pos++;
		op = static_cast<OpCode>(opcode);
		assert(static_cast<value_type>(op) == opcode);
	}
	else
	{
		op = *pos++;
	}

	// Handling the invalid opcode case is optional.
	if constexpr(std::is_invocable_v<Visitor, OpCode, std::nullopt_t>)
	{
		if(not opcode_exists(static_cast<underlying_type>(op)))
		{
			return visitor(
				first,
				last,
				pos,
				op,
				BadOpcodeError(op, "Given op is not a valid WASM opcode.")
			);
		}
	}
	else
	{
		if(not opcode_exists(static_cast<underlying_type>(op)))
			assert(false);
	}

	if(op >= OpCode::I32_LOAD and op <= OpCode::I64_STORE32)
	{
		wasm_uint32_t flags, offset;
		std::tie(flags, offset, pos) = detail::read_memory_immediate(pos, last);
		return visitor(first, last, pos, op, flags, offset);
	}
	else if(
		(op >= OpCode::GET_LOCAL and op <= OpCode::SET_GLOBAL)
		or (op == OpCode::CALL or op == OpCode::CALL_INDIRECT)
		or (op == OpCode::BR or op == OpCode::BR_IF)
		or (op == OpCode::ELSE)
	)
	{
		wasm_uint32_t value;
		std::tie(value, pos) = detail::read_serialized_immediate<wasm_uint32_t>(pos, last);
		return visitor(first, last, pos, op, value);
	}
	else if(op == OpCode::BLOCK or op == OpCode::IF)
	{
		assert(first != last);
		LanguageType tp = static_cast<LanguageType>(*first++);
		wasm_uint32_t label;
		std::tie(label, pos) = detail::read_serialized_immediate<wasm_uint32_t>(pos, last);
		return visitor(first, last, pos, op, tp, label);
	}
	else if(op == OpCode::LOOP)
	{
		assert(first != last);
		LanguageType tp = static_cast<LanguageType>(*first++);
		return visitor(first, last, pos, op, tp);
	}
	else if(op == OpCode::BR_TABLE)
	{
		wasm_uint32_t len;
		std::tie(len, pos) = detail::read_serialized_immediate<wasm_uint32_t>(pos, last);
		auto base = pos;
		std::advance(pos, (1 + len) * sizeof(wasm_uint32_t));
		return visitor(first, last, pos, op, base, len);
	}

	switch(op)
	{
	case OpCode::I32_CONST: {
		wasm_sint32_t v;
		std::tie(v, pos) = detail::read_serialized_immediate<wasm_sint32_t>(pos, last);
		return visitor(first, last, pos, op, v);
		break;
	}
	case OpCode::I64_CONST: {
		wasm_sint64_t v;
		std::tie(v, pos) = detail::read_serialized_immediate<wasm_sint64_t>(pos, last);
		return visitor(first, last, pos, op, v);
		break;
	}
	case OpCode::F32_CONST: {
		wasm_float32_t v;
		std::tie(v, pos) = detail::read_serialized_immediate<wasm_float32_t>(pos, last);
		return visitor(first, last, pos, op, v);
		break;
	}
	case OpCode::F64_CONST: {
		wasm_float64_t v;
		std::tie(v, pos) = detail::read_serialized_immediate<wasm_float64_t>(pos, last);
		return visitor(first, last, pos, op, v);
		break;
	}
	default:
		return visitor(op, first, pos, last);
	}
	assert(false and "Internal Error: All cases should have been handled by this point.");
}

struct MemoryImmediate:
	public std::pair<const wasm_uint32_t, const wasm_uint32_t>
{
	using std::pair<const wasm_uint32_t, const wasm_uint32_t>::pair;
};

wasm_uint32_t flags(const MemoryImmediate& immed)
{ return immed.first; }

wasm_uint32_t offset(const MemoryImmediate& immed)
{ return immed.second; }

struct BlockImmediate:
	public std::pair<const LanguageType, const wasm_uint32_t>
{
	using std::pair<const LanguageType, const wasm_uint32_t>::pair;
};

gsl::span<const LanguageType> signature(const BlockImmediate& immed)
{ return gsl::span<const LanguageType>(&(immed.first), 1u); }

std::size_t arity(const BlockImmediate& immed)
{ return signature(immed).size(); }

wasm_uint32_t offset(const BlockImmediate& immed)
{ return immed.second; }

struct BranchTableImmediate
{
	template <class ... T>
	BranchTableImmediate(T&& ... args):
		table_(std::forward<T>(args)...)
	{

	}

	wasm_uint32_t at(wasm_uint32_t idx) const
	{
		idx = std::min(std::size_t(table_.size() - 1u), std::size_t(idx));
		wasm_uint32_t depth;
		std::memcpy(&depth, table_.data() + idx, sizeof(depth));
		return depth;
	}

private:
	const gsl::span<const char[sizeof(wasm_uint32_t)]> table_;
};

struct CodeView
{
	struct Iterator;
	using value_type = WasmInstruction;
	using pointer = WasmInstruction*;
	using const_pointer = const WasmInstruction*;
	using reference = WasmInstruction&;
	using const_reference = WasmInstruction&;
	using size_type = std::string_view::size_type;
	using difference_type = std::string_view::difference_type;
	using iterator = Iterator;
	using const_iterator = iterator;

private:

	struct InstructionVisitor {

		template <class ... T>
		WasmInstruction operator()(
			const char* first,
			const char* last,
			const char* pos,
			OpCode op,
			T&& ... args
		)
		{
			assert(first < pos);
			assert(pos <= last);
			return make_instr(
				std::string_view(first, pos - first),
				op,
				last,
				std::forward<T>(args)...
			);
		}

	private:
		// Overload for instructions with immediate operands
		WasmInstruction make_instr(std::string_view view, OpCode op, const char* last)
		{ return WasmInstruction(view, op, last, std::monostate()); }

		template <
			class T,
			/* enable overloads for the simple alternatives. */
			class = std::enable_if_t<
				std::disjunction_v<
					std::is_same_v<std::decay_t<T>, wasm_uint32_t>,
					std::is_same_v<std::decay_t<T>, wasm_sint32_t>,
					std::is_same_v<std::decay_t<T>, wasm_sint64_t>,
					std::is_same_v<std::decay_t<T>, wasm_float32_t>,
					std::is_same_v<std::decay_t<T>, wasm_float64_t>
				>
			>
		>
		WasmInstruction make_instr(std::string_view view, OpCode op, const char* last, T&& arg)
		{ return WasmInstruction(view, op, last, std::forward<T>(arg)); }

		/// Loop overload
		WasmInstruction make_instr(std::string_view view, OpCode op, const char* last, LanguageType tp)
		{ return WasmInstruction(view, op, last, tp); }

		/// Branch table overload
		WasmInstruction make_instr(std::string_view view, OpCode op, const char* last, const char* base, wasm_uint32_t len)
		{
			assert(last > base);
			std::size_t byte_count = base - last;
			std::size_t bytes_needed = sizeof(wasm_uint32_t) * (len + 1u);
			assert(byte_count >= bytes_needed);
			// the table and 'view' should end at the same byte address
			assert(view.data() + view.size() == (base + (len + 1u) * sizeof(wasm_uint32_t)));
			using buffer_type = const char[sizeof(wasm_uint32_t)];
			auto table = gsl::span<buffer_type>(reinterpret_cast<buffer_type*>(base), len + 1u);
			return WasmInstruction(view, op, last, table);
		}

		/// Block overload
		WasmInstruction make_instr(std::string_view view, OpCode op, const char* last, LanguageType tp, wasm_uint32_t label)
		{ return WasmInstruction(view, op, last, BlockImmediate(tp, label)); }

		/// Memory overload
		WasmInstruction make_instr(std::string_view view, OpCode op, const char* last, wasm_uint32_t flags, wasm_uint32_t offset)
		{ return WasmInstruction(view, op, last, MemoryImmediate(flags, offset)); }

		/// Invalid opcode overload
		[[noreturn]]
		WasmInstruction make_instr(std::string_view, OpCode, const char*, const BadOpCodeError& err)
		{ throw err; }

	};

public:

	struct Iterator {
		using value_type = CodeView::value_type;
		using difference_type = CodeView::difference_type;
		using pointer = CodeView::pointer;
		using reference = CodeView::value_type;
		using iterator_category = std::input_iterator_tag;

		friend bool operator==(const Iterator& left, const Iterator& right)
		{ return left.code_ == right.code_; }

		friend bool operator!=(const Iterator& left, const Iterator& right)
		{ return not (left == right); }
	private:
		gsl::not_null<CodeView*> code_;
	};

	CodeView(const WasmFunction& func):
		function_(&func),
		code_(code(func))
	{

	}

	const WasmFunction* function() const
	{ return function_; }

	OpCode current_op() const
	{
		assert(ready());
		return static_cast<WasmInstruction>(code_.front());
	}

	bool done() const
	{ return not code_.empty(); }

	std::optional<WasmInstruction> next_instruction() const
	{
		assert(ready());
		if(code_.size() == 0)
			return std::nullopt;
		return visit_opcode(InstructionVisitor{}, code_.data(), code_.data() + code_.size());
	}

	const char* pos() const
	{ return code_.data(); }

	void advance(const CodeView& other)
	{
		assert(function() == other.function());
		assert(ready());
		assert(code_.data() + code_.size() == other.code_.data() + other.code_.size());
		assert(code_.data() < other.code_.data());
		code_ = other.code_;
	}

private:

	bool ready() const
	{
		if(done())
			return false;
		assert(opcode_exists(code_.front()));
		return true;
	}
	const WasmFunction* function_;
	std::string_view code_;
};

} /* namespace opc */
} /* namespace wasm */

#endif /* VM_CODE_CODE_VIEW_H */

max_stars_repo_name: tvanslyke/wasm-cpp

max_issues_repo_name: tvanslyke/wasm-cpp

max_forks_repo_name: tvanslyke/wasm-cpp


text:
// Copyright (c) 2005 - 2009 Marc de Kamps, Dave Harrison
// 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 the copyright holder 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.
//
//      If you use this software in work leading to a scientific publication, you should cite
//      the 'currently valid reference', which can be found at http://miind.sourceforge.net
#ifndef _CODE_LIBS_NUMTOOLSLIB_INTERPOLATION_INCLUDE_GUARD
#define _CODE_LIBS_NUMTOOLSLIB_INTERPOLATION_INCLUDE_GUARD
// Date:   9-02-2009
// Author: Dave Harrison

#include <valarray>
#include <utility>
#include <gsl/gsl_interp.h>

namespace NumtoolsLib
{
    //! Define types of interpolation
    enum InterpType
    {
        INTERP_LINEAR,
        INTERP_CSPLINE,
        INTERP_AKIMA
    };

    class Interpolator
    {
        public:
            /*!
             Interpolator object to wrap GSL interpolation functions
             interp_type Type of Interpolation to use: INTERP_LINEAR/_CSPLINE/_AKIMA
             xs Array of x values
             ys array of y values
             Throws NumtoolsException on GSL error
             */
            Interpolator(const InterpType interp_type, const std::valarray<double>& xs,
                    const std::valarray<double>& ys);

            ~Interpolator();

            /*!
             Retrieve an interpolated y value for position x
             Throws NumtoolsLibException on error
             */
            double InterpValue(const double x);

        private:
            //! The accelerator object
            gsl_interp_accel*		_acc;
            //! The interpolation object
            gsl_interp*				_interp;
            //! The x values
            std::valarray<double>& _x;
            //! The y values
            std::valarray<double>& _y;
            //! The size of the x & y arrays
            const size_t			_size;
    };
}

#endif // include guard

max_stars_repo_name: dekamps/miind

max_issues_repo_name: dekamps/miind

max_forks_repo_name: dekamps/miind