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

Copyright 2011-2019, Javier Burguete Tolosa.

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

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

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

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

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

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

/**
 * Function to obtain the coefficients of a 5 steps 3rd order Runge-Kutta
 * method.
 */
int
rk_tb_5_3 (Optimize * optimize) ///< Optimize struct.
{
  long double A[3], B[3], C[3], D[3];
  long double *tb, *r;
#if DEBUG_RK_5_3
  fprintf (stderr, "rk_tb_5_3: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t5 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  b21 (tb) = r[2];
  t3 (tb) = r[3];
  b31 (tb) = r[4];
  b32 (tb) = r[5];
  t4 (tb) = r[6];
  b41 (tb) = r[7];
  b42 (tb) = r[8];
  b43 (tb) = r[9];
  b54 (tb) = r[10];
  A[0] = t1 (tb);
  B[0] = t2 (tb);
  C[0] = t3 (tb);
  D[0] = 0.5L - b54 (tb) * t4 (tb);
  A[1] = A[0] * t1 (tb);
  B[1] = B[0] * t2 (tb);
  C[1] = C[0] * t3 (tb);
  D[1] = 1.L / 3.L - b54 (tb) * sqr (t4 (tb));
  A[2] = 0.L;
  B[2] = b21 (tb) * t1 (tb);
  C[2] = b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb);
  D[2] = 1.L / 6.L - b54 (tb) * (b41 (tb) * t1 (tb) + b42 (tb) * t2 (tb)
                                 + b43 (tb) * t3 (tb));
  solve_3 (A, B, C, D);
  if (isnan (D[0]) || isnan (D[1]) || isnan (D[2]))
    return 0;
  b53 (tb) = D[2];
  b52 (tb) = D[1];
  b51 (tb) = D[0];
  rk_b_5 (tb);
#if DEBUG_RK_5_3
  fprintf (stderr, "rk_tb_5_3: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 5 steps 3rd order, 4th order in
 * equations depending only in time, Runge-Kutta method.
 */
int
rk_tb_5_3t (Optimize * optimize)        ///< Optimize struct.
{
  long double A[4], B[4], C[4], D[4], E[4];
  long double *tb, *r;
#if DEBUG_RK_5_3
  fprintf (stderr, "rk_tb_5_3t: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t5 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  b21 (tb) = r[2];
  t3 (tb) = r[3];
  b31 (tb) = r[4];
  b32 (tb) = r[5];
  t4 (tb) = r[6];
  b41 (tb) = r[7];
  b42 (tb) = r[8];
  b43 (tb) = r[9];
  A[0] = t1 (tb);
  B[0] = t2 (tb);
  C[0] = t3 (tb);
  D[0] = t4 (tb);
  E[0] = 0.5L;
  A[1] = A[0] * t1 (tb);
  B[1] = B[0] * t2 (tb);
  C[1] = C[0] * t3 (tb);
  D[1] = D[0] * t4 (tb);
  E[1] = 1.L / 3.L;
  A[2] = A[1] * t1 (tb);
  B[2] = B[1] * t2 (tb);
  C[2] = C[1] * t3 (tb);
  D[2] = D[1] * t4 (tb);
  E[2] = 0.25L;
  A[3] = 0.L;
  B[3] = b21 (tb) * t1 (tb);
  C[3] = b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb);
  D[3] = b41 (tb) * t1 (tb) + b42 (tb) * t2 (tb) + b43 (tb) * t3 (tb);
  E[3] = 1.L / 6.L;
  solve_4 (A, B, C, D, E);
  if (isnan (E[0]) || isnan (E[1]) || isnan (E[2]) || isnan (E[3]))
    return 0;
  b54 (tb) = E[3];
  b53 (tb) = E[2];
  b52 (tb) = E[1];
  b51 (tb) = E[0];
  rk_b_5 (tb);
#if DEBUG_RK_5_3
  fprintf (stderr, "rk_tb_5_3t: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 5 steps 2nd-3rd order Runge-Kutta
 * pair.
 */
int
rk_tb_5_3p (Optimize * optimize)        ///< Optimize struct.
{
  long double *tb;
#if DEBUG_RK_5_3
  fprintf (stderr, "rk_tb_5_3p: start\n");
#endif
  if (!rk_tb_5_3 (optimize))
    return 0;
  tb = optimize->coefficient;
  e51 (tb) = 0.5L / t1 (tb);
  e52 (tb) = e53 (tb) = 0.L;
  rk_e_5 (tb);
#if DEBUG_RK_5_3
  fprintf (stderr, "rk_tb_5_3p: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 5 steps 2nd-3rd order, 3rd-4th order
 * in equations depending only in time, Runge-Kutta pair.
 */
int
rk_tb_5_3tp (Optimize * optimize)       ///< Optimize struct.
{
  long double *tb;
#if DEBUG_RK_5_3
  fprintf (stderr, "rk_tb_5_3tp: start\n");
#endif
  if (!rk_tb_5_3t (optimize))
    return 0;
  tb = optimize->coefficient;
  e53 (tb) = 0.L;
  e52 (tb) = (1.L / 3.L - 0.5L * t1 (tb)) / (t2 (tb) * (t2 (tb) - t1 (tb)));
  if (isnan (e52 (tb)))
    return 0;
  e51 (tb) = (0.5L - e52 (tb) * t2 (tb)) / t1 (tb);
  if (isnan (e51 (tb)))
    return 0;
  rk_e_5 (tb);
#if DEBUG_RK_5_3
  fprintf (stderr, "rk_tb_5_3tp: end\n");
#endif
  return 1;
}

/**
 * Function to calculate the objective function of a 5 steps 3rd order
 * Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_5_3 (RK * rk)   ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_5_3
  fprintf (stderr, "rk_objective_tb_5_3: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b50 (tb) < 0.L)
    o += b50 (tb);
  if (b51 (tb) < 0.L)
    o += b51 (tb);
  if (b52 (tb) < 0.L)
    o += b52 (tb);
  if (b53 (tb) < 0.L)
    o += b53 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L
    + fmaxl (1.L, fmaxl (t1 (tb), fmaxl (t2 (tb), fmaxl (t3 (tb), t4 (tb)))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_5_3
  fprintf (stderr, "rk_objective_tb_5_3: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_5_3: end\n");
#endif
  return o;
}

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

/**
 * Function to calculate the objective function of a 5 steps 2nd-3rd order
 * Runge-Kutta pair.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_5_3p (RK * rk)  ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_5_3
  fprintf (stderr, "rk_objective_tb_5_3p: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b50 (tb) < 0.L)
    o += b50 (tb);
  if (b51 (tb) < 0.L)
    o += b51 (tb);
  if (b52 (tb) < 0.L)
    o += b52 (tb);
  if (b53 (tb) < 0.L)
    o += b53 (tb);
  if (e50 (tb) < 0.L)
    o += e50 (tb);
  if (e51 (tb) < 0.L)
    o += e51 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L
    + fmaxl (1.L, fmaxl (t1 (tb), fmaxl (t2 (tb), fmaxl (t3 (tb), t4 (tb)))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_5_3
  fprintf (stderr, "rk_objective_tb_5_3p: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_5_3p: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 5 steps 2nd-3rd order,
 * 3rd-4th order in equations depending only in time, Runge-Kutta pair.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_5_3tp (RK * rk) ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_5_3
  fprintf (stderr, "rk_objective_tb_5_3tp: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b50 (tb) < 0.L)
    o += b50 (tb);
  if (b51 (tb) < 0.L)
    o += b51 (tb);
  if (b52 (tb) < 0.L)
    o += b52 (tb);
  if (b53 (tb) < 0.L)
    o += b53 (tb);
  if (b54 (tb) < 0.L)
    o += b54 (tb);
  if (e50 (tb) < 0.L)
    o += e50 (tb);
  if (e51 (tb) < 0.L)
    o += e51 (tb);
  if (e52 (tb) < 0.L)
    o += e52 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L
    + fmaxl (1.L, fmaxl (t1 (tb), fmaxl (t2 (tb), fmaxl (t3 (tb), t4 (tb)))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_5_3
  fprintf (stderr, "rk_objective_tb_5_3tp: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_5_3tp: end\n");
#endif
  return o;
}

max_stars_repo_name: jburguete/ode

max_issues_repo_name: jburguete/ode

max_forks_repo_name: jburguete/ode


text:
#ifndef UTILS_H
#define UTILS_H

#include <gsl/gsl_vector.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_sf_psi.h>
#include <math.h>
#include <assert.h>
#include <time.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <unistd.h>
#include <dirent.h>
#include <memory.h>

#define outlog(format, args...) \
    fprintf(stderr, format, args); \
    fprintf(stderr, "\n");

int compare (const void * a, const void * b);

inline double safe_log(double x) {
  if (x <= 0)
    return(-10000);
  else
    return(log(x));
}
double log_sum(double, double);

inline double vget(const gsl_vector* v, int i) { return(gsl_vector_get(v, i)); }

inline void vset(gsl_vector* v, int i, double x) { gsl_vector_set(v, i, x); }

// Increment a vector element by a double.
inline void vinc(gsl_vector* v, int i, double x) {
  vset(v, i, vget(v, i) + x);
}

inline double mget(const gsl_matrix* m, int i, int j)
{ return(gsl_matrix_get(m, i, j)); }

inline void mset(gsl_matrix* m, int i, int j, double x)
{ gsl_matrix_set(m, i, j, x); }

// Increment a matrix element by a double.
void minc(gsl_matrix*, int, int, double);

void col_sum(const gsl_matrix*, gsl_vector*);
void row_sum(const gsl_matrix*, gsl_vector*);

void vct_fprintf(FILE* file, const gsl_vector* v);
void mtx_fprintf(FILE* file, const gsl_matrix* m);
void mtx_fscanf(FILE* file, gsl_matrix* m);

inline bool check_sym(const gsl_matrix *m) {
  for (size_t i = 0; i < m->size1-1; i ++)
    for (size_t j=i; j < m->size2; j ++)
      if (mget(m, i, j) != mget(m, j, i)) {
        printf("not sym\n");
        return false;
      }
  return true;
}

double log_det(const gsl_matrix*);

void matrix_inverse(const gsl_matrix*, gsl_matrix*);
void matrix_vector_solve(const gsl_matrix* m, const gsl_vector* b, gsl_vector* v);

void sym_eigen(gsl_matrix*, gsl_vector*, gsl_matrix*);

inline double vsum(const gsl_vector* v) {
  double val = 0;
  int i, size = v->size;
  for (i = 0; i < size; i++)
    val += vget(v, i);
  return(val);
}

double vnorm(const gsl_vector * v);

void gsl_vector_apply(gsl_vector* x, double(*fun)(double));
void vct_log(gsl_vector* v);
void mtx_log(gsl_matrix* x);
void vct_exp(gsl_vector* x);
void mtx_exp(gsl_matrix* x);

double mahalanobis_distance(const gsl_matrix * m, const gsl_vector* u, const gsl_vector* v);
double mahalanobis_prod(const gsl_matrix * m, const gsl_vector* u, const gsl_vector* v);
double matrix_dot_prod(const gsl_matrix * m1, const gsl_matrix* m2);

void choose_k_from_n(int k, int n, int* result, int* src);

double log_normalize(gsl_vector* x);
double vnormalize(gsl_vector* x);

int dir_exists(const char *dname);
bool file_exists(const char * filename);
void make_directory(const char* name);

double  digamma(double x);
unsigned int rmultinomial(const gsl_vector* v);
double  rgamma(double a, double b);
double  rbeta(double a, double b);
unsigned int rbernoulli(double p);
double runiform();
void rshuffle (void* base, size_t n, size_t size);
unsigned long int runiform_int(unsigned long int n);

// new and free random number generator
gsl_rng* new_random_number_generator(long seed);
void free_random_number_generator(gsl_rng * random_number_generator);

#endif

ext: h


text:
/* randist/binomial_tpe.c
 *
 * Copyright (C) 1996, 2003, 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 <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_pow_int.h>
#include <gsl/gsl_sf_gamma.h>

/* The binomial distribution has the form,

   f(x) =  n!/(x!(n-x)!) * p^x (1-p)^(n-x) for integer 0 <= x <= n
        =  0                               otherwise

   This implementation follows the public domain ranlib function
   "ignbin", the bulk of which is the BTPE (Binomial Triangle
   Parallelogram Exponential) algorithm introduced in
   Kachitvichyanukul and Schmeiser[1].  It has been translated to use
   modern C coding standards.

   If n is small and/or p is near 0 or near 1 (specifically, if
   n*min(p,1-p) < SMALL_MEAN), then a different algorithm, called
   BINV, is used which has an average runtime that scales linearly
   with n*min(p,1-p).

   But for larger problems, the BTPE algorithm takes the form of two
   functions b(x) and t(x) -- "bottom" and "top" -- for which b(x) <
   f(x)/f(M) < t(x), with M = floor(n*p+p).  b(x) defines a triangular
   region, and t(x) includes a parallelogram and two tails.  Details
   (including a nice drawing) are in the paper.

   [1] Kachitvichyanukul, V. and Schmeiser, B. W.  Binomial Random
   Variate Generation.  Communications of the ACM, 31, 2 (February,
   1988) 216.

   Note, Bruce Schmeiser (personal communication) points out that if
   you want very fast binomial deviates, and you are happy with
   approximate results, and/or n and n*p are both large, then you can
   just use gaussian estimates: mean=n*p, variance=n*p*(1-p).

   This implementation by James Theiler, April 2003, after obtaining
   permission -- and some good advice -- from Drs. Kachitvichyanukul
   and Schmeiser to use their code as a starting point, and then doing
   a little bit of tweaking.

   Additional polishing for GSL coding standards by Brian Gough.  */

#define SMALL_MEAN 14           /* If n*p < SMALL_MEAN then use BINV
                                   algorithm. The ranlib
                                   implementation used cutoff=30; but
                                   on my computer 14 works better */

#define BINV_CUTOFF 110         /* In BINV, do not permit ix too large */

#define FAR_FROM_MEAN 20        /* If ix-n*p is larger than this, then
                                   use the "squeeze" algorithm.
                                   Ranlib used 20, and this seems to
                                   be the best choice on my machine as
                                   well */

#define LNFACT(x) gsl_sf_lnfact(x)

inline static double
Stirling (double y1)
{
  double y2 = y1 * y1;
  double s =
    (13860.0 -
     (462.0 - (132.0 - (99.0 - 140.0 / y2) / y2) / y2) / y2) / y1 / 166320.0;
  return s;
}

unsigned int
gsl_ran_binomial_tpe (const gsl_rng * rng, double p, unsigned int n)
{
  return gsl_ran_binomial (rng, p, n);
}

unsigned int
gsl_ran_binomial (const gsl_rng * rng, double p, unsigned int n)
{
  int ix;                       /* return value */
  int flipped = 0;
  double q, s, np;

  if (n == 0)
    return 0;

  if (p > 0.5)
    {
      p = 1.0 - p;              /* work with small p */
      flipped = 1;
    }

  q = 1 - p;
  s = p / q;
  np = n * p;

  /* Inverse cdf logic for small mean (BINV in K+S) */

  if (np < SMALL_MEAN)
    {
      double f0 = gsl_pow_uint (q, n);   /* f(x), starting with x=0 */

      while (1)
        {
          /* This while(1) loop will almost certainly only loop once; but
           * if u=1 to within a few epsilons of machine precision, then it
           * is possible for roundoff to prevent the main loop over ix to
           * achieve its proper value.  following the ranlib implementation,
           * we introduce a check for that situation, and when it occurs,
           * we just try again.
           */

          double f = f0;
          double u = gsl_rng_uniform (rng);

          for (ix = 0; ix <= BINV_CUTOFF; ++ix)
            {
              if (u < f)
                goto Finish;
              u -= f;
              /* Use recursion f(x+1) = f(x)*[(n-x)/(x+1)]*[p/(1-p)] */
              f *= s * (n - ix) / (ix + 1);
            }

          /* It should be the case that the 'goto Finish' was encountered
           * before this point was ever reached.  But if we have reached
           * this point, then roundoff has prevented u from decreasing
           * all the way to zero.  This can happen only if the initial u
           * was very nearly equal to 1, which is a rare situation.  In
           * that rare situation, we just try again.
           *
           * Note, following the ranlib implementation, we loop ix only to
           * a hardcoded value of SMALL_MEAN_LARGE_N=110; we could have
           * looped to n, and 99.99...% of the time it won't matter.  This
           * choice, I think is a little more robust against the rare
           * roundoff error.  If n>LARGE_N, then it is technically
           * possible for ix>LARGE_N, but it is astronomically rare, and
           * if ix is that large, it is more likely due to roundoff than
           * probability, so better to nip it at LARGE_N than to take a
           * chance that roundoff will somehow conspire to produce an even
           * larger (and more improbable) ix.  If n<LARGE_N, then once
           * ix=n, f=0, and the loop will continue until ix=LARGE_N.
           */
        }
    }
  else
    {
      /* For n >= SMALL_MEAN, we invoke the BTPE algorithm */

      int k;

      double ffm = np + p;      /* ffm = n*p+p             */
      int m = (int) ffm;        /* m = int floor[n*p+p]    */
      double fm = m;            /* fm = double m;          */
      double xm = fm + 0.5;     /* xm = half integer mean (tip of triangle)  */
      double npq = np * q;      /* npq = n*p*q            */

      /* Compute cumulative area of tri, para, exp tails */

      /* p1: radius of triangle region; since height=1, also: area of region */
      /* p2: p1 + area of parallelogram region */
      /* p3: p2 + area of left tail */
      /* p4: p3 + area of right tail */
      /* pi/p4: probability of i'th area (i=1,2,3,4) */

      /* Note: magic numbers 2.195, 4.6, 0.134, 20.5, 15.3 */
      /* These magic numbers are not adjustable...at least not easily! */

      double p1 = floor (2.195 * sqrt (npq) - 4.6 * q) + 0.5;

      /* xl, xr: left and right edges of triangle */
      double xl = xm - p1;
      double xr = xm + p1;

      /* Parameter of exponential tails */
      /* Left tail:  t(x) = c*exp(-lambda_l*[xl - (x+0.5)]) */
      /* Right tail: t(x) = c*exp(-lambda_r*[(x+0.5) - xr]) */

      double c = 0.134 + 20.5 / (15.3 + fm);
      double p2 = p1 * (1.0 + c + c);

      double al = (ffm - xl) / (ffm - xl * p);
      double lambda_l = al * (1.0 + 0.5 * al);
      double ar = (xr - ffm) / (xr * q);
      double lambda_r = ar * (1.0 + 0.5 * ar);
      double p3 = p2 + c / lambda_l;
      double p4 = p3 + c / lambda_r;

      double var, accept;
      double u, v;              /* random variates */

    TryAgain:

      /* generate random variates, u specifies which region: Tri, Par, Tail */
      u = gsl_rng_uniform (rng) * p4;
      v = gsl_rng_uniform (rng);

      if (u <= p1)
        {
          /* Triangular region */
          ix = (int) (xm - p1 * v + u);
          goto Finish;
        }
      else if (u <= p2)
        {
          /* Parallelogram region */
          double x = xl + (u - p1) / c;
          v = v * c + 1.0 - fabs (x - xm) / p1;
          if (v > 1.0 || v <= 0.0)
            goto TryAgain;
          ix = (int) x;
        }
      else if (u <= p3)
        {
          /* Left tail */
          ix = (int) (xl + log (v) / lambda_l);
          if (ix < 0)
            goto TryAgain;
          v *= ((u - p2) * lambda_l);
        }
      else
        {
          /* Right tail */
          ix = (int) (xr - log (v) / lambda_r);
          if (ix > (double) n)
            goto TryAgain;
          v *= ((u - p3) * lambda_r);
        }

      /* At this point, the goal is to test whether v <= f(x)/f(m)
       *
       *  v <= f(x)/f(m) = (m!(n-m)! / (x!(n-x)!)) * (p/q)^{x-m}
       *
       */

      /* Here is a direct test using logarithms.  It is a little
       * slower than the various "squeezing" computations below, but
       * if things are working, it should give exactly the same answer
       * (given the same random number seed).  */

#ifdef DIRECT
      var = log (v);

      accept =
        LNFACT (m) + LNFACT (n - m) - LNFACT (ix) - LNFACT (n - ix)
        + (ix - m) * log (p / q);

#else /* SQUEEZE METHOD */

      /* More efficient determination of whether v < f(x)/f(M) */

      k = abs (ix - m);

      if (k <= FAR_FROM_MEAN)
        {
          /*
           * If ix near m (ie, |ix-m|<FAR_FROM_MEAN), then do
           * explicit evaluation using recursion relation for f(x)
           */
          double g = (n + 1) * s;
          double f = 1.0;

          var = v;

          if (m < ix)
            {
              int i;
              for (i = m + 1; i <= ix; i++)
                {
                  f *= (g / i - s);
                }
            }
          else if (m > ix)
            {
              int i;
              for (i = ix + 1; i <= m; i++)
                {
                  f /= (g / i - s);
                }
            }

          accept = f;
        }
      else
        {
          /* If ix is far from the mean m: k=ABS(ix-m) large */

          var = log (v);

          if (k < npq / 2 - 1)
            {
              /* "Squeeze" using upper and lower bounds on
               * log(f(x)) The squeeze condition was derived
               * under the condition k < npq/2-1 */
              double amaxp =
                k / npq * ((k * (k / 3.0 + 0.625) + (1.0 / 6.0)) / npq + 0.5);
              double ynorm = -(k * k / (2.0 * npq));
              if (var < ynorm - amaxp)
                goto Finish;
              if (var > ynorm + amaxp)
                goto TryAgain;
            }

          /* Now, again: do the test log(v) vs. log f(x)/f(M) */

#if USE_EXACT
          /* This is equivalent to the above, but is a little (~20%) slower */
          /* There are five log's vs three above, maybe that's it? */

          accept = LNFACT (m) + LNFACT (n - m)
            - LNFACT (ix) - LNFACT (n - ix) + (ix - m) * log (p / q);

#else /* USE STIRLING */
          /* The "#define Stirling" above corresponds to the first five
           * terms in asymptoic formula for
           * log Gamma (y) - (y-0.5)log(y) + y - 0.5 log(2*pi);
           * See Abramowitz and Stegun, eq 6.1.40
           */

          /* Note below: two Stirling's are added, and two are
           * subtracted.  In both K+S, and in the ranlib
           * implementation, all four are added.  I (jt) believe that
           * is a mistake -- this has been confirmed by personal
           * correspondence w/ Dr. Kachitvichyanukul.  Note, however,
           * the corrections are so small, that I couldn't find an
           * example where it made a difference that could be
           * observed, let alone tested.  In fact, define'ing Stirling
           * to be zero gave identical results!!  In practice, alv is
           * O(1), ranging 0 to -10 or so, while the Stirling
           * correction is typically O(10^{-5}) ...setting the
           * correction to zero gives about a 2% performance boost;
           * might as well keep it just to be pendantic.  */

          {
            double x1 = ix + 1.0;
            double w1 = n - ix + 1.0;
            double f1 = fm + 1.0;
            double z1 = n + 1.0 - fm;

            accept = xm * log (f1 / x1) + (n - m + 0.5) * log (z1 / w1)
              + (ix - m) * log (w1 * p / (x1 * q))
              + Stirling (f1) + Stirling (z1) - Stirling (x1) - Stirling (w1);
          }
#endif
#endif
        }

      if (var <= accept)
        {
          goto Finish;
        }
      else
        {
          goto TryAgain;
        }
    }

Finish:

  return (flipped) ? (n - ix) : (unsigned int)ix;
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* monte/gsl_monte_plain.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 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.
 */

/* Plain Monte-Carlo. */

/* Author: MJB */

#ifndef __GSL_MONTE_PLAIN_H__
#define __GSL_MONTE_PLAIN_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 <stdio.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_rng.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 dim;
  double *x;
} gsl_monte_plain_state;

GSL_FUN int
gsl_monte_plain_integrate (const gsl_monte_function * f,
                           const double xl[], const double xu[],
                           const size_t dim,
                           const size_t calls,
                           gsl_rng * r,
                           gsl_monte_plain_state * state,
                           double *result, double *abserr);

GSL_FUN gsl_monte_plain_state* gsl_monte_plain_alloc(size_t dim);

GSL_FUN int gsl_monte_plain_init(gsl_monte_plain_state* state);

GSL_FUN void gsl_monte_plain_free (gsl_monte_plain_state* state);

__END_DECLS

#endif /* __GSL_MONTE_PLAIN_H__ */

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:
#pragma once

#include <tuple>

#include "arcana/threading/dispatcher.h"

#include "router.h"

#include <gsl/gsl>

namespace arcana
{
    //
    // A mediator is an event pool that ensures all events are sent
    // through the right dispatcher (execution context).
    // This is usefull when you want anyone to be able to send
    // events but only want to process the events from one thread.
    //
    template<typename DispatcherT, typename... EventTs>
    class mediator
    {
        using router_t = router<EventTs...>;

    public:
        using dispatcher_t = DispatcherT;

        explicit mediator(dispatcher_t& dispatcher)
            : m_dispatcher{ dispatcher }
        {}

        template<typename T>
        void send(T&& evt)
        {
            m_dispatcher.queue([ this, evt = std::forward<T>(evt) ]() { m_router.fire(evt); });
        }

        template<typename EventT, typename T>
        ticket add_listener(T&& listener)
        {
            GSL_CONTRACT_CHECK("thread affinity", m_dispatcher.get_affinity().check());
            return m_router.template add_listener<EventT>(std::forward<T>(listener));
        }

        dispatcher_t& dispatcher()
        {
            return m_dispatcher;
        }

    private:
        dispatcher_t& m_dispatcher;
        router_t m_router;
    };
}


text:
/**
 *
 * @file core_dlatro.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Azzam Haidar
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:49 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_double
 *
 *  CORE_dlatro transposes a m-by-n matrix out of place.
 *
 *******************************************************************************
 *
 * @param[in] uplo
 *          Specifies whether the matrix A is upper triangular or lower
 *          triangular:
 *          = PlasmaUpper: the upper triangle of A and the lower triangle of B
 *          are referenced.
 *          = PlasmaLower: the lower triangle of A and the upper triangle of B
 *          are referenced.
 *          = PlasmaUpperLower: All A and B are referenced.
 *
 * @param[in] trans
 *          Specifies whether the matrix A is transposed, not transposed or
 *          ugate transposed:
 *          = PlasmaNoTrans:   B is a copy of A (equivalent to dlacpy);
 *          = PlasmaTrans:     B is the transpose of A;
 *          = PlasmaTrans: B is the ugate transpose of A.
 *
 * @param[in] M
 *         Number of rows of the matrix A and number of columns of the matrix B, if trans == Pasma[Conj]Trans.
 *         Number of rows of the matrix A and the matrix B, if trans == PasmaNoTrans.
 *
 * @param[in] N
 *         Number of columns of the matrix A and number of rows of the matrix B, if trans == Pasma[Conj]Trans.
 *         Number of columns of the matrix A and of the matrix B, if trans == PlasmaNoTrans.
 *
 * @param[in] A
 *         Matrix of size LDA-by-N, if trans == Pasma[Conj]Trans.
 *         Matrix of size LDA-by-M, if trans == PasmaNoTrans.
 *
 * @param[in] LDA
 *         The leading dimension of the array A.
 *         LDA >= max(1,M), if trans == Pasma[Conj]Trans.
 *         LDA >= max(1,N), if trans == PasmaNoTrans.
 *
 * @param[out] B
 *         Matrix of size LDB-by-M, if trans == Pasma[Conj]Trans.
 *         Matrix of size LDB-by-N, if trans == PasmaNoTrans.
 *
 * @param[in] LDB
 *         The leading dimension of the array B.
 *         LDB >= max(1,N), if trans == Pasma[Conj]Trans.
 *         LDB >= max(1,M), if trans == PasmaNoTrans.
 *
 *
 *******************************************************************************
 *
 * @return
 *         \retval PLASMA_SUCCESS successful exit
 *         \retval <0 if INFO = -k, the k-th argument had an illegal value
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dlatro = PCORE_dlatro
#define CORE_dlatro PCORE_dlatro
#endif
int CORE_dlatro(PLASMA_enum uplo, PLASMA_enum trans,
                int M, int N,
                const double *A, int LDA,
                      double *B, int LDB)
{
    int i, j;

    /* Check input arguments */
    if ((uplo != PlasmaUpper) && (uplo != PlasmaLower) && (uplo != PlasmaUpperLower) ) {
        coreblas_error(1, "Illegal value of uplo");
        return -1;
    }
    if ((trans != PlasmaTrans) && (trans != PlasmaNoTrans) && (trans != PlasmaTrans) ) {
        coreblas_error(2, "Illegal value of trans");
        return -2;
    }
    if (M < 0) {
        coreblas_error(3, "Illegal value of M");
        return -3;
    }
    if (N < 0) {
        coreblas_error(4, "Illegal value of N");
        return -4;
    }
    if ( (LDA < max(1,M)) && (M > 0) ) {
        coreblas_error(6, "Illegal value of LDA");
        return -6;
    }
    if ( (LDB < max(1,N)) && (N > 0) ) {
        coreblas_error(8, "Illegal value of LDB");
        return -8;
    }

    if (trans == PlasmaNoTrans) {
        CORE_dlacpy(uplo, M, N, A, LDA, B, LDB);
    }
    else {
        if (trans == PlasmaTrans) {
            if(uplo == PlasmaUpper) {
                for(j=0; j<N; j++)
                    for(i=0; i<min(j+1,M); i++)
                        B[j+i*LDB] = (A[i+j*LDA]);
            }
            else if(uplo == PlasmaLower) {
                for(j=0;j<N;j++)
                    for(i=j;i<M;i++)
                        B[j+i*LDB] = (A[i+j*LDA]);
            }
            else {
                for(j=0;j<N;j++)
                    for(i=0;i<M;i++)
                        B[j+i*LDB] = (A[i+j*LDA]);
            }
        }
        else {
            if(uplo==PlasmaUpper) {
                for(j=0;j<N;j++)
                    for(i=0;i<min(j+1,M);i++)
                        B[j+i*LDB] = A[i+j*LDA];
            }
            else if(uplo==PlasmaLower) {
                for(j=0;j<N;j++)
                    for(i=j;i<M;i++)
                        B[j+i*LDB] = A[i+j*LDA];
            }
            else {
                for(j=0;j<N;j++)
                    for(i=0;i<M;i++)
                        B[j+i*LDB] = A[i+j*LDA];
            }
        }
    }

    return PLASMA_SUCCESS;
}


text:
/*
 * TableFunction.h
 *
 * Author:
 *       Oleg Kalashev
 *
 * Copyright (c) 2020 Institute for Nuclear Research, RAS
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

#ifndef TABLEFUNCTION_H_
#define TABLEFUNCTION_H_

#include <iostream>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

#include "MathUtils.h"
#include "TableReader.h"

namespace Utils {

template<typename X = double >
class IScale
{
public:
	virtual X ToScale(X aX) = 0;
	virtual X FromScale(X aXprime) = 0;
	virtual IScale<X>* Clone() const = 0;
	virtual ~IScale(){}
};

template<typename X = double >
class LogScaleX : public IScale<X>
{
public:
	X ToScale(X aX);
	X FromScale(X aXprime);
	IScale<X>* Clone() const;
};

typedef LogScaleX<double> LogScale;

template<typename X = double >
class TableFunctionX : public SmartReferencedObj, virtual public FunctionX<X>{
public:
	static int FindLeftX(const std::vector<X>&, X xValue);
	TableFunctionX(const std::vector<X>& _x, const std::vector<X>& _y, X _leftVal=0., X _rightVal=0.);

	//takes ownership of aXscale and aYscale
	TableFunctionX(const std::vector<X>& _x, const std::vector<X>& _y, IScale<X>* aXscale, IScale<X>* aYscale, X _leftVal=0., X _rightVal=0.);

	//stores smart ref to aReader
	TableFunctionX(TableReaderX<X>* aReader, X _leftVal=0., X _rightVal=0.);

	bool InTableRange(X aArg) const;

	virtual X f(X _x) const;
	virtual X f_scaled(X _x) const = 0;
	virtual X Xmin() const;
	virtual X Xmax() const;

	void SetAutoLimits();
	void Print(std::ostream& aOut, X aXcoef=1., X aYcoef=1.);
private:
	//common part of constructors
	void Init();
protected:
	int FindLeftX(X xValue) const;
	//used for cloning
	TableFunctionX(const TableFunctionX<X>& aTableFunction);

	std::vector<X>			fCloneX;
	std::vector<X>			fCloneY;
	const std::vector<X>&	m_x;
	const std::vector<X>&	m_y;
	X						m_leftVal;
	X						m_rightVal;
	SafePtr<IScale<X> >		fXscale;
	SafePtr<IScale<X> >		fYscale;
	SmartPtr<TableReaderX<X> >	fTable;
	X						fXmin;
	X						fXmax;
};

typedef TableFunctionX<double> TableFunction;

template<typename X = double >
class LinearFuncX  : public TableFunctionX<X>
{
public:
	LinearFuncX(const std::vector<X>& _x,const  std::vector<X>& _y, X _leftVal=0., X _rightVal=0.):
	  TableFunctionX<X>(_x, _y, _leftVal, _rightVal)  {};

	//stores smart ref to aReader
	LinearFuncX(TableReaderX<X>* aReader, X _leftVal=0., X _rightVal=0.):
		TableFunctionX<X>(aReader,_leftVal,_rightVal) {};

	//takes ownership of aXscale and aYscale
	LinearFuncX(const std::vector<X>& _x,const std::vector<X>& _y, IScale<X>* aXscale, IScale<X>* aYscale, X _leftVal=0., X _rightVal=0.):
		TableFunctionX<X>(_x, _y, aXscale, aYscale, _leftVal, _rightVal) {}

	X f_scaled(X _x) const;
	FunctionX<X>* Clone() const { return new LinearFuncX<X>(*this); }
protected:
	LinearFuncX(const TableFunctionX<X>& aTableFunction):
		TableFunctionX<X>(aTableFunction) {};
};

typedef LinearFuncX<double> LinearFunc;

class GSLTableFunc  : public TableFunction
{
public:
	GSLTableFunc(const Vector& _x, const Vector& _y, const gsl_interp_type * aInterpType=gsl_interp_linear,
			double _leftVal=0., double _rightVal=0.):
	  TableFunction(_x, _y, _leftVal, _rightVal),fAcc(0),fSpline(0) {Init(aInterpType);}

	//stores smart ref to aReader
	GSLTableFunc(TableReader *aReader, const gsl_interp_type * aInterpType=gsl_interp_linear, double _leftVal=0., double _rightVal=0.):
		TableFunction(aReader,_leftVal,_rightVal),fAcc(0),fSpline(0) {Init(aInterpType);}

	//takes ownership of aXscale and aYscale
	//_x and _y are values are modified with scale functions given by aXscale and aYscale parameters
	GSLTableFunc(const Vector& _x, const Vector& _y, IScale<double>* aXscale, IScale<double>* aYscale, const gsl_interp_type * aInterpType=gsl_interp_linear, double _leftVal=0., double _rightVal=0.):
		TableFunction(_x, _y, aXscale, aYscale, _leftVal, _rightVal) {Init(aInterpType);}

	virtual ~GSLTableFunc();

	double f_scaled(double _x) const;
	Function* Clone() const { return new GSLTableFunc(*this); }
protected:
	GSLTableFunc(const GSLTableFunc& aTableFunction):
		TableFunction(aTableFunction) {Init(aTableFunction.fSpline->interp->type);};
private:
	void Init(const gsl_interp_type * aInterpType);
    gsl_interp_accel *fAcc;
    gsl_spline *fSpline;
    const gsl_interp_type *fInterpType;
};

template<typename X = double >
class MatrixFunctionX : public Function2X<X>
{
public:
	MatrixFunctionX(const std::string& aFile)
	{
		std::ifstream dataFile(aFile.c_str());
		if(!dataFile)
			Exception::Throw("Failed to open " + aFile);
		std::string header;
		std::getline(dataFile, header);
		int logscX,logscY;
		const char* headerFormat = "# minX=%lg stepX=%lg logscaleX=%d minY=%lg stepY=%lg logscaleY=%d";
		if(sscanf(header.c_str(),headerFormat,
				  &xMin,&xStep,&logscX,&yMin,&yStep,&logscY)!=6)
			Exception::Throw("Invalid header format in " + aFile + "\nExpected format: " + headerFormat);
		logScaleX = (bool)logscX;
		logScaleY = (bool)logscY;
		if((logScaleX && (xMin<=0 || xStep<=1.))||xStep<=0.)
			Exception::Throw("Invalid X scale in " + aFile);
		if((logScaleY && (yMin<=0 || yStep<=1.))||yStep<=0.)
			Exception::Throw("Invalid Y scale in " + aFile);
		fData = new TableReaderX<X>(dataFile, TableReader::Auto);
		nCols = fData->numberOfColumns();
		if(logScaleX) {
			xMax = xMin*pow(xStep,nCols-1);
			xStep = log(xStep);//convert multiplier to log step
		}
		else{
			xMax = xMin + xStep*(nCols-1);
		}
		nRows = fData->getColumn(0).size();
		if(logScaleY){
			yMax = yMin*pow(yStep,nRows-1);
			yStep = log(yStep);//convert multiplier to log step
		}
		else{
			yMax = yMin + yStep*(nRows-1);
		}
	}

	X f(double x, double y) const
	{
		double binX, binY;

		if(logScaleX)
		{//log scale
			binX = log(x/xMin)/xStep;
		}
		else
		{//linear scale
			binX = (x-xMin)/xStep;
		}
		if(binX<0. || binX>nCols-1)
			return 0.;
		if(logScaleY)
		{//log scale
			binY = log(y/yMin)/yStep;
		}
		else
		{//linear scale
			binY = (y-yMin)/yStep;
		}
		if(binY<0. || binY>nRows-1)
			return 0.;
		int iX = floor(binX);
		int iY = floor(binY);
		double weightX = 1.-binX+iX;
		double weightY = 1.-binY+iY;

		//linear f interpolation //TODO implement logscale f interpolation
		double f_11 = fData->getColumn(iX)[iY]*weightX*weightY;
		double f_21 = (weightX<1.)? (fData->getColumn(iX+1)[iY]*(1.-weightX)*weightY) : 0.;
		double f_12 = (weightY<1.)? (fData->getColumn(iX)[iY+1]*weightX*(1.-weightY)) : 0.;
		double f_22 = (weightX<1.&&weightY<1.)?(fData->getColumn(iX+1)[binY+1]*(1.-weightX)*(1.-weightY)) : 0.;
		double result = f_11 + f_21 + f_12 + f_22;
		return result;
	}

	virtual X MinArg(int aArgNo) const {return aArgNo==1?xMin:yMin;};
	virtual X MaxArg(int aArgNo) const {return aArgNo==1?xMax:yMax;}

private:
	SafePtr<TableReaderX<X> >  fData;
	X xMin;
	X yMin;
	X xMax;
	X yMax;
	X xStep;
	X yStep;
	int nCols;
	int nRows;
	bool logScaleX;
	bool logScaleY;
};

typedef MatrixFunctionX<double> MatrixFunction;

} /* namespace Utils */

#endif /* TABLEFUNCTION_H_ */

max_stars_repo_name: alexkorochkin/mcray

max_issues_repo_name: alexkorochkin/mcray

max_forks_repo_name: alexkorochkin/mcray


text:
#ifndef VMICORE_VMIINITDATA_H
#define VMICORE_VMIINITDATA_H

#include <filesystem>
#include <gsl/pointers>
#include <libvmi/libvmi.h>

class VmiInitData
{
  public:
    gsl::owner<vmi_init_data_t*> data{};

    explicit VmiInitData(const std::filesystem::path& socketPath);

    ~VmiInitData();
};

#endif // VMICORE_VMIINITDATA_H

max_stars_repo_name: GDATASoftwareAG/smartvmi

max_issues_repo_name: GDATASoftwareAG/smartvmi

max_forks_repo_name: GDATASoftwareAG/smartvmi


text:
/*
 * Copyright 2020 Makani Technologies LLC
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include "common/c_math/linalg.h"

#include <assert.h>
#include <float.h>
#include <math.h>
#include <stddef.h>
#include <stdint.h>

#include <gsl/gsl_cblas.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_vector.h>

#include "common/c_math/gsl_linalg_extra.h"
#include "common/c_math/linalg_common.h"

#if defined(__arm__)
// Turn off any error handling by the GSL library.  The standard way
// of doing this is with gsl_set_error_handler_off; however this
// method does not require bringing in more components of GSL.
void gsl_error(const char *reason, const char *file, int line, int gsl_errno) {
  (void)reason;
  (void)file;
  (void)line;
  (void)gsl_errno;
}
#endif  // defined(__arm__)

// Operations on vectors with any number of elements.

bool VecIsSize(const Vec *v, int32_t l) { return v != NULL && l == v->length; }

double *VecPtr(const Vec *v, int32_t i) {
  assert(0 <= i && i < v->length);
  return &v->d[i];
}

double VecGet(const Vec *v, int32_t i) {
  assert(0 <= i && i < v->length);
  return v->d[i];
}

// v = [0, 0, ... ]
const Vec *VecZero(Vec *v) {
  for (int32_t i = 0; i < v->length; ++i) {
    v->d[i] = 0.0;
  }
  return v;
}

// v_out = a*v + v_out
const Vec *VecAxpy(double a, const Vec *v, Vec *v_out) {
  assert(v->length == v_out->length);
  for (int32_t i = 0; i < v->length; ++i) {
    v_out->d[i] += a * v->d[i];
  }
  return v_out;
}

// Changes the length of a vector.  Returns true on success, false
// otherwise.
bool VecResize(int32_t length, Vec *v) {
  if (v->length == length) {
    return true;
  } else if (v->variable_len && 0 <= length && length <= v->max_length) {
    v->length = length;
    return true;
  }
  assert(false);
  return false;
}

// v_out = v_in
const Vec *VecCopy(const Vec *v_in, Vec *v_out) {
  VecResize(v_in->length, v_out);
  for (int32_t i = 0; i < v_out->length; ++i) {
    v_out->d[i] = v_in->d[i];
  }
  return v_out;
}

// Copy the contents of an array in a Vec.
const Vec *VecInit(const double a_in[], int32_t length, Vec *v_out) {
  VecResize(length, v_out);
  for (int32_t i = 0; i < v_out->length; ++i) {
    v_out->d[i] = a_in[i];
  }
  return v_out;
}

// v_out = scale * v_in
const Vec *VecScale(const Vec *v_in, double scale, Vec *v_out) {
  VecResize(v_in->length, v_out);
  for (int32_t i = 0; i < v_in->length; ++i) {
    v_out->d[i] = scale * v_in->d[i];
  }
  return v_out;
}

// v_out = v0 + v1
const Vec *VecAdd(const Vec *v0, const Vec *v1, Vec *v_out) {
  assert(v0->length == v1->length);
  VecResize(v0->length, v_out);
  for (int32_t i = 0; i < v_out->length; ++i) {
    v_out->d[i] = v0->d[i] + v1->d[i];
  }
  return v_out;
}

// v_out = v0 + v1 + v2
const Vec *VecAdd3(const Vec *v0, const Vec *v1, const Vec *v2, Vec *v_out) {
  assert(v0->length == v1->length && v0->length == v2->length);
  VecResize(v0->length, v_out);
  for (int32_t i = 0; i < v_out->length; ++i) {
    v_out->d[i] = v0->d[i] + v1->d[i] + v2->d[i];
  }
  return v_out;
}

// v_out = v0 - v1
const Vec *VecSub(const Vec *v0, const Vec *v1, Vec *v_out) {
  assert(v0->length == v1->length);
  VecResize(v0->length, v_out);
  for (int32_t i = 0; i < v_out->length; ++i) {
    v_out->d[i] = v0->d[i] - v1->d[i];
  }
  return v_out;
}

// v_out = c0 * v0 + c1 * v1
const Vec *VecLinComb(double c0, const Vec *v0, double c1, const Vec *v1,
                      Vec *v_out) {
  assert(v0->length == v1->length);
  VecResize(v0->length, v_out);
  for (int32_t i = 0; i < v_out->length; ++i) {
    v_out->d[i] = c0 * v0->d[i] + c1 * v1->d[i];
  }
  return v_out;
}

// v_out = c0 * v0 + c1 * v1 + c2 * v2
const Vec *VecLinComb3(double c0, const Vec *v0, double c1, const Vec *v1,
                       double c2, const Vec *v2, Vec *v_out) {
  assert(v0->length == v1->length && v0->length == v2->length);
  VecResize(v0->length, v_out);
  for (int32_t i = 0; i < v_out->length; ++i) {
    v_out->d[i] = c0 * v0->d[i] + c1 * v1->d[i] + c2 * v2->d[i];
  }
  return v_out;
}

// v_out = v0 .* v1
const Vec *VecMult(const Vec *v0, const Vec *v1, Vec *v_out) {
  assert(v0->length == v1->length);
  VecResize(v0->length, v_out);
  for (int32_t i = 0; i < v_out->length; ++i) {
    v_out->d[i] = v0->d[i] * v1->d[i];
  }
  return v_out;
}

// v0^T * v1
double VecDot(const Vec *v0, const Vec *v1) {
  assert(v0->length == v1->length);
  double tmp = 0.0;
  for (int32_t i = 0; i < v0->length; ++i) {
    tmp += v0->d[i] * v1->d[i];
  }
  return tmp;
}

// Computes the 2-norm of a vector.
//
// TODO: This does not yet handle overflow and underflow
// gracefully.
double VecNorm(const Vec *v) { return sqrt(VecNormSquared(v)); }

double VecNormBound(const Vec *v, double low) { return fmax(low, VecNorm(v)); }

double VecNormSquared(const Vec *v) { return VecDot(v, v); }

const Vec *VecNormalize(const Vec *v_in, Vec *v_out) {
  return VecScale(v_in, 1.0 / VecNormBound(v_in, 1e-9), v_out);
}

// Return the elements of v_in where inds != 0. The MATLAB logical indexing
// equivalent is: v_out = v_in(inds).
const Vec *VecSlice(const Vec *v_in, const int32_t inds[], Vec *v_out) {
  int32_t i_out = 0;
  for (int32_t i = 0; i < v_in->length; ++i) {
    if (inds[i]) {
      assert(i_out < v_out->max_length);
      // Don't use VecPtr here because we haven't set v_out->length yet.
      v_out->d[i_out] = VecGet(v_in, i);
      i_out++;
    }
  }

  // Do not resize v_out until the end of the function if you want to
  // be able to reuse and input as an output.
  VecResize(i_out, v_out);
  return v_out;
}

// Assigns elements of a larger vector with those from a smaller one. The
// MATLAB logical indexing equivalent is: v_out(inds) = v_in. The sum(inds != 0)
// must be equal to length(v_in). Warning: there is no check that
// length(inds) == length(v_out).
const Vec *VecSliceSet(const Vec *v_in, const int32_t inds[], Vec *v_out) {
  assert(v_in != v_out);

  int32_t i_in = 0;
  for (int32_t i = 0; i < v_out->length; ++i) {
    if (inds[i]) {
      *VecPtr(v_out, i) = VecGet(v_in, i_in);
      i_in++;
    }
  }

  // Assert that the sum(inds != 0) is equal to length(v_in).
  assert(v_in->length == i_in);

  return v_out;
}

// Operations on arbitrary matrices.

bool MatIsSize(const Mat *m, int32_t nr, int32_t nc) {
  return m != NULL && m->nr == nr && m->nc == nc;
}

double *MatPtr(const Mat *m, int32_t i, int32_t j) {
  assert(i < m->nr && j < m->nc);
  return &m->d[i * m->nc + j];
}

double MatGet(const Mat *m, int32_t i, int32_t j) {
  assert(i < m->nr && j < m->nc);
  return m->d[i * m->nc + j];
}

bool MatResize(int32_t nr, int32_t nc, Mat *m) {
  if (m->nr == nr && m->nc == nc) {
    return true;
  } else if (m->variable_dim && nr * nc <= m->max_size) {
    m->nr = nr;
    m->nc = nc;
    return true;
  }
  assert(false);
  return false;
}

// Copy the contents of an array into a Mat.
const Mat *MatInit(const double *d, int32_t nr, int32_t nc, Mat *m) {
  MatResize(nr, nc, m);
  for (int32_t i = 0; i < m->nr; ++i) {
    for (int32_t j = 0; j < m->nc; ++j) {
      *MatPtr(m, i, j) = d[i * m->nc + j];
    }
  }
  return m;
}

// m_out = scale * m_in
const Mat *MatScale(const Mat *m_in, double scale, Mat *m_out) {
  MatResize(m_in->nr, m_in->nc, m_out);
  for (int32_t i = 0; i < m_out->nr; ++i) {
    for (int32_t j = 0; j < m_out->nc; ++j) {
      *MatPtr(m_out, i, j) = scale * MatGet(m_in, i, j);
    }
  }
  return m_out;
}

const Mat *MatZero(Mat *m) {
  MatArrZero(m->nr, m->nc, m->d);
  return m;
}

const Mat *MatCopy(const Mat *a, Mat *b) {
  if (a != b) {  // No work to do if the Mat pointers are equal.
    MatResize(a->nr, a->nc, b);
    // This is still safe if the data pointers are equal.
    MatArrCopy(a->d, a->nr, a->nc, b->d);
  }
  return b;
}

// Copies the sub-matrix src(start_row + (1:num_rows), start_col + (1:num_cols))
// into dest(dest_row + (1:num_rows), dest_col + (1:num_cols)).
//
// It is not safe for src and dest to share any memory.
//
// Returns:
//   A pointer to dest.
const Mat *MatSubmatSet(const Mat *src, int32_t start_row, int32_t start_col,
                        int32_t num_rows, int32_t num_cols, int32_t dest_row,
                        int32_t dest_col, Mat *dest) {
  assert(src != dest && src->d != dest->d);
  assert(start_row >= 0 && start_col >= 0);
  assert(num_rows >= 0 && num_cols >= 0);
  assert(dest_row >= 0 && dest_col >= 0);
  // Make sure there is sufficient data and sufficient room.
  assert(src->nr >= start_row + num_rows && src->nc >= start_col + num_cols);
  assert(dest->nr >= dest_row + num_rows && dest->nc >= dest_col + num_cols);

  for (int32_t i = 0; i < num_rows; ++i) {
    for (int32_t j = 0; j < num_cols; ++j) {
      *MatPtr(dest, dest_row + i, dest_col + j) =
          MatGet(src, i + start_row, j + start_col);
    }
  }
  return dest;
}

const Mat *MatI(int32_t n, Mat *m) {
  MatResize(n, n, m);
  MatArrI(n, m->d);
  return m;
}

const Mat *MatMult(const Mat *m1, const Mat *m2, Mat *m_out) {
  assert(m1->nc == m2->nr);
  MatResize(m1->nr, m2->nc, m_out);
  MatArrMult(m1->d, m1->nr, m1->nc, m2->d, m2->nc, m_out->d);
  return m_out;
}

// Generalized matrix-vector multiply.
//
//   z = alpha*op(a)*x + beta*y.
//
// Args:
//   ta: Indicates if a should be transposed.
//   alpha, a, x, beta: See formula above.
//   y: Can be NULL if beta = 0.0.
//   z: Resulting vector.  It is safe to reuse y as z.
const Vec *MatVecGenMult(TransposeType ta, double alpha, const Mat *a,
                         const Vec *x, double beta, const Vec *y, Vec *z) {
  assert(a != NULL && x != NULL && z != NULL);
  assert((ta == kNoTrans ? a->nc : a->nr) == x->length);
  const int32_t zl = (ta == kNoTrans) ? a->nr : a->nc;
  assert(y == NULL || y->length == zl);
  assert(y != NULL || beta == 0.0);
  if (y == NULL) beta = 0.0;

  if (y == NULL || beta == 0.0) {
    VecResize(zl, z);
  } else {
    VecCopy(y, z);
  }

  MatArrGemv(ta, alpha, a->d, a->nr, a->nc, x->d, beta, z->d);
  return z;
}

// Generalized matrix-matrix multiply.
//
//   d = alpha*op(a)*op(b) + beta*c.
//
// Args:
//   ta: Indicates if a should be transposed.
//   tb: Indicates if b should be transposed.
//   alpha, a, b, beta: See above formula.
//   c: Can be NULL if beta = 0.0.
//   d: Resulting matrix.  It is safe to reuse c as d.
const Mat *MatGenMult(TransposeType ta, TransposeType tb, double alpha,
                      const Mat *a, const Mat *b, double beta, const Mat *c,
                      Mat *d) {
  assert(a != NULL && b != NULL && d != NULL);
  assert((ta == kNoTrans ? a->nc : a->nr) == (tb == kNoTrans ? b->nr : b->nc));
  const int32_t nr_d = (ta == kNoTrans ? a->nr : a->nc);
  const int32_t nc_d = (tb == kNoTrans ? b->nc : b->nr);
  assert(c == NULL || (c->nr == nr_d && c->nc == nc_d));
  assert(c != NULL || beta == 0.0);
  if (c == NULL) beta = 0.0;

  if (c == NULL || beta == 0.0)
    MatResize(nr_d, nc_d, d);
  else
    MatCopy(c, d);

  MatArrGemm(ta, tb, alpha, a->d, a->nr, a->nc, b->d, nc_d, beta, d->d);
  return d;
}

const Mat *MatMult3(const Mat *m0, const Mat *m1, const Mat *m2, Mat *m_out) {
  assert(m0->nc == m1->nr && m1->nc == m2->nr);
  MatResize(m0->nr, m2->nc, m_out);
  MAT(m0->nr, m1->nc, m0m1);
  return MatMult(MatMult(m0, m1, &m0m1), m2, m_out);
}

const Vec *MatVecMult(const Mat *m, const Vec *v_in, Vec *v_out) {
  assert(m->nc == v_in->length);
  VecResize(m->nr, v_out);
  MatArrMult(m->d, m->nr, m->nc, v_in->d, 1, v_out->d);
  return v_out;
}

const Vec *MatTransVecMult(const Mat *m, const Vec *v_in, Vec *v_out) {
  assert(m->nr == v_in->length);
  VecResize(m->nc, v_out);
  MatArrGemm(kTrans, kNoTrans, 1.0, m->d, m->nr, m->nc, v_in->d, 1, 0.0,
             v_out->d);
  return v_out;
}

const Mat *MatTrans(const Mat *m, Mat *m_out) {
  MatResize(m->nc, m->nr, m_out);
  MatArrTrans(m->d, m->nr, m->nc, m_out->d);
  return m_out;
}

const Mat *MatAdd(const Mat *m0, const Mat *m1, Mat *m_out) {
  assert(m0->nr == m1->nr && m0->nc == m1->nc);
  MatResize(m0->nr, m0->nc, m_out);
  for (int32_t i = 0; i < m0->nr; ++i) {
    for (int32_t j = 0; j < m0->nc; ++j) {
      *MatPtr(m_out, i, j) = MatGet(m0, i, j) + MatGet(m1, i, j);
    }
  }
  return m_out;
}

const Mat *MatSub(const Mat *m0, const Mat *m1, Mat *m_out) {
  assert(m0->nr == m1->nr && m0->nc == m1->nc);
  MatResize(m0->nr, m0->nc, m_out);
  for (int32_t i = 0; i < m0->nr; ++i) {
    for (int32_t j = 0; j < m0->nc; ++j) {
      *MatPtr(m_out, i, j) = MatGet(m0, i, j) - MatGet(m1, i, j);
    }
  }
  return m_out;
}

void MatQrDecomp(const Mat *A, Mat *Q, Mat *R) {
  MatResize(A->nr, A->nc, R);
  double *Qd = NULL;
  if (Q != NULL) {
    MatResize(A->nr, A->nr, Q);
    Qd = Q->d;
  }
  MatArrQrDecomp(A->d, A->nr, A->nc, Qd, R->d);
}

// Calculates the "thin" singular value decomposition of an m x n
// matrix A, defined by:
//
//   A = U * S * V'
//
// where U is an m x n column-orthogonal matrix (i.e. U'*U = I), S is
// an n x n diagonal matrix with non-decreasing elements, and V is an
// n x n orthogonal matrix (i.e. V'*V = V*V' = I).  The "thin" aspect
// is equivalent to MATLAB's "economy" SVD.
//
// Args:
//   A: Input m x n matrix with m >= n.
//   U: Output m x n column-orthogonal matrix.
//   s: Output n length vector of the diagonal elements of S.
//   V: Output n x n orthogonal matrix.
//
// Returns:
//   kLinalgErrorNone on success and kLinalgErrorMaxIter if GSL's
//   solver reached its maximum number of iterations.
int32_t MatThinSvDecomp(const Mat *A, Mat *U, Vec *s, Mat *V) {
  assert(A != NULL && U != NULL && s != NULL && V != NULL);
  assert(A->nc <= VEC_MAX_ELEMENTS);
  assert(A->nr >= A->nc);

  MatCopy(A, U);
  VecResize(A->nc, s);
  MatResize(A->nc, A->nc, V);

  double work_data[VEC_MAX_ELEMENTS];
  gsl_matrix_view U_gsl =
      gsl_matrix_view_array(U->d, (size_t)U->nr, (size_t)U->nc);
  gsl_matrix_view V_gsl =
      gsl_matrix_view_array(V->d, (size_t)V->nr, (size_t)V->nc);
  gsl_vector_view s_gsl = gsl_vector_view_array(s->d, (size_t)s->length);
  gsl_vector_view work = gsl_vector_view_array(work_data, (size_t)s->length);
  int32_t ret = gsl_linalg_SV_decomp(&U_gsl.matrix, &V_gsl.matrix,
                                     &s_gsl.vector, &work.vector);
  // These are the only errors that gsl_linalg_SV_decomp should ever
  // return, assuming we're passing in matrices and vectors of the
  // appropriate sizes.
  assert(ret == GSL_SUCCESS || ret == GSL_EMAXITER);

  if (ret == GSL_EMAXITER) {
    assert(false);
    return kLinalgErrorMaxIter;
  }
  return kLinalgErrorNone;
}

// Calculates the rank of a matrix.  Performs a QR decomposition (with
// pivoting) of the input matrix and finds the number of diagonal
// elements of R with magnitude greater than a specified tolerance.
//
// Args:
//   A: Input matrix.
//   tol: Tolerance for checking linear independence of vectors.
//
// Returns:
//   Rank of the matrix.
int32_t MatRank(const Mat *A, double tol) {
  assert(A != NULL);
  assert(A->nr * A->nc <= MAT_MAX_ELEMENTS);
  assert(A->nr <= VEC_MAX_ELEMENTS && A->nc <= VEC_MAX_ELEMENTS);
  assert(tol >= DBL_EPSILON);

  int32_t min_dim = (A->nr < A->nc) ? A->nr : A->nc;
  if (min_dim == 0) return 0;

  double qr_data[MAT_MAX_ELEMENTS];
  double tau_data[VEC_MAX_ELEMENTS];
  double norm_data[VEC_MAX_ELEMENTS];
  size_t perm_data[VEC_MAX_ELEMENTS];
  int signum;
  gsl_matrix_const_view A_gsl =
      gsl_matrix_const_view_array(A->d, (size_t)A->nr, (size_t)A->nc);
  gsl_matrix_view QR =
      gsl_matrix_view_array(qr_data, (size_t)A->nr, (size_t)A->nc);
  gsl_vector_view tau = gsl_vector_view_array(tau_data, (size_t)min_dim);
  gsl_permutation p = {(size_t)A->nc, perm_data};
  gsl_vector_view norm = gsl_vector_view_array(norm_data, (size_t)A->nc);

  gsl_matrix_memcpy(&QR.matrix, &A_gsl.matrix);
  int32_t ret = gsl_linalg_QRPT_decomp(&QR.matrix, &tau.vector, &p, &signum,
                                       &norm.vector);
  assert(ret == GSL_SUCCESS);
  (void)ret;

  int32_t rank = 0;
  while (rank < min_dim) {
    if (fabs(gsl_matrix_get(&QR.matrix, (size_t)rank, (size_t)rank)) <= tol) {
      break;
    }
    rank++;
  }

  return rank;
}

// Solves A*x = b for x similar to x = A \ b.
int32_t MatVecLeftDivide(const Mat *A, const Vec *b, Vec *x) {
  assert(A != NULL && b != NULL && x != NULL);
  assert(b != x);
  assert(A->nr == b->length && A->nr > 0);
  VecResize(A->nc, x);
  if (A->nc == 0) return kLinalgErrorNone;

  size_t m = (size_t)A->nr, n = (size_t)A->nc;
  gsl_matrix_const_view A_gsl = gsl_matrix_const_view_array(A->d, m, n);
  gsl_matrix_const_view b_gsl = gsl_matrix_const_view_array(b->d, m, 1U);
  gsl_matrix_view x_gsl = gsl_matrix_view_array(x->d, n, 1U);
  int32_t ret =
      GslMatrixDivide(CblasLeft, &A_gsl.matrix, &b_gsl.matrix, &x_gsl.matrix);
  assert(ret == GSL_SUCCESS || ret == GSL_ESING);
  return (ret == GSL_ESING) ? kLinalgErrorSingularMat : kLinalgErrorNone;
}

// Solves x*A = b for x similar to x = b / A.
int32_t MatVecRightDivide(const Mat *A, const Vec *b, Vec *x) {
  assert(A != NULL && b != NULL && x != NULL);
  assert(b != x);
  assert(A->nc == b->length && A->nc > 0);
  VecResize(A->nr, x);
  if (A->nr == 0) return kLinalgErrorNone;

  size_t m = (size_t)A->nr, n = (size_t)A->nc;
  gsl_matrix_const_view A_gsl = gsl_matrix_const_view_array(A->d, m, n);
  gsl_matrix_const_view b_gsl = gsl_matrix_const_view_array(b->d, 1U, n);
  gsl_matrix_view x_gsl = gsl_matrix_view_array(x->d, 1U, m);
  int32_t ret =
      GslMatrixDivide(CblasRight, &A_gsl.matrix, &b_gsl.matrix, &x_gsl.matrix);
  assert(ret == GSL_SUCCESS || ret == GSL_ESING);
  return (ret == GSL_ESING) ? kLinalgErrorSingularMat : kLinalgErrorNone;
}

// Solves A*X = b for X similar to X = A \ B.
int32_t MatMatLeftDivide(const Mat *A, const Mat *B, Mat *X) {
  assert(A != NULL && B != NULL && X != NULL);
  assert(A != X && B != X);
  assert(A->nr == B->nr && A->nr > 0);
  MatResize(A->nc, B->nc, X);
  if (A->nc == 0) return kLinalgErrorNone;

  size_t m = (size_t)A->nr, n = (size_t)A->nc, k = (size_t)B->nc;
  gsl_matrix_const_view A_gsl = gsl_matrix_const_view_array(A->d, m, n);
  gsl_matrix_const_view B_gsl = gsl_matrix_const_view_array(B->d, m, k);
  gsl_matrix_view X_gsl = gsl_matrix_view_array(X->d, n, k);
  int32_t ret =
      GslMatrixDivide(CblasLeft, &A_gsl.matrix, &B_gsl.matrix, &X_gsl.matrix);
  assert(ret == GSL_SUCCESS || ret == GSL_ESING);
  return (ret == GSL_ESING) ? kLinalgErrorSingularMat : kLinalgErrorNone;
}

// Solves X*A = B for X similar to X = B / A.
int32_t MatMatRightDivide(const Mat *A, const Mat *B, Mat *X) {
  assert(A != NULL && B != NULL && X != NULL);
  assert(A != X && B != X);
  assert(A->nc == B->nc && A->nc > 0);
  MatResize(B->nr, A->nr, X);
  if (A->nr == 0) return kLinalgErrorNone;

  size_t m = (size_t)A->nr, n = (size_t)A->nc, k = (size_t)B->nr;
  gsl_matrix_const_view A_gsl = gsl_matrix_const_view_array(A->d, m, n);
  gsl_matrix_const_view B_gsl = gsl_matrix_const_view_array(B->d, k, n);
  gsl_matrix_view X_gsl = gsl_matrix_view_array(X->d, k, m);
  int32_t ret =
      GslMatrixDivide(CblasRight, &A_gsl.matrix, &B_gsl.matrix, &X_gsl.matrix);
  assert(ret == GSL_SUCCESS || ret == GSL_ESING);
  return (ret == GSL_ESING) ? kLinalgErrorSingularMat : kLinalgErrorNone;
}

bool MatIsUpperTriangular(const Mat *r) {
  return MatArrIsUpperTriangular(r->d, r->nr, r->nc);
}

bool MatIsLowerTriangular(const Mat *l) {
  return MatArrIsLowerTriangular(l->d, l->nr, l->nc);
}

bool MatHasNonnegDiag(const Mat *m) {
  assert(m != NULL);
  for (int32_t i = 0; i < m->nc && i < m->nr; ++i) {
    if (MatGet(m, i, i) < 0.0) return false;
  }
  return true;
}

// Calculates the matrix square root of P = A'*A + B'*B by calculating a QR
// decomposition, Q*R = [A ; B].  It is assumed that the matrix [A ; B] is full
// row rank (i.e. has linearly independent columns).
//
// Q is an orthogonal matrix and R "trapezoidal matrix", i.e. R = [sqrt_P ; 0]
// where sqrt_P is upper triangular.  Thus:
//   P = A'*A + B'*B = [A ; B]'*[A ; B] = [sqrT_P ; 0]'*Q'*Q*[sqrt_P ; 0].
//
// Args:
//   A: na-by-nc matrix.
//   B: nb-by-nc matrix.
//   sqrt_P: Upper-triangular matrix with positive diagonal elements such that
//     sqrt_P'*sqrt_P = P.
//
// Note: it is safe to reuse A or B as sqrt_P (assuming shapes match).
void MatSqrtSum(const Mat *A, const Mat *B, Mat *sqrt_P) {
  assert(A != NULL && B != NULL && sqrt_P != NULL);
  assert(A->nc == B->nc);
  assert(A->nr + B->nr >= A->nc);  // This is required for full row rank.
  MatResize(A->nc, A->nc, sqrt_P);

  // Define a temporary matrix, stack = [A ; B].
  MAT(A->nr + B->nr, A->nc, stack);
  MatSubmatSet(A, 0, 0, A->nr, A->nc, 0, 0, &stack);
  MatSubmatSet(B, 0, 0, B->nr, B->nc, A->nr, 0, &stack);
  MatQrDecomp(&stack, NULL, &stack);
  MatSubmatSet(&stack, 0, 0, A->nc, A->nc, 0, 0, sqrt_P);
  // Next, enforce a non-negative diagonal.
  for (int32_t i = 0; i < sqrt_P->nr; ++i) {
    if (MatGet(sqrt_P, i, i) < 0.0) {
      for (int32_t j = i; j < sqrt_P->nc; ++j) {
        *MatPtr(sqrt_P, i, j) *= -1.0;
      }
    }
  }
}

int32_t MatMatBackSub(const Mat *R, const Mat *b, Mat *x) {
  assert(b->nr == R->nr && b->nc == x->nc);
  assert(x->nr == R->nc);
  return MatArrBackSub(R->d, R->nr, R->nc, b->d, b->nc, x->d);
}

int32_t MatMatForwardSub(const Mat *L, const Mat *b, Mat *x) {
  assert(b->nr == L->nr && b->nc == x->nc);
  assert(x->nr == L->nc);
  return MatArrForwardSub(L->d, L->nr, L->nc, b->d, b->nc, x->d);
}

int32_t MatVecBackSub(const Mat *R, const Vec *b, Vec *x) {
  assert(b->length == R->nr);
  assert(x->length == R->nc);
  return MatArrBackSub(R->d, R->nr, R->nc, b->d, 1, x->d);
}

int32_t MatVecForwardSub(const Mat *L, const Vec *b, Vec *x) {
  assert(b->length == L->nr);
  assert(x->length == L->nc);
  return MatArrForwardSub(L->d, L->nr, L->nc, b->d, 1, x->d);
}

// Returns the rows and columns of m where rows != 0 and cols != 0
const Mat *MatSlice(const Mat *m, const int32_t rows[], const int32_t cols[],
                    Mat *m_out) {
  int32_t nr_out = 0, nc_out = 0;
  if (rows == NULL) {
    nr_out = m->nr;
  } else {
    for (int32_t i = 0; i < m->nr; ++i)
      if (rows[i]) nr_out++;
  }
  if (cols == NULL) {
    nc_out = m->nc;
  } else {
    for (int32_t i = 0; i < m->nc; ++i)
      if (cols[i]) nc_out++;
  }

  if (m_out->variable_dim) {
    assert(nr_out * nc_out <= m_out->max_size);
  } else {
    assert(m_out->nr == nr_out && m_out->nc == nc_out);
  }

  int32_t i_out = 0;
  for (int32_t i = 0; i < m->nr; ++i) {
    if (rows == NULL || rows[i]) {
      int32_t j_out = 0;
      for (int32_t j = 0; j < m->nc; ++j) {
        if (cols == NULL || cols[j]) {
          // Don't use MatPtr here because we haven't set m_out->nc yet.
          m_out->d[i_out * nc_out + j_out] = MatGet(m, i, j);
          j_out++;
        }
      }
      i_out++;
    }
  }

  // Do not resize m until the end of the function if you want to be
  // able to reuse the input as the output.  Otherwise, the matrix
  // indexing in MatGet will be off.
  MatResize(nr_out, nc_out, m_out);
  return m_out;
}

const double *MatArrMult(const double *m1, int32_t nr1, int32_t nc1,
                         const double *m2, int32_t nc2, double *m_out) {
  assert(m_out != m1 && m_out != m2);
  return MatArrGemm(kNoTrans, kNoTrans, 1.0, m1, nr1, nc1, m2, nc2, 0.0, m_out);
}

// This is modeled after the BLAS function GEMV, which calculates:
//     y <-- alpha*op(A)*x + beta*y,
// where op either returns the matrix or its transpose.
//
// Args:
//   ta: Whether to transpose A.
//   alpha: See above formula.
//   a: Data for A.
//   nra: Number of rows of A without transpose applied.
//   nca: Number of columns of A without transpose applied.
//   x: Data for x.
//   beta: See above formula.
//   c: Data for y.
const double *MatArrGemv(TransposeType ta, double alpha, const double *a,
                         int32_t nra, int32_t nca, const double *x, double beta,
                         double *y) {
  assert(x != y);
  int32_t ly = ((ta == kNoTrans) ? nra : nca);
  int32_t lx = ((ta == kNoTrans) ? nca : nra);
  double aij, axj;
  for (int32_t i = 0; i < ly; ++i) {  // For each entry of y.
    axj = 0.0;
    for (int32_t j = 0; j < lx; ++j) {  // Walk down the row of op(A).
      aij = (ta == kNoTrans) ? a[i * nca + j] : a[j * nca + i];
      axj += aij * x[j];
    }
    if (beta == 0.0) {
      y[i] = alpha * axj;
    } else {
      y[i] = alpha * axj + beta * y[i];
    }
  }
  return y;
}

// This is modeled after the BLAS function GEMM, which calculates:
//     C <-- alpha*op(A)*op(B) + beta*C,
// where op either returns the matrix or its transpose.
//
// Args:
//   ta: Whether to transpose A.
//   tb: Whether to transpose B.
//   alpha: See above formula.
//   a: Data for A.
//   nra: Number of rows of A without transpose applied.
//   nca: Number of columns of A without transpose applied.
//   b: Data for B.
//   ncc: Number of columns of C.
//   beta: See above formula.
//   c: Data for C.
const double *MatArrGemm(TransposeType ta, TransposeType tb, double alpha,
                         const double *a, int32_t nra, int32_t nca,
                         const double *b, int32_t ncc, double beta, double *c) {
  assert(a != c && b != c);
  int32_t nrc = (ta == kNoTrans) ? nra : nca;
  int32_t nk = (ta == kNoTrans) ? nca : nra;
  double ABij, aik, bkj;
  for (int32_t i = 0; i < nrc; ++i) {
    for (int32_t j = 0; j < ncc; ++j) {
      ABij = 0.0;
      for (int32_t k = 0; k < nk; ++k) {
        aik = (ta == kNoTrans) ? a[i * nca + k] : a[k * nca + i];
        bkj = (tb == kNoTrans) ? b[k * ncc + j] : b[j * nk + k];
        ABij += aik * bkj;
      }
      if (beta == 0.0) {
        c[i * ncc + j] = alpha * ABij;
      } else {
        c[i * ncc + j] = alpha * ABij + beta * c[i * ncc + j];
      }
    }
  }
  return c;
}

const double *MatArrZero(int32_t nr, int32_t nc, double *m_out) {
  for (int32_t i = 0; i < nr; ++i) {
    for (int32_t j = 0; j < nc; ++j) {
      m_out[i * nc + j] = 0.0;
    }
  }
  return m_out;
}

void MatArrCopy(const double *d_in, int32_t nr, int32_t nc, double *d_out) {
  for (int32_t i = 0; i < nr; ++i) {
    for (int32_t j = 0; j < nc; ++j) {
      d_out[i * nc + j] = d_in[i * nc + j];
    }
  }
}

void MatArrTrans(const double *m_in, int32_t nr, int32_t nc, double *m_out) {
  assert(m_out != m_in && m_out != NULL);
  for (int32_t i = 0; i < nr; ++i) {
    for (int32_t j = 0; j < nc; ++j) {
      m_out[j * nr + i] = m_in[i * nc + j];
    }
  }
}

const double *MatArrI(int32_t nr, double *m_out) {
  MatArrZero(nr, nr, m_out);
  for (int32_t i = 0; i < nr; ++i) {
    m_out[i * nr + i] = 1.0;
  }
  return m_out;
}

// QR decomposition.  Decompose a matrix into orthogonal and
// upper-right-triangular components.
//
// You may reuse a as r, and you may skip calculating q by passing a
// NULL pointer.
void MatArrQrDecomp(const double *a, int32_t nr, int32_t nc, double *q,
                    double *r) {
  assert(q != a && r != q && r != NULL);
  assert(nr <= VEC_MAX_ELEMENTS && nc <= VEC_MAX_ELEMENTS);

  double tau_data[VEC_MAX_ELEMENTS];
  gsl_matrix_const_view A =
      gsl_matrix_const_view_array(a, (size_t)nr, (size_t)nc);
  gsl_matrix_view QR = gsl_matrix_view_array(r, (size_t)nr, (size_t)nc);
  gsl_vector_view tau =
      gsl_vector_view_array(tau_data, (nr < nc) ? (size_t)nr : (size_t)nc);
  gsl_matrix_memcpy(&QR.matrix, &A.matrix);
  gsl_linalg_QR_decomp(&QR.matrix, &tau.vector);

  if (q != NULL) {
    gsl_matrix_view Q = gsl_matrix_view_array(q, (size_t)nr, (size_t)nr);
    gsl_matrix_view R = gsl_matrix_view_array(r, (size_t)nr, (size_t)nc);
    gsl_linalg_QR_unpack(&QR.matrix, &tau.vector, &Q.matrix, &R.matrix);
  } else {
    for (int32_t i = 1; i < nr; ++i) {
      for (int32_t j = 0; j < i && j < nc; ++j) {
        r[i * nc + j] = 0.0;
      }
    }
  }
}

bool MatArrIsUpperTriangular(const double *r, int32_t nr, int32_t nc) {
  for (int32_t i = 1; i < nr; ++i) {
    for (int32_t j = 0; j < i && j < nc; ++j) {
      if (fabs(r[i * nc + j]) >= DBL_EPSILON) {
        return false;
      }
    }
  }
  return true;
}

bool MatArrIsLowerTriangular(const double *l, int32_t nr, int32_t nc) {
  for (int32_t i = 0; i < nr; ++i) {
    for (int32_t j = i + 1; j < nc; ++j) {
      if (fabs(l[i * nc + j]) >= DBL_EPSILON) {
        return false;
      }
    }
  }
  return true;
}

// Solves R*X = B using back substitution where R is an upper
// triangular matrix.  If the triangular matrix contains a zero
// diagonal element (i.e. the matrix is singular), we set the
// corresponding x element to zero.
//
// Returns kLinalgErrorNone on success and kLinalgErrorSingularMat
// if the matrix was singular.
//
// Note: Confirmed safe to reuse b as x.
int32_t MatArrBackSub(const double *r, int32_t nr, int32_t nc, const double *b,
                      int32_t nc_b, double *x) {
  assert(r != x);
  assert(nr >= nc);
  assert(nr >= 0 && nc >= 0 && nc_b >= 0);
  assert(MatArrIsUpperTriangular(r, nr, nc));
  (void)nr;

  int32_t err = kLinalgErrorNone;
  // Solve for the k-th column of x.
  for (int32_t k = 0; k < nc_b; ++k) {
    // Iterate over the rows of r.
    for (int32_t i = nc - 1; i >= 0; --i) {
      if (fabs(r[i * nc + i]) <= DBL_EPSILON) {
        err = kLinalgErrorSingularMat;
        x[i * nc_b + k] = 0.0;
      } else {
        double sum = 0.0;
        for (int32_t j = nc - 1; j > i; --j) {
          sum += x[j * nc_b + k] * r[i * nc + j];
        }
        x[i * nc_b + k] = (b[i * nc_b + k] - sum) / r[i * nc + i];
      }
    }
    // TODO: If nr > nc, then the bottom of r is all zeros.
    // Should we test if b has non-zero elements for these rows?
  }
  return err;
}

// Solves L*X = B using forward substitution where L is a lower
// triangular matrix.  If the triangular matrix contains a zero
// diagonal element (i.e. the matrix is singular), we set the
// corresponding x element to zero.
//
// Returns kLinalgErrorNone on success and kLinalgErrorSingularMat
// if the matrix is singular.
//
// Note: Confirmed safe to reuse b as x.
int32_t MatArrForwardSub(const double *l, int32_t nr, int32_t nc,
                         const double *b, int32_t nc_b, double *x) {
  assert(nr <= nc);
  assert(MatArrIsLowerTriangular(l, nr, nc));
  int32_t err = kLinalgErrorNone;
  for (int32_t k = 0; k < nc_b; ++k) {  // Solve for the k-th column of x.
    for (int32_t i = 0; i < nr; ++i) {  // Iterate over the rows of l.
      if (fabs(l[i * nc + i]) <= DBL_EPSILON) {
        err = kLinalgErrorSingularMat;
        x[i * nc_b + k] = 0.0;
      } else {
        double sum = 0.0;
        for (int32_t j = 0; j < i; ++j) {
          sum += x[j * nc_b + k] * l[i * nc + j];
        }
        x[i * nc_b + k] = (b[i * nc_b + k] - sum) / l[i * nc + i];
      }
    }
    for (int32_t i = nr; i < nc; ++i) {
      x[i * nc_b + k] = 0.0;
    }
  }
  return err;
}


text:
/**
 *
 * @file testing_ssygv.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Azzam Haidar
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:45:18 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>
#include "testing_smain.h"

#undef COMPLEX
#define REAL

static int check_orthogonality(int itype, int uplo, int N,
                               float *Z, int LDZ,
                               float *B, float *CHOLB, int LDB,
                               float eps);
static int check_reduction(int itype, int uplo, int N, float *D,
                           float *A, int LDA,
                           float *B, int LDB,
                           float *Z, int LDZ,
                           float eps );
static int check_solution(int N, float *E1, float *E2, float eps);

int testing_ssygv(int argc, char **argv)
{
    /* Check for number of arguments*/
    if (argc != 3) {
        USAGE("HEGV", "N LDA LDB",
              "   - N    : size of the matrices A and B\n"
              "   - LDA  : leading dimension of the matrix A\n"
              "   - LDB  : leading dimension of the matrix B\n");
        return -1;
    }

    float      eps = LAPACKE_slamch_work('e');
    PLASMA_enum vec = PlasmaVec;
    int    N        = atoi(argv[0]);
    int    LDA      = atoi(argv[1]);
    int    LDB      = atoi(argv[2]);
    int    LDQ      = LDA;
    int    LDAxN    = LDA*N;
    int    LDBxN    = LDB*N;
    int    LDQxN    = LDQ*N;

    int info_ortho     = 0;
    int info_solution  = 0;
    int info_reduction = 0;
    int i, u;

    float *A1    = (float *)malloc(LDAxN*sizeof(float));
    float *A2    = (float *)malloc(LDAxN*sizeof(float));
    float *B1    = (float *)malloc(LDBxN*sizeof(float));
    float *B2    = (float *)malloc(LDBxN*sizeof(float));
    float *Q     = (float *)malloc(LDQxN*sizeof(float));
    float *Ainit = (float *)malloc(LDAxN*sizeof(float));
    float *Binit = (float *)malloc(LDBxN*sizeof(float));
    float *W1 = (float *)malloc(N*sizeof(float));
    float *W2 = (float *)malloc(N*sizeof(float));
    float *work = (float *)malloc(3*N* sizeof(float));
    PLASMA_desc *T;

    /* Check if unable to allocate memory */
    if ((!A1)||(!A2)||(!B1)||(!B2)||(!Q)||(!Ainit)||(!Binit)){
        printf("Out of Memory \n ");
        return -2;
    }

    PLASMA_Disable(PLASMA_AUTOTUNING);
    PLASMA_Set(PLASMA_TILE_SIZE, 120);
    PLASMA_Set(PLASMA_INNER_BLOCK_SIZE, 20);

    PLASMA_Enable(PLASMA_WARNINGS);
    PLASMA_Enable(PLASMA_ERRORS);
    PLASMA_Alloc_Workspace_ssygv(N, N, &T);

    /*----------------------------------------------------------
    *  TESTING SSYGV
    */
    /* Initialize A1 and Ainit */
    PLASMA_splgsy(0., N, A1, LDA, 5198);
    PLASMA_slacpy(PlasmaUpperLower, N, N, A1, LDA, Ainit, LDA);

    /* Initialize B1 and Binit */
    PLASMA_splgsy((float)N, N, B1, LDB, 4321 );
    PLASMA_slacpy(PlasmaUpperLower, N, N, B1, LDB, Binit, LDB);

    printf("\n");
    printf("------ TESTS FOR PLASMA SSYGV ROUTINE -------  \n");
    printf("        Size of the Matrix %d by %d\n", N, N);
    printf("\n");
    printf(" The matrix A is randomly generated for each test.\n");
    printf("============\n");
    printf(" The relative machine precision (eps) is to be %e \n",eps);
    printf(" Computational tests pass if scaled residuals are less than 60.\n");

    /*----------------------------------------------------------
     *  TESTING SSYGV
     */

    for (i=0; i<3; i++) {
        for (u=0; u<2; u++) {
            memcpy(A2, Ainit, LDAxN*sizeof(float));
            memcpy(B2, Binit, LDBxN*sizeof(float));
            /* CALL SSYGV with itype= 1, 2, 3 and uplo= L,U */
            PLASMA_ssygv(itype[i], vec, uplo[u], N, A2, LDA, B2, LDB, W2, T, Q, LDQ);

            /* CALL LAPACK to compute eigenvalues. */
            memcpy(A1, Ainit, LDAxN*sizeof(float));
            memcpy(B1, Binit, LDBxN*sizeof(float));
            LAPACKE_ssygv( LAPACK_COL_MAJOR,
                     itype[i], lapack_const(vec), lapack_const(uplo[u]),
                     N, A1, LDA, B1, LDB, W1);

            /* CHECK */
            info_solution = check_solution(N, W1, W2, eps);
            /* Check the orthogonality, reduction and the eigen solutions */
            if (vec == PlasmaVec){
                info_ortho = check_orthogonality(itype[i], uplo[u], N, Q, LDQ, Binit, B1, LDB, eps);
                info_reduction = check_reduction(itype[i], uplo[u], N, W2, Ainit, LDA, Binit, LDB, Q, LDQ, eps );
            }

            if ( (info_ortho == 0) & (info_reduction == 0) & (info_solution == 0)) {
                printf("***************************************************\n");
                printf(" ---- TESTING SSYGV (%s, %s) ...................... PASSED !\n", itypestr[i], uplostr[u]);
                printf("***************************************************\n");
            }
            else {
                printf("************************************************\n");
                printf(" - TESTING SSYGV (%s, %s) ... FAILED !\n", itypestr[i], uplostr[u]);
                printf("************************************************\n");
            }
        }
    }

    PLASMA_Dealloc_Handle_Tile(&T);
    free(A1);
    free(A2);
    free(B1);
    free(B2);
    free(Q);
    free(Ainit);
    free(Binit);
    free(W1);
    free(W2);
    free(work);

    return 0;
}

/*-------------------------------------------------------------------
 * Check the orthogonality of Z
 */
static int check_orthogonality(int itype, int uplo, int N,
                               float *Z, int LDZ,
                               float *B, float *CHOLB, int LDB,
                               float eps)
{
    static float zone  =  1.0;
    static float mzone = -1.0;
    static float zzero =  0.0;
    static float             done  =  1.0;
    static float             mdone = -1.0;
    PLASMA_enum trans;
    float  normQ, result;
    int     info_ortho;
    float *work = (float *)malloc(N*sizeof(float));
    float *TEMP    = (float *)malloc(N*N*sizeof(float));
    float *Id = (float *) malloc(N*N*sizeof(float));
    char *str;

    /* Build the idendity matrix */
    LAPACKE_slaset_work(LAPACK_COL_MAJOR, 'A', N, N, 0., 1., Id, N);

    /* Perform orth test*/
    if ( itype == 3 ) {
        /*
         * if ITYPE = 3, Z**T*inv(B)*Z = I.
         * inv(B) = inv(L)**T * inv(L) or inv(U)*inv(U)**T
         */
        str = "Z**T*inv(B)*Z";

        if( uplo==PlasmaUpper ) {
            trans = PlasmaTrans;
        } else {
            trans = PlasmaNoTrans;
        }

        /* Compute inv(L)*Z or inv(U)**T*Z */
        LAPACKE_slacpy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaUpperLower),
                            N, N, Z, LDZ, TEMP, N);

        cblas_strsm(CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo, (CBLAS_TRANSPOSE)trans,
                    CblasNonUnit, N, N,
                    (zone), CHOLB, LDB,
                                       TEMP,  N);

        /*
         * Compute Z**T*inv(B)*Z-Id = (Z**T*inv(L)**T) * (inv(L)*Z) - Id
         *
         * Note: Z**T*inv(L)**T is the ConjTranspose of the previous result, so we use SSYRK
         */
        cblas_ssyrk(CblasColMajor, (CBLAS_UPLO)uplo, CblasTrans,
                    N, N, done, TEMP, N, mdone, Id, N);
    } else {
        /*
         * if ITYPE = 1 or 2, Z**T*B*Z = I;
         */
        str = "Z**T*B*Z";
        cblas_ssymm(CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo, N, N,
                    (zone),  B,    LDB,
                                        Z,    LDZ,
                    (zzero), TEMP, N);
        cblas_sgemm(CblasColMajor, CblasTrans, CblasNoTrans, N, N, N,
                    (zone),  Z,    LDZ,
                                        TEMP, N,
                    (mzone), Id  , N);
    }

    normQ = LAPACKE_slansy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), 'U', N, Id, N, work);
    result = normQ / (N * eps);
    printf(" ======================================================\n");
    printf(" | Id-%s |_oo / (N*eps)          : %15.3E \n",str, result );
    printf(" ======================================================\n");

    if ( isnan(result) || isinf(result) || (result > 60.0) ) {
        printf("-- Orthogonality is suspicious ! \n");
        info_ortho=1;
    }
    else {
        printf("-- Orthogonality is CORRECT ! \n");
        info_ortho=0;
    }

    free(work); free(Id); free(TEMP);
    return info_ortho;
}

/*------------------------------------------------------------
 *  Check the reduction
 */
static int check_reduction(int itype, int uplo, int N, float *D,
                           float *A, int LDA,
                           float *B, int LDB,
                           float *Z, int LDZ,
                           float eps )
{
    float zone  =  1.0;
    float zzero =  0.0;
    float mzone = -1.0;
    float Anorm, Znorm, Rnorm, result;
    int info_reduction;
    int i;
    char *str;
    float *TEMP     = (float *)malloc(N*N*sizeof(float));
    float *Residual = (float *)malloc(N*N*sizeof(float));
    float *work = (float *)malloc(N*sizeof(float));

    if( itype == 1 ) {
        /*
         * | A Z - B Z D | / ( |A| |Z| n ulp )
         */
        str = " A Z - B Z D ";

        /* Compute TEMP = Z*D */
        LAPACKE_slacpy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaUpperLower),
                            N, N, Z, LDZ, TEMP, N);

        for (i=0; i<N; i++) {
            cblas_sscal(N, D[i], TEMP + i*N, 1);
        }

        /* Compute Residual = B*Z*D = B*TEMP */
        cblas_ssymm(CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo, N, N,
                    (zone),  B,        LDB,
                                        TEMP,     N,
                    (zzero), Residual, N);

        /* Compute A*Z - B*Z*D */
        cblas_ssymm(CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo, N, N,
                    (zone),  A,        LDA,
                                        Z,        LDZ,
                    (mzone), Residual, N);
    }
    else if( itype == 2 ) {
        /*
         * | A B Z - Z D | / ( |A| |Z| n ulp )
         */
        str = " A B Z - Z D ";

        /* Compute Residual = Z*D */
        LAPACKE_slacpy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaUpperLower),
                            N, N, Z, LDZ, Residual, N);

        for (i=0; i<N; i++) {
            cblas_sscal(N, D[i], Residual + i*N, 1);
        }

        /* Compute TEMP = B*Z */
        cblas_ssymm(CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo, N, N,
                    (zone),  B,    LDB,
                                        Z,    LDZ,
                    (zzero), TEMP, N);

        /* Compute A*B*Z - Z*D = A*TEMP-Residual */
        cblas_ssymm(CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo, N, N,
                    (zone),  A,        LDA,
                                        TEMP,     N,
                    (mzone), Residual, N);
    } else {
        /*
         * | B A Z - Z D | / ( |A| |Z| n ulp )
         */
        str = " B A Z - Z D ";

        /* Compute Residual = Z*D */
        LAPACKE_slacpy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaUpperLower),
                            N, N, Z, LDZ, Residual, N);

        for (i=0; i<N; i++) {
            cblas_sscal(N, D[i], Residual + i*N, 1);
        }

        /* Compute TEMP = A*Z */
        cblas_ssymm(CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo, N, N,
                    (zone),  A,    LDA,
                                        Z,    LDZ,
                    (zzero), TEMP, N);

        /* Compute B*A*Z - Z*D = B*TEMP-Residual */
        cblas_ssymm(CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo, N, N,
                    (zone),  B,        LDB,
                                        TEMP,     N,
                    (mzone), Residual, N);
    }

    Rnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaOneNorm), N, N, Residual, N,   work);
    Anorm = LAPACKE_slansy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaOneNorm), lapack_const(uplo), N, A,        LDA, work);
    Znorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaOneNorm), N, N, Z,        LDZ, work);

    result = Rnorm / ( Anorm * Znorm * N * eps);
    printf(" ======================================================\n");
    printf(" | %s |_oo/(|A|_oo.|Z|_oo.N.eps) : %15.3E \n", str, result );
    printf(" ======================================================\n");

    if ( isnan(result) || isinf(result) || (result > 60.0) ) {
        printf("-- Reduction is suspicious ! \n");
        info_reduction = 1;
    }
    else {
        printf("-- Reduction is CORRECT ! \n");
        info_reduction = 0;
    }

    free(TEMP); free(Residual); free(work);
    return info_reduction;
}
/*--------------------------------------------------------------
 * Check the solution
 */
static int check_solution(int N, float *E1, float *E2, float eps)
{
    int info_solution, i;
    float resid;
    float maxtmp;
    float maxel = fabs( fabs(E1[0]) - fabs(E2[0]) );
    float maxeig = max( fabs(E1[0]), fabs(E2[0]) );

    for (i=1; i<N; i++) {
        resid  = fabs(fabs(E1[i])-fabs(E2[i]));
        maxtmp = max(fabs(E1[i]), fabs(E2[i]));

        /* Update */
        maxeig = max(maxtmp, maxeig);
        maxel  = max(resid,  maxel );
    }

    maxel = maxel / (maxeig * N * eps);
    printf(" ======================================================\n");
    printf(" | D - eigcomputed | / (|D| * N * eps) : %15.3E \n",  maxel );
    printf(" ======================================================\n");

    printf("============\n");
    printf("Checking the eigenvalues of A\n");
    if ( isnan(maxel) || isinf(maxel) || (maxel > 100) ) {
        printf("-- The eigenvalues are suspicious ! \n");
        info_solution = 1;
    }
    else{
        printf("-- The eigenvalues are CORRECT ! \n");
        info_solution = 0;
    }
    return info_solution;
}


text:
#ifndef Two_Group_Diffusion_h
#define Two_Group_Diffusion_h

#include <vector>

// Trilinos includes
#include "Kokkos_DefaultNode.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_RCP.hpp"
#include "Tpetra_CrsGraph.hpp"
#include "Tpetra_CrsMatrix.hpp"
#include "Tpetra_MultiVector.hpp"
#include "Tpetra_Operator.hpp"
#include "Tpetra_RowMatrix.hpp"
#include "Tpetra_Vector.hpp"

// SCALE includes
#include "Nemesis/comm/global.hh"

// vendored includes
#include <gsl/gsl>

// enrico includes
#include "Assembly_Model.h"
#include "Neutronics_Solver.h"
#include "Two_Group_Cross_Sections.h"

namespace enrico {

//===========================================================================//
/*!
 * \class Two_Group_Diffusion
 * \brief Two-group, 3D neutron diffusion solver.
 */
//===========================================================================//

class Two_Group_Diffusion : public Neutronics_Solver {
public:
  enum BC_TYPE { VACUUM, REFLECT };
  enum FACE { LO_X, HI_X, LO_Y, HI_Y, LO_Z, HI_Z };

  //@{
  //! Typedefs
  using XS = Two_Group_Cross_Sections;
  using XS_Data = XS::XS_Data;
  using Vec_Dbl = std::vector<double>;
  using Vec_Int = std::vector<int>;
  using Vec_BC = std::vector<BC_TYPE>;
  using SP_Assembly = std::shared_ptr<Assembly_Model>;
  using Pin_Map = std::vector<Assembly_Model::PIN_TYPE>;
  using ST = double;
  using LO = int;
  using GO = int;
  using NODE = KokkosClassic::DefaultNode::DefaultNodeType;
  using MV = Tpetra::MultiVector<ST, LO, GO, NODE>;
  using VECTOR = Tpetra::Vector<ST, LO, GO, NODE>;
  using OP = Tpetra::Operator<ST, LO, GO, NODE>;
  using MAP = Tpetra::Map<LO, GO, NODE>;
  using MATRIX = Tpetra::CrsMatrix<ST, LO, GO, NODE>;
  using GRAPH = Tpetra::CrsGraph<LO, GO, NODE>;
  using RCP_PL = Teuchos::RCP<Teuchos::ParameterList>;
  //@}

private:
  // >>> DATA
  SP_Assembly d_assembly;

  Vec_Dbl d_dx;
  Vec_Dbl d_dy;
  Vec_Dbl d_dz;
  Pin_Map d_pin_map;
  XS d_xs;
  size_t d_Nx, d_Ny, d_Nz, d_num_cells;

  // Solver parameters
  double d_tol;
  int d_max_iters;

  // Boundary conditions
  Vec_BC d_bcs;

  // Data stored over single group
  Teuchos::RCP<MAP> d_map;
  Teuchos::RCP<GRAPH> d_graph;
  Teuchos::RCP<MATRIX> d_A_therm, d_A_fast;
  Teuchos::RCP<VECTOR> d_x_therm, d_x_fast;
  Teuchos::RCP<VECTOR> d_b_therm, d_b_fast;

public:
  // Constructor
  Two_Group_Diffusion(SP_Assembly assembly, RCP_PL params, const Vec_Dbl& dz);

  // Solve
  virtual void solve(const Vec_Dbl& temperatures,
                     const Vec_Dbl& densities,
                     Vec_Dbl& power) override;

private:
  int cellid(int ix, int iy, int iz)
  {
    Expects(ix >= 0);
    Expects(ix < d_Nx);
    Expects(iy >= 0);
    Expects(iy < d_Ny);
    Expects(iz >= 0);
    Expects(iz < d_Nz);
    int cell = ix + d_Nx * (iy + d_Ny * iz);
    Ensures(cell >= 0);
    Ensures(cell < d_num_cells);
    return cell;
  }

  void build_matrices(const std::vector<XS_Data>& xs_data);

  // Compute 2-norm of vector
  double vec_norm(const Teuchos::ArrayRCP<double>& vec)
  {
    double nrm = 0.0;
    for (auto val : vec)
      nrm += val * val;
    return std::sqrt(nrm);
  }

  // Apply scale factor to vector data
  void scale_vec(Teuchos::ArrayRCP<double>& vec, double factor)
  {
    for (auto& val : vec)
      val *= factor;
  }
};

//---------------------------------------------------------------------------//
} // end namespace enrico

//---------------------------------------------------------------------------//
#endif // Two_Group_Diffusion_h

//---------------------------------------------------------------------------//
// end of Two_Group_Diffusion.h
//---------------------------------------------------------------------------//

max_stars_repo_name: pshriwise/enrico

max_issues_repo_name: pshriwise/enrico

max_forks_repo_name: pshriwise/enrico


text:
/*
 * Licensed to the Apache Software Foundation (ASF) under one
 * or more contributor license agreements.  See the NOTICE file
 * distributed with this work for additional information
 * regarding copyright ownership.  The ASF licenses this file
 * to you under the Apache License, Version 2.0 (the
 * "License"); you may not use this file except in compliance
 * with the License.  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing,
 * software distributed under the License is distributed on an
 * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
 * KIND, either express or implied.  See the License for the
 * specific language governing permissions and limitations
 * under the License.
 */

#ifndef LIBHE_H
#define LIBHE_H

#include <cassert>
#include <algorithm>
#include <optional>
#include <gsl/span>

#include "seal/seal.h"
#include "seal/util/common.h"
#include "seal/util/rlwe.h"
#include "seal/util/polyarithsmallmod.h"

using namespace std;
using namespace seal;

class RawPolynomData {
    vector<uint64_t > _data;
    size_t _size;

public:
    explicit RawPolynomData(const SEALContext& context);

    SEAL_NODISCARD inline const size_t& size() const { return _size; };
    SEAL_NODISCARD inline uint64_t* data() { return _data.data(); };
    SEAL_NODISCARD inline gsl::span<uint64_t > data_span() { return { data(), size() }; };

    void set_data(vector<uint64_t >& data);
};

gsl::span<Ciphertext::ct_coeff_type > data_span(Ciphertext& c, size_t n);

RawPolynomData generate_a(const SEALContext& context);

EncryptionParameters generateParameters();

size_t get_slot_count(const SEALContext& ctx);

// returns a vector filled with random double values between 0 and 1
vector<double> random_plaintext_data(size_t count);

struct GlobalState {
    SEALContext context;
    RawPolynomData a;
    double scale;

    explicit GlobalState(double _scale);
};

class Client {
    GlobalState _gs;
    CKKSEncoder _encoder;
    SecretKey _partial_secret_key;
    PublicKey _partial_public_key;
    std::optional<PublicKey> _public_key = std::nullopt;
    std::unique_ptr<Encryptor> _encryptor = nullptr;

    SEAL_NODISCARD static PublicKey generate_partial_public_key(const SecretKey &secret_key, const SEALContext &context, RawPolynomData& a);

public:
    explicit Client(GlobalState global_state);

    SEAL_NODISCARD inline const SEALContext& context() const { return _gs.context; };
    SEAL_NODISCARD inline const PublicKey& partial_public_key() const { return _partial_public_key; };
    SEAL_NODISCARD inline const CKKSEncoder& encoder() const { return _encoder; };
    SEAL_NODISCARD inline CKKSEncoder& encoder() { return _encoder; };
    SEAL_NODISCARD inline const Encryptor& encryptor() const { assert(_encryptor != nullptr); return *_encryptor; };
    SEAL_NODISCARD inline const PublicKey& public_key() { return *_public_key; };
    inline void set_public_key(const PublicKey& pk) { _public_key = make_optional(pk); };

    Ciphertext encrypted_data(gsl::span<const double> plain_data);

    Plaintext partial_decryption(const Ciphertext& encrypted);
};

// adds b to a in place
template<typename T> void sum_first_poly_inplace(const SEALContext& context, T& a, const T& b) {
    auto &context_data = *context.get_context_data(a.parms_id());
    auto &parms = context_data.parms();
    auto &coeff_modulus = parms.coeff_modulus();
    size_t coeff_count = parms.poly_modulus_degree();
    size_t coeff_modulus_size = coeff_modulus.size();

    // by dereferencing we get only the first poly
    auto summand_iter = *util::ConstPolyIter(b.data(), coeff_count, coeff_modulus_size);
    auto sum_iter = *util::ConstPolyIter(a.data(), coeff_count, coeff_modulus_size);
    auto result_iter = *util::PolyIter(a.data(), coeff_count, coeff_modulus_size);
    // see Evaluator::add_inplace
    util::add_poly_coeffmod(sum_iter, summand_iter, coeff_modulus_size, coeff_modulus, result_iter);
}

// This function adds the first polys in summands to sum (either Ciphertext or Plaintext).
template<typename T> T sum_first_polys_inplace(const SEALContext& context, T& sum, gsl::span<const T> summands) {
    for (size_t i = 0; i < summands.size(); i++) {
        sum_first_poly_inplace(context, sum, summands[i]);
    }
    return sum;
}

// This function sums the first polys in summands (either Ciphertext or Plaintext).
template<typename T> T sum_first_polys(const SEALContext& context, gsl::span<const T> summands) {
    T sum = summands[0];
    sum_first_polys_inplace(context, sum, gsl::span(&summands.data()[1], summands.size() - 1));
    return sum;
}

class Server {
    GlobalState _gs;
    PublicKey _public_key;

public:
    explicit Server(GlobalState global_state);

    SEAL_NODISCARD inline RawPolynomData& a() { return _gs.a; };
    SEAL_NODISCARD inline const SEALContext& context() const { return _gs.context; };
    SEAL_NODISCARD inline const PublicKey& public_key() const { return _public_key; };

    void accumulate_partial_public_keys(gsl::span<const Ciphertext> partial_pub_keys);

    Ciphertext sum_data(vector<Ciphertext>&& data) const;

    vector<double> average(const Ciphertext& encrypted_sum, gsl::span<const Plaintext> partial_decryptions) const;
};

#endif //LIBHE_H

max_stars_repo_name: apache/systemml

max_issues_repo_name: Baunsgaard/systemml

max_forks_repo_name: Baunsgaard/systemml