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

#include <gsl/gsl_blas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_vector.h>

#include <argtable2.h>

#include "qdm.h"

void
qdm_ispline_vector(
    gsl_vector *result,
    const double tau,
    const size_t spline_df,
    const gsl_vector *knots
)
{
  size_t bin = qdm_vector_search(knots, tau);

  gsl_vector_set(result, 0, 1);

  for (size_t m = 0; m < result->size - 1; m++) {
    double v;

    if (bin < m) {
      v = 0;
    } else if (bin - spline_df + 1 > m) {
      v = 1;
    } else if (bin == m) {
      double n = gsl_sf_pow_int(tau - gsl_vector_get(knots, m), 3);
      double d1 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m);
      double d2 = gsl_vector_get(knots, m + 2) - gsl_vector_get(knots, m);
      double d3 = gsl_vector_get(knots, m + 1) - gsl_vector_get(knots, m);

      v = n / (d1 * d2 * d3);
    } else if (bin == m + 1) {
      double i1n = 3 * tau;
      double i1d = gsl_vector_get(knots, m + 2) - gsl_vector_get(knots, m + 1);
      double i1 = i1n / i1d;

      double i2n = -gsl_sf_pow_int(tau - gsl_vector_get(knots, m), 3);
      double i2d1 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m);
      double i2d2 = gsl_vector_get(knots, m + 2) - gsl_vector_get(knots, m);
      double i2d3 = gsl_vector_get(knots, m + 2) - gsl_vector_get(knots, m + 1);
      double i2 = i2n / (i2d1 * i2d2 * i2d3);

      double i3n = gsl_sf_pow_int(gsl_vector_get(knots, m + 3) - tau, 3);
      double i3d1 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m + 1);
      double i3d2 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m);
      double i3d3 = gsl_vector_get(knots, m + 2) - gsl_vector_get(knots, m + 1);
      double i3 = i3n / (i3d1 * i3d2 * i3d3);

      double i4an = -3 * gsl_vector_get(knots, m + 1);
      double i4ad = gsl_vector_get(knots, m + 2) - gsl_vector_get(knots, m + 1);
      double i4a = i4an / i4ad;
      double i4bn = gsl_sf_pow_int(gsl_vector_get(knots, m + 1) - gsl_vector_get(knots, m), 3);
      double i4bd1 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m);
      double i4bd2 = gsl_vector_get(knots, m + 2) - gsl_vector_get(knots, m);
      double i4bd3 = gsl_vector_get(knots, m + 2) - gsl_vector_get(knots, m + 1);
      double i4b = i4bn / (i4bd1 * i4bd2 * i4bd3);
      double i4cn = gsl_sf_pow_int(gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m + 1), 3);
      double i4cd1 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m + 1);
      double i4cd2 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m);
      double i4cd3 = gsl_vector_get(knots, m + 2) - gsl_vector_get(knots, m + 1);
      double i4c = i4cn / (i4cd1 * i4cd2 * i4cd3);
      double i4 = i4a + i4b - i4c;

      double i5 = 0;
      if (m != 1) {
        double i5n = gsl_sf_pow_int(gsl_vector_get(knots, m + 1) - gsl_vector_get(knots, m), 3);
        double i5d1 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m);
        double i5d2 = gsl_vector_get(knots, m + 2) - gsl_vector_get(knots, m);
        double i5d3 = gsl_vector_get(knots, m + 1) - gsl_vector_get(knots, m);
        i5 = i5n / (i5d1 * i5d2 * i5d3);
      }

      v = i1 + i2 + i3 + i4 + i5;
    } else {
      double n = gsl_sf_pow_int(gsl_vector_get(knots, m + 3) - tau, 3);
      double d1 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m);
      double d2 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m + 1);
      double d3 = gsl_vector_get(knots, m + 3) - gsl_vector_get(knots, m + 2);
      double d = d1 * d2 * d3;

      v = 1 - n / d;
    }

    gsl_vector_set(result, m + 1, v);
  }
}

void
qdm_ispline_matrix(
    gsl_matrix *result,
    const gsl_vector *taus,
    const size_t spline_df,
    const gsl_vector *knots
)
{
  for (size_t i = 0; i < result->size1; i++) {
    gsl_vector_view row = gsl_matrix_row(result, i);
    qdm_ispline_vector(&row.vector, gsl_vector_get(taus, i), spline_df, knots);
  }
}

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
#ifndef __RNG_H__
#define __RNG_H__

#include <gsl/gsl_rng.h>
#include <ctime>
#include<iostream>
#include <cstring>
#include <fstream>

#define ABS(x) ((x) >= 0 ? (x) : -(x))
#define min(a,b) ((a) <= (b) ? (a) : (b))
#define max(a,b) ((a) >= (b) ? (a) : (b))
#define OVERFLO 1e100
#define UNDERFLO 1e-100

namespace weakarg
{

extern gsl_rng * rng;

unsigned long makerng(bool fast=false);
unsigned long seedrng(unsigned long seed=0);
int saverng(std::string fname);
int loadrng(std::string fname);

double RandomReal(double low, double high);
int RandomInteger(int low, int high);
double rnd();
double RGamma(double n,double lambda);
void RDirichlet(const double * a, const int k, double * b);
long RPoisson(double mu);
double RExpon(double av);
double RNormal(double mu,double sd) ;
double fsign( double num, double sign );
double sexpo(void);
 double snorm();
 double genexp(double av);
 long ignpoi(double mean);
 long ignuin(int low, int high);
 double genunf(double low, double high);
 long   Binomial(int n, double p);
 long   Binomial1(int n, double p);
 double BinoProb(int n, double p,int i);
 void LogRDirichlet (const double *a, const int k, double *b,double *c);

} // end namespace weakarg
#endif // __RNG_H__

ext: h


text:
/*
 * SPDX-License-Identifier:  BSD-3-Clause
 *
 * example_04-ArrayOfStructs-Naive-Omp-SIMD-Tiled.c :
 *      Example of SPH Density Calculation using a
 *      naive implementation of the main density loop,
 *      no neighbours earch, and Array of Structs (AoS)
 *      data layout, OpenMP parallelization and SIMD
 *      directives on the kernel and density calculation.
 *      This incorporates strip mining and exchange to
 *      implement cache blocking and support performance
 *      for large number of particles that would otherwise
 *      be lost.
 *
 * (C) Copyright 2021 José Hugo Elsas
 * Author: José Hugo Elsas <jhelsas@gmail.com>
 *
 * Command Line Options:
 *   -runs  <int>   : Set the number of repetitions (runs) for
 *                      calculating the density. The value of
 *                      the density is based on the last
 *                      iteration.
 *                    Default value: 1
 *   -run_seed <int>: Flag to set an alternative seed use for
 *                      for the PRNG. Instead of feeding seed
 *                      to the PRNG directly, it feeds
 *                      seed + iteration, as to generate different
 *                      configurations for each iteration.
 *                    Default value: 0 - (possible 0/1)
 *   -seed     <int>: Set the seed to use for the SPH particles
 *                      uniform position generation in the box
 *                    Default value: 123123123
 *
 *   -N        <int>: Set the number of SPH particles to be used
 *                    Default value: 1e5 = 100,000
 *   -h      <float>: Set the value of the smoothing kernel
 *                      parameter h, which corresponds to half
 *                      of the support of the kernel.
 *                    Default value: 0.05
 *
 *   -Nx       <int>: Set the number of Cells in the X direction
 *                    Default value: 10
 *   -Ny       <int>: Set the number of Cells in the Y direction
 *                    Default value: 10
 *   -Nz       <int>: Set the number of Cells in the Z direction
 *                    Default value: 10
 *
 *   -Xmin   <float>: Set the lower bound in the X direction for
 *                      the Cell Linked List box
 *                    Default value: 0.0
 *   -Ymin   <float>: Set the lower bound in the Y direction for
 *                    the Cell Linked List box
 *                      Default value: 0.0
 *   -Ymin   <float>: Set the lower bound in the Z direction for
 *                      the Cell Linked List box
 *                    Default value: 0.0
 *
 *   -Xmax   <float>: Set the lower bound in the X direction for
 *                      the Cell Linked List box
 *                    Default value: 1.0
 *   -Ymax   <float>: Set the lower bound in the Y direction for
 *                      the Cell Linked List box
 *                    Default value: 1.0
 *   -Zmax   <float>: Set the lower bound in the Z direction for
 *                      the Cell Linked List box
 *                    Default value: 1.0
 */

#include <math.h>
#include <ctype.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <limits.h>
#include <unistd.h>
#include <stdbool.h>
#include <sys/time.h>
#include <inttypes.h>

#include <omp.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_heapsort.h>

#include "sph_data_types.h"
#include "sph_linked_list.h"
#include "sph_utils.h"

#ifndef M_PI
#define M_PI (3.14159265358979323846)
#endif

#define COMPUTE_BLOCKS 1

int main_loop(int run, bool run_seed, int64_t N, double h, long int seed,
              linkedListBox *box, SPHparticle *lsph, double *times);

int compute_density_3d_naive_omp_simd_tiled(int N,double h,SPHparticle *lsph);

double w_bspline_3d_constant(double h);

#pragma omp declare simd
double w_bspline_3d_simd(double q);

int main(int argc, char **argv){
  bool run_seed = false;       // By default the behavior is is to use the same seed
  int err, runs = 1;           // By default the main loop only runs once
  long int seed = 123123123;   // The default seed is 123123123
  int64_t N = 100000;          // The default number of particles is N = 1e5 = 100,000
  double h=0.05;               // The default kernel smoothing length is h = 0.05
  linkedListBox *box;          // Uninitialized Box containing the cells for the cell linked list method
  SPHparticle *lsph;           // Uninitialized array of SPH particles

  box = (linkedListBox*)malloc(1*sizeof(linkedListBox)); // Create a box representing the entire 3d domain

  // allow for command line customization of the run
  arg_parse(argc,argv,&N,&h,&seed,&runs,&run_seed,box);  // Parse the command line options
                                                         // line arguments and override default values
  lsph = (SPHparticle*)malloc(N*sizeof(SPHparticle));    // Create an array of N particles

  double times[runs*COMPUTE_BLOCKS];

  for(int run=0;run<runs;run+=1)
    main_loop(run,run_seed,N,h,seed,box,lsph,times);

  bool is_cll = false;
  const char *prefix = "ex04,naive,AoS,omp,simd,tiled";
  print_time_stats(prefix,is_cll,N,h,seed,runs,lsph,box,times);
  print_sph_particles_density(prefix,is_cll,N,h,seed,runs,lsph,box);

  free(lsph);
  safe_free_box(box);

  return 0;
}

/*
 *  Function main_loop:
 *    Runs the main loop of the program, including the particle array generation,
 *    density calculation and the timings annotations.
 *
 *    Arguments:
 *       run <int>            : index (or value) or the present iteration
 *       run_seed <bool>      : boolean defining whether to use run index for seed or not
 *       N <int>              : Number of SPH particles to be used in the run
 *       h <double>           : Smoothing Length for the Smoothing Kernel w_bspline
 *       seed <long int>      : seed for GSL PRNG generator to generate particle positions
 *       box  <linkedListBox> : Box of linked list cells, encapsulating the 3d domain
 *       lsph <SPHparticle>   : Array (pointer) of SPH particles to be updated
 *       times <double>       : Array to store the computation timings to be updated
 *    Returns:
 *       0                    : error code returned
 *       lsph <SPHparticle>   : SPH particle array is updated in the rho field by reference
 *       times <double>       : Times is updated by reference
 */
int main_loop(int run, bool run_seed, int64_t N, double h, long int seed,
              linkedListBox *box, SPHparticle *lsph, double *times)
{
  int err;

  if(run_seed)
    err = gen_unif_rdn_pos_box(N,seed+run,box,lsph);
  else
    err = gen_unif_rdn_pos_box(N,seed,box,lsph);

  if(err)
    fprintf(stderr,"error in gen_unif_rdn_pos\n");

  // ------------------------------------------------------ //

  double t0,t1;

  t0 = omp_get_wtime();

  compute_density_3d_naive_omp_simd_tiled(N,h,lsph);        // Compute the density for all particles

  t1 = omp_get_wtime();

  // ------------------------------------------------------ //

  times[COMPUTE_BLOCKS*run+0] = t1-t0;                      // Only one component to measure time

  return 0;
}

/*
 *  Function compute_density_3d_naive_omp_simd_tiled:
 *    Computes the SPH density from the particles implementing a strip mine and exchange
 *    strategy to re-use data in cache over the direct loop. It executes calculations
 *    in parallel for the outer-most loop using openMP and SIMD in inner-most loop,
 *    though SIMD only for limited success.
 *
 *    Reference: https://en.wikipedia.org/wiki/Loop_nest_optimization
 *
 *    Arguments:
 *       N <int>              : Number of SPH particles to be used in the run
 *       h <double>           : Smoothing Length for the Smoothing Kernel w_bspline
 *       lsph <SPHparticle>   : Array (pointer) of SPH particles to be updated
 *    Returns:
 *       0                    : error code returned
 *       lsph <SPHparticle>   : SPH particle array is updated in the rho field by reference
 */
int compute_density_3d_naive_omp_simd_tiled(int N,double h,SPHparticle *lsph){
  const double inv_h = 1./h;                                       // Pre-invert the smoothing distance
  const double kernel_constant = w_bspline_3d_constant(h);         // Pre-compute the 3d normalization constant
  const int64_t STRIP = 500;                                       // Setting the size of the strip or block

  #pragma omp parallel for                                         // Run the iteration in Parallel
  for(int64_t ii=0;ii<N;ii+=1)                                     // Iterate
    lsph[ii].rho = 0.;                                             // Pre-initialize the density to zero

  #pragma omp parallel for                                         // Run the iteration in i in parallel
  for(int64_t i=0;i<N;i+=STRIP){                                   // Breaking up the i and j iterations in blocks
    for(int64_t j=0;j<N;j+=STRIP){                                 // of size STRIP to do data re-use and cache blocking
      for(int64_t ii=i;ii < ((i+STRIP<N)?(i+STRIP):N); ii+=1){     // Iterate a block over ii
        double xii = lsph[ii].r.x;                                 // Load the position in X for ii
        double yii = lsph[ii].r.y;                                 // Load the position in Y for ii
        double zii = lsph[ii].r.z;                                 // Load the position in Z for ii
        double rhoii = 0.0;                                        // Initialize partial density ii density to zero

        #pragma omp simd                                           // Hint at the compiler to vectorize this loop
        for(int64_t jj=j;jj < ((j+STRIP<N)?(j+STRIP):N); jj+=1 ){  // and iterate over the jj part of the block
          double q = 0.;                                           // initialize the distance variable

          double xij = xii-lsph[jj].r.x;                           // Load and subtract jj particle's X position component
          double yij = yii-lsph[jj].r.y;                           // Load and subtract jj particle's Y position component
          double zij = zii-lsph[jj].r.z;                           // Load and subtract jj particle's Z position component

          q += xij*xij;                                            // Add the jj contribution to the ii distance in X
          q += yij*yij;                                            // Add the jj contribution to the ii distance in Y
          q += zij*zij;                                            // Add the jj contribution to the ii distance in Z

          q = sqrt(q)*inv_h;                                       // Sqrt and normalizing the distance by the smoothing lengh

          rhoii += lsph[jj].nu*w_bspline_3d_simd(q);               // Add up the contribution from the jj particle
        }                                                          // to the intermediary density and then
        lsph[ii].rho += kernel_constant*rhoii;                     // add the intermediary density to the full density
      }
    }
  }

  return 0;
}

/*
 *  Function w_bspline_3d_constant:
 *    Returns the 3d normalization constant for the cubic b-spline SPH smoothing kernel
 *
 *    Arguments:
 *       h <double>           : Smoothing Length for the Smoothing Kernel w_bspline
 *    Returns:
 *       3d bspline normalization density <double>
 */
double w_bspline_3d_constant(double h){
  return 3./(2.*M_PI*h*h*h);  // 3d normalization value for the b-spline kernel
}

/*
 *  Function w_bspline_3d_simd:
 *    Returns the un-normalized value of the cubic b-spline SPH smoothing kernel
 *
 *    Arguments:
 *       q <double>           : Distance between particles normalized by the smoothing length h
 *    Returns:
 *       wq <double>          : Unnormalized value of the kernel
 *
 *    Observation:
 *       Why not else if(q<2.)?
 *       Because if you use "else if", the compiler refuses to vectorize,
 *       This results in a large slowdown, as of 2.5x slower for example_04
 */
#pragma omp declare simd
double w_bspline_3d_simd(double q){                                // Use as input the normalized distance
  double wq = 0.0;
  double wq1 = (0.6666666666666666 - q*q + 0.5*q*q*q);             // The first polynomial of the spline
  double wq2 = 0.16666666666666666*(2.-q)*(2.-q)*(2.-q);           // The second polynomial of the spline

  if(q<2.)                                                         // If the distance is below 2
    wq = wq2;                                                      // Use the 2nd polynomial for the spline

  if(q<1.)                                                         // If the distance is below 1
    wq = wq1;                                                      // Use the 1nd polynomial for the spline

  return wq;                                                       // return which ever value corresponds to the distance
}

ext: c

max_stars_repo_name: jhelsas/sphalerite

max_issues_repo_name: jhelsas/sphalerite

max_forks_repo_name: jhelsas/sphalerite


text:
/*
 *  Player - One Hell of a Robot Server
 *  Copyright (C) 2000  Brian Gerkey   &  Kasper Stoy
 *                      gerkey@usc.edu    kaspers@robotics.usc.edu
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This library 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
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */
/**************************************************************************
 * Desc: Vector functions
 * Author: Andrew Howard
 * Date: 10 Dec 2002
 * CVS: $Id: pf_vector.c 6345 2008-04-17 01:36:39Z gerkey $
 *************************************************************************/

#include <math.h>
//#include <gsl/gsl_matrix.h>
//#include <gsl/gsl_eigen.h>
//#include <gsl/gsl_linalg.h>

#include "amcl/pf/pf_vector.h"
#include "amcl/pf/eig3.h"

// Return a zero vector
pf_vector_t pf_vector_zero()
{
  pf_vector_t c;

  c.v[0] = 0.0;
  c.v[1] = 0.0;
  c.v[2] = 0.0;

  return c;
}

// Check for NAN or INF in any component
int pf_vector_finite(pf_vector_t a)
{
  int i;

  for (i = 0; i < 3; i++)
    if (!isfinite(a.v[i]))
      return 0;

  return 1;
}

// Print a vector
void pf_vector_fprintf(pf_vector_t a, FILE *file, const char *fmt)
{
  int i;

  for (i = 0; i < 3; i++)
  {
    fprintf(file, fmt, a.v[i]);
    fprintf(file, " ");
  }
  fprintf(file, "\n");

  return;
}

// Simple vector addition
pf_vector_t pf_vector_add(pf_vector_t a, pf_vector_t b)
{
  pf_vector_t c;

  c.v[0] = a.v[0] + b.v[0];
  c.v[1] = a.v[1] + b.v[1];
  c.v[2] = a.v[2] + b.v[2];

  return c;
}

// Simple vector subtraction
pf_vector_t pf_vector_sub(pf_vector_t a, pf_vector_t b)
{
  pf_vector_t c;

  c.v[0] = a.v[0] - b.v[0];
  c.v[1] = a.v[1] - b.v[1];
  c.v[2] = a.v[2] - b.v[2];

  return c;
}

// Transform from local to global coords (a + b)
pf_vector_t pf_vector_coord_add(pf_vector_t a, pf_vector_t b)
{
  pf_vector_t c;

  c.v[0] = b.v[0] + a.v[0] * cos(b.v[2]) - a.v[1] * sin(b.v[2]);
  c.v[1] = b.v[1] + a.v[0] * sin(b.v[2]) + a.v[1] * cos(b.v[2]);
  c.v[2] = b.v[2] + a.v[2];
  c.v[2] = atan2(sin(c.v[2]), cos(c.v[2]));

  return c;
}

// Transform from global to local coords (a - b)
pf_vector_t pf_vector_coord_sub(pf_vector_t a, pf_vector_t b)
{
  pf_vector_t c;

  c.v[0] = +(a.v[0] - b.v[0]) * cos(b.v[2]) + (a.v[1] - b.v[1]) * sin(b.v[2]);
  c.v[1] = -(a.v[0] - b.v[0]) * sin(b.v[2]) + (a.v[1] - b.v[1]) * cos(b.v[2]);
  c.v[2] = a.v[2] - b.v[2];
  c.v[2] = atan2(sin(c.v[2]), cos(c.v[2]));

  return c;
}

// Return a zero matrix
pf_matrix_t pf_matrix_zero()
{
  int i, j;
  pf_matrix_t c;

  for (i = 0; i < 3; i++)
    for (j = 0; j < 3; j++)
      c.m[i][j] = 0.0;

  return c;
}

// Check for NAN or INF in any component
int pf_matrix_finite(pf_matrix_t a)
{
  int i, j;

  for (i = 0; i < 3; i++)
    for (j = 0; j < 3; j++)
      if (!isfinite(a.m[i][j]))
        return 0;

  return 1;
}

// Print a matrix
void pf_matrix_fprintf(pf_matrix_t a, FILE *file, const char *fmt)
{
  int i, j;

  for (i = 0; i < 3; i++)
  {
    for (j = 0; j < 3; j++)
    {
      fprintf(file, fmt, a.m[i][j]);
      fprintf(file, " ");
    }
    fprintf(file, "\n");
  }
  return;
}

/*
// Compute the matrix inverse
pf_matrix_t pf_matrix_inverse(pf_matrix_t a, double *det)
{
  double lndet;
  int signum;
  gsl_permutation *p;
  gsl_matrix_view A, Ai;

  pf_matrix_t ai;

  A = gsl_matrix_view_array((double*) a.m, 3, 3);
  Ai = gsl_matrix_view_array((double*) ai.m, 3, 3);

  // Do LU decomposition
  p = gsl_permutation_alloc(3);
  gsl_linalg_LU_decomp(&A.matrix, p, &signum);

  // Check for underflow
  lndet = gsl_linalg_LU_lndet(&A.matrix);
  if (lndet < -1000)
  {
    //printf("underflow in matrix inverse lndet = %f", lndet);
    gsl_matrix_set_zero(&Ai.matrix);
  }
  else
  {
    // Compute inverse
    gsl_linalg_LU_invert(&A.matrix, p, &Ai.matrix);
  }

  gsl_permutation_free(p);

  if (det)
    *det = exp(lndet);

  return ai;
}
*/

// Decompose a covariance matrix [a] into a rotation matrix [r] and a diagonal
// matrix [d] such that a = r d r^T.
void pf_matrix_unitary(pf_matrix_t *r, pf_matrix_t *d, pf_matrix_t a)
{
  int i, j;
  /*
  gsl_matrix *aa;
  gsl_vector *eval;
  gsl_matrix *evec;
  gsl_eigen_symmv_workspace *w;

  aa = gsl_matrix_alloc(3, 3);
  eval = gsl_vector_alloc(3);
  evec = gsl_matrix_alloc(3, 3);
  */

  double aa[3][3];
  double eval[3];
  double evec[3][3];

  for (i = 0; i < 3; i++)
  {
    for (j = 0; j < 3; j++)
    {
      //gsl_matrix_set(aa, i, j, a.m[i][j]);
      aa[i][j] = a.m[i][j];
    }
  }

  // Compute eigenvectors/values
  /*
  w = gsl_eigen_symmv_alloc(3);
  gsl_eigen_symmv(aa, eval, evec, w);
  gsl_eigen_symmv_free(w);
  */

  eigen_decomposition(aa,evec,eval);

  *d = pf_matrix_zero();
  for (i = 0; i < 3; i++)
  {
    //d->m[i][i] = gsl_vector_get(eval, i);
    d->m[i][i] = eval[i];
    for (j = 0; j < 3; j++)
    {
      //r->m[i][j] = gsl_matrix_get(evec, i, j);
      r->m[i][j] = evec[i][j];
    }
  }

  //gsl_matrix_free(evec);
  //gsl_vector_free(eval);
  //gsl_matrix_free(aa);

  return;
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_sf_zeta.h>

/* Calculate the Hurwitz zeta function by invoking the GNU scientific
   library */
/* Only intended to be used with R; not the most elegant integration but it
   may do for now */

/* Two arguments, the power used in the sum (s) and the baseline value (q) */

double s, q, zeta; /* arguments to zeta, value */
char *program_name; /* name program is invoked under, for errors */

int main(int argc, char* argv[]) {
  void usage(void);	/* Warn users about proper usage */

  program_name = argv[0];
  if (argc != 3) {
    usage();
  }
  s = atof(&argv[1][0]);
  if (s <= 1.0) {
    usage();
  }
  q = atof(&argv[2][0]);
  if (q <= 0.0) {
    usage();
  }
  zeta = gsl_sf_hzeta(s,q);
  printf("%.18e\n",zeta);
  return(0);
}

void usage(void) {
  (void) fprintf(stderr, "Usage is %s [floating-point exponent > 1] [floating-point additive constant > 0]\n", program_name);
  exit(8);
}

max_stars_repo_name: jguerber/spatialwarnings

max_stars_repo_licenses: Unlicense MIT

max_issues_repo_name: jguerber/spatialwarnings

max_issues_repo_licenses: Unlicense MIT

max_forks_repo_name: jguerber/spatialwarnings

max_forks_repo_licenses: Unlicense MIT


text:
/*
 * This file is part of MXE. See LICENSE.md for licensing information.
 */

/*  taken from http://www.netlib.org/lapack/lapacke.html */

/*  Calling CGEQRF and CUNGQR to compute Q with workspace querying */

#include <stdio.h>
#include <stdlib.h>
#include <lapacke_utils.h>
#include <cblas.h>

int main (int argc, const char * argv[])
{
   (void)argc;
   (void)argv;

   lapack_complex_float *a,*tau,*r,*work,one,zero,query;
   lapack_int info,m,n,lda,lwork;
   int i,j;
   float err;
   m = 10;   n = 5;   lda = m;
   one = lapack_make_complex_float(1.0,0.0);
   zero= lapack_make_complex_float(0.0,0.0);
   a = calloc(m*n,sizeof(lapack_complex_float));
   r = calloc(n*n,sizeof(lapack_complex_float));
   tau = calloc(m,sizeof(lapack_complex_float));
   for(j=0;j<n;j++)
      for(i=0;i<m;i++)
         a[i+j*m] = lapack_make_complex_float(i+1,j+1);
   info = LAPACKE_cgeqrf_work(LAPACK_COL_MAJOR,m,n,a,lda,tau,&query,-1);
   lwork = (lapack_int)query;
   info = LAPACKE_cungqr_work(LAPACK_COL_MAJOR,m,n,n,a,lda,tau,&query,-1);
   lwork = MAX(lwork,(lapack_int)query);
   work = calloc(lwork,sizeof(lapack_complex_float));
   info = LAPACKE_cgeqrf_work(LAPACK_COL_MAJOR,m,n,a,lda,tau,work,lwork);
   info = LAPACKE_cungqr_work(LAPACK_COL_MAJOR,m,n,n,a,lda,tau,work,lwork);
   for(j=0;j<n;j++)
      for(i=0;i<n;i++)
         r[i+j*n]=(i==j)?-one:zero;
   cblas_cgemm(CblasColMajor,CblasConjTrans,CblasNoTrans,
               n,n,m,&one,a,lda,a,lda,&one,r,n);
   err=0.0;
   for(i=0;i<n;i++)
      for(j=0;j<n;j++)
         err=MAX(err,cabs(r[i+j*n]));
   printf("error=%e\n",err);
   free(work);
   free(tau);
   free(r);
   free(a);
   return(info);
}

max_stars_repo_name: Lord-Kamina/mxe

max_issues_repo_name: Lord-Kamina/mxe

max_forks_repo_name: Lord-Kamina/mxe


text:
/* avl.c
 *
 * Copyright (C) 1998-2002, 2004 Free Software Foundation, Inc.
 * Copyright (C) 2018 Patrick Alken
 *
 * 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.
 */

/* This code is originally from GNU libavl, with some modifications */

#include <config.h>
#include <stdlib.h>
#include <stddef.h>
#include <string.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_bst.h>
#include <gsl/gsl_errno.h>

typedef struct gsl_bst_avl_node avl_node;
typedef gsl_bst_avl_table avl_table;
typedef gsl_bst_avl_traverser avl_traverser;

#ifndef AVL_MAX_HEIGHT
#define AVL_MAX_HEIGHT GSL_BST_AVL_MAX_HEIGHT
#endif

/* Function types. */
typedef void avl_item_func (void *avl_item, void *avl_param);
typedef void *avl_copy_func (void *avl_item, void *avl_param);

/* tree functions */
static int avl_init(const gsl_bst_allocator * allocator,
                    gsl_bst_cmp_function * compare, void * params, void * vtable);
static size_t avl_nodes (const void * vtable);
static int avl_empty (void * vtable);
static void ** avl_probe (void * item, avl_table * table);
static void * avl_insert (void * item, void * vtable);
static void * avl_find (const void *item, const void * vtable);
static void * avl_remove (const void *item, void * vtable);

/* traverser functions */
static int avl_t_init (void * vtrav, const void * vtable);
static void * avl_t_first (void * vtrav, const void * vtable);
static void * avl_t_last (void * vtrav, const void * vtable);
static void * avl_t_find (const void * item, void * vtrav, const void * vtable);
static void * avl_t_insert (void * item, void * vtrav, void * vtable);
static void * avl_t_copy (void * vtrav, const void * vsrc);
static void * avl_t_next (void * vtrav);
static void * avl_t_prev (void * vtrav);
static void * avl_t_cur (const void * vtrav);
static void * avl_t_replace (void * vtrav, void * new_item);
static void avl_trav_refresh (avl_traverser * trav);

#if 0
static avl_table * avl_copy (const avl_table *, avl_copy_func *,
                             avl_item_func *);
static void *avl_replace (avl_table *, void *);
#endif

static int
avl_init(const gsl_bst_allocator * allocator,
         gsl_bst_cmp_function * compare, void * params, void * vtable)
{
  avl_table * table = (avl_table *) vtable;

  table->avl_alloc = allocator;
  table->avl_compare = compare;
  table->avl_param = params;

  table->avl_root = NULL;
  table->avl_count = 0;
  table->avl_generation = 0;

  return GSL_SUCCESS;
}

static size_t
avl_nodes (const void * vtable)
{
  const avl_table * table = (const avl_table *) vtable;
  return table->avl_count;
}

/* empty tree (delete all nodes) but do not free the tree itself */
static int
avl_empty (void * vtable)
{
  avl_table * table = (avl_table *) vtable;
  avl_node *p, *q;

  for (p = table->avl_root; p != NULL; p = q)
    {
      if (p->avl_link[0] == NULL)
        {
          q = p->avl_link[1];
          table->avl_alloc->free (p, table->avl_param);
        }
      else
        {
          q = p->avl_link[0];
          p->avl_link[0] = q->avl_link[1];
          q->avl_link[1] = p;
        }
    }

  table->avl_root = NULL;
  table->avl_count = 0;
  table->avl_generation = 0;

  return GSL_SUCCESS;
}

/*
avl_probe()
  Inserts |item| into |tree| and returns a pointer to |item|'s address.
If a duplicate item is found in the tree, returns a pointer to the existing
item without inserting |item|.
Returns |NULL| in case of memory allocation failure.
*/

static void **
avl_probe (void * item, avl_table * table)
{
  avl_node *y, *z; /* top node to update balance factor, and parent */
  avl_node *p, *q; /* iterator, and parent */
  avl_node *n;     /* newly inserted node */
  avl_node *w;     /* new root of rebalanced subtree */
  int dir;         /* direction to descend */

  unsigned char da[AVL_MAX_HEIGHT]; /* cached comparison results */
  int k = 0;                        /* number of cached results */

  z = (avl_node *) &table->avl_root;
  y = table->avl_root;
  dir = 0;
  for (q = z, p = y; p != NULL; q = p, p = p->avl_link[dir])
    {
      int cmp = table->avl_compare (item, p->avl_data, table->avl_param);

      if (cmp == 0)
        return &p->avl_data;

      if (p->avl_balance != 0)
        z = q, y = p, k = 0;

      da[k++] = dir = cmp > 0;
    }

  /* allocate a new node */
  n = q->avl_link[dir] = table->avl_alloc->alloc (sizeof *n, table->avl_param);
  if (n == NULL)
    return NULL;

  table->avl_count++;

  n->avl_data = item;
  n->avl_link[0] = n->avl_link[1] = NULL;
  n->avl_balance = 0;

  if (y == NULL)
    return &n->avl_data;

  for (p = y, k = 0; p != n; p = p->avl_link[da[k]], k++)
    if (da[k] == 0)
      p->avl_balance--;
    else
      p->avl_balance++;

  if (y->avl_balance == -2)
    {
      avl_node *x = y->avl_link[0];
      if (x->avl_balance == -1)
        {
          w = x;
          y->avl_link[0] = x->avl_link[1];
          x->avl_link[1] = y;
          x->avl_balance = y->avl_balance = 0;
        }
      else
        {
          w = x->avl_link[1];
          x->avl_link[1] = w->avl_link[0];
          w->avl_link[0] = x;
          y->avl_link[0] = w->avl_link[1];
          w->avl_link[1] = y;
          if (w->avl_balance == -1)
            x->avl_balance = 0, y->avl_balance = +1;
          else if (w->avl_balance == 0)
            x->avl_balance = y->avl_balance = 0;
          else /* |w->avl_balance == +1| */
            x->avl_balance = -1, y->avl_balance = 0;
          w->avl_balance = 0;
        }
    }
  else if (y->avl_balance == +2)
    {
      avl_node *x = y->avl_link[1];
      if (x->avl_balance == +1)
        {
          w = x;
          y->avl_link[1] = x->avl_link[0];
          x->avl_link[0] = y;
          x->avl_balance = y->avl_balance = 0;
        }
      else
        {
          w = x->avl_link[0];
          x->avl_link[0] = w->avl_link[1];
          w->avl_link[1] = x;
          y->avl_link[1] = w->avl_link[0];
          w->avl_link[0] = y;
          if (w->avl_balance == +1)
            x->avl_balance = 0, y->avl_balance = -1;
          else if (w->avl_balance == 0)
            x->avl_balance = y->avl_balance = 0;
          else /* |w->avl_balance == -1| */
            x->avl_balance = +1, y->avl_balance = 0;
          w->avl_balance = 0;
        }
    }
  else
    return &n->avl_data;

  z->avl_link[y != z->avl_link[0]] = w;
  table->avl_generation++;

  return &n->avl_data;
}

/*
avl_insert()
  Inserts |item| into |table|. Returns NULL if item successfully inserted,
or if a memory allocation error occurred. Otherwise, return duplicate
item.
*/

static void *
avl_insert (void * item, void * vtable)
{
  void **p = avl_probe (item, vtable);
  return p == NULL || *p == item ? NULL : *p;
}

/*
avl_find()
  Search for |item| in |table| and return a pointer to item if found.
Return NULL if not found.
*/

static void *
avl_find (const void * item, const void * vtable)
{
  const avl_table * table = (const avl_table *) vtable;
  avl_node *p;

  for (p = table->avl_root; p != NULL; )
    {
      int cmp = table->avl_compare (item, p->avl_data, table->avl_param);

      if (cmp < 0)
        p = p->avl_link[0];
      else if (cmp > 0)
        p = p->avl_link[1];
      else /* |cmp == 0| */
        return p->avl_data;
    }

  return NULL;
}

/*
avl_remove()
  Deletes from |table| and returns an item matching |item|.
Returns a null pointer if no matching item found.
*/

static void *
avl_remove (const void * item, void * vtable)
{
  avl_table * table = (avl_table *) vtable;

  /* stack of nodes */
  avl_node *pa[AVL_MAX_HEIGHT];      /* nodes */
  unsigned char da[AVL_MAX_HEIGHT];  /* |link[]| indexes */
  int k;                             /* stack pointer */

  avl_node *p;                       /* traverses tree to find node to delete */
  int cmp;                           /* result of comparison between |item| and |p| */

  k = 0;
  p = (avl_node *) &table->avl_root;
  for (cmp = -1; cmp != 0;
       cmp = table->avl_compare (item, p->avl_data, table->avl_param))
    {
      int dir = cmp > 0;

      pa[k] = p;
      da[k++] = dir;

      p = p->avl_link[dir];
      if (p == NULL)
        return NULL;
    }

  item = p->avl_data;

  if (p->avl_link[1] == NULL)
    pa[k - 1]->avl_link[da[k - 1]] = p->avl_link[0];
  else
    {
      avl_node *r = p->avl_link[1];
      if (r->avl_link[0] == NULL)
        {
          r->avl_link[0] = p->avl_link[0];
          r->avl_balance = p->avl_balance;
          pa[k - 1]->avl_link[da[k - 1]] = r;
          da[k] = 1;
          pa[k++] = r;
        }
      else
        {
          avl_node *s;
          int j = k++;

          for (;;)
            {
              da[k] = 0;
              pa[k++] = r;
              s = r->avl_link[0];
              if (s->avl_link[0] == NULL)
                break;

              r = s;
            }

          s->avl_link[0] = p->avl_link[0];
          r->avl_link[0] = s->avl_link[1];
          s->avl_link[1] = p->avl_link[1];
          s->avl_balance = p->avl_balance;

          pa[j - 1]->avl_link[da[j - 1]] = s;
          da[j] = 1;
          pa[j] = s;
        }
    }

  table->avl_alloc->free (p, table->avl_param);

  while (--k > 0)
    {
      avl_node *y = pa[k];

      if (da[k] == 0)
        {
          y->avl_balance++;
          if (y->avl_balance == +1)
            break;
          else if (y->avl_balance == +2)
            {
              avl_node *x = y->avl_link[1];
              if (x->avl_balance == -1)
                {
                  avl_node *w;
                  w = x->avl_link[0];
                  x->avl_link[0] = w->avl_link[1];
                  w->avl_link[1] = x;
                  y->avl_link[1] = w->avl_link[0];
                  w->avl_link[0] = y;
                  if (w->avl_balance == +1)
                    x->avl_balance = 0, y->avl_balance = -1;
                  else if (w->avl_balance == 0)
                    x->avl_balance = y->avl_balance = 0;
                  else /* |w->avl_balance == -1| */
                    x->avl_balance = +1, y->avl_balance = 0;
                  w->avl_balance = 0;
                  pa[k - 1]->avl_link[da[k - 1]] = w;
                }
              else
                {
                  y->avl_link[1] = x->avl_link[0];
                  x->avl_link[0] = y;
                  pa[k - 1]->avl_link[da[k - 1]] = x;
                  if (x->avl_balance == 0)
                    {
                      x->avl_balance = -1;
                      y->avl_balance = +1;
                      break;
                    }
                  else
                    x->avl_balance = y->avl_balance = 0;
                }
            }
        }
      else
        {
          y->avl_balance--;
          if (y->avl_balance == -1)
            break;
          else if (y->avl_balance == -2)
            {
              avl_node *x = y->avl_link[0];
              if (x->avl_balance == +1)
                {
                  avl_node *w;
                  w = x->avl_link[1];
                  x->avl_link[1] = w->avl_link[0];
                  w->avl_link[0] = x;
                  y->avl_link[0] = w->avl_link[1];
                  w->avl_link[1] = y;
                  if (w->avl_balance == -1)
                    x->avl_balance = 0, y->avl_balance = +1;
                  else if (w->avl_balance == 0)
                    x->avl_balance = y->avl_balance = 0;
                  else /* |w->avl_balance == +1| */
                    x->avl_balance = -1, y->avl_balance = 0;
                  w->avl_balance = 0;
                  pa[k - 1]->avl_link[da[k - 1]] = w;
                }
              else
                {
                  y->avl_link[0] = x->avl_link[1];
                  x->avl_link[1] = y;
                  pa[k - 1]->avl_link[da[k - 1]] = x;
                  if (x->avl_balance == 0)
                    {
                      x->avl_balance = +1;
                      y->avl_balance = -1;
                      break;
                    }
                  else
                    x->avl_balance = y->avl_balance = 0;
                }
            }
        }
    }

  table->avl_count--;
  table->avl_generation++;

  return (void *) item;
}

#if 0
/* Inserts |item| into |table|, replacing any duplicate item.
   Returns |NULL| if |item| was inserted without replacing a duplicate,
   or if a memory allocation error occurred.
   Otherwise, returns the item that was replaced. */
static void *
avl_replace (avl_table *table, void *item)
{
  void **p = avl_probe (table, item);
  if (p == NULL || *p == item)
    return NULL;
  else
    {
      void *r = *p;
      *p = item;
      return r;
    }
}

/* Destroys |new| with |avl_destroy (new, destroy)|,
   first setting right links of nodes in |stack| within |new|
   to null pointers to avoid touching uninitialized data. */
static void
copy_error_recovery (avl_node **stack, int height,
                     avl_table *new, avl_item_func *destroy)
{
  for (; height > 2; height -= 2)
    stack[height - 1]->avl_link[1] = NULL;
  avl_destroy (new, destroy);
}

/* Copies |org| to a newly created tree, which is returned.
   If |copy != NULL|, each data item in |org| is first passed to |copy|,
   and the return values are inserted into the tree,
   with |NULL| return values taken as indications of failure.
   On failure, destroys the partially created new tree,
   applying |destroy|, if non-null, to each item in the new tree so far,
   and returns |NULL|.
   If |allocator != NULL|, it is used for allocation in the new tree.
   Otherwise, the same allocator used for |org| is used. */
static avl_table *
avl_copy (const avl_table *org, avl_copy_func *copy,
          avl_item_func *destroy, struct libavl_allocator *allocator)
{
  avl_node *stack[2 * (AVL_MAX_HEIGHT + 1)];
  int height = 0;

  avl_table *new;
  const avl_node *x;
  avl_node *y;

  new = avl_alloc (org->avl_compare, org->avl_param,
                   allocator != NULL ? allocator : org->avl_alloc);
  if (new == NULL)
    return NULL;
  new->avl_count = org->avl_count;
  if (new->avl_count == 0)
    return new;

  x = (const avl_node *) &org->avl_root;
  y = (avl_node *) &new->avl_root;
  for (;;)
    {
      while (x->avl_link[0] != NULL)
        {
          y->avl_link[0] =
            new->allocator->alloc (sizeof *y->avl_link[0],
                                   new->avl_param);
          if (y->avl_link[0] == NULL)
            {
              if (y != (avl_node *) &new->avl_root)
                {
                  y->avl_data = NULL;
                  y->avl_link[1] = NULL;
                }

              copy_error_recovery (stack, height, new, destroy);
              return NULL;
            }

          stack[height++] = (avl_node *) x;
          stack[height++] = y;
          x = x->avl_link[0];
          y = y->avl_link[0];
        }
      y->avl_link[0] = NULL;

      for (;;)
        {
          y->avl_balance = x->avl_balance;
          if (copy == NULL)
            y->avl_data = x->avl_data;
          else
            {
              y->avl_data = copy (x->avl_data, org->avl_param);
              if (y->avl_data == NULL)
                {
                  y->avl_link[1] = NULL;
                  copy_error_recovery (stack, height, new, destroy);
                  return NULL;
                }
            }

          if (x->avl_link[1] != NULL)
            {
              y->avl_link[1] =
                new->allocator->alloc (sizeof *y->avl_link[1],
                                       new->avl_param);
              if (y->avl_link[1] == NULL)
                {
                  copy_error_recovery (stack, height, new, destroy);
                  return NULL;
                }

              x = x->avl_link[1];
              y = y->avl_link[1];
              break;
            }
          else
            y->avl_link[1] = NULL;

          if (height <= 2)
            return new;

          y = stack[--height];
          x = stack[--height];
        }
    }
}

#endif

/*
avl_t_init()
  Initializes |trav| for use with |tree| and selects the null node.
*/

static int
avl_t_init (void * vtrav, const void * vtable)
{
  avl_traverser * trav = (avl_traverser *) vtrav;
  const avl_table * table = (const avl_table *) vtable;

  trav->avl_table = table;
  trav->avl_node = NULL;
  trav->avl_height = 0;
  trav->avl_generation = table->avl_generation;

  return GSL_SUCCESS;
}

/*
avl_t_first()
  Initializes |trav| for |tree| and selects and returns a pointer to its least-valued item.
  Returns |NULL| if |tree| contains no nodes.
*/

static void *
avl_t_first (void * vtrav, const void * vtable)
{
  const avl_table * table = (const avl_table *) vtable;
  avl_traverser * trav = (avl_traverser *) vtrav;
  avl_node *x;

  trav->avl_table = table;
  trav->avl_height = 0;
  trav->avl_generation = table->avl_generation;

  x = table->avl_root;
  if (x != NULL)
    {
      while (x->avl_link[0] != NULL)
        {
          if (trav->avl_height >= AVL_MAX_HEIGHT)
            {
              GSL_ERROR_NULL("traverser height exceeds maximum", GSL_ETABLE);
            }

          trav->avl_stack[trav->avl_height++] = x;
          x = x->avl_link[0];
        }
    }

  trav->avl_node = x;

  return x != NULL ? x->avl_data : NULL;
}

/*
avl_t_last()
  Initializes |trav| for |tree| and selects and returns a pointer to its greatest-valued item.
  Returns |NULL| if |tree| contains no nodes.
*/

static void *
avl_t_last (void * vtrav, const void * vtable)
{
  const avl_table * table = (const avl_table *) vtable;
  avl_traverser * trav = (avl_traverser *) vtrav;
  avl_node *x;

  trav->avl_table = table;
  trav->avl_height = 0;
  trav->avl_generation = table->avl_generation;

  x = table->avl_root;
  if (x != NULL)
    {
      while (x->avl_link[1] != NULL)
        {
          if (trav->avl_height >= AVL_MAX_HEIGHT)
            {
              GSL_ERROR_NULL("traverser height exceeds maximum", GSL_ETABLE);
            }

          trav->avl_stack[trav->avl_height++] = x;
          x = x->avl_link[1];
        }
    }

  trav->avl_node = x;

  return x != NULL ? x->avl_data : NULL;
}

/*
avl_t_find()
  Searches for |item| in |table|. If found, initializes |trav| to
the item found and returns the item as well.  If there is no matching
item, initializes |trav| to the null item and returns |NULL|.
*/

static void *
avl_t_find (const void * item, void * vtrav, const void * vtable)
{
  const avl_table * table = (const avl_table *) vtable;
  avl_traverser * trav = (avl_traverser *) vtrav;
  avl_node *p, *q;

  trav->avl_table = table;
  trav->avl_height = 0;
  trav->avl_generation = table->avl_generation;

  for (p = table->avl_root; p != NULL; p = q)
    {
      int cmp = table->avl_compare (item, p->avl_data, table->avl_param);

      if (cmp < 0)
        q = p->avl_link[0];
      else if (cmp > 0)
        q = p->avl_link[1];
      else /* |cmp == 0| */
        {
          trav->avl_node = p;
          return p->avl_data;
        }

      if (trav->avl_height >= AVL_MAX_HEIGHT)
        {
          GSL_ERROR_NULL("traverser height exceeds maximum", GSL_ETABLE);
        }

      trav->avl_stack[trav->avl_height++] = p;
    }

  trav->avl_height = 0;
  trav->avl_node = NULL;

  return NULL;
}

/*
avl_t_insert()
  Attempts to insert |item| into |table|.  If |item| is inserted
successfully, it is returned and |trav| is initialized to its location.
If a duplicate is found, it is returned and |trav| is initialized to
its location.  No replacement of the item occurs.
If a memory allocation failure occurs, |NULL| is returned and |trav|
is initialized to the null item.
*/

static void *
avl_t_insert (void * item, void * vtrav, void * vtable)
{
  avl_table * table = (avl_table *) vtable;
  avl_traverser * trav = (avl_traverser *) vtrav;
  void **p;

  p = avl_probe (item, table);
  if (p != NULL)
    {
      trav->avl_table = table;
      trav->avl_node = ((avl_node *) ((char *) p - offsetof (avl_node, avl_data)));
      trav->avl_generation = table->avl_generation - 1;
      return *p;
    }
  else
    {
      avl_t_init (vtrav, vtable);
      return NULL;
    }
}

/*
avl_t_copy()
  Initializes |trav| to have the same current node as |src|.
*/

static void *
avl_t_copy (void * vtrav, const void * vsrc)
{
  const avl_traverser * src = (const avl_traverser *) vsrc;
  avl_traverser * trav = (avl_traverser *) vtrav;

  if (trav != src)
    {
      trav->avl_table = src->avl_table;
      trav->avl_node = src->avl_node;
      trav->avl_generation = src->avl_generation;
      if (trav->avl_generation == trav->avl_table->avl_generation)
        {
          trav->avl_height = src->avl_height;
          memcpy (trav->avl_stack, (const void *) src->avl_stack,
                  sizeof *trav->avl_stack * trav->avl_height);
        }
    }

  return trav->avl_node != NULL ? trav->avl_node->avl_data : NULL;
}

/*
avl_t_next()
  Returns the next data item in in-order within the tree
being traversed with |trav|, or if there are no more data items returns NULL.
*/

static void *
avl_t_next (void * vtrav)
{
  avl_traverser * trav = (avl_traverser *) vtrav;
  avl_node *x;

  if (trav->avl_generation != trav->avl_table->avl_generation)
    avl_trav_refresh (trav);

  x = trav->avl_node;
  if (x == NULL)
    {
      return avl_t_first (vtrav, trav->avl_table);
    }
  else if (x->avl_link[1] != NULL)
    {
      if (trav->avl_height >= AVL_MAX_HEIGHT)
        {
          GSL_ERROR_NULL("traverser height exceeds maximum", GSL_ETABLE);
        }

      trav->avl_stack[trav->avl_height++] = x;
      x = x->avl_link[1];

      while (x->avl_link[0] != NULL)
        {
          if (trav->avl_height >= AVL_MAX_HEIGHT)
            {
              GSL_ERROR_NULL("traverser height exceeds maximum", GSL_ETABLE);
            }

          trav->avl_stack[trav->avl_height++] = x;
          x = x->avl_link[0];
        }
    }
  else
    {
      avl_node *y;

      do
        {
          if (trav->avl_height == 0)
            {
              trav->avl_node = NULL;
              return NULL;
            }

          y = x;
          x = trav->avl_stack[--trav->avl_height];
        }
      while (y == x->avl_link[1]);
    }

  trav->avl_node = x;

  return x->avl_data;
}

/*
avl_t_prev()
  Returns the previous data item in inorder within the tree being
traversed with |trav|, or if there are no more data items returns NULL.
*/

static void *
avl_t_prev (void * vtrav)
{
  avl_traverser * trav = (avl_traverser *) vtrav;
  avl_node *x;

  if (trav->avl_generation != trav->avl_table->avl_generation)
    avl_trav_refresh (trav);

  x = trav->avl_node;
  if (x == NULL)
    {
      return avl_t_last (vtrav, trav->avl_table);
    }
  else if (x->avl_link[0] != NULL)
    {
      if (trav->avl_height >= AVL_MAX_HEIGHT)
        {
          GSL_ERROR_NULL("traverser height exceeds maximum", GSL_ETABLE);
        }

      trav->avl_stack[trav->avl_height++] = x;
      x = x->avl_link[0];

      while (x->avl_link[1] != NULL)
        {
          if (trav->avl_height >= AVL_MAX_HEIGHT)
            {
              GSL_ERROR_NULL("traverser height exceeds maximum", GSL_ETABLE);
            }

          trav->avl_stack[trav->avl_height++] = x;
          x = x->avl_link[1];
        }
    }
  else
    {
      avl_node *y;

      do
        {
          if (trav->avl_height == 0)
            {
              trav->avl_node = NULL;
              return NULL;
            }

          y = x;
          x = trav->avl_stack[--trav->avl_height];
        }
      while (y == x->avl_link[0]);
    }

  trav->avl_node = x;

  return x->avl_data;
}

static void *
avl_t_cur (const void * vtrav)
{
  const avl_traverser * trav = (const avl_traverser *) vtrav;
  return trav->avl_node != NULL ? trav->avl_node->avl_data : NULL;
}

/*
avl_t_replace()
  Replaces the current item in |trav| by |new| and returns the item replaced.
|trav| must not have the null item selected. The new item must not upset the
ordering of the tree.
*/

static void *
avl_t_replace (void * vtrav, void * new_item)
{
  avl_traverser * trav = (avl_traverser *) vtrav;
  void *old;

  old = trav->avl_node->avl_data;
  trav->avl_node->avl_data = new_item;

  return old;
}

/*
avl_trav_refresh()
  Refreshes the stack of parent pointers in |trav| and updates its generation number
*/

static void
avl_trav_refresh (avl_traverser * trav)
{
  trav->avl_generation = trav->avl_table->avl_generation;

  if (trav->avl_node != NULL)
    {
      gsl_bst_cmp_function *cmp = trav->avl_table->avl_compare;
      void *param = trav->avl_table->avl_param;
      avl_node *node = trav->avl_node;
      avl_node *i;

      trav->avl_height = 0;
      for (i = trav->avl_table->avl_root; i != node; )
        {
          if (trav->avl_height >= AVL_MAX_HEIGHT)
            {
              GSL_ERROR_VOID("traverser height exceeds maximum", GSL_ETABLE);
            }

          trav->avl_stack[trav->avl_height++] = i;
          i = i->avl_link[cmp (node->avl_data, i->avl_data, param) > 0];
        }
    }
}

static const gsl_bst_type avl_tree_type =
{
  "AVL",
  sizeof(avl_node),
  avl_init,
  avl_nodes,
  avl_insert,
  avl_find,
  avl_remove,
  avl_empty,

  avl_t_init,
  avl_t_first,
  avl_t_last,
  avl_t_find,
  avl_t_insert,
  avl_t_copy,
  avl_t_next,
  avl_t_prev,
  avl_t_cur,
  avl_t_replace
};

const gsl_bst_type * gsl_bst_avl = &avl_tree_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/**
 * @brief handle cosmological functions like power spectrum, growth, etc.
 *
 * @file core_power.h
 * @author Michal Vrastil
 * @date 2018-06-21
 */

#pragma once
#include "precision.hpp"
#include <gsl/gsl_spline.h>

/*******************//**
* FORWARD DECLARATIONS *
************************/

class Cosmo_Param; ///< declaration in params.hpp
class Sim_Param; ///< declaration in params.hpp
template <typename T, size_t N> class Data_Vec; ///< declaration in class_data_vec.hpp

/**
 * @brief Initialize CCL power spectrum
 *
 * @param cosmo cosmological parameters
 */
void norm_pwr(Cosmo_Param& cosmo);

/**
 * @brief when computing growth factor outside CCL range we need to normalize the growth factor; \f$D(a=1)\equiv1\f$
 *
 * @param cosmo cosmological parameters
 * @return FTYPE_t unnormalized growth factor at 'a = 1'
 */
FTYPE_t norm_growth_factor(const Cosmo_Param& cosmo);

/**
 * @brief
 *
 * @param a
 * @param cosmo cosmological parameters
 * @return FTYPE_t
 */
FTYPE_t growth_factor(FTYPE_t a, const Cosmo_Param& cosmo);

/**
 * @brief
 *
 * @param a
 * @param cosmo cosmological parameters
 * @return FTYPE_t
 */
FTYPE_t growth_rate(FTYPE_t a, const Cosmo_Param& cosmo);

/**
 * @brief
 *
 * @param a
 * @param cosmo cosmological parameters
 * @return FTYPE_t
 */
FTYPE_t growth_change(FTYPE_t a, const Cosmo_Param& cosmo);

/**
 * @brief
 *
 * @param a
 * @param cosmo cosmological parameters
 * @return FTYPE_t
 */
FTYPE_t Omega_lambda(FTYPE_t a, const Cosmo_Param& cosmo);

/**
 * @brief
 *
 * @param a
 * @param k
 * @param cosmo cosmological parameters
 * @return FTYPE_t
 */
FTYPE_t lin_pow_spec(FTYPE_t a, FTYPE_t k, const Cosmo_Param& cosmo);

/**
 * @brief
 *
 * @param a
 * @param k
 * @param cosmo cosmological parameters
 * @return FTYPE_t
 */
FTYPE_t non_lin_pow_spec(FTYPE_t a, FTYPE_t k, const Cosmo_Param& cosmo);

/**
 * @class:	Interp_obj
 * @brief:	linear interpolation (Steffen) of data [x, y]
 */
class Interp_obj
{
public:
    Interp_obj(): is_init(false) {}
    ~Interp_obj();
    double operator()(double x) const;
    template <typename T, size_t N>
    void init(const Data_Vec<T, N>& data);
    double x_min, x_max;

private:
    bool is_init;
    gsl_spline* spline;
    gsl_interp_accel* acc;
};

/**
 * @class:	Extrap_Pk
 * @brief:	linear interpolation of data [k, P(k)] within 'useful' range
 *          fit to primordial P_i(k) below the 'useful' range
 *          fit to Padé approximant R [0/3] above the 'useful' range
 *
 *  Steffen interpolation of data [k, P(k)] within range k_min, k_max
 *  fit to primordial P_i(k) below this range, fit A*k^ns above
 */
template <typename T, size_t N>
class Extrap_Pk : public Interp_obj
{
public:
    Extrap_Pk(const Data_Vec<T, N>& data, const Sim_Param& sim);
    Extrap_Pk(const Data_Vec<T, N>& data, const Sim_Param& sim, const size_t m_l, const size_t n_u);
    Extrap_Pk(const Data_Vec<T, N>& data, const Sim_Param& sim, const size_t m_l, const size_t n_l,
              const size_t m_u, const size_t n_u);
    double operator()(double k) const;

    void fit_lin(const Data_Vec<T, N>& data, const size_t m, const size_t n, double& A);
    void fit_power_law(const Data_Vec<T, N>& data, const size_t m, const size_t n, double& A, double& n_s);

    double A_low; ///< amplitude of linear power in lower range
    const Cosmo_Param& cosmo;
    double A_up, n_s; ///< scale-free power spectrum in upper range
    T k_min, k_max; ///< interpolation range
};

/**
 * @class:	Extrap_Pk_Nl
 * @brief:	creates Extrapolate object (linear power spectrum) from data and store non-linear parameters
            call 'operator()(k)' based on k_split (upper range of the linear)
 */
template <typename T, size_t N>
class Extrap_Pk_Nl : public Extrap_Pk<T, N>
{
public:
    Extrap_Pk_Nl(const Data_Vec<T, N>& data, const Sim_Param &sim, T A_nl, T a_eff);
    const T A_nl, a_eff, k_split;
    double operator()(double k) const;
};

/**
 * @brief compute correlation function and store results
 *
 * @tparam P callable object with 'operator()(double)'
 * @param sim simulation parameters
 * @param P_k power spectrum (callable)
 * @param corr_func_binned object to store binned correlation function
 */
template<class P>
void gen_corr_func_binned_gsl_qawf(const Sim_Param &sim, const P& P_k, Data_Vec<FTYPE_t, 2>& corr_func_binned);

/**
 * @brief compute linear correlation function and store results
 *
 * @param sim simulation parameters
 * @param a scale factor
 * @param corr_func_binned object to store binned correlation function
 */
void gen_corr_func_binned_gsl_qawf_lin(const Sim_Param &sim, FTYPE_t a, Data_Vec<FTYPE_t, 2>& corr_func_binned);

/**
 * @brief compute non-linear correlation function and store results
 *
 * @param sim simulation parameters
 * @param a scale factor
 * @param corr_func_binned object to store binned correlation function
 */
void gen_corr_func_binned_gsl_qawf_nl(const Sim_Param &sim, FTYPE_t a, Data_Vec<FTYPE_t, 2>& corr_func_binned);

/**
 * @brief compute amplitude of density fluctuation and store results
 *
 * @tparam P callable object with 'operator()(double)'
 * @param sim simulation parameters
 * @param P_k power spectrum (callable)
 * @param sigma_binned object to store binned correlation function
 */
template<class P>
void gen_sigma_binned_gsl_qawf(const Sim_Param &sim, const P& P_k, Data_Vec<FTYPE_t, 2>& sigma_binned);

/**
 * @brief compute linear amplitude of density fluctuation and store results
 *
 * @param sim simulation parameters
 * @param a scale factor
 * @param sigma_binned object to store binned correlation function
 */
void gen_sigma_func_binned_gsl_qawf_lin(const Sim_Param &sim, FTYPE_t a, Data_Vec<FTYPE_t, 2>& sigma_binned);

/**
 * @brief compute non-linear amplitude of density fluctuation and store results
 *
 * @param sim simulation parameters
 * @param a scale factor
 * @param sigma_binned object to store binned correlation function
 */
void gen_sigma_func_binned_gsl_qawf_nl(const Sim_Param &sim, FTYPE_t a, Data_Vec<FTYPE_t, 2>& sigma_binned);

max_stars_repo_name: sjvs/FastSim

max_issues_repo_name: sjvs/FastSim

max_forks_repo_name: sjvs/FastSim


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

#ifndef __GSL_MATRIX_UINT_H__
#define __GSL_MATRIX_UINT_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_uint.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 size1;
  size_t size2;
  size_t tda;
  unsigned int * data;
  gsl_block_uint * block;
  int owner;
} gsl_matrix_uint;

typedef struct
{
  gsl_matrix_uint matrix;
} _gsl_matrix_uint_view;

typedef _gsl_matrix_uint_view gsl_matrix_uint_view;

typedef struct
{
  gsl_matrix_uint matrix;
} _gsl_matrix_uint_const_view;

typedef const _gsl_matrix_uint_const_view gsl_matrix_uint_const_view;

/* Allocation */

gsl_matrix_uint *
gsl_matrix_uint_alloc (const size_t n1, const size_t n2);

gsl_matrix_uint *
gsl_matrix_uint_calloc (const size_t n1, const size_t n2);

gsl_matrix_uint *
gsl_matrix_uint_alloc_from_block (gsl_block_uint * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

gsl_matrix_uint *
gsl_matrix_uint_alloc_from_matrix (gsl_matrix_uint * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

gsl_vector_uint *
gsl_vector_uint_alloc_row_from_matrix (gsl_matrix_uint * m,
                                        const size_t i);

gsl_vector_uint *
gsl_vector_uint_alloc_col_from_matrix (gsl_matrix_uint * m,
                                        const size_t j);

void gsl_matrix_uint_free (gsl_matrix_uint * m);

/* Views */

_gsl_matrix_uint_view
gsl_matrix_uint_submatrix (gsl_matrix_uint * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

_gsl_vector_uint_view
gsl_matrix_uint_row (gsl_matrix_uint * m, const size_t i);

_gsl_vector_uint_view
gsl_matrix_uint_column (gsl_matrix_uint * m, const size_t j);

_gsl_vector_uint_view
gsl_matrix_uint_diagonal (gsl_matrix_uint * m);

_gsl_vector_uint_view
gsl_matrix_uint_subdiagonal (gsl_matrix_uint * m, const size_t k);

_gsl_vector_uint_view
gsl_matrix_uint_superdiagonal (gsl_matrix_uint * m, const size_t k);

_gsl_matrix_uint_view
gsl_matrix_uint_view_array (unsigned int * base,
                             const size_t n1,
                             const size_t n2);

_gsl_matrix_uint_view
gsl_matrix_uint_view_array_with_tda (unsigned int * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

_gsl_matrix_uint_view
gsl_matrix_uint_view_vector (gsl_vector_uint * v,
                              const size_t n1,
                              const size_t n2);

_gsl_matrix_uint_view
gsl_matrix_uint_view_vector_with_tda (gsl_vector_uint * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

_gsl_matrix_uint_const_view
gsl_matrix_uint_const_submatrix (const gsl_matrix_uint * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

_gsl_vector_uint_const_view
gsl_matrix_uint_const_row (const gsl_matrix_uint * m,
                            const size_t i);

_gsl_vector_uint_const_view
gsl_matrix_uint_const_column (const gsl_matrix_uint * m,
                               const size_t j);

_gsl_vector_uint_const_view
gsl_matrix_uint_const_diagonal (const gsl_matrix_uint * m);

_gsl_vector_uint_const_view
gsl_matrix_uint_const_subdiagonal (const gsl_matrix_uint * m,
                                    const size_t k);

_gsl_vector_uint_const_view
gsl_matrix_uint_const_superdiagonal (const gsl_matrix_uint * m,
                                      const size_t k);

_gsl_matrix_uint_const_view
gsl_matrix_uint_const_view_array (const unsigned int * base,
                                   const size_t n1,
                                   const size_t n2);

_gsl_matrix_uint_const_view
gsl_matrix_uint_const_view_array_with_tda (const unsigned int * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

_gsl_matrix_uint_const_view
gsl_matrix_uint_const_view_vector (const gsl_vector_uint * v,
                                    const size_t n1,
                                    const size_t n2);

_gsl_matrix_uint_const_view
gsl_matrix_uint_const_view_vector_with_tda (const gsl_vector_uint * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

unsigned int   gsl_matrix_uint_get(const gsl_matrix_uint * m, const size_t i, const size_t j);
void    gsl_matrix_uint_set(gsl_matrix_uint * m, const size_t i, const size_t j, const unsigned int x);

unsigned int * gsl_matrix_uint_ptr(gsl_matrix_uint * m, const size_t i, const size_t j);
const unsigned int * gsl_matrix_uint_const_ptr(const gsl_matrix_uint * m, const size_t i, const size_t j);

void gsl_matrix_uint_set_zero (gsl_matrix_uint * m);
void gsl_matrix_uint_set_identity (gsl_matrix_uint * m);
void gsl_matrix_uint_set_all (gsl_matrix_uint * m, unsigned int x);

int gsl_matrix_uint_fread (FILE * stream, gsl_matrix_uint * m) ;
int gsl_matrix_uint_fwrite (FILE * stream, const gsl_matrix_uint * m) ;
int gsl_matrix_uint_fscanf (FILE * stream, gsl_matrix_uint * m);
int gsl_matrix_uint_fprintf (FILE * stream, const gsl_matrix_uint * m, const char * format);

int gsl_matrix_uint_memcpy(gsl_matrix_uint * dest, const gsl_matrix_uint * src);
int gsl_matrix_uint_swap(gsl_matrix_uint * m1, gsl_matrix_uint * m2);

int gsl_matrix_uint_swap_rows(gsl_matrix_uint * m, const size_t i, const size_t j);
int gsl_matrix_uint_swap_columns(gsl_matrix_uint * m, const size_t i, const size_t j);
int gsl_matrix_uint_swap_rowcol(gsl_matrix_uint * m, const size_t i, const size_t j);
int gsl_matrix_uint_transpose (gsl_matrix_uint * m);
int gsl_matrix_uint_transpose_memcpy (gsl_matrix_uint * dest, const gsl_matrix_uint * src);

unsigned int gsl_matrix_uint_max (const gsl_matrix_uint * m);
unsigned int gsl_matrix_uint_min (const gsl_matrix_uint * m);
void gsl_matrix_uint_minmax (const gsl_matrix_uint * m, unsigned int * min_out, unsigned int * max_out);

void gsl_matrix_uint_max_index (const gsl_matrix_uint * m, size_t * imax, size_t *jmax);
void gsl_matrix_uint_min_index (const gsl_matrix_uint * m, size_t * imin, size_t *jmin);
void gsl_matrix_uint_minmax_index (const gsl_matrix_uint * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

int gsl_matrix_uint_isnull (const gsl_matrix_uint * m);

int gsl_matrix_uint_add (gsl_matrix_uint * a, const gsl_matrix_uint * b);
int gsl_matrix_uint_sub (gsl_matrix_uint * a, const gsl_matrix_uint * b);
int gsl_matrix_uint_mul_elements (gsl_matrix_uint * a, const gsl_matrix_uint * b);
int gsl_matrix_uint_div_elements (gsl_matrix_uint * a, const gsl_matrix_uint * b);
int gsl_matrix_uint_scale (gsl_matrix_uint * a, const double x);
int gsl_matrix_uint_add_constant (gsl_matrix_uint * a, const double x);
int gsl_matrix_uint_add_diagonal (gsl_matrix_uint * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
int gsl_matrix_uint_get_row(gsl_vector_uint * v, const gsl_matrix_uint * m, const size_t i);
int gsl_matrix_uint_get_col(gsl_vector_uint * v, const gsl_matrix_uint * m, const size_t j);
int gsl_matrix_uint_set_row(gsl_matrix_uint * m, const size_t i, const gsl_vector_uint * v);
int gsl_matrix_uint_set_col(gsl_matrix_uint * m, const size_t j, const gsl_vector_uint * v);

/* inline functions if you are using GCC */

#ifdef HAVE_INLINE
extern inline
unsigned int
gsl_matrix_uint_get(const gsl_matrix_uint * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
    }
#endif
  return m->data[i * m->tda + j] ;
}

extern inline
void
gsl_matrix_uint_set(gsl_matrix_uint * m, const size_t i, const size_t j, const unsigned int x)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
    }
#endif
  m->data[i * m->tda + j] = x ;
}

extern inline
unsigned int *
gsl_matrix_uint_ptr(gsl_matrix_uint * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (unsigned int *) (m->data + (i * m->tda + j)) ;
}

extern inline
const unsigned int *
gsl_matrix_uint_const_ptr(const gsl_matrix_uint * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (const unsigned int *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_UINT_H__ */

max_stars_repo_name: andrewkern/segSiteHMM

max_issues_repo_name: andrewkern/segSiteHMM

max_forks_repo_name: andrewkern/segSiteHMM


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

#ifndef __GSL_STATISTICS_USHORT_H__
#define __GSL_STATISTICS_USHORT_H__

#include <stddef.h>
#include <gsl/gsl_types.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

GSL_EXPORT double gsl_stats_ushort_mean (const unsigned short data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_ushort_variance (const unsigned short data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_ushort_sd (const unsigned short data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_ushort_variance_with_fixed_mean (const unsigned short data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_ushort_sd_with_fixed_mean (const unsigned short data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_ushort_absdev (const unsigned short data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_ushort_skew (const unsigned short data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_ushort_kurtosis (const unsigned short data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_ushort_lag1_autocorrelation (const unsigned short data[], const size_t stride, const size_t n);

GSL_EXPORT double gsl_stats_ushort_covariance (const unsigned short data1[], const size_t stride1,const unsigned short data2[], const size_t stride2, const size_t n);

GSL_EXPORT double gsl_stats_ushort_variance_m (const unsigned short data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_ushort_sd_m (const unsigned short data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_ushort_absdev_m (const unsigned short data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_ushort_skew_m_sd (const unsigned short data[], const size_t stride, const size_t n, const double mean, const double sd);
GSL_EXPORT double gsl_stats_ushort_kurtosis_m_sd (const unsigned short data[], const size_t stride, const size_t n, const double mean, const double sd);
GSL_EXPORT double gsl_stats_ushort_lag1_autocorrelation_m (const unsigned short data[], const size_t stride, const size_t n, const double mean);

GSL_EXPORT double gsl_stats_ushort_covariance_m (const unsigned short data1[], const size_t stride1,const unsigned short data2[], const size_t stride2, const size_t n, const double mean1, const double mean2);

GSL_EXPORT double gsl_stats_ushort_pvariance (const unsigned short data1[], const size_t stride1, const size_t n1, const unsigned short data2[], const size_t stride2, const size_t n2);
GSL_EXPORT double gsl_stats_ushort_ttest (const unsigned short data1[], const size_t stride1, const size_t n1, const unsigned short data2[], const size_t stride2, const size_t n2);

GSL_EXPORT unsigned short gsl_stats_ushort_max (const unsigned short data[], const size_t stride, const size_t n);
GSL_EXPORT unsigned short gsl_stats_ushort_min (const unsigned short data[], const size_t stride, const size_t n);
GSL_EXPORT void gsl_stats_ushort_minmax (unsigned short * min, unsigned short * max, const unsigned short data[], const size_t stride, const size_t n);

GSL_EXPORT size_t gsl_stats_ushort_max_index (const unsigned short data[], const size_t stride, const size_t n);
GSL_EXPORT size_t gsl_stats_ushort_min_index (const unsigned short data[], const size_t stride, const size_t n);
GSL_EXPORT void gsl_stats_ushort_minmax_index (size_t * min_index, size_t * max_index, const unsigned short data[], const size_t stride, const size_t n);

GSL_EXPORT double gsl_stats_ushort_median_from_sorted_data (const unsigned short sorted_data[], const size_t stride, const size_t n) ;
GSL_EXPORT double gsl_stats_ushort_quantile_from_sorted_data (const unsigned short sorted_data[], const size_t stride, const size_t n, const double f) ;

__END_DECLS

#endif /* __GSL_STATISTICS_USHORT_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
/**
 * @file batchg_zgemv.c
 *
 * Part of API test for Batched BLAS routines.
 *
 * @author  Samuel  D. Relton
 * @author  Pedro   V. Lara
 * @author  Mawussi Zounon
 * @date 2016-06-01
 *
 * @precisions normal z -> c d s
 *
 **/

#include <cblas.h>
#include "bblas.h"

#define COMPLEX

void batchg_zgemv(
	const enum BBLAS_TRANS *trans,
	const int *m, const int *n,
	const BBLAS_Complex64_t *alpha,
	const BBLAS_Complex64_t **arrayA, const int *lda,
	const BBLAS_Complex64_t **arrayx, const int *incx,
	const BBLAS_Complex64_t *beta,
	BBLAS_Complex64_t **arrayy, const int *incy,
	const int group_count, const int *group_size,
	int* info)
{
	/* Local variables */
	char func_name[15] = "batchg_zgemv";
	int group_iter = 0;
	int offset = 0; // How many subproblems solved so far

	/* Check group_count */
	if (group_count < 0)
	{
		xerbla_batch(func_name, BBLAS_ERR_GROUP_COUNT, -1);
		return;
	}

	/* Check group_size and call fixed batch computation */
	for (group_iter = 0; group_iter < group_count; group_iter++)
	{
		if (group_size[group_iter] < 0)
		{
			xerbla_batch(func_name, BBLAS_ERR_GROUP_SIZE, group_iter);
			info[group_iter] = BBLAS_ERR_GROUP_SIZE;
			continue;
		}

		/* Call fixed batch computation on the group */
		batchf_zgemv(
			trans[group_iter],
			m[group_iter], n[group_iter],
			alpha[group_iter],
			arrayA+offset, lda[group_iter],
			arrayx+offset, incx[group_iter],
			beta[group_iter],
			arrayy+offset, incy[group_iter],
			group_size[group_iter], info[group_iter]);

		/* Update offset */
		offset += group_size[group_iter];
	} // End of group loop
}
#undef COMPLEX

max_stars_repo_name: mawussi/BBLAS-group


text:
#ifndef libceed_solids_examples_cl_options_h
#define libceed_solids_examples_cl_options_h

#include <petsc.h>
#include "../include/structs.h"

// Process general command line options
PetscErrorCode ProcessCommandLineOptions(MPI_Comm comm, AppCtx app_ctx);

#endif // libceed_solids_examples_cl_options_h

max_stars_repo_name: wence-/libCEED

max_issues_repo_name: wence-/libCEED

max_forks_repo_name: wence-/libCEED


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * 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.
 */

#pragma once

#ifndef dirty_08c91b95_68e8_437c_b721_3ad67fdc25ac_h
#define dirty_08c91b95_68e8_437c_b721_3ad67fdc25ac_h

#include <gslib/std.h>
#include <ariel/type.h>

__ariel_begin__

class dirty_list
{
public:
    typedef list<rect> rect_list;
    typedef rect_list::iterator iterator;
    typedef rect_list::const_iterator const_iterator;

protected:
    int         _width, _height;
    int         _cap;
    bool        _whole;

private:
    rect_list   _rclist;

public:
    dirty_list();
    dirty_list(int w, int h);
    int size() const { return (int)_rclist.size(); }
    void clear();
    void set_dimension(int w, int h);
    void add(rect rc);
    void set_whole() { _whole = true;}
    bool is_whole() const { return _whole; }
    bool is_dirty(const rect& rc) const;
    iterator begin() { return _rclist.begin(); }
    const_iterator begin() const { return _rclist.begin(); }
    iterator end() { return _rclist.end(); }
    const_iterator end() const { return _rclist.end(); }
};

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#ifndef QDM_PARAMETERS_H
#define QDM_PARAMETERS_H 1

#include <gsl/gsl_vector.h>

typedef struct {
  unsigned long int rng_seed;

  int acc_check;
  int burn;
  int iter;
  int knot_try;
  int spline_df;
  int thin;

  double month;
  double knot;

  double tau_high;
  double tau_low;

  double theta_min;
  double theta_tune_sd;

  double xi_high;
  double xi_low;
  double xi_prior_mean;
  double xi_prior_var;
  double xi_tune_sd;

  double bound;

  int debug;
  int tau_table;

  bool truncate;
} qdm_parameters;

void
qdm_parameters_fprint(
    FILE *f,
    const qdm_parameters *p
);

int
qdm_parameters_write(
    hid_t id,
    const qdm_parameters *p
);

int
qdm_parameters_read(
    hid_t id,
    qdm_parameters *p
);

#endif /* QDM_PARAMETERS_H */


text:
#ifndef CTETRA_TETRA_H
#define CTETRA_TETRA_H

#include <stdlib.h>
#include <gsl/gsl_vector.h>
#include "ecache.h"

typedef double (*tetra_SumFn)(double E, double E1, double E2, double E3, double E4, double num_tetra);

double tetra_SumTetra(tetra_SumFn F, double E, EnergyCache *Ecache);

void sortEs(double Es[4]);

#endif // CTETRA_TETRA_H

max_stars_repo_name: tflovorn/ctetra

max_issues_repo_name: tflovorn/ctetra

max_forks_repo_name: tflovorn/ctetra


text:
// sammlung verwendeter Funktionen
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_gamma.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <omp.h>

#ifndef HEADER_FILE_UTIL
#define HEADER_FILE_UTIL

double Norm_Diff(double *vec1,double *vec2, int dim, int p); //berechnet norm von || vec1-vec2 ||_p
void Lattice_Setup(double **positions,int  numbertoinf,int dim, double latticespacing); // Legt Koordianten auf mxn Gittermatrix der dimension numbertoinf**dim, dim
int ipow(int base, int exp);
int Kahan_Sum(int N, double input[], double *ans);
int Sign_Sum(int N, double input[], double *ans);
int Neumaier_Sum(int N, double input[N], double *ans);
void deleteSpaces(char src[], char dst[]);

void vec_zero(double *v,int  N);
void vec_zero_i(int *v,int  N);

int Update_Lattice_Position(int * pos, const double *vec, const double lattice_spacing);

// begin stepper with different external forces -------//

void VVerlet_Step(int N,  double  x[N/2] , double v[N/2], double a[N/2], double *t,
					void (*deriv) (double *y, double *ans, double t,int N));

void VVerlet_Step_Target_Square(const int N,  double  x[N/2] , double v[N/2], double a[N/2], double *t,
					void (*deriv) (double *y, double *ans, double t,int N));

void VVerlet_Step_deriv(const int N,  double  x[N/2] , double v[N/2], double a[N/2], double *t, //mach einzelnen update Schritt nach Velocity-Verlet ohne Hinderniss
					void (*derivmethod) (double *y, double *ans, double t,int N));

void VVerlet_Step_Yukawa(const int N,  double  x[N/2] , double v[N/2], double a[N/2], double *t, //mach einzelnen update Schritt nach Velocity-Verlet ohne Hinderniss
					void (*derivmethod) (double *y, double *ans, double t,int N));

void VVerlet_Step_Pore_Rectangle(const int N,  double  x[N/2] , double v[N/2], double a[N/2], double *t, // wie Verlet step, aber Teilchen reflektiert stehen falls in Target (Quadrat)
					void (*derivmethod) (double *y, double *ans, double t,int N));

void VVerlet_Step_Pore_Yukawa(const int N,  double  x[N/2] , double v[N/2], double a[N/2], double *t, //mach einzelnen update Schritt nach Velocity-Verlet ohne Hinderniss
					void (*derivmethod) (double *y, double *ans, double t,int N));

// end stepper with different external forces -------//

void VVerlet(int N, double *y, double **ans, double *t,
			void (*deriv) (double *y, double *ans, double t, int N));

void VVerlet_parallel(const int N, const double *y, double **ans, double *t,
				void (*derivmethod) (double *y, double *ans, double t, int N));

void VVerlet_parallel_burn(const int N, const  double *y, double **ans, double *t,							// Velocity Verlet für Start y, Ausgabe ans[LengthT][ORDER] zu Zeiten T
				void (*derivmethod) (double *y, double *ans, double t, int N));

void Bath_Setup(double *y, char *label, gsl_rng * r);

void Ommega_Setup(gsl_rng * r);

void Gamma_Setup(double * gamma);

void deriv(double *y, double *ans, double t, int N);
void deriv_parallel(const double *yin, double *ans, double t, int N);

void deriv_Yukawa(const double *y, double *ans, double t, int N);

#endif


text:
#include <stdio.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_matrix_double.h>
// #include "covar.h"

void print_mat(gsl_matrix* C)
{
    int i, j;
    for (i = 0; i < C->size1; i++)
    {
        printf("\n");
        for (j = 0; j < C->size2; j++)
        {
            printf("%f", C->data[C->size2*i+j]);
            if(j < C->size2 - 1) printf(", ");
        }
    }
    printf("\n");
}

int main()
{
    int i, j;

/*
    Okay, here's what we want to do
    1)  Collect *all* the data - 10 seconds at standstill should be sufficient, with
            x, y, z, yaw, pitch, roll; for both pose and twist
    2)  twist gives us euler angles (yaw pitch roll), but pose gives us quaternions,
            so we'll need to convert
    2)  This leads to a 6 by N matrix of values
    3)  We propose a cov_gen() function
            input: 6 by N matrix of observation values
            output: 6 by 6 matrix of covariance values
    4) We can put a ROS wrapper around this function as a service routine,
            i.e. the output of cov_gen() will serve as our estimate of
            the covariance matrix for loam_velodyne's odometry topic
 */

    // const int rsize = 5;
    // const int csize = 3;
    #define rsize 8
    #define csize 6

    double dataA[rsize][csize] =
        {{9, 6, 9, 3, 2, 1},
         {9, 9, 3, 4, 3, 5},
         {6, 6, 6, 6, 7, 3},
         {6, 6, 9, 4, 5, 1},
         {9, 9, 3, 4, 3, 5},
         {6, 6, 6, 6, 7, 3},
         {6, 6, 9, 4, 5, 1},
         {3, 3, 3, 6, 4, 1}};

    // double dataA[rsize][csize] =
    //     {{90, 60, 90},
    //      {90, 90, 30},
    //      {60, 60, 60},
    //      {60, 60, 90},
    //      {30, 30, 30}};

    // double dataA[5][3] =   {{ 24,   0,  30},
    //                         { 24,  30, -30},
    //                         { -6,   0,   0},
    //                         { -6,   0,  30},
    //                         {-36, -30, -30}};

    // double dataA[6][6] = {{1,1,1,1,1,1},
    //                       {2,2,2,2,2,2},
    //                       {3,3,3,3,3,3},
    //                       {4,4,4,4,4,4},
    //                       {5,5,5,5,5,5},
    //                       {6,6,6,6,6,6}};

    gsl_vector_view a, b;

    gsl_matrix *A, *C;
    A = gsl_matrix_alloc(rsize, csize);
    C = gsl_matrix_alloc(csize, csize);

    for ( i = 0; i < rsize; i++)
        for ( j = 0; j < csize; j++)
            gsl_matrix_set (A, i, j, dataA[i][j]);

    print_mat(A);

    for (i = 0; i < A->size2; i++) {
        for (j = 0; j < A->size2; j++) {
          a = gsl_matrix_column (A, i);
          b = gsl_matrix_column (A, j);
          double cov = gsl_stats_covariance(a.vector.data,
                                            a.vector.stride,
                                            b.vector.data,
                                            b.vector.stride,
                                            rsize);
          gsl_matrix_set (C, i, j, cov);
        }
      }

    print_mat(C);

    return 0;
}


text:
#ifndef ANGLE_CONVERSION_H_1QTGVG9X
#define ANGLE_CONVERSION_H_1QTGVG9X

#include <gsl/gsl>
#include <sens_loc/math/constants.h>
#include <type_traits>

namespace sens_loc {

/// Implement necessary math functionality in this namespace, which is mostly
/// used by conversions.
namespace math {

/// This function converts degree into radians
/// \tparam T precision of the type (e.g. \c float or \c double)
/// \param degree angle in degree
/// \returns the same angle in radians
/// \sa rad_to_deg
/// \pre \f$0. <= degree <= 360.\f$
/// \post \f$0. <= result <= 2\pi\f$
/// \post \c rad_to_deg(result) == \p degree
template <typename T>  // requires(Float(T))
constexpr inline T deg_to_rad(T degree) {
    static_assert(std::is_floating_point_v<T>);

    Expects(degree <= T(360.));
    Expects(degree >= T(0.));
    const T radians = degree / T(180.) * pi<T>;
    Ensures(radians >= T(0.));
    Ensures(radians <= T(2. * pi<T>));

    return radians;
}

/// This function converts radians into degrees
/// \tparam T precision of the type (e.g. \c float or \c double)
/// \param radians angle in radians
/// \returns the same angle in degree
/// \sa deg_to_rad
/// \pre \f$0. <= radians <= 2\pi\f$
/// \post \f$0. <= result <= 360.\f$
/// \post \c deg_to_rad(result) == \p radians
template <typename T>  // requires(Float(T))
constexpr inline T rad_to_deg(T radians) {
    static_assert(std::is_floating_point_v<T>);

    Expects(radians >= T(0.));
    Expects(radians <= T(2. * pi<T>));
    const T degree = radians * T(180.) / pi<T>;
    Ensures(degree <= T(360.));
    Ensures(degree >= T(0.));

    return degree;
}
}  // namespace math
}  // namespace sens_loc

#endif /* end of include guard: ANGLE_CONVERSION_H_1QTGVG9X */

max_stars_repo_name: JonasToth/depth-conversions

max_issues_repo_name: JonasToth/depth-conversions

max_forks_repo_name: JonasToth/depth-conversions


text:
/* vector/gsl_vector_char.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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.
 */

#ifndef __GSL_VECTOR_CHAR_H__
#define __GSL_VECTOR_CHAR_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 <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_char.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  char *data;
  gsl_block_char *block;
  int owner;
}
gsl_vector_char;

typedef struct
{
  gsl_vector_char vector;
} _gsl_vector_char_view;

typedef _gsl_vector_char_view gsl_vector_char_view;

typedef struct
{
  gsl_vector_char vector;
} _gsl_vector_char_const_view;

typedef const _gsl_vector_char_const_view gsl_vector_char_const_view;

/* Allocation */

GSL_FUN gsl_vector_char *gsl_vector_char_alloc (const size_t n);
GSL_FUN gsl_vector_char *gsl_vector_char_calloc (const size_t n);

GSL_FUN gsl_vector_char *gsl_vector_char_alloc_from_block (gsl_block_char * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

GSL_FUN gsl_vector_char *gsl_vector_char_alloc_from_vector (gsl_vector_char * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

GSL_FUN void gsl_vector_char_free (gsl_vector_char * v);

/* Views */

GSL_FUN _gsl_vector_char_view
gsl_vector_char_view_array (char *v, size_t n);

GSL_FUN _gsl_vector_char_view
gsl_vector_char_view_array_with_stride (char *base,
                                         size_t stride,
                                         size_t n);

GSL_FUN _gsl_vector_char_const_view
gsl_vector_char_const_view_array (const char *v, size_t n);

GSL_FUN _gsl_vector_char_const_view
gsl_vector_char_const_view_array_with_stride (const char *base,
                                               size_t stride,
                                               size_t n);

GSL_FUN _gsl_vector_char_view
gsl_vector_char_subvector (gsl_vector_char *v,
                            size_t i,
                            size_t n);

GSL_FUN _gsl_vector_char_view
gsl_vector_char_subvector_with_stride (gsl_vector_char *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_FUN _gsl_vector_char_const_view
gsl_vector_char_const_subvector (const gsl_vector_char *v,
                                  size_t i,
                                  size_t n);

GSL_FUN _gsl_vector_char_const_view
gsl_vector_char_const_subvector_with_stride (const gsl_vector_char *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_FUN void gsl_vector_char_set_zero (gsl_vector_char * v);
GSL_FUN void gsl_vector_char_set_all (gsl_vector_char * v, char x);
GSL_FUN int gsl_vector_char_set_basis (gsl_vector_char * v, size_t i);

GSL_FUN int gsl_vector_char_fread (FILE * stream, gsl_vector_char * v);
GSL_FUN int gsl_vector_char_fwrite (FILE * stream, const gsl_vector_char * v);
GSL_FUN int gsl_vector_char_fscanf (FILE * stream, gsl_vector_char * v);
GSL_FUN int gsl_vector_char_fprintf (FILE * stream, const gsl_vector_char * v,
                              const char *format);

GSL_FUN int gsl_vector_char_memcpy (gsl_vector_char * dest, const gsl_vector_char * src);

GSL_FUN int gsl_vector_char_reverse (gsl_vector_char * v);

GSL_FUN int gsl_vector_char_swap (gsl_vector_char * v, gsl_vector_char * w);
GSL_FUN int gsl_vector_char_swap_elements (gsl_vector_char * v, const size_t i, const size_t j);

GSL_FUN char gsl_vector_char_max (const gsl_vector_char * v);
GSL_FUN char gsl_vector_char_min (const gsl_vector_char * v);
GSL_FUN void gsl_vector_char_minmax (const gsl_vector_char * v, char * min_out, char * max_out);

GSL_FUN size_t gsl_vector_char_max_index (const gsl_vector_char * v);
GSL_FUN size_t gsl_vector_char_min_index (const gsl_vector_char * v);
GSL_FUN void gsl_vector_char_minmax_index (const gsl_vector_char * v, size_t * imin, size_t * imax);

GSL_FUN int gsl_vector_char_add (gsl_vector_char * a, const gsl_vector_char * b);
GSL_FUN int gsl_vector_char_sub (gsl_vector_char * a, const gsl_vector_char * b);
GSL_FUN int gsl_vector_char_mul (gsl_vector_char * a, const gsl_vector_char * b);
GSL_FUN int gsl_vector_char_div (gsl_vector_char * a, const gsl_vector_char * b);
GSL_FUN int gsl_vector_char_scale (gsl_vector_char * a, const double x);
GSL_FUN int gsl_vector_char_add_constant (gsl_vector_char * a, const double x);

GSL_FUN int gsl_vector_char_isnull (const gsl_vector_char * v);
GSL_FUN int gsl_vector_char_ispos (const gsl_vector_char * v);
GSL_FUN int gsl_vector_char_isneg (const gsl_vector_char * v);
GSL_FUN int gsl_vector_char_isnonneg (const gsl_vector_char * v);

GSL_FUN INLINE_DECL char gsl_vector_char_get (const gsl_vector_char * v, const size_t i);
GSL_FUN INLINE_DECL void gsl_vector_char_set (gsl_vector_char * v, const size_t i, char x);
GSL_FUN INLINE_DECL char * gsl_vector_char_ptr (gsl_vector_char * v, const size_t i);
GSL_FUN INLINE_DECL const char * gsl_vector_char_const_ptr (const gsl_vector_char * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
char
gsl_vector_char_get (const gsl_vector_char * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

INLINE_FUN
void
gsl_vector_char_set (gsl_vector_char * v, const size_t i, char x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

INLINE_FUN
char *
gsl_vector_char_ptr (gsl_vector_char * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (char *) (v->data + i * v->stride);
}

INLINE_FUN
const char *
gsl_vector_char_const_ptr (const gsl_vector_char * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const char *) (v->data + i * v->stride);
}
#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_CHAR_H__ */

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs