text:
/* eigen/gsl_eigen.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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_EIGEN_H__
#define __GSL_EIGEN_H__

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.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 size;
  double * d;
  double * sd;
} gsl_eigen_symm_workspace;

gsl_eigen_symm_workspace * gsl_eigen_symm_alloc (const size_t n);
void gsl_eigen_symm_free (gsl_eigen_symm_workspace * w);
int gsl_eigen_symm (gsl_matrix * A, gsl_vector * eval, gsl_eigen_symm_workspace * w);

typedef struct {
  size_t size;
  double * d;
  double * sd;
  double * gc;
  double * gs;
} gsl_eigen_symmv_workspace;

gsl_eigen_symmv_workspace * gsl_eigen_symmv_alloc (const size_t n);
void gsl_eigen_symmv_free (gsl_eigen_symmv_workspace * w);
int gsl_eigen_symmv (gsl_matrix * A, gsl_vector * eval, gsl_matrix * evec, gsl_eigen_symmv_workspace * w);

typedef struct {
  size_t size;           /* matrix size */
  size_t max_iterations; /* max iterations since last eigenvalue found */
  size_t n_iter;         /* number of iterations since last eigenvalue found */
  size_t n_evals;        /* number of eigenvalues found so far */

  gsl_vector *hv2;       /* temporary 2-by-1 householder vector */
  gsl_vector *hv3;       /* temporary 3-by-1 householder vector */

  int compute_t;         /* compute Schur form T = Z^t A Z */

  gsl_matrix *H;         /* pointer to Hessenberg matrix */

  gsl_matrix *Z;         /* pointer to Schur vector matrix */
} gsl_eigen_francis_workspace;

gsl_eigen_francis_workspace * gsl_eigen_francis_alloc (void);
void gsl_eigen_francis_free (gsl_eigen_francis_workspace * w);
void gsl_eigen_francis_T (const int compute_t,
                          gsl_eigen_francis_workspace * w);
int gsl_eigen_francis (gsl_matrix * H, gsl_vector_complex * eval,
                       gsl_eigen_francis_workspace * w);
int gsl_eigen_francis_Z (gsl_matrix * H, gsl_vector_complex * eval,
                         gsl_matrix * Z,
                         gsl_eigen_francis_workspace * w);

typedef struct {
  size_t size;                 /* size of matrices */
  gsl_vector *diag;            /* diagonal matrix elements from balancing */
  gsl_vector *tau;             /* Householder coefficients */
  gsl_matrix *Z;               /* pointer to Z matrix */
  int do_balance;              /* perform balancing transformation? */
  size_t n_evals;              /* number of eigenvalues found */

  gsl_eigen_francis_workspace *francis_workspace_p;
} gsl_eigen_nonsymm_workspace;

gsl_eigen_nonsymm_workspace * gsl_eigen_nonsymm_alloc (const size_t n);
void gsl_eigen_nonsymm_free (gsl_eigen_nonsymm_workspace * w);
void gsl_eigen_nonsymm_params (const int compute_t, const int balance,
                               gsl_eigen_nonsymm_workspace *w);
int gsl_eigen_nonsymm (gsl_matrix * A, gsl_vector_complex * eval,
                       gsl_eigen_nonsymm_workspace * w);
int gsl_eigen_nonsymm_Z (gsl_matrix * A, gsl_vector_complex * eval,
                         gsl_matrix * Z, gsl_eigen_nonsymm_workspace * w);

typedef struct {
  size_t size;                 /* size of matrices */
  gsl_vector *work;            /* scratch workspace */
  gsl_vector *work2;           /* scratch workspace */
  gsl_vector *work3;           /* scratch workspace */

  gsl_matrix *Z;               /* pointer to Schur vectors */

  gsl_eigen_nonsymm_workspace *nonsymm_workspace_p;
} gsl_eigen_nonsymmv_workspace;

gsl_eigen_nonsymmv_workspace * gsl_eigen_nonsymmv_alloc (const size_t n);
void gsl_eigen_nonsymmv_free (gsl_eigen_nonsymmv_workspace * w);
int gsl_eigen_nonsymmv (gsl_matrix * A, gsl_vector_complex * eval,
                        gsl_matrix_complex * evec,
                        gsl_eigen_nonsymmv_workspace * w);
int gsl_eigen_nonsymmv_Z (gsl_matrix * A, gsl_vector_complex * eval,
                          gsl_matrix_complex * evec, gsl_matrix * Z,
                          gsl_eigen_nonsymmv_workspace * w);

typedef struct {
  size_t size;
  double * d;
  double * sd;
  double * tau;
} gsl_eigen_herm_workspace;

gsl_eigen_herm_workspace * gsl_eigen_herm_alloc (const size_t n);
void gsl_eigen_herm_free (gsl_eigen_herm_workspace * w);
int gsl_eigen_herm (gsl_matrix_complex * A, gsl_vector * eval,
                         gsl_eigen_herm_workspace * w);

typedef struct {
  size_t size;
  double * d;
  double * sd;
  double * tau;
  double * gc;
  double * gs;
} gsl_eigen_hermv_workspace;

gsl_eigen_hermv_workspace * gsl_eigen_hermv_alloc (const size_t n);
void gsl_eigen_hermv_free (gsl_eigen_hermv_workspace * w);
int gsl_eigen_hermv (gsl_matrix_complex * A, gsl_vector * eval,
                           gsl_matrix_complex * evec,
                           gsl_eigen_hermv_workspace * w);

typedef enum {
  GSL_EIGEN_SORT_VAL_ASC,
  GSL_EIGEN_SORT_VAL_DESC,
  GSL_EIGEN_SORT_ABS_ASC,
  GSL_EIGEN_SORT_ABS_DESC
}
gsl_eigen_sort_t;

/* Sort eigensystem results based on eigenvalues.
 * Sorts in order of increasing value or increasing
 * absolute value.
 *
 * exceptions: GSL_EBADLEN
 */

int gsl_eigen_symmv_sort(gsl_vector * eval, gsl_matrix * evec,
                         gsl_eigen_sort_t sort_type);

int gsl_eigen_hermv_sort(gsl_vector * eval, gsl_matrix_complex * evec,
                         gsl_eigen_sort_t sort_type);

int gsl_eigen_nonsymmv_sort(gsl_vector_complex * eval,
                            gsl_matrix_complex * evec,
                            gsl_eigen_sort_t sort_type);

/* The following functions are obsolete: */

/* Eigensolve by Jacobi Method
 *
 * The data in the matrix input is destroyed.
 *
 * exceptions:
 */
int
gsl_eigen_jacobi(gsl_matrix * matrix,
                      gsl_vector * eval,
                      gsl_matrix * evec,
                      unsigned int max_rot,
                      unsigned int * nrot);

/* Invert by Jacobi Method
 *
 * exceptions:
 */
int
gsl_eigen_invert_jacobi(const gsl_matrix * matrix,
                             gsl_matrix * ainv,
                             unsigned int max_rot);

__END_DECLS

#endif /* __GSL_EIGEN_H__ */

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#include <stdio.h>
#include <string.h>
#include "erl_nif.h"
#include <cblas.h>

#define MAX(a,b) (((a)>(b))?(a):(b))

static ERL_NIF_TERM atom_ok;
static ERL_NIF_TERM atom_true;

static ERL_NIF_TERM atom_rowmaj;
static ERL_NIF_TERM atom_colmaj;

static ERL_NIF_TERM atom_notransp;
static ERL_NIF_TERM atom_transpose;
static ERL_NIF_TERM atom_conjugatet;

static ERL_NIF_TERM atom_upper;
static ERL_NIF_TERM atom_lower;

static ERL_NIF_TERM atom_left;
static ERL_NIF_TERM atom_right;

static ERL_NIF_TERM atom_nonunit;
static ERL_NIF_TERM atom_unit;

static ERL_NIF_TERM make_cont(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[]);
static ERL_NIF_TERM from_list(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[]);
static ERL_NIF_TERM to_values(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[]);
static ERL_NIF_TERM to_idx_list(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[]);
static ERL_NIF_TERM cont_size(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[]);
static ERL_NIF_TERM cont_update(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[]);

static ERL_NIF_TERM drotg(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[]);
static ERL_NIF_TERM drot(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[]);
static ERL_NIF_TERM dcopy(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[]);
static ERL_NIF_TERM l1d_1cont(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[]);
static ERL_NIF_TERM l1d_2cont(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[]);
static ERL_NIF_TERM dgemv(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[]);
static ERL_NIF_TERM dtrmv(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[]);
static ERL_NIF_TERM dgemm(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[]);

static ErlNifFunc nif_funcs[] = {
    {"make_cont", 2, make_cont},
    {"from_list", 1, from_list},
    {"to_tuple_list_impl", 4, to_values},
    {"cont_size", 1, cont_size},
    {"values", 2, to_idx_list},
    {"update", 2, cont_update},
    {"update", 3, cont_update},

    {"rotg", 2, drotg},
    {"rot",  9, drot},
    {"copy", 4, dcopy},
    {"copy", 7, dcopy},
    {"one_vec", 6, l1d_1cont},
    {"two_vec", 9, l1d_2cont},

    {"gemv_impl", 15, dgemv},
    {"trmv_impl", 11, dtrmv},
    {"trmv_impl", 12, dtrmv},

    {"gemm_impl", 15, dgemm},
    {"trmm_impl", 15, dgemm},
};

static ErlNifResourceType *avec_r;

typedef struct {
    double v[1];
} Avec;

static ERL_NIF_TERM mk_avec(ErlNifEnv *env, unsigned int n, Avec **avec) {
    ERL_NIF_TERM term;
    *avec = enif_alloc_resource(avec_r, n*sizeof(double));
    term = enif_make_resource(env, *avec);
    enif_release_resource(*avec);
    return term;
}

#define AVEC_SIZE(AV) 	(enif_sizeof_resource(AV) / 8)

/* ---------------------------------------------------*/

static ERL_NIF_TERM make_cont(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[])
{
    unsigned int n;
    ERL_NIF_TERM res;
    Avec *avec = NULL;
    double *data;

    if(!enif_get_uint(env, argv[0], &n)) return enif_make_badarg(env);
    res = mk_avec(env, n, &avec);

    if(enif_is_binary(env, argv[1])) {
	ErlNifBinary bin;
	enif_inspect_binary(env, argv[1], &bin);
	if(n*sizeof(double) != bin.size)
	    return enif_make_badarg(env);
	memcpy(avec->v, bin.data, sizeof(double)*n);
    } else if(enif_is_identical(argv[1], atom_true)) {
	data = avec->v;
	memset(data, 0, sizeof(double)*n);
    }

    return res;
}

static ERL_NIF_TERM from_list(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[])
{
    unsigned int len;
    int dim, i,j, tmp;
    ERL_NIF_TERM hd, tail, res;
    const ERL_NIF_TERM *current;
    Avec *avec = NULL;
    double *data;

    if(!enif_get_list_length(env, argv[0], &len) || len == 0)
	return enif_make_badarg(env);

    enif_get_list_cell(env, argv[0], &hd, &tail);
    if(!enif_get_tuple(env, hd, &dim, &current))
	dim = 1;
    res = mk_avec(env, len*dim, &avec);

    enif_consume_timeslice(env, (len*dim)/1000);

    data = avec->v;
    if(dim == 1) {
	for(i=0; i<len; i++) {
	    if(!enif_get_double(env, hd, data))
		return enif_make_badarg(env);
	    data++;
	    enif_get_list_cell(env, tail, &hd, &tail);
	}
    } else {
	for(i=0; i<len; i++) {
	    if(!enif_get_tuple(env, hd, &tmp, &current) || tmp != dim)
		return enif_make_badarg(env);
	    for(j=0; j<dim; j++) {
		if(!enif_get_double(env, current[j], data))
		    return enif_make_badarg(env);
		data++;
	    }
	    enif_get_list_cell(env, tail, &hd, &tail);
	}
    }
    return res;
}

static ERL_NIF_TERM to_values(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[])
{
    ERL_NIF_TERM tail, *tmp;
    Avec *avec = NULL;
    double *arr;
    unsigned int idx, n, dim, i, j, max;

    if(!enif_get_uint(env, argv[0], &idx)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[1], &n))   return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[2], &dim))  return enif_make_badarg(env);

    if(!enif_get_resource(env, argv[3], avec_r, (void **) &avec))
	return enif_make_badarg(env);
    max = AVEC_SIZE(avec);

    enif_consume_timeslice(env, max/1000);

    if(idx+n > max) return enif_make_badarg(env);

    if(n == 1 && dim == 0)  /* get_value() */
	return enif_make_double(env, avec->v[idx]);

    if(dim == 0) { /* to_list */
	tail = enif_make_list(env, 0);
	arr = avec->v + idx + n-1;
	for(i=0; i < max; i++) {
	    tail = enif_make_list_cell(env, enif_make_double(env, *arr), tail);
	    arr -= 1;
	}
	return tail;
    }
    if(dim == n) { /* to_tuple */
	tmp = (ERL_NIF_TERM*) malloc(sizeof(ERL_NIF_TERM)*n);
	arr = avec->v;
	for(i=0; i < n; i++) {
	    tmp[i] = enif_make_double(env, arr[idx+i]);
	}
	tail = enif_make_tuple_from_array(env, tmp, n);
	free(tmp);
	return tail;
    }

    /* list of tuples */
    if(n % dim != 0)
	return enif_make_badarg(env);

    arr = avec->v + idx + n-dim;
    tmp = (ERL_NIF_TERM*) malloc(sizeof(ERL_NIF_TERM)*dim);
    tail = enif_make_list(env, 0);

    for(i=0; i < max / dim; i++) {
	for(j=0; j < dim; j++, arr++) {
	    tmp[j] = enif_make_double(env, *arr);
	}
	tail = enif_make_list_cell(env, enif_make_tuple_from_array(env, tmp, dim), tail);
	arr -= 2*dim;
    }
    free(tmp);
    return tail;
}

static ERL_NIF_TERM to_idx_list(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[])
{
    ERL_NIF_TERM hd, tail, *tmp;
    Avec *avec = NULL;
    double *arr;
    unsigned int idx, n, i;

    if(!enif_get_list_length(env, argv[0], &n) || n == 0)
	return enif_make_badarg(env);
    if(!enif_get_resource(env, argv[1], avec_r, (void **) &avec))
	return enif_make_badarg(env);

    enif_consume_timeslice(env, n/1000);
    tmp = (ERL_NIF_TERM*) malloc(sizeof(ERL_NIF_TERM)*n);
    arr = avec->v;
    tail = argv[0];
    for(i=0; i < n; i++) {
	enif_get_list_cell(env, tail, &hd, &tail);
	if(!enif_get_uint(env, hd, &idx) || idx >= AVEC_SIZE(avec))
	    return enif_make_badarg(env);
	tmp[i] = enif_make_tuple2(env, hd, enif_make_double(env, arr[idx]));
    }
    return enif_make_list_from_array(env, tmp, n);
}

static ERL_NIF_TERM cont_size(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[])
{
    Avec *avec = NULL;

    if(!enif_get_resource(env, argv[0], avec_r, (void **) &avec)) {
	return enif_make_badarg(env);
    }

    return enif_make_uint(env, AVEC_SIZE(avec));
}

static ERL_NIF_TERM cont_update(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[])
{
    ERL_NIF_TERM hd, tail;
    const ERL_NIF_TERM *curr;
    Avec *avec = NULL;
    double *arr, temp;
    int i=0, n=0;
    unsigned int idx, max;

    if(argc == 2) { /* tuple list of values */
	if(!enif_is_list(env, argv[0]))
	    return enif_make_badarg(env);
	if(!enif_get_resource(env, argv[1], avec_r, (void **) &avec))
	    return enif_make_badarg(env);
	arr = avec->v;
	max = AVEC_SIZE(avec);
	tail = argv[0];
	while(!enif_is_empty_list(env, tail)) {
	    if(!enif_get_list_cell(env, tail, &hd, &tail))
		return enif_make_badarg(env);
	    if(!enif_get_tuple(env, hd, &i, &curr) || i != 2)
		return enif_make_badarg(env);
	    if(!enif_get_uint(env, curr[0], &idx) || idx >= max)
		return enif_make_badarg(env);
	    if(!enif_get_double(env, curr[1], &temp))
		return enif_make_badarg(env);
	    arr[idx] = temp;
	    n++;
	}
	enif_consume_timeslice(env, n/1000);
	return atom_ok;
    }
    /* update(Idx, V|[Vs], Vec) */
    if(!enif_get_uint(env, argv[0], &idx))
	return enif_make_badarg(env);
    if(!enif_is_list(env, argv[1])) {
	/* update(Idx, Double, Vec) */
	if(!enif_get_double(env, argv[1], &temp))
	    return enif_make_badarg(env);
	if(!enif_get_resource(env, argv[2], avec_r, (void **) &avec))
	    return enif_make_badarg(env);
	max = AVEC_SIZE(avec);
	if(idx >= max)
	    return enif_make_badarg(env);
	arr = avec->v;
	arr[idx] = temp;
	return atom_ok;
    }
    /* update(Idx, [Values], Vec) */
    tail = argv[1];
    if(!enif_get_resource(env, argv[2], avec_r, (void **) &avec))
	return enif_make_badarg(env);
    max = AVEC_SIZE(avec);
    arr = avec->v+idx;
    while(!enif_is_empty_list(env, tail)) {
	if(!enif_get_list_cell(env, tail, &hd, &tail))
	    return enif_make_badarg(env);
	if(!enif_get_double(env, hd, &temp))
	    return enif_make_badarg(env);
	*arr = temp;
	arr++;
	idx++;
	n++;
	if(idx > max)
	    return enif_make_badarg(env);
    }
    enif_consume_timeslice(env, n/1000);
    return atom_ok;
}

/* ---------------------------------------------------*/

static ERL_NIF_TERM drotg(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[])
{/* (a, b) */
    double a,b,c,s;

    if(!enif_get_double(env, argv[0], &a) || !enif_get_double(env, argv[1], &b)) {
	return enif_make_badarg(env);
    }

    cblas_drotg(&a,&b,&c,&s);
    return enif_make_tuple4(env,
			    enif_make_double(env, a), enif_make_double(env, b),
			    enif_make_double(env, c), enif_make_double(env, s));

}

static ERL_NIF_TERM drot(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[])
{
    unsigned int n, xs, ys;
    int xi, yi;
    Avec *ax, *ay;
    double *x, *y, c, s;

    if(!enif_get_uint(env, argv[0], &n)) return enif_make_badarg(env);
    if(!enif_get_resource(env, argv[1], avec_r, (void **) &ax)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[2], &xs)) return enif_make_badarg(env);
    if(!enif_get_int(env, argv[3], &xi)) return enif_make_badarg(env);
    if(!enif_get_resource(env, argv[4], avec_r, (void **) &ay)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[5], &ys)) return enif_make_badarg(env);
    if(!enif_get_int(env, argv[6], &yi)) return enif_make_badarg(env);
    if(!enif_get_double(env, argv[7], &c)) return enif_make_badarg(env);
    if(!enif_get_double(env, argv[8], &s)) return enif_make_badarg(env);

    x = ax->v + xs;
    y = ay->v + ys;
    /* array limit checks */
    if((xs+n*abs(xi)-1) > AVEC_SIZE(ax)) return enif_make_badarg(env);
    if((ys+n*abs(yi)-1) > AVEC_SIZE(ax)) return enif_make_badarg(env);

    enif_consume_timeslice(env, n/1000);
    cblas_drot(n, x, xi, y, xi, c, s);

    return atom_ok;
}

static ERL_NIF_TERM l1d_1cont(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[])
{
    unsigned int op, n, xs;
    int xi, tmp;
    Avec *ax;
    double *x, alpha;
    ERL_NIF_TERM res = atom_ok;

    if(!enif_get_uint(env, argv[0], &n)) return enif_make_badarg(env);
    if(!enif_get_double(env, argv[1], &alpha)) return enif_make_badarg(env);
    if(!enif_get_resource(env, argv[2], avec_r, (void **) &ax)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[3], &xs)) return enif_make_badarg(env);
    if(!enif_get_int(env, argv[4], &xi)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[5], &op)) return enif_make_badarg(env);
    x = ax->v + xs;
    /* array limit checks */
    if((xs+n*abs(xi)-1) > AVEC_SIZE(ax)) return enif_make_badarg(env);

    switch(op) {
    case 0:
	cblas_dscal(n, alpha, x, xi);
	break;
    case 1:
	res = enif_make_double(env, cblas_dnrm2(n, x, xi));
	break;
    case 2:
	res = enif_make_double(env, cblas_dasum(n, x, xi));
	break;
    case 3:
	tmp = cblas_idamax(n, x, xi);
	res = enif_make_tuple2(env,
			       enif_make_uint(env, tmp),
			       enif_make_double(env, x[tmp]));
	break;
    }
    enif_consume_timeslice(env, n/1000);
    return res;
}

static ERL_NIF_TERM dcopy(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[])
{
    unsigned int n, xs, ys = 0;
    int xi, yi=1;
    Avec *ax, *ay;
    double *x, *y;
    ERL_NIF_TERM res;

    if(!enif_get_uint(env, argv[0], &n)) return enif_make_badarg(env);
    if(!enif_get_resource(env, argv[1], avec_r, (void **) &ax)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[2], &xs)) return enif_make_badarg(env);
    if(!enif_get_int(env, argv[3], &xi)) return enif_make_badarg(env);
    /* array limit checks */
    if((xs+n*abs(xi)-1) > AVEC_SIZE(ax)) return enif_make_badarg(env);

    if(argc == 4) {  /* copy to new vector */
	res = mk_avec(env, n, &ay);
    } else {
	if(!enif_get_resource(env, argv[4], avec_r, (void **) &ay)) return enif_make_badarg(env);
	if(!enif_get_uint(env, argv[5], &ys)) return enif_make_badarg(env);
	if(!enif_get_int(env, argv[6], &yi)) return enif_make_badarg(env);
	if((ys+n*abs(yi)-1) > AVEC_SIZE(ay)) return enif_make_badarg(env);
	res = atom_ok;
    }

    x = ax->v + xs;
    y = ay->v + ys;
    cblas_dcopy(n, x, xi, y, yi);
    enif_consume_timeslice(env, n/1000);

    return res;
}

static ERL_NIF_TERM l1d_2cont(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[])
{
    unsigned int op, n, xs, ys;
    int xi, yi;
    Avec *ax, *ay;
    double *x, *y, alpha, res;

    if(!enif_get_uint(env, argv[0], &n)) return enif_make_badarg(env);
    if(!enif_get_double(env, argv[1], &alpha)) return enif_make_badarg(env);
    if(!enif_get_resource(env, argv[2], avec_r, (void **) &ax)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[3], &xs)) return enif_make_badarg(env);
    if(!enif_get_int(env, argv[4], &xi)) return enif_make_badarg(env);
    if(!enif_get_resource(env, argv[5], avec_r, (void **) &ay)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[6], &ys)) return enif_make_badarg(env);
    if(!enif_get_int(env, argv[7], &yi)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[8], &op)) return enif_make_badarg(env);
    x = ax->v + xs;
    y = ay->v + ys;
    /* array limit checks */
    if((xs+n*abs(xi)-1) > AVEC_SIZE(ax)) return enif_make_badarg(env);
    if((ys+n*abs(yi)-1) > AVEC_SIZE(ay)) return enif_make_badarg(env);

    enif_consume_timeslice(env, n/1000);
    switch(op) {
    case 0:
	cblas_daxpy(n, alpha, x, xi, y, yi);
	break;
    case 1:
	cblas_dswap(n, x, xi, y, yi);
	break;
    case 2:
	res = cblas_ddot(n, x, xi, y, yi);
	return enif_make_double(env, res);
    }
    return atom_ok;
}

/* ---------------------------------------------------*/

static ERL_NIF_TERM dgemv(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[])
{
    unsigned int m, n, lda, xs, ys, op;
    int xi, yi;
    Avec *aa, *ax, *ay;
    double *a, *x, *y, alpha, beta;
    enum CBLAS_ORDER order = 0;
    enum CBLAS_TRANSPOSE trans = 0;
    enum CBLAS_UPLO uplo = 0;

    if(enif_is_identical(argv[0], atom_rowmaj)) order = CblasRowMajor;
    else if(enif_is_identical(argv[0], atom_colmaj)) order = CblasColMajor;
    else return enif_make_badarg(env);

    if(enif_is_identical(argv[1], atom_notransp)) trans = CblasNoTrans;
    else if(enif_is_identical(argv[1], atom_transpose)) trans = CblasTrans;
    else if(enif_is_identical(argv[1], atom_conjugatet)) trans = CblasConjTrans;
    else if(enif_is_identical(argv[1], atom_lower)) uplo = CblasLower;
    else if(enif_is_identical(argv[1], atom_upper)) uplo = CblasUpper;
    else return enif_make_badarg(env);

    if(!enif_get_uint(env, argv[2], &m)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[3], &n)) return enif_make_badarg(env);
    if(!enif_get_double(env, argv[4], &alpha)) return enif_make_badarg(env);
    if(!enif_get_resource(env, argv[5], avec_r, (void **) &aa)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[6], &lda)) return enif_make_badarg(env);
    if(!enif_get_resource(env, argv[7], avec_r, (void **) &ax)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[8], &xs)) return enif_make_badarg(env);
    if(!enif_get_int(env, argv[9], &xi)) return enif_make_badarg(env);
    if(!enif_get_double(env, argv[10], &beta)) return enif_make_badarg(env);
    if(!enif_get_resource(env, argv[11], avec_r, (void **) &ay)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[12], &ys)) return enif_make_badarg(env);
    if(!enif_get_int(env, argv[13], &yi)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[14], &op)) return enif_make_badarg(env);

    a = aa->v;
    x = ax->v + xs;
    y = ay->v + ys;

    /* array limit checks */
    if((xs+n*abs(xi)-1) > AVEC_SIZE(ax)) return enif_make_badarg(env);
    if((ys+n*abs(yi)-1) > AVEC_SIZE(ay)) return enif_make_badarg(env);

    switch(op) {
    case 0:
	cblas_dgemv(order, trans, m, n, alpha, a, lda, x, xi, beta, y, yi); break;
    case 1:
	cblas_dsymv(order, uplo, n, alpha, a, lda, x, xi, beta, y, yi); break;
    case 2:
	cblas_dspmv(order, uplo, n, alpha, a, x, xi, beta, y, yi); break;
    case 3:
	cblas_dger(order, m, n, alpha, x, xi, y, yi, a, lda); break;
    case 4:
	cblas_dsyr2(order, uplo, n, alpha, x, xi, y, yi, a, lda); break;
    case 5:
	cblas_dspr2(order, uplo, n, alpha, x, xi, y, yi, a); break;
    }
    enif_consume_timeslice(env, MAX(1,m)*n/1000);
    return atom_ok;
}

static ERL_NIF_TERM dtrmv(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[])
{
    unsigned int n, lda, xs, op;
    int xi;
    Avec *aa, *ax;
    double *a, *x, alpha;
    enum CBLAS_ORDER order = 0;
    enum CBLAS_TRANSPOSE trans = 0;
    enum CBLAS_UPLO uplo = 0;
    enum CBLAS_DIAG diag = 0;

    if(enif_is_identical(argv[0], atom_rowmaj)) order = CblasRowMajor;
    else if(enif_is_identical(argv[0], atom_colmaj)) order = CblasColMajor;
    else return enif_make_badarg(env);

    if(enif_is_identical(argv[1], atom_lower)) uplo = CblasLower;
    else if(enif_is_identical(argv[1], atom_upper)) uplo = CblasUpper;
    else return enif_make_badarg(env);

    if(!enif_get_uint(env, argv[10], &op)) return enif_make_badarg(env);
    if(op < 4) {
	if(enif_is_identical(argv[2], atom_notransp)) trans = CblasNoTrans;
	else if(enif_is_identical(argv[2], atom_transpose)) trans = CblasTrans;
	else if(enif_is_identical(argv[2], atom_conjugatet)) trans = CblasConjTrans;
	else return enif_make_badarg(env);

	if(enif_is_identical(argv[3], atom_nonunit)) diag = CblasNonUnit;
	else if(enif_is_identical(argv[3], atom_unit)) diag = CblasUnit;
	else return enif_make_badarg(env);
    }

    if(!enif_get_uint(env, argv[4], &n)) return enif_make_badarg(env);
    if(!enif_get_resource(env, argv[5], avec_r, (void **) &aa)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[6], &lda)) return enif_make_badarg(env);
    if(!enif_get_resource(env, argv[7], avec_r, (void **) &ax)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[8], &xs)) return enif_make_badarg(env);
    if(!enif_get_int(env, argv[9], &xi)) return enif_make_badarg(env);
    if(argc > 11)
	if(!enif_get_double(env, argv[11], &alpha)) return enif_make_badarg(env);

    a = aa->v;
    x = ax->v + xs;

    /* array limit checks */
    if((xs+n*abs(xi)-1) > AVEC_SIZE(ax)) return enif_make_badarg(env);

    switch(op) {
    case 0:
	cblas_dtrmv(order, uplo, trans, diag, n, a, lda, x, xi);break;
    case 1:
	cblas_dtpmv(order, uplo, trans, diag, n, a, x, xi); break;
    case 2:
	cblas_dtrsv(order, uplo, trans, diag, n, a, lda, x, xi);break;
    case 3:
	cblas_dtpsv(order, uplo, trans, diag, n, a, x, xi);break;
    case 4:
	cblas_dsyr(order, uplo, n, alpha, x, xi, a, lda); break;
    case 5:
	cblas_dspr(order, uplo, n, alpha, x, xi, a); break;
    }

    enif_consume_timeslice(env, n*n/1000);
    return atom_ok;
}

/* ---------------------------------------------------*/

static ERL_NIF_TERM dgemm(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[])
{
    unsigned int m, n, k, lda, ldb, ldc, op;
    Avec *aa, *ab, *ac;
    double alpha, beta, *a, *b, *c = NULL;
    enum CBLAS_ORDER order = 0;
    enum CBLAS_TRANSPOSE ta = 0, tb = 0;
    enum CBLAS_UPLO uplo = 0;
    enum CBLAS_SIDE side = 0;
    enum CBLAS_DIAG diag = 0;

    if(enif_is_identical(argv[0], atom_rowmaj)) order = CblasRowMajor;
    else if(enif_is_identical(argv[0], atom_colmaj)) order = CblasColMajor;
    else return enif_make_badarg(env);

    if(!enif_get_uint(env, argv[3], &m)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[4], &n)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[5], &k)) return enif_make_badarg(env);

    if(!enif_get_double(env, argv[6], &alpha)) return enif_make_badarg(env);

    if(!enif_get_resource(env, argv[7], avec_r, (void **) &aa)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[8], &lda)) return enif_make_badarg(env);

    if(!enif_get_resource(env, argv[9], avec_r, (void **) &ab)) return enif_make_badarg(env);
    if(!enif_get_uint(env, argv[10], &ldb)) return enif_make_badarg(env);

    if(!enif_get_double(env, argv[11], &beta)) return enif_make_badarg(env);

    if(!enif_get_uint(env, argv[12], &op)) return enif_make_badarg(env);
    if(op < 3) {
	if(!enif_get_resource(env, argv[13], avec_r, (void **) &ac))
	    return enif_make_badarg(env);
	if(!enif_get_uint(env, argv[14], &ldc)) return enif_make_badarg(env);
	c = ac->v;
    }
    a = aa->v;
    b = ab->v;
    switch(op) {
    case 0: {
	if(enif_is_identical(argv[1], atom_notransp)) ta = CblasNoTrans;
	else if(enif_is_identical(argv[1], atom_transpose)) ta = CblasTrans;
	else if(enif_is_identical(argv[1], atom_conjugatet)) ta = CblasConjTrans;
	else return enif_make_badarg(env);

	if(enif_is_identical(argv[2], atom_notransp)) tb = CblasNoTrans;
	else if(enif_is_identical(argv[2], atom_transpose)) tb = CblasTrans;
	else if(enif_is_identical(argv[2], atom_conjugatet)) tb = CblasConjTrans;
	else return enif_make_badarg(env);

	cblas_dgemm(order, ta,tb, m,n,k, alpha, a,lda, b,ldb, beta, c,ldc);
	break;}
    case 1: {
	if(enif_is_identical(argv[1], atom_left)) side = CblasLeft;
	else if(enif_is_identical(argv[1], atom_right)) side = CblasRight;
	else return enif_make_badarg(env);

	if(enif_is_identical(argv[2], atom_lower)) uplo = CblasLower;
	else if(enif_is_identical(argv[2], atom_upper)) uplo = CblasUpper;
	else return enif_make_badarg(env);

	cblas_dsymm(order, side, uplo, m,n, alpha, a,lda, b,ldb, beta, c,ldc);
	break; }
    case 2: {
	if(enif_is_identical(argv[1], atom_lower)) uplo = CblasLower;
	else if(enif_is_identical(argv[1], atom_upper)) uplo = CblasUpper;
	else return enif_make_badarg(env);

	if(enif_is_identical(argv[2], atom_notransp)) tb = CblasNoTrans;
	else if(enif_is_identical(argv[2], atom_transpose)) tb = CblasTrans;
	else if(enif_is_identical(argv[2], atom_conjugatet)) tb = CblasConjTrans;

	// fprintf(stderr, "%d %d,%d,%d %.2f %.2f %d,%d,%d\r\n", __LINE__, m,n,k, alpha, beta, lda,ldb,ldc);
	cblas_dsyr2k(order, uplo, tb, n, k, alpha, a,lda, b,ldb, beta, c,ldc);
	break; }
    case 3: {
	if(enif_is_identical(argv[1], atom_lower)) uplo = CblasLower;
	else if(enif_is_identical(argv[1], atom_upper)) uplo = CblasUpper;
	else return enif_make_badarg(env);

	if(enif_is_identical(argv[2], atom_notransp)) tb = CblasNoTrans;
	else if(enif_is_identical(argv[2], atom_transpose)) tb = CblasTrans;
	else if(enif_is_identical(argv[2], atom_conjugatet)) tb = CblasConjTrans;

	if(enif_is_identical(argv[13], atom_left)) side = CblasLeft;
	else if(enif_is_identical(argv[13], atom_right)) side = CblasRight;
	else return enif_make_badarg(env);

	if(enif_is_identical(argv[14], atom_nonunit)) diag = CblasNonUnit;
	else if(enif_is_identical(argv[14], atom_unit)) diag = CblasUnit;
	else return enif_make_badarg(env);

	cblas_dtrmm(order,side,uplo,tb,diag, m,n, alpha, a,lda, b, ldb);
	break; }
    case 4: {
	if(enif_is_identical(argv[1], atom_lower)) uplo = CblasLower;
	else if(enif_is_identical(argv[1], atom_upper)) uplo = CblasUpper;
	else return enif_make_badarg(env);

	if(enif_is_identical(argv[2], atom_notransp)) tb = CblasNoTrans;
	else if(enif_is_identical(argv[2], atom_transpose)) tb = CblasTrans;
	else if(enif_is_identical(argv[2], atom_conjugatet)) tb = CblasConjTrans;

	if(enif_is_identical(argv[13], atom_left)) side = CblasLeft;
	else if(enif_is_identical(argv[13], atom_right)) side = CblasRight;
	else return enif_make_badarg(env);

	if(enif_is_identical(argv[14], atom_nonunit)) diag = CblasNonUnit;
	else if(enif_is_identical(argv[14], atom_unit)) diag = CblasUnit;
	else return enif_make_badarg(env);

	cblas_dtrsm(order,side,uplo,tb,diag, m,n, alpha, a,lda, b, ldb);
	break; }
    case 5: {
	if(enif_is_identical(argv[1], atom_lower)) uplo = CblasLower;
	else if(enif_is_identical(argv[1], atom_upper)) uplo = CblasUpper;
	else return enif_make_badarg(env);

	if(enif_is_identical(argv[2], atom_notransp)) tb = CblasNoTrans;
	else if(enif_is_identical(argv[2], atom_transpose)) tb = CblasTrans;
	else if(enif_is_identical(argv[2], atom_conjugatet)) tb = CblasConjTrans;

	cblas_dsyrk(order,uplo,tb, n,k, alpha, a,lda, beta, b,ldb);
	break; }
    }

    enif_consume_timeslice(env, MAX(m,k)*n/1000);
    return atom_ok;
}

/* ---------------------------------------------------*/

static int load(ErlNifEnv* env, void** priv_data, ERL_NIF_TERM load_info)
{
    atom_ok = enif_make_atom(env,"ok");
    atom_true = enif_make_atom(env,"true");

    atom_notransp = enif_make_atom(env,"no_transp");
    atom_transpose = enif_make_atom(env,"transp");
    atom_conjugatet = enif_make_atom(env,"conj_transp");

    atom_rowmaj = enif_make_atom(env,"row_maj");
    atom_colmaj = enif_make_atom(env,"col_maj");

    atom_upper = enif_make_atom(env,"upper");
    atom_lower = enif_make_atom(env,"lower");

    atom_left = enif_make_atom(env,"left");
    atom_right = enif_make_atom(env,"right");

    atom_nonunit = enif_make_atom(env,"non_unit");
    atom_unit = enif_make_atom(env,"unit");

    avec_r  = enif_open_resource_type(env, "eblas", "avec", NULL, ERL_NIF_RT_CREATE, NULL);
    return 0;
}

static int upgrade(ErlNifEnv* env, void** priv_data, void** old_priv_data,
		   ERL_NIF_TERM load_info)
{
    return 0;
}

static void unload(ErlNifEnv* env, void* priv_data)
{

}

ERL_NIF_INIT(blasd_raw,nif_funcs,load,NULL,upgrade,unload)

ext: c


text:
/* cheb/test.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.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include "gsl_chebyshev.h"

double f_sin (double x, void * p) {
  p = 0;
  return sin(x);
}

int
main(void)
{
  double tol = 100.0 * GSL_DBL_EPSILON;
  double x;

  gsl_cheb_series * cs = gsl_cheb_alloc(40);
  gsl_cheb_series * csd = gsl_cheb_alloc(40);
  gsl_cheb_series * csi = gsl_cheb_alloc(40);

  gsl_function F_sin;

  F_sin.function = f_sin;
  F_sin.params = 0;

  gsl_ieee_env_setup();

  gsl_cheb_init(cs, &F_sin, -M_PI, M_PI);

  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r = gsl_cheb_eval(cs, x);
    gsl_test_abs(r, sin(x), tol, "gsl_cheb_eval, sin(%.3g)", x);
  }

  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r, e;
    gsl_cheb_eval_err(cs, x, &r, &e);
    gsl_test_abs(r, sin(x), tol, "gsl_cheb_eval_err, sin(%.3g)", x);
    gsl_test_factor(fabs(r-sin(x)) + GSL_DBL_EPSILON, e, 10.0,
                    "gsl_cheb_eval_err, error sin(%.3g)", x);
  }

  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r = gsl_cheb_eval_n(cs, 25, x);
    gsl_test_abs(r, sin(x), tol, "gsl_cheb_eval_n, sin(%.3g)", x);
  }

  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r, e;
    gsl_cheb_eval_n_err(cs, 25, x, &r, &e);
    gsl_test_abs(r, sin(x), 100.0 * tol, "gsl_cheb_eval_n_err, deriv sin(%.3g)", x);
    gsl_test_factor(fabs(r-sin(x)) + GSL_DBL_EPSILON, e, 10.0,
                    "gsl_cheb_eval_n_err, error sin(%.3g)", x);
  }

  /* Test derivative */

  gsl_cheb_calc_deriv(csd, cs);

  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r = gsl_cheb_eval(csd, x);
    gsl_test_abs(r, cos(x), 1600 * tol, "gsl_cheb_eval, deriv sin(%.3g)", x);
  }

#ifdef TEST_DERIVATIVE_ERR
  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r, e;
    gsl_cheb_eval_err(csd, x, &r, &e);
    gsl_test_abs(r, cos(x), tol, "gsl_cheb_eval_err, deriv sin(%.3g)", x);
    gsl_test_factor(fabs(r-cos(x)) + GSL_DBL_EPSILON, e, 10.0,
                    "gsl_cheb_eval_err, deriv error sin(%.3g)", x);
  }
#endif

  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r = gsl_cheb_eval_n(csd, 25, x);
    gsl_test_abs(r, cos(x), 1600 * tol, "gsl_cheb_eval_n, deriv sin(%.3g)", x);
  }

#ifdef TEST_DERIVATIVE_ERR
  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r, e;
    gsl_cheb_eval_n_err(csd, 25, x, &r, &e);
    gsl_test_abs(r, cos(x), 100.0 * tol, "gsl_cheb_eval_n_err, deriv sin(%.3g)", x);
    gsl_test_factor(fabs(r-cos(x)) + GSL_DBL_EPSILON, e, 10.0,
                    "gsl_cheb_eval_n_err, deriv error sin(%.3g)", x);
  }
#endif

  /* Test integral */

  gsl_cheb_calc_integ(csi, cs);

  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r = gsl_cheb_eval(csi, x);
    gsl_test_abs(r, -(1+cos(x)), tol, "gsl_cheb_eval, integ sin(%.3g)", x);
  }

#ifdef TEST_INTEGRAL_ERR
  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r, e;
    gsl_cheb_eval_err(csi, x, &r, &e);
    gsl_test_abs(r, -(1+cos(x)), tol, "gsl_cheb_eval_err, integ sin(%.3g)", x);
    gsl_test_factor(fabs(r-(-1-cos(x))) + GSL_DBL_EPSILON, e, 10.0,
                    "gsl_cheb_eval_err, integ error sin(%.3g)", x);
  }
#endif

  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r = gsl_cheb_eval_n(csi, 25, x);
    gsl_test_abs(r, -(1+cos(x)), tol, "gsl_cheb_eval_n, integ sin(%.3g)", x);
  }

#ifdef TEST_INTEGRAL_ERR
  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r, e;
    gsl_cheb_eval_n_err(csi, 25, x, &r, &e);
    gsl_test_abs(r, -(1+cos(x)), 100.0 * tol, "gsl_cheb_eval_n_err, integ sin(%.3g)", x);
    gsl_test_factor(fabs(r-(-1-cos(x))) + GSL_DBL_EPSILON, e, 10.0,
                    "gsl_cheb_eval_n_err, integ error sin(%.3g)", x);
  }
#endif

  gsl_cheb_free(csi);
  gsl_cheb_free(csd);
  gsl_cheb_free(cs);

  exit (gsl_test_summary());
}


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>

double f  (double x,void *params);
double f  (double x,void *params)
{
        double n=* (double *) params;
        double f=exp (-x) * pow (x,n);
        return f;
}

void integral_gen (int nmin, int nmax, double vals[]);
void integral_gen (int nmin, int nmax, double vals[])
{
	double a=0.,b=1;
        double abserr=1.e-9,relerr=1.e-9;
        double result;
        double error;
        double n=10;
        size_t np=1000;
	 gsl_integration_workspace *w = gsl_integration_workspace_alloc (np);
                gsl_function F;
                F.function = &f;

	for (int i=nmin;i<=nmax;i++)
	{
        	F.params =&n;
        	gsl_integration_qag (&F, a, b, abserr, relerr,np,GSL_INTEG_GAUSS15, w, &result, &error);
		vals[i]=result;
	}
	gsl_integration_workspace_free (w);
}


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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_chebyshev.h>

double f_T0 (double x, void * p) {
  p = 0;
  return 1.0;
}

double f_T1 (double x, void * p) {
  p = 0;
  return x;
}

double f_T2 (double x, void * p) {
  p = 0;
  return 2*x*x - 1;
}

double f_sin (double x, void * p) {
  p = 0;
  return sin(x);
}

int
main(void)
{
  double tol = 100.0 * GSL_DBL_EPSILON;
  double ftol = 20.0;
  double x;
  size_t i;

  gsl_cheb_series * cs = gsl_cheb_alloc(40);
  gsl_cheb_series * csd = gsl_cheb_alloc(40);
  gsl_cheb_series * csi = gsl_cheb_alloc(40);

  gsl_function F_sin, F_T0, F_T1, F_T2;

  F_sin.function = f_sin;
  F_sin.params = 0;

  F_T0.function = f_T0;
  F_T0.params = 0;

  F_T1.function = f_T1;
  F_T1.params = 0;

  F_T2.function = f_T2;
  F_T2.params = 0;

  gsl_ieee_env_setup();

  gsl_cheb_init(cs, &F_T0, -1.0, 1.0);

  {
    size_t expected = 40;
    size_t order = gsl_cheb_order (cs);
    size_t size = gsl_cheb_size (cs);
    double * p = gsl_cheb_coeffs (cs);
    gsl_test(order != expected, "gsl_cheb_order");
    gsl_test(size != expected + 1, "gsl_cheb_size");
    gsl_test(p != cs->c, "gsl_cheb_coeffs");
  }

  for (i = 0; i<cs->order; i++)
    {
      double c_exp = (i == 0) ? 2.0 : 0.0;
      gsl_test_abs (cs->c[i], c_exp, tol, "c[%d] for T_0(x)", i);
    }

  gsl_cheb_init(cs, &F_T1, -1.0, 1.0);

  for (i = 0; i<cs->order; i++)
    {
      double c_exp = (i == 1) ? 1.0 : 0.0;
      gsl_test_abs (cs->c[i], c_exp, tol, "c[%d] for T_1(x)", i);
    }

  gsl_cheb_init(cs, &F_T2, -1.0, 1.0);

  for (i = 0; i<cs->order; i++)
    {
      double c_exp = (i == 2) ? 1.0 : 0.0;
      gsl_test_abs (cs->c[i], c_exp, tol, "c[%d] for T_2(x)", i);
    }

  gsl_cheb_init(cs, &F_sin, -M_PI, M_PI);

  gsl_test_abs (cs->c[0], 0.0, tol, "c[0] for F_sin(x)");
  gsl_test_abs (cs->c[1], 5.69230686359506e-01, tol, "c[1] for F_sin(x)");
  gsl_test_abs (cs->c[2], 0.0, tol, "c[2] for F_sin(x)");
  gsl_test_abs (cs->c[3], -6.66916672405979e-01, tol, "c[3] for F_sin(x)");
  gsl_test_abs (cs->c[4], 0.0, tol, "c[4] for F_sin(x)");
  gsl_test_abs (cs->c[5], 1.04282368734237e-01, tol, "c[5] for F_sin(x)");

  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r = gsl_cheb_eval(cs, x);
    gsl_test_abs(r, sin(x), tol, "gsl_cheb_eval, sin(%.3g)", x);
  }

  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r, e;
    gsl_cheb_eval_err(cs, x, &r, &e);
    gsl_test_abs(r, sin(x), tol, "gsl_cheb_eval_err, sin(%.3g)", x);
    gsl_test_factor(fabs(r-sin(x)) + GSL_DBL_EPSILON, e, ftol,
                    "gsl_cheb_eval_err, error sin(%.3g)", x);
  }

  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r = gsl_cheb_eval_n(cs, 25, x);
    gsl_test_abs(r, sin(x), tol, "gsl_cheb_eval_n, sin(%.3g)", x);
  }

  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r, e;
    gsl_cheb_eval_n_err(cs, 25, x, &r, &e);
    gsl_test_abs(r, sin(x), 100.0 * tol, "gsl_cheb_eval_n_err, deriv sin(%.3g)", x);
    gsl_test_factor(fabs(r-sin(x)) + GSL_DBL_EPSILON, e, ftol,
                    "gsl_cheb_eval_n_err, error sin(%.3g)", x);
  }

  /* Test derivative */

  gsl_cheb_calc_deriv(csd, cs);

  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r = gsl_cheb_eval(csd, x);
    gsl_test_abs(r, cos(x), 1600 * tol, "gsl_cheb_eval, deriv sin(%.3g)", x);
  }

#ifdef TEST_DERIVATIVE_ERR
  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r, e;
    gsl_cheb_eval_err(csd, x, &r, &e);
    gsl_test_abs(r, cos(x), tol, "gsl_cheb_eval_err, deriv sin(%.3g)", x);
    gsl_test_factor(fabs(r-cos(x)) + GSL_DBL_EPSILON, e, ftol,
                    "gsl_cheb_eval_err, deriv error sin(%.3g)", x);
  }
#endif

  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r = gsl_cheb_eval_n(csd, 25, x);
    gsl_test_abs(r, cos(x), 1600 * tol, "gsl_cheb_eval_n, deriv sin(%.3g)", x);
  }

#ifdef TEST_DERIVATIVE_ERR
  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r, e;
    gsl_cheb_eval_n_err(csd, 25, x, &r, &e);
    gsl_test_abs(r, cos(x), 100.0 * tol, "gsl_cheb_eval_n_err, deriv sin(%.3g)", x);
    gsl_test_factor(fabs(r-cos(x)) + GSL_DBL_EPSILON, e, ftol,
                    "gsl_cheb_eval_n_err, deriv error sin(%.3g)", x);
  }
#endif

  /* Test integral */

  gsl_cheb_calc_integ(csi, cs);

  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r = gsl_cheb_eval(csi, x);
    gsl_test_abs(r, -(1+cos(x)), tol, "gsl_cheb_eval, integ sin(%.3g)", x);
  }

#ifdef TEST_INTEGRAL_ERR
  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r, e;
    gsl_cheb_eval_err(csi, x, &r, &e);
    gsl_test_abs(r, -(1+cos(x)), tol, "gsl_cheb_eval_err, integ sin(%.3g)", x);
    gsl_test_factor(fabs(r-(-1-cos(x))) + GSL_DBL_EPSILON, e, ftol,
                    "gsl_cheb_eval_err, integ error sin(%.3g)", x);
  }
#endif

  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r = gsl_cheb_eval_n(csi, 25, x);
    gsl_test_abs(r, -(1+cos(x)), tol, "gsl_cheb_eval_n, integ sin(%.3g)", x);
  }

#ifdef TEST_INTEGRAL_ERR
  for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
    double r, e;
    gsl_cheb_eval_n_err(csi, 25, x, &r, &e);
    gsl_test_abs(r, -(1+cos(x)), 100.0 * tol, "gsl_cheb_eval_n_err, integ sin(%.3g)", x);
    gsl_test_factor(fabs(r-(-1-cos(x))) + GSL_DBL_EPSILON, e, ftol,
                    "gsl_cheb_eval_n_err, integ error sin(%.3g)", x);
  }
#endif

  gsl_cheb_free(csi);
  gsl_cheb_free(csd);
  gsl_cheb_free(cs);

  exit (gsl_test_summary());
}


text:
#include <petsc.h>
#include <petsc/private/petscimpl.h>

#if PETSC_VERSION_LT(3,17,0)
#undef  SETERRQ
#define SETERRQ(comm,ierr,...) return PetscError(comm,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr,PETSC_ERROR_INITIAL,__VA_ARGS__)
#endif

PETSC_EXTERN PetscErrorCode IGA_Partition(PetscInt,PetscInt,PetscInt,const PetscInt[],PetscInt[],PetscInt[]);
PETSC_EXTERN PetscErrorCode IGA_Distribute(PetscInt,const PetscInt[],const PetscInt[],const PetscInt[],PetscInt[],PetscInt[]);

static inline
PetscInt IGA_CUT2D(PetscInt M,PetscInt N,
                   PetscInt m,PetscInt n)
{return M*(n-1) + N*(m-1);}

static inline
PetscInt IGA_CUT3D(PetscInt M,PetscInt N,PetscInt P,
                   PetscInt m,PetscInt n,PetscInt p)
{return (N*P*(m-1) + M*P*(n-1) + M*N*(p-1));}

static inline
PetscInt IGA_PART2D_INNER(PetscInt size,
                          PetscInt   M,PetscInt   N,
                          PetscInt *_m,PetscInt *_n)
{
  PetscInt m,n;
  m = (PetscInt)(0.5 + sqrt(((double)M)/((double)N)*((double)size)));
  if (m == 0) {m = 1;} while (m > 0 && size % m) m--;
  n = size / m;
  *_m = m; *_n = n;
  return IGA_CUT2D(M,N,m,n);
}

static inline
void IGA_PART2D(PetscInt size,
                PetscInt   M,PetscInt   N,
                PetscInt *_m,PetscInt *_n)
{
  PetscInt m,n;
  PetscInt m1,n1,a;
  PetscInt m2,n2,b;
  a = IGA_PART2D_INNER(size,M,N,&m1,&n1);
  b = IGA_PART2D_INNER(size,N,M,&n2,&m2);
  if (a<b) {m = m1; n = n1;}
  else     {m = m2; n = n2;}
  if (M == N && n < m) {PetscInt t = m; m = n; n = t;}
  *_m = m; *_n = n;
}

static inline
PetscInt IGA_PART3D_INNER(PetscInt size,
                          PetscInt   M,PetscInt   N,PetscInt   P,
                          PetscInt *_m,PetscInt *_n,PetscInt *_p)
{
  PetscInt m,n,p,C;
  PetscInt mm,nn,pp,CC;
  /**/
  m = (PetscInt)(0.5 + pow(((double)M*(double)M)/((double)N*(double)P)*(double)size,1./3.));
  if (m == 0) {m = 1;} while (m > 0 && size % m) m--;
  /**/
  IGA_PART2D(size/m,N,P,&n,&p);
  C = IGA_CUT3D(M,N,P,m,n,p);
  /**/
  for (mm=m; mm>=1; mm--) {
    if (size % mm) continue;
    IGA_PART2D(size/mm,N,P,&nn,&pp);
    CC = IGA_CUT3D(M,N,P,mm,nn,pp);
    if (CC < C) {m = mm; n = nn; p = pp; C = CC;}
  }
  /**/
  for (nn=n; nn>=1; nn--) {
    if (size % nn) continue;
    IGA_PART2D(size/nn,M,P,&mm,&pp);
    CC = IGA_CUT3D(M,N,P,mm,nn,pp);
    if (CC < C) {m = mm; n = nn; p = pp; C = CC;}
  }
  /**/
  for (pp=p; pp>=1; pp--) {
    if (size % pp) continue;
    IGA_PART2D(size/pp,M,N,&mm,&nn);
    CC = IGA_CUT3D(M,N,P,mm,nn,pp);
    if (CC < C) {m = mm; n = nn; p = pp; C = CC;}
  }
  /**/
  *_m = m; *_n = n; *_p = p;
  return IGA_CUT3D(M,N,P,m,n,p);
}

static inline
void IGA_PART3D(PetscInt size,
                PetscInt   M,PetscInt   N,PetscInt   P,
                PetscInt *_m,PetscInt *_n,PetscInt *_p)
{
  PetscInt m[3],n[3],p[3],C[3],k,i=0,Cmin=PETSC_MAX_INT,t;
  C[0] = IGA_PART3D_INNER(size,M,N,P,&m[0],&n[0],&p[0]);
  C[1] = IGA_PART3D_INNER(size,N,M,P,&n[1],&m[1],&p[1]);
  C[2] = IGA_PART3D_INNER(size,P,M,N,&p[2],&m[2],&n[2]);
  for (k=0; k<3; k++) if (C[k]<Cmin) {Cmin=C[k]; i=k;}
  if (M == N && n[i] < m[i]) {t = m[i]; m[i] = n[i]; n[i] = t;}
  if (M == P && p[i] < m[i]) {t = m[i]; m[i] = p[i]; p[i] = t;}
  if (N == P && p[i] < n[i]) {t = n[i]; n[i] = p[i]; p[i] = t;}
  *_m = m[i]; *_n = n[i]; *_p = p[i];
}

static inline
void IGA_Part2D(PetscInt size,
                PetscInt   M,PetscInt   N,
                PetscInt *_m,PetscInt *_n)
{
  PetscInt m,n;
  m = *_m; n = *_n;
  if (m < 1 && n < 1) IGA_PART2D(size,M,N,&m,&n);
  else if (m < 1)     m = size / n;
  else if (n < 1)     n = size / m;
  *_m = m; *_n = n;
}

static inline
void IGA_Part3D(PetscInt size,
                PetscInt   M,PetscInt   N,PetscInt   P,
                PetscInt *_m,PetscInt *_n,PetscInt *_p)
{
  PetscInt m,n,p;
  m = *_m; n = *_n; p = *_p;
  if (m < 1 && n < 1 && p < 1) IGA_PART3D(size,M,N,P,&m,&n,&p);
  else if (m < 1 && n < 1)     IGA_PART2D(size/p,M,N,&m,&n);
  else if (m < 1 && p < 1)     IGA_PART2D(size/n,M,P,&m,&p);
  else if (n < 1 && p < 1)     IGA_PART2D(size/m,N,P,&n,&p);
  else if (m < 1)              m = size / (n*p);
  else if (n < 1)              n = size / (m*p);
  else if (p < 1)              p = size / (m*n);
  *_m = m; *_n = n; *_p = p;
}

PetscErrorCode IGA_Partition(PetscInt size,PetscInt rank,
                             PetscInt dim,const PetscInt N[],
                             PetscInt n[],PetscInt i[])
{
  PetscInt k,p=1,m[3]={1,1,1};
  PetscFunctionBegin;
  PetscValidIntPointer(N,4);
  PetscValidIntPointer(n,5);
  if (i) PetscValidIntPointer(i,6);
  if (size < 1)
    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,
             "Number of partitions %D must be positive",size);
  if (i && (rank < 0 || rank >= size))
    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Partition index %D must be in range [0,%D]",rank,size-1);

  switch (dim) {
  case 3:  IGA_Part3D(size,N[0],N[1],N[2],&n[0],&n[1],&n[2]); break;
  case 2:  IGA_Part2D(size,N[0],N[1],&n[0],&n[1]); break;
  case 1:  if (n[0] < 1) n[0] = size; break;
  default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of dimensions %D must be in range [1,3]",dim);
  }
  for (k=0; k<dim; k++) p *= (m[k] = n[k]);
  if (p != size) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Bad partition, prod(%D,%D,%D) != %D",m[0],m[1],m[2],size);
  for (k=0; k<dim; k++)
    if (N[k] < n[k]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Partition %D is too fine, %D elements in %D parts",k,N[k],n[k]);
  if (i)
    for (k=0; k<dim; k++) {
      i[k] = rank % n[k];
      rank -= i[k];
      rank /= n[k];
    }
  PetscFunctionReturn(0);
}

static inline
void IGA_Dist1D(PetscInt size,PetscInt rank,
                PetscInt N,PetscInt *n,PetscInt *s)
{
  *n = N/size + ((N % size) > rank);
  *s = rank * (N/size) + (((N % size) > rank) ? rank : (N % size));
}

PetscErrorCode IGA_Distribute(PetscInt dim,
                              const PetscInt size[],const PetscInt rank[],
                              const PetscInt N[],PetscInt n[],PetscInt s[])
{
  PetscInt k;
  PetscFunctionBegin;
  PetscValidIntPointer(size,2);
  PetscValidIntPointer(rank,3);
  PetscValidIntPointer(N,4);
  PetscValidIntPointer(n,5);
  PetscValidIntPointer(s,6);
  if (dim < 1 || dim > 3)
    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of dimensions %D must be in range [1,3]",dim);
  for (k=0; k<dim; k++) {
    if (size[k] < 1)
      SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of partitions %D must be positive",size[k]);
    if (rank[k] < 0 || rank[k] >= size[k])
      SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Partition index %D must be in range [0,%D]",rank[k],size[k]-1);
    if (N[k] < 0)
      SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of elements %D must be non-negative",N[k]);
  }
  for (k=0; k<dim; k++)
    IGA_Dist1D(size[k],rank[k],N[k],&n[k],&s[k]);
  PetscFunctionReturn(0);
}

max_stars_repo_name: dalcinl/PetIGA

max_issues_repo_name: dalcinl/PetIGA

max_forks_repo_name: dalcinl/PetIGA


text:
#pragma once

#include "rev/Utilities.h"
#include "rev/geometry/Tools.h"
#include <array>
#include <glm/glm.hpp>
#include <gsl/gsl_assert>
#include <iostream>
#include <limits>
#include <set>
#include <unordered_set>
#include <variant>
#include <vector>

namespace rev {

template <typename SurfaceData>
struct Triangle {
    std::array<glm::vec3, 3> vertices;
    SurfaceData data;

    AxisAlignedBoundingBox getBoundingBox() const
    {
        AxisAlignedBoundingBox box;
        for (const auto& vertex : vertices) {
            box.expandToVertex(vertex);
        }
        return box;
    }

    glm::vec3 getNormal() const
    {
        auto e1 = vertices[1] - vertices[0];
        auto e2 = vertices[2] - vertices[0];
        return glm::normalize(glm::cross(e1, e2));
    }

    bool isLeftOfPlane(const AxisAlignedPlane& plane) const
    {
        for (const auto& vertex : vertices) {
            if (vertex[plane.dimensionIndex] < plane.boundary) {
                return true;
            }
        }
        return false;
    }

    bool isRightOfPlane(const AxisAlignedPlane& plane) const
    {
        for (const auto& vertex : vertices) {
            if (vertex[plane.dimensionIndex] > plane.boundary) {
                return true;
            }
        }
        return false;
    }

    std::optional<AxisAlignedBoundingBox> clippedBoundingBox(const AxisAlignedBoundingBox& box)
    {
        // Alternate between these two buffers for input and output
        std::array<glm::vec3, 9> buf1{};
        std::array<glm::vec3, 9> buf2{};

        // Fill the input buffer
        for (int i = 0; i < 3; i++) {
            buf1[i] = vertices[i];
        }

        // Clip the triangle by all sides of our bounding box
        auto inputIter = buf1.data();
        auto inputEnd = inputIter + 3;
        auto outputIter = buf2.data();

        // First the min planes in each dimension
        for (uint8_t dim = 0; dim < 3; dim++) {
            gsl::span<glm::vec3> inputRange(inputIter, inputEnd);
            OutputIteratorPolygonBuilder builder(outputIter);
            AxisAlignedPlane boundary = { dim, box.minimum[dim] };
            boundary.splitConvexPolygon(inputRange, NullPolygonBuilder(), builder);

            if (std::distance(outputIter, builder.getIterator()) < 3) {
                return {};
            }

            std::swap(inputIter, outputIter);
            inputEnd = builder.getIterator();
        }

        // Now the max planes in each dimension
        for (uint8_t dim = 0; dim < 3; dim++) {
            gsl::span<glm::vec3> inputRange(inputIter, inputEnd);
            OutputIteratorPolygonBuilder builder(outputIter);
            AxisAlignedPlane boundary = { dim, box.maximum[dim] };
            boundary.splitConvexPolygon(inputRange, builder, NullPolygonBuilder());

            if (std::distance(outputIter, builder.getIterator()) < 3) {
                return {};
            }

            std::swap(inputIter, outputIter);
            inputEnd = builder.getIterator();
        }

        gsl::span<glm::vec3> inputRange(inputIter, inputEnd);
        auto clippedBox = smallestBoxContainingVertices(inputRange);
        clippedBox.reduceToBox(box);
        return clippedBox;
    }

    struct Hit {
        glm::vec2 uv;
        float t;
    };

    std::optional<Hit> castRay(const Ray& ray) const
    {
        glm::vec3 edge1 = vertices[1] - vertices[0];
        glm::vec3 edge2 = vertices[2] - vertices[0];

        glm::vec3 p = glm::cross(ray.direction, edge2);
        float determinant = glm::dot(edge1, p);
        if (abs(determinant) < std::numeric_limits<float>::epsilon()) {
            // The ray is parallel.
            return std::nullopt;
        }

        glm::vec3 fromV0 = ray.origin - vertices[0];
        float u = glm::dot(fromV0, p) / determinant;
        if (u < 0.0f || u > 1.0f) {
            return std::nullopt;
        }

        glm::vec3 q = glm::cross(fromV0, edge1);
        float v = glm::dot(ray.direction, q) / determinant;
        if (v < 0.0f || v > 1.0f) {
            return std::nullopt;
        }

        float t = glm::dot(edge2, q) / determinant;
        if (t < 0.0f) {
            return std::nullopt;
        }

        return Hit{ { u, v }, t };
    }

    glm::vec3 baryCentricToCartesian(glm::vec2 uv)
    {
        float u = uv[0];
        float v = uv[1];
        float w = 1.0f - (u + v);
        return w * vertices[0] + u * vertices[1] + v * vertices[2];
    }

    glm::vec2 cartesianToBarycentric(glm::vec3 position)
    {
        glm::vec3 v0 = vertices[1] - vertices[0];
        glm::vec3 v1 = vertices[2] - vertices[0];
        glm::vec3 v2 = position - vertices[0];

        float d00 = glm::dot(v0, v0);
        float d01 = glm::dot(v0, v1);
        float d11 = glm::dot(v1, v1);
        float d20 = glm::dot(v2, v0);
        float d21 = glm::dot(v2, v1);
        float denom = d00 * d11 - d01 * d01;

        float v = (d11 * d20 - d01 * d21) / denom;
        float w = (d00 * d21 - d01 * d20) / denom;
        float u = 1.0f - v - w;

        return { u, v };
    }

    std::optional<float> intersectsSphere(const Sphere& sphere)
    {
        glm::vec3 normal = getNormal();
        float t = glm::dot(normal, (vertices[0] - sphere.center));
        if (sphere.radius < abs(t)) {
            // Closest point to the triangle is outside the sphere.
            return std::nullopt;
        }

        glm::vec3 closestPoint = sphere.center - (normal * t);
        glm::vec2 cpBary = cartesianToBarycentric(closestPoint);
        size_t outsideCount = 0;
        glm::bvec3 isOutside{};
        if (cpBary[0] < 0.0f) {
            outsideCount++;
            isOutside[0] = true;
        }

        if (cpBary[1] < 0.0f) {
            outsideCount++;
            isOutside[1] = true;
        }

        if ((1.0f - (cpBary[0] + cpBary[1])) < 0.0f) {
            outsideCount++;
            isOutside[2] = true;
        }

        if (outsideCount == 0) {
            // Closest point is inside the triangle and inside the sphere.
            return t;
        }

        // Calculate the square of the radius of the cross section of the sphere passing through
        // the triangle's plane.
        float r2 = (sphere.radius * sphere.radius) - (t * t);

        if (outsideCount == 1) {
            // We need to find the distance to the edge (squared) here.
            glm::vec3 edge;
            glm::vec3 toVertex;
            if (isOutside[0]) {
                edge = vertices[2] - vertices[1];
                toVertex = closestPoint - vertices[1];
            } else if (isOutside[1]) {
                edge = vertices[2] - vertices[0];
                toVertex = closestPoint - vertices[0];
            } else {
                edge = vertices[1] - vertices[0];
                toVertex = closestPoint - vertices[0];
            }

            float toVertexLength = glm::length(toVertex);
            float legLength = glm::dot(edge, toVertex) / toVertexLength;

            return ((toVertexLength * toVertexLength) - (legLength * legLength)) < r2;
        }

        // Otherwise, one of the vertices is actually the closest point on the triangle.
        glm::vec3 toVertex;
        if (!isOutside[0]) {
            toVertex = closestPoint - vertices[0];
        } else if (!isOutside[1]) {
            toVertex = closestPoint - vertices[1];
        } else {
            toVertex = closestPoint - vertices[2];
        }

        return glm::dot(toVertex, toVertex) < r2;
    }
};

template <typename SurfaceData>
struct LeafNode {
    std::unordered_set<Triangle<SurfaceData>*> triangles;
};

template <typename SurfaceData>
struct BranchNode;

template <typename SurfaceData>
using MapNode = std::variant<LeafNode<SurfaceData>, BranchNode<SurfaceData>>;

template <typename SurfaceData>
struct BranchNode {
    AxisAlignedPlane split;

    std::unique_ptr<MapNode<SurfaceData>> left;
    std::unique_ptr<MapNode<SurfaceData>> right;
};

template <typename SurfaceData>
class KDTree {
public:
    KDTree(std::vector<std::unique_ptr<Triangle<SurfaceData>>> triangles,
        std::unique_ptr<MapNode<SurfaceData>> root, const AxisAlignedBoundingBox boundingBox)
        : _triangles(std::move(triangles))
        , _root(std::move(root))
        , _boundingBox(boundingBox)
    {
    }

    struct Hit {
        Triangle<SurfaceData>* triangle;
        glm::vec2 uv;
        float t;
    };

    std::optional<Hit> castRay(
        const Ray& ray, float maxDistance = std::numeric_limits<float>::infinity()) const
    {
        if (_boundingBox.containsPoint(ray.origin)) {
            AxisAlignedBoundingBox::Hit entryPoint{ ray.origin, 0.0f, 0 };
            return getHit(*_root, ray, _boundingBox, entryPoint, maxDistance);
        } else {
            auto hit = _boundingBox.castExternalRay(ray);
            if (!hit || (hit->t > maxDistance)) {
                return std::nullopt;
            }
            return getHit(*_root, ray, _boundingBox, *hit, maxDistance);
        }
    }

    template <typename Visitor>
    void visitTrianglesIntersectingSphere(const Sphere& sphere, Visitor&& visitor) const
    {
        std::visit(
            [this, &sphere, &visitor](const auto& node) {
                visitTrianglesIntersectingSphereInNode(
                    sphere, std::forward<Visitor>(visitor), _boundingBox, node);
            },
            *_root);
    }

    void dump() const
    {
        std::cout << "[KDTree]{" << std::endl;
        printBoundingBox(_boundingBox, 2);
        std::visit([this](const auto& root) { printNode(_boundingBox, root, 2); }, *_root);
        std::cout << "}" << std::endl;
    }

private:
    std::optional<Hit> getHit(const MapNode<SurfaceData>& node, const Ray& ray,
        const AxisAlignedBoundingBox& box, AxisAlignedBoundingBox::Hit& entryPoint,
        float maxDistance) const
    {
        return std::visit(
            [this, &ray, &box, &entryPoint, maxDistance](
                const auto& node) { return getHit(node, ray, box, entryPoint, maxDistance); },
            node);
    }

    std::optional<Hit> getHit(const BranchNode<SurfaceData>& node, const Ray& ray,
        const AxisAlignedBoundingBox& box, AxisAlignedBoundingBox::Hit& entryPoint,
        float maxDistance) const
    {
        auto [leftBox, rightBox] = box.split(node.split);
        float position = entryPoint.intersectionPoint[node.split.dimensionIndex];
        bool goLeft;
        if (position < node.split.boundary) {
            goLeft = true;
        } else if (position > node.split.boundary) {
            goLeft = false;
        } else if (ray.direction[node.split.dimensionIndex] < 0.0f) {
            goLeft = false;
        } else if (ray.direction[node.split.dimensionIndex] > 0.0f) {
            goLeft = true;
        } else {
            return std::nullopt;
        }

        MapNode<SurfaceData>* childNode = goLeft ? node.left.get() : node.right.get();
        auto hit = getHit(*childNode, ray, (goLeft ? leftBox : rightBox), entryPoint, maxDistance);
        if (hit) {
            return hit;
        }
        if (entryPoint.t > maxDistance) {
            return std::nullopt;
        }
        if (entryPoint.planeDimension != node.split.dimensionIndex) {
            return std::nullopt;
        }

        bool shouldTryOtherChild;
        if (goLeft) {
            shouldTryOtherChild = (ray.direction[node.split.dimensionIndex] > 0.0f);
        } else {
            shouldTryOtherChild = (ray.direction[node.split.dimensionIndex] < 0.0f);
        }

        if (shouldTryOtherChild) {
            childNode = goLeft ? node.right.get() : node.left.get();
            return getHit(*childNode, ray, (goLeft ? rightBox : leftBox), entryPoint, maxDistance);
        }
        return std::nullopt;
    }

    std::optional<Hit> getHit(const LeafNode<SurfaceData>& node, const Ray& ray,
        const AxisAlignedBoundingBox& box, AxisAlignedBoundingBox::Hit& entryPoint,
        float maxDistance) const
    {
        float t = std::numeric_limits<float>::infinity();
        std::optional<Hit> bestHit;
        for (const auto& triangle : node.triangles) {
            auto hit = triangle->castRay(ray);
            if (hit) {
                if (!bestHit || (hit->t < bestHit->t)) {
                    bestHit = Hit{ triangle, hit->uv, hit->t };
                }
            }
        }

        if (bestHit && !(bestHit->t > maxDistance)) {
            return bestHit;
        }

        auto boxHit = box.castInternalRay(ray);
        entryPoint = boxHit;

        return std::nullopt;
    }

    template <typename Visitor>
    void visitTrianglesIntersectingSphereInNode(const Sphere& sphere, Visitor&& visitor,
        const AxisAlignedBoundingBox& box, const LeafNode<SurfaceData>& node) const
    {
        if (!box.intersectsSphere(sphere)) {
            return;
        }

        for (const auto& triangle : node.triangles) {
            auto hit = triangle->intersectsSphere(sphere);
            if (hit) {
                visitor(*triangle, *hit);
            }
        }
    }

    template <typename Visitor>
    void visitTrianglesIntersectingSphereInNode(const Sphere& sphere, Visitor&& visitor,
        const AxisAlignedBoundingBox& box, const BranchNode<SurfaceData>& node) const
    {
        if (!box.intersectsSphere(sphere)) {
            return;
        }

        auto [leftBox, rightBox] = box.split(node.split);
        std::visit(
            [this, leftBox = leftBox, &sphere, &visitor](const auto& node) {
                visitTrianglesIntersectingSphereInNode(
                    sphere, std::forward<Visitor>(visitor), leftBox, node);
            },
            *node.left);
        std::visit(
            [this, rightBox = rightBox, &sphere, &visitor](const auto& node) {
                visitTrianglesIntersectingSphereInNode(
                    sphere, std::forward<Visitor>(visitor), rightBox, node);
            },
            *node.right);
    }

    void printIndent(size_t indent = 0) const
    {
        for (size_t i = 0; i < indent; i++) {
            std::cout << " ";
        }
    }

    void printVertex(const glm::vec3& vertex) const
    {
        std::cout << "(" << vertex.x << ", " << vertex.y << ", " << vertex.z << ")";
    }

    void printBoundingBox(const AxisAlignedBoundingBox& box, size_t indent = 0) const
    {
        printIndent(indent);
        std::cout << "[Box]{ Min: ";
        printVertex(box.minimum);
        std::cout << " Max: ";
        printVertex(box.maximum);
        std::cout << " }" << std::endl;
    }

    void printTriangle(const Triangle<SurfaceData>& triangle, size_t indent = 0) const
    {
        printIndent(indent);
        std::cout << "[Triangle]{" << std::endl;
        for (const auto& vertex : triangle.vertices) {
            printIndent(indent + 2);
            printVertex(vertex);
            std::cout << std::endl;
        }
        printIndent(indent);
        std::cout << "}" << std::endl;
    }

    void printNode(const AxisAlignedBoundingBox& box, const LeafNode<SurfaceData>& node,
        size_t indent = 0) const
    {
        printIndent(indent);
        std::cout << "[Leaf]{" << std::endl;
        printBoundingBox(box, indent + 2);
        for (const auto& triangle : node.triangles) {
            printTriangle(*triangle, indent + 2);
        }
        printIndent(indent);
        std::cout << "}" << std::endl;
    }

    void printNode(const AxisAlignedBoundingBox& box, const BranchNode<SurfaceData>& node,
        size_t indent = 0) const
    {
        printIndent(indent);
        std::cout << "[Branch]{" << std::endl;
        printBoundingBox(box, indent + 2);
        auto [leftBox, rightBox] = box.split(node.split);
        std::visit([this, indent, leftBox = leftBox](
                       const auto& childNode) { printNode(leftBox, childNode, indent + 2); },
            *node.left);
        std::visit([this, indent, rightBox = rightBox](
                       const auto& childNode) { printNode(rightBox, childNode, indent + 2); },
            *node.right);
        printIndent(indent);
        std::cout << "}" << std::endl;
    }

    std::vector<std::unique_ptr<Triangle<SurfaceData>>> _triangles;
    std::unique_ptr<MapNode<SurfaceData>> _root;
    AxisAlignedBoundingBox _boundingBox;
};

template <typename SurfaceData>
class KDTreeBuilder {
public:
    void addTriangle(std::array<glm::vec3, 3> vertices, SurfaceData data)
    {
        Expects(!glm::any(glm::isnan(vertices[0])));
        Expects(!glm::any(glm::isnan(vertices[1])));
        Expects(!glm::any(glm::isnan(vertices[2])));

        auto triangle
            = std::make_unique<Triangle<SurfaceData>>(Triangle<SurfaceData>{ vertices, data });
        auto boundingBox = smallestBoxContainingVertices(vertices);
        _boundingBox.expandToBox(boundingBox);

        for (const auto& event : buildTriangleEvents(boundingBox, triangle.get())) {
            _events.insert(event);
        }
        _triangles.push_back(std::move(triangle));
    }

    KDTree<SurfaceData> build()
    {
        Expects(!_triangles.empty());
        Expects(!_events.empty());

        std::unordered_set<Triangle<SurfaceData>*> triangleSet;
        for (const auto& triangle : _triangles) {
            triangleSet.insert(triangle.get());
        }
        auto rootNode = createNode(_boundingBox, std::move(_events), triangleSet);
        return KDTree<SurfaceData>{ std::move(_triangles), std::move(rootNode), _boundingBox };
    }

private:
    static constexpr float kTraversalCost = 15.0f;
    static constexpr float kIntersectionCost = 20.0f;
    static constexpr float kEmptySplitDiscount = 0.8f;

    enum class Side {
        Left,
        Right,
    };

    struct Event {
        enum class Type {
            Ending = 0,
            Planar = 1,
            Starting = 2,
        };

        bool operator<(const Event& other) const
        {
            auto tuplify = [](const Event& event) {
                return std::tie(event.separationPlane.boundary,
                    event.separationPlane.dimensionIndex, event.type, event.triangle);
            };
            return tuplify(*this) < tuplify(other);
        }

        Triangle<SurfaceData>* triangle;
        AxisAlignedPlane separationPlane;
        Type type;
    };

    std::unique_ptr<MapNode<SurfaceData>> createNode(const AxisAlignedBoundingBox& box,
        std::set<Event> events, std::unordered_set<Triangle<SurfaceData>*> triangles)
    {
        size_t triangleCount = triangles.size();
        if (triangleCount > 0) {
            auto [plane, side, cost] = findBestSplit(box, events, triangleCount);
            float terminateCost = static_cast<float>(triangleCount) * kIntersectionCost;
            if (terminateCost > cost) {
                std::unordered_set<Triangle<SurfaceData>*> leftTriangles;
                std::unordered_set<Triangle<SurfaceData>*> rightTriangles;

                for (const auto& event : events) {
                    const auto& eventPlane = event.separationPlane;
                    if (eventPlane.dimensionIndex != plane.dimensionIndex) {
                        continue;
                    }

                    if (eventPlane.boundary < plane.boundary) {
                        if (event.type != Event::Type::Starting) {
                            size_t erased = triangles.erase(event.triangle);
                            Expects(erased == 1);
                            leftTriangles.insert(event.triangle);

                            Expects(event.triangle->isLeftOfPlane(plane));
                        }
                    } else if (eventPlane.boundary > plane.boundary) {
                        if (event.type != Event::Type::Ending) {
                            size_t erased = triangles.erase(event.triangle);
                            Expects(erased == 1);
                            rightTriangles.insert(event.triangle);

                            Expects(event.triangle->isRightOfPlane(plane));
                        }
                    } else {
                        // On the boundary
                        size_t erased = triangles.erase(event.triangle);
                        Expects(erased == 1);
                        switch (event.type) {
                        case Event::Type::Starting:
                            rightTriangles.insert(event.triangle);
                            Expects(event.triangle->isRightOfPlane(plane));
                            break;
                        case Event::Type::Ending:
                            leftTriangles.insert(event.triangle);
                            Expects(event.triangle->isLeftOfPlane(plane));
                            break;
                        case Event::Type::Planar:
                            if (side == Side::Left) {
                                leftTriangles.insert(event.triangle);
                            } else {
                                rightTriangles.insert(event.triangle);
                            }
                            break;
                        }
                    }
                }

                std::set<Event> leftEvents;
                std::set<Event> rightEvents;
                for (const auto& event : events) {
                    if (leftTriangles.count(event.triangle)) {
                        leftEvents.insert(event);
                    } else if (rightTriangles.count(event.triangle)) {
                        rightEvents.insert(event);
                    }
                }

                auto [leftBox, rightBox] = box.split(plane);
                // Overlapping triangles remain
                for (const auto& triangle : triangles) {
                    Expects(triangle->isLeftOfPlane(plane));
                    Expects(triangle->isRightOfPlane(plane));

                    auto leftClipBox = triangle->clippedBoundingBox(leftBox);
                    if (leftClipBox) {
                        for (const auto& event : buildTriangleEvents(*leftClipBox, triangle))
                        {
                            leftEvents.insert(event);
                        }
                        leftTriangles.insert(triangle);
                    }

                    auto rightClipBox = triangle->clippedBoundingBox(rightBox);
                    if (rightClipBox) {
                        for (const auto& event : buildTriangleEvents(*rightClipBox, triangle))
                        {
                            rightEvents.insert(event);
                        }
                        rightTriangles.insert(triangle);
                    }
                }
                triangles.clear();
                events.clear();

                return std::make_unique<MapNode<SurfaceData>>(BranchNode<SurfaceData>{
                    plane,
                    createNode(leftBox, std::move(leftEvents), std::move(leftTriangles)),
                    createNode(rightBox, std::move(rightEvents), std::move(rightTriangles)),
                });
            }
        }
        return std::make_unique<MapNode<SurfaceData>>(
            LeafNode<SurfaceData>{ std::move(triangles) });
    }

    std::set<Event> buildTriangleEvents(
        const AxisAlignedBoundingBox& boundingBox, Triangle<SurfaceData>* triangle)
    {
        AxisAlignedBoundingBox actualBox = triangle->getBoundingBox();
        std::set<Event> events;
        for (uint8_t k = 0; k < 3; k++) {
            float minimum = boundingBox.minimum[k];
            float maximum = boundingBox.maximum[k];
            Expects(!(minimum < actualBox.minimum[k]));
            Expects(!(maximum > actualBox.maximum[k]));
            if (minimum < maximum) {
                events.insert(
                    Event{ triangle, AxisAlignedPlane{ k, maximum }, Event::Type::Ending });
                events.insert(
                    Event{ triangle, AxisAlignedPlane{ k, minimum }, Event::Type::Starting });
            } else {
                events.insert(
                    Event{ triangle, AxisAlignedPlane{ k, minimum }, Event::Type::Planar });
            }
        }
        return events;
    }

    std::tuple<AxisAlignedPlane, Side, float> findBestSplit(
        const AxisAlignedBoundingBox& box, const std::set<Event>& events, size_t triangleCount)
    {
        float bestCost = std::numeric_limits<float>::infinity();
        AxisAlignedPlane bestPlane;
        Side bestSide;

        std::array<size_t, 3> leftTriangleCounts{ 0 };
        std::array<size_t, 3> rightTriangleCounts{ triangleCount, triangleCount, triangleCount };
        auto iter = events.begin();
        auto eventsEnd = events.end();
        while (iter != eventsEnd) {
            const AxisAlignedPlane& candidatePlane = iter->separationPlane;
            uint8_t dimension = candidatePlane.dimensionIndex;
            size_t endingTriangleCount = 0;
            size_t planarTriangleCount = 0;
            size_t startingTriangleCount = 0;
            auto validate = [&eventsEnd, &candidatePlane, dimension](const auto& iter) {
                return (iter != eventsEnd)
                    && !(iter->separationPlane.boundary > candidatePlane.boundary)
                    && (iter->separationPlane.dimensionIndex == dimension);
            };
            while (validate(iter) && (iter->type == Event::Type::Ending)) {
                endingTriangleCount++;
                iter++;
            }
            while (validate(iter) && (iter->type == Event::Type::Planar)) {
                planarTriangleCount++;
                iter++;
            }
            while (validate(iter) && (iter->type == Event::Type::Starting)) {
                startingTriangleCount++;
                iter++;
            }
            rightTriangleCounts[dimension] -= planarTriangleCount;
            rightTriangleCounts[dimension] -= endingTriangleCount;
            auto [side, cost] = surfaceAreaCostHeuristic(box, candidatePlane,
                leftTriangleCounts[dimension], planarTriangleCount, rightTriangleCounts[dimension]);
            if (cost < bestCost) {
                bestCost = cost;
                bestPlane = candidatePlane;
                bestSide = side;
            }
            leftTriangleCounts[dimension] += planarTriangleCount;
            leftTriangleCounts[dimension] += startingTriangleCount;
        }
        return { bestPlane, bestSide, bestCost };
    }

    std::pair<Side, float> surfaceAreaCostHeuristic(const AxisAlignedBoundingBox& box,
        const AxisAlignedPlane& splitPlane, size_t leftTriangleCount, size_t planarTriangleCount,
        size_t rightTriangleCount)
    {
        auto [leftBox, rightBox] = box.split(splitPlane);
        float boxSurfaceArea = box.getSurfaceArea();

        float leftArea = leftBox.getSurfaceArea() / boxSurfaceArea;
        float rightArea = rightBox.getSurfaceArea() / boxSurfaceArea;
        bool removesVolume = (leftBox.getVolume() > 0.0f) && (rightBox.getVolume() > 0.0f);

        auto cost = [leftArea, rightArea, removesVolume](
                        size_t leftTriangleCount, size_t rightTriangleCount) {
            float cost = static_cast<float>(leftTriangleCount) * leftArea
                + static_cast<float>(rightTriangleCount) * rightArea;
            if (removesVolume && (!leftTriangleCount || !rightTriangleCount)) {
                // We slightly prefer splits that remove empty space.
                cost *= kEmptySplitDiscount;
            }
            return cost;
        };

        float leftCost = cost(leftTriangleCount + planarTriangleCount, rightTriangleCount);
        float rightCost = cost(leftTriangleCount, planarTriangleCount + rightTriangleCount);
        if (leftCost < rightCost) {
            return { Side::Left, (kIntersectionCost * leftCost) + kTraversalCost };
        } else {
            return { Side::Right, (kIntersectionCost * rightCost) + kTraversalCost };
        }
    }

    AxisAlignedBoundingBox _boundingBox;
    std::vector<std::unique_ptr<Triangle<SurfaceData>>> _triangles;
    std::set<Event> _events;
};
}

ext: h

max_stars_repo_name: eyebrowsoffire/rev

max_issues_repo_name: eyebrowsoffire/rev

max_forks_repo_name: eyebrowsoffire/rev


text:
// Copyright András Vukics 2006–2020. Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE.txt)
/// \briefFile{Random & related}
#ifndef CPPQEDCORE_UTILS_RANDOM_H_REENTRANT

#ifndef CPPQEDCORE_UTILS_RANDOM_H_INCLUDED
#define CPPQEDCORE_UTILS_RANDOM_H_INCLUDED

#include "Pars.h"

#include "pcg_random.hpp"
#include "XoshiroCpp.hpp"

#ifdef CPPQED_HAS_GSL
#include <gsl/gsl_rng.h>
#endif // CPPQED_HAS_GSL

#include <boost/range/algorithm/generate.hpp>

#ifdef BZ_HAVE_BOOST_SERIALIZATION
#include <boost/serialization/string.hpp>
#include <boost/serialization/array.hpp>
#endif // BZ_HAVE_BOOST_SERIALIZATION

#include <memory>
#include <optional>
#include <random>
#include <stdexcept>
#include <sstream>
#include <optional>

namespace randomutils {

template<typename Engine, typename BASE>
struct Pars;

#ifdef CPPQED_HAS_GSL

/// Wraps GSL random number generators into the concept [UniformRandomBitGenerator](http://en.cppreference.com/w/cpp/named_req/UniformRandomBitGenerator)
class GSL_Engine
{
public:
  using result_type=unsigned long;

  explicit GSL_Engine(result_type s/*, const gsl_rng_type* ran_gen_type=gsl_rng_taus2*/);

  void seed(result_type value);

  result_type operator()();

  result_type min() const;
  result_type max() const;

  void write(std::ostream& os) const;

  void read(std::istream& is);

private:
  const std::shared_ptr<gsl_rng> ranGen_;
  const std::string name_;

};

inline std::ostream& operator<<(std::ostream& os, const GSL_Engine& r) {r.write(os); return os;}
inline std::istream& operator>>(std::istream& is, GSL_Engine& r) {r.read(is); return is;}

#define CPPQEDCORE_UTILS_RANDOM_H_REENTRANT
#define CPPQEDCORE_UTILS_RANDOM_H_RANDOMENGINE GSL_Engine
#define CPPQEDCORE_UTILS_RANDOM_H_STRING "GSL_Engine state"

#include "Random.h"

} // randomutils

#endif // CPPQED_HAS_GSL

namespace randomutils {

template<typename Distribution, typename Engine, typename Array, typename... DCP>
Engine& fill(Array a, Engine& re, DCP&&... dcp)
{
  Distribution d{std::forward<DCP>(dcp)...};
  std::generate(a.begin(),a.end(),[&]() {return d(re);});
  return re;
}

template<typename Distribution, typename Engine, typename Array, typename... DCP>
void fill(Array a, unsigned long seed, DCP&&... dcp)
{
  Engine re{seed};
  fill<Distribution>(a,re,std::forward<DCP>(dcp)...);
}

template<typename Engine>
constexpr auto EngineID_v = std::nullopt;

template<>
inline const std::string EngineID_v<GSL_Engine> = "GSL_Taus2";

template<>
inline const std::string EngineID_v<std::mt19937_64> = "STD_MT19937_64";

template<>
inline const std::string EngineID_v<pcg64> = "PCG64";

template<>
inline const std::string EngineID_v<XoshiroCpp::Xoshiro256PlusPlus> = "Xoshiro256pp";

template<typename Engine>
struct EngineWithParameters
{
  EngineWithParameters(unsigned long s, unsigned long p) : engine{s,p}, seed{s}, prngStream{p} {}

  std::ostream& stream(std::ostream& os) const {return os << "Random engine: "<<EngineID_v<Engine><<". Parameters: seed="<<seed<<", streamOrdo=" << prngStream;}

  Engine engine;
  const unsigned long seed, prngStream;

};

// template<>
// inline std::list<pcg64> independentPRNG_Streams<pcg64>(size_t n, pcg64::result_type seed, pcg64::result_type ordoStream)
// {
//   std::list<pcg64> res;
//   for (size_t i=0; i<n; ++i) res.emplace_back(seed,1001+ordoStream+i);
//   return res;
// }

template<typename BASE>
struct Pars<pcg64,BASE> : BASE
{
  unsigned long &seed, &prngStream; ///< PCG64 generator seed and ordinal of independent stream

  Pars(parameters::Table& p, const std::string& mod="")
    : BASE{p,mod},
      seed(p.addTitle("StochasticTrajectory",mod).add("seed",mod,"Random number generator seed",1001ul)),
      prngStream(p.add("prngStream",mod,"Random number generator independent stream ordinal",1ul))
    {}
};

#ifndef NDEBUG
#pragma GCC warning "TODO: implement nextStream optionally from previousState"
#endif // NDEBUG
template<typename BASE>
pcg64 streamOfOrdo(const Pars<pcg64,BASE>& p/*, std::optional<EngineState> previousState = std::nullopt*/)
{
  return pcg64{p.seed,1001+p.prngStream};
}

template<typename BASE>
void incrementForNextStream(const Pars<pcg64,BASE>& p) {++p.prngStream;}

} // randomutils

namespace boost { namespace serialization {

template<typename Ar>
void load(Ar& ar, XoshiroCpp::Xoshiro256PlusPlus& xpp, unsigned) {
  XoshiroCpp::Xoshiro256PlusPlus::state_type state;
  ar & state;
  xpp.deserialize(state);
}

template<typename Ar>
void save(Ar& ar, XoshiroCpp::Xoshiro256PlusPlus const& xpp, unsigned) {
  auto state{xpp.serialize()};
  ar & state;
}

template<typename Ar>
void serialize(Ar& ar, XoshiroCpp::Xoshiro256PlusPlus& xpp, unsigned version) {
  if (typename Ar::is_saving())
    save(ar, xpp, version);
  else
    load(ar, xpp, version);
}

#define CPPQEDCORE_UTILS_RANDOM_H_REENTRANT
#define CPPQEDCORE_UTILS_RANDOM_H_RANDOMENGINE std::mt19937_64
#define CPPQEDCORE_UTILS_RANDOM_H_STRING "mersenne_twister_engine state"

#include "Random.h"

#define CPPQEDCORE_UTILS_RANDOM_H_REENTRANT
#define CPPQEDCORE_UTILS_RANDOM_H_RANDOMENGINE pcg64
#define CPPQEDCORE_UTILS_RANDOM_H_STRING "pcg64 state"

#include "Random.h"

} } // boost::serialization

#endif // CPPQEDCORE_UTILS_RANDOM_H_INCLUDED

#else // CPPQEDCORE_UTILS_RANDOM_H_REENTRANT

template<typename Ar>
void load(Ar& ar, CPPQEDCORE_UTILS_RANDOM_H_RANDOMENGINE& randomEngine, unsigned) {
  std::string text;
  ar & text;
  std::istringstream iss(text);

  if (!(iss >> randomEngine))
    throw std::invalid_argument(CPPQEDCORE_UTILS_RANDOM_H_STRING);
}

template<typename Ar>
void save(Ar& ar, CPPQEDCORE_UTILS_RANDOM_H_RANDOMENGINE const& randomEngine, unsigned) {
  std::ostringstream oss;
  if (!(oss << randomEngine))
    throw std::invalid_argument(CPPQEDCORE_UTILS_RANDOM_H_STRING);
  std::string text = oss.str();
  ar & text;
}

template<typename Ar>
void serialize(Ar& ar, CPPQEDCORE_UTILS_RANDOM_H_RANDOMENGINE& randomEngine, unsigned version) {
  if (typename Ar::is_saving())
    save(ar, randomEngine, version);
  else
    load(ar, randomEngine, version);
}

#undef CPPQEDCORE_UTILS_RANDOM_H_STRING
#undef CPPQEDCORE_UTILS_RANDOM_H_RANDOMENGINE
#undef CPPQEDCORE_UTILS_RANDOM_H_REENTRANT

#endif // CPPQEDCORE_UTILS_RANDOM_H_REENTRANT

max_stars_repo_name: vukics/cppqed

max_issues_repo_name: vukics/cppqed

max_forks_repo_name: vukics/cppqed