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

#include <config.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

#define N 100000
void testMoments (double (*f) (void), const char *name,
		   double a, double b, double p);
void testPDF (double (*f) (void), double (*pdf)(double), const char *name);
void testDiscretePDF (double (*f) (void), double (*pdf)(unsigned int),
			const char *name);

void test_shuffle (void);
void test_choose (void);
double test_beta (void);
double test_beta_pdf (double x);
double test_bernoulli (void);
double test_bernoulli_pdf (unsigned int n);
double test_binomial (void);
double test_binomial_pdf (unsigned int n);
double test_binomial_large (void);
double test_binomial_large_pdf (unsigned int n);
double test_cauchy (void);
double test_cauchy_pdf (double x);
double test_chisq (void);
double test_chisq_pdf (double x);
double test_discrete1 (void);
double test_discrete1_pdf (unsigned int n);
double test_discrete2 (void);
double test_discrete2_pdf (unsigned int n);
double test_erlang (void);
double test_erlang_pdf (double x);
double test_exponential (void);
double test_exponential_pdf (double x);
double test_exppow0 (void);
double test_exppow0_pdf (double x);
double test_exppow1 (void);
double test_exppow1_pdf (double x);
double test_exppow1a (void);
double test_exppow1a_pdf (double x);
double test_exppow2 (void);
double test_exppow2_pdf (double x);
double test_exppow2a (void);
double test_exppow2a_pdf (double x);
double test_fdist (void);
double test_fdist_pdf (double x);
double test_flat (void);
double test_flat_pdf (double x);
double test_gamma (void);
double test_gamma_pdf (double x);
double test_gamma1 (void);
double test_gamma1_pdf (double x);
double test_gamma_int (void);
double test_gamma_int_pdf (double x);
double test_gamma_large (void);
double test_gamma_large_pdf (double x);
double test_gaussian (void);
double test_gaussian_pdf (double x);
double test_gaussian_ratio_method (void);
double test_gaussian_ratio_method_pdf (double x);
double test_gaussian_tail (void);
double test_gaussian_tail_pdf (double x);
double test_gaussian_tail1 (void);
double test_gaussian_tail1_pdf (double x);
double test_gaussian_tail2 (void);
double test_gaussian_tail2_pdf (double x);
double test_ugaussian (void);
double test_ugaussian_pdf (double x);
double test_ugaussian_ratio_method (void);
double test_ugaussian_ratio_method_pdf (double x);
double test_ugaussian_tail (void);
double test_ugaussian_tail_pdf (double x);
double test_bivariate_gaussian1 (void);
double test_bivariate_gaussian1_pdf (double x);
double test_bivariate_gaussian2 (void);
double test_bivariate_gaussian2_pdf (double x);
double test_bivariate_gaussian3 (void);
double test_bivariate_gaussian3_pdf (double x);
double test_bivariate_gaussian4 (void);
double test_bivariate_gaussian4_pdf (double x);
double test_gumbel1 (void);
double test_gumbel1_pdf (double x);
double test_gumbel2 (void);
double test_gumbel2_pdf (double x);
double test_geometric (void);
double test_geometric_pdf (unsigned int x);
double test_geometric1 (void);
double test_geometric1_pdf (unsigned int x);
double test_hypergeometric1 (void);
double test_hypergeometric1_pdf (unsigned int x);
double test_hypergeometric2 (void);
double test_hypergeometric2_pdf (unsigned int x);
double test_hypergeometric3 (void);
double test_hypergeometric3_pdf (unsigned int x);
double test_hypergeometric4 (void);
double test_hypergeometric4_pdf (unsigned int x);
double test_hypergeometric5 (void);
double test_hypergeometric5_pdf (unsigned int x);
double test_hypergeometric6 (void);
double test_hypergeometric6_pdf (unsigned int x);
double test_landau (void);
double test_landau_pdf (double x);
double test_levy1 (void);
double test_levy1_pdf (double x);
double test_levy2 (void);
double test_levy2_pdf (double x);
double test_levy1a (void);
double test_levy1a_pdf (double x);
double test_levy2a (void);
double test_levy2a_pdf (double x);
double test_levy_skew1 (void);
double test_levy_skew1_pdf (double x);
double test_levy_skew2 (void);
double test_levy_skew2_pdf (double x);
double test_levy_skew1a (void);
double test_levy_skew1a_pdf (double x);
double test_levy_skew2a (void);
double test_levy_skew2a_pdf (double x);
double test_levy_skew1b (void);
double test_levy_skew1b_pdf (double x);
double test_levy_skew2b (void);
double test_levy_skew2b_pdf (double x);
double test_logistic (void);
double test_logistic_pdf (double x);
double test_lognormal (void);
double test_lognormal_pdf (double x);
double test_logarithmic (void);
double test_logarithmic_pdf (unsigned int n);
double test_negative_binomial (void);
double test_negative_binomial_pdf (unsigned int n);
double test_pascal (void);
double test_pascal_pdf (unsigned int n);
double test_pareto (void);
double test_pareto_pdf (double x);
double test_poisson (void);
double test_poisson_pdf (unsigned int x);
double test_poisson_large (void);
double test_poisson_large_pdf (unsigned int x);
double test_dir2d (void);
double test_dir2d_pdf (double x);
double test_dir2d_trig_method (void);
double test_dir2d_trig_method_pdf (double x);
double test_dir3dxy (void);
double test_dir3dxy_pdf (double x);
double test_dir3dyz (void);
double test_dir3dyz_pdf (double x);
double test_dir3dzx (void);
double test_dir3dzx_pdf (double x);
double test_rayleigh (void);
double test_rayleigh_pdf (double x);
double test_rayleigh_tail (void);
double test_rayleigh_tail_pdf (double x);
double test_tdist1 (void);
double test_tdist1_pdf (double x);
double test_tdist2 (void);
double test_tdist2_pdf (double x);
double test_laplace (void);
double test_laplace_pdf (double x);
double test_weibull (void);
double test_weibull_pdf (double x);
double test_weibull1 (void);
double test_weibull1_pdf (double x);

gsl_rng *r_global;

int
main (void)
{
  gsl_ieee_env_setup ();

  gsl_rng_env_setup() ;
  r_global = gsl_rng_alloc (gsl_rng_default);

#define FUNC(x)  test_ ## x,                     "test gsl_ran_" #x
#define FUNC2(x) test_ ## x, test_ ## x ## _pdf, "test gsl_ran_" #x

  test_shuffle() ;
  test_choose() ;

  testMoments (FUNC (ugaussian), 0.0, 100.0, 0.5);
  testMoments (FUNC (ugaussian), -1.0, 1.0, 0.6826895);
  testMoments (FUNC (ugaussian), 3.0, 3.5, 0.0011172689);
  testMoments (FUNC (ugaussian_tail), 3.0, 3.5, 0.0011172689/0.0013498981);
  testMoments (FUNC (exponential), 0.0, 1.0, 1- exp(-0.5));
  testMoments (FUNC (cauchy), 0.0, 10000.0, 0.5);

  testMoments (FUNC (discrete1), -0.5, 0.5, 0.59 );
  testMoments (FUNC (discrete1), 0.5, 1.5, 0.40 );
  testMoments (FUNC (discrete1), 1.5, 3.5, 0.01 );

  testPDF (FUNC2(beta));
  testPDF (FUNC2(cauchy));
  testPDF (FUNC2(chisq));
  testPDF (FUNC2(erlang));
  testPDF (FUNC2(exponential));

  testPDF (FUNC2(exppow0));
  testPDF (FUNC2(exppow1));
  testPDF (FUNC2(exppow1a));
  testPDF (FUNC2(exppow2));
  testPDF (FUNC2(exppow2a));

  testPDF (FUNC2(fdist));
  testPDF (FUNC2(flat));
  testPDF (FUNC2(gamma));
  testPDF (FUNC2(gamma1));
  testPDF (FUNC2(gamma_int));
  testPDF (FUNC2(gamma_large));
  testPDF (FUNC2(gaussian));
  testPDF (FUNC2(gaussian_ratio_method));
  testPDF (FUNC2(ugaussian));
  testPDF (FUNC2(ugaussian_ratio_method));
  testPDF (FUNC2(gaussian_tail));
  testPDF (FUNC2(gaussian_tail1));
  testPDF (FUNC2(gaussian_tail2));
  testPDF (FUNC2(ugaussian_tail));

  testPDF (FUNC2(bivariate_gaussian1));
  testPDF (FUNC2(bivariate_gaussian2));
  testPDF (FUNC2(bivariate_gaussian3));
  testPDF (FUNC2(bivariate_gaussian4));

  testPDF (FUNC2(gumbel1));
  testPDF (FUNC2(gumbel2));
  testPDF (FUNC2(landau));
  testPDF (FUNC2(levy1));
  testPDF (FUNC2(levy2));
  testPDF (FUNC2(levy1a));
  testPDF (FUNC2(levy2a));
  testPDF (FUNC2(levy_skew1));
  testPDF (FUNC2(levy_skew2));
  testPDF (FUNC2(levy_skew1a));
  testPDF (FUNC2(levy_skew2a));
  testPDF (FUNC2(levy_skew1b));
  testPDF (FUNC2(levy_skew2b));
  testPDF (FUNC2(logistic));
  testPDF (FUNC2(lognormal));
  testPDF (FUNC2(pareto));
  testPDF (FUNC2(rayleigh));
  testPDF (FUNC2(rayleigh_tail));
  testPDF (FUNC2(tdist1));
  testPDF (FUNC2(tdist2));
  testPDF (FUNC2(laplace));
  testPDF (FUNC2(weibull));
  testPDF (FUNC2(weibull1));

  testPDF (FUNC2(dir2d));
  testPDF (FUNC2(dir2d_trig_method));
  testPDF (FUNC2(dir3dxy));
  testPDF (FUNC2(dir3dyz));
  testPDF (FUNC2(dir3dzx));

  testDiscretePDF (FUNC2(discrete1));
  testDiscretePDF (FUNC2(discrete2));
  testDiscretePDF (FUNC2(poisson));
  testDiscretePDF (FUNC2(poisson_large));
  testDiscretePDF (FUNC2(bernoulli));
  testDiscretePDF (FUNC2(binomial));
  testDiscretePDF (FUNC2(binomial_large));
  testDiscretePDF (FUNC2(geometric));
  testDiscretePDF (FUNC2(geometric1));
  testDiscretePDF (FUNC2(hypergeometric1));
  testDiscretePDF (FUNC2(hypergeometric2));
  testDiscretePDF (FUNC2(hypergeometric3));
  testDiscretePDF (FUNC2(hypergeometric4));
  testDiscretePDF (FUNC2(hypergeometric5));
  testDiscretePDF (FUNC2(hypergeometric6));
  testDiscretePDF (FUNC2(logarithmic));
  testDiscretePDF (FUNC2(negative_binomial));
  testDiscretePDF (FUNC2(pascal));

  exit (gsl_test_summary());
}

void
test_shuffle (void)
{
  double count[10][10] ;
  int x[10] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9} ;
  int i, j, status = 0;

  for (i = 0; i < 10; i++)
    {
      for (j = 0; j < 10; j++)
	{
	  count[i][j] = 0 ;
	}
    }

  for (i = 0 ; i < N; i++)
    {
      for (j = 0; j < 10; j++)
	x[j] = j ;

      gsl_ran_shuffle (r_global, x, 10, sizeof(int)) ;

      for (j = 0; j < 10; j++)
	count[x[j]][j] ++ ;
    }

  for (i = 0; i < 10; i++)
    {
      for (j = 0; j < 10; j++)
	{
	  double expected = N / 10.0 ;
	  double d = fabs(count[i][j] - expected);
	  double sigma = d / sqrt(expected) ;
	  if (sigma > 5 && d > 1)
	    {
	      status = 1 ;
	      gsl_test (status,
			"gsl_ran_shuffle %d,%d (%g observed vs %g expected)",
			i, j, count[i][j]/N, 0.1) ;
	    }
	}
    }

  gsl_test (status, "gsl_ran_shuffle on {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}") ;

}

void
test_choose (void)
{
  double count[10] ;
  int x[10] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9} ;
  int y[3] = {0, 1, 2} ;
  int i, j, status = 0;

  for (i = 0; i < 10; i++)
    {
      count[i] = 0 ;
    }

  for (i = 0 ; i < N; i++)
    {
      for (j = 0; j < 10; j++)
	x[j] = j ;

      gsl_ran_choose (r_global, y, 3, x, 10, sizeof(int)) ;

      for (j = 0; j < 3; j++)
	count[y[j]]++ ;
    }

  for (i = 0; i < 10; i++)
    {
      double expected = 3.0 * N / 10.0 ;
      double d = fabs(count[i] - expected);
      double sigma = d / sqrt(expected) ;
      if (sigma > 5 && d > 1)
	{
	  status = 1 ;
	  gsl_test (status,
		    "gsl_ran_choose %d (%g observed vs %g expected)",
		    i, count[i]/N, 0.1) ;
	}
    }

  gsl_test (status, "gsl_ran_choose (3) on {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}") ;

}

void
testMoments (double (*f) (void), const char *name,
	      double a, double b, double p)
{
  int i;
  double count = 0, expected, sigma;
  int status;

  for (i = 0; i < N; i++)
    {
      double r = f ();
      if (r < b && r > a)
	count++;
    }

  expected = p * N;
  sigma = fabs (count - expected) / sqrt (expected);

  status = (sigma > 3);

  gsl_test (status, "%s [%g,%g] (%g observed vs %g expected)",
	    name, a, b, count / N, p);
}

#define BINS 100

void
testPDF (double (*f) (void), double (*pdf)(double), const char *name)
{
  double count[BINS], p[BINS];
  double a = -5.0, b = +5.0 ;
  double dx = (b - a) / BINS ;
  int i,j,status = 0, status_i =0 ;

  for (i = 0; i < BINS; i++)
    count[i] = 0 ;

  for (i = 0; i < N; i++)
    {
      double r = f ();
      if (r < b && r > a)
	{
	  j =  (int)((r - a)/dx) ;
	  count[j]++;
	}
    }

  for (i = 0; i < BINS; i++)
    {
      /* Compute an approximation to the integral of p(x) from x to
         x+dx using Simpson's rule */

      double x = a + i * dx ;
#define STEPS 100
      double sum = 0 ;

      if (fabs(x) < 1e-10) /* hit the origin exactly */
	x = 0.0 ;

      for (j = 1; j < STEPS; j++)
	sum += pdf(x + j * dx / STEPS) ;

      p[i] =  0.5 * (pdf(x) + 2*sum + pdf(x + dx - 1e-7)) * dx / STEPS ;
    }

  for (i = 0; i < BINS; i++)
    {
      double x = a + i * dx ;
      double d = fabs(count[i] - N*p[i]) ;
      if (p[i] != 0)
	{
	  double s = d / sqrt(N*p[i]) ;
	  status_i = (s > 5) && (d > 1) ;
	}
      else
	{
	  status_i = (count[i] != 0) ;
	}
      status |= status_i ;
      if (status_i)
	gsl_test (status_i, "%s [%g,%g) (%g/%d=%g observed vs %g expected)",
		  name, x, x+dx, count[i],N,count[i]/N, p[i]) ;
    }

  if (status == 0)
    gsl_test (status, "%s, sampling against pdf over range [%g,%g) ",
	      name, a, b) ;
}

void
testDiscretePDF (double (*f) (void), double (*pdf)(unsigned int), const char *name)
{
  double count[BINS], p[BINS];
  unsigned int i ;
  int status = 0, status_i =0 ;

  for (i = 0; i < BINS; i++)
    count[i] = 0 ;

  for (i = 0; i < N; i++)
    {
      int r = (int)(f ());
      if (r>= 0 && r < BINS)
	count[r]++;
    }

  for (i = 0; i < BINS; i++)
    p[i] =  pdf(i) ;

  for (i = 0; i < BINS; i++)
    {
      double d = fabs(count[i] - N*p[i]) ;
      if (p[i] != 0)
	{
	  double s = d/sqrt(N*p[i]) ;
	  status_i = (s > 5) && (d > 1);
	}
      else
	{
	  status_i = (count[i] != 0) ;
	}
      status |= status_i ;
      if (status_i)
	gsl_test (status_i, "%s i=%d (%g observed vs %g expected)",
		  name, i, count[i]/N, p[i]) ;
    }

  if (status == 0)
    gsl_test (status, "%s, sampling against pdf over range [%d,%d) ",
	      name, 0, BINS) ;
}

double
test_beta (void)
{
  return gsl_ran_beta (r_global, 2.0, 3.0);
}

double
test_beta_pdf (double x)
{
  return gsl_ran_beta_pdf (x, 2.0, 3.0);
}

double
test_bernoulli (void)
{
  return gsl_ran_bernoulli (r_global, 0.3);
}

double
test_bernoulli_pdf (unsigned int n)
{
  return gsl_ran_bernoulli_pdf (n, 0.3);
}

double
test_binomial (void)
{
  return gsl_ran_binomial (r_global, 0.3, 5);
}

double
test_binomial_pdf (unsigned int n)
{
  return gsl_ran_binomial_pdf (n, 0.3, 5);
}

double
test_binomial_large (void)
{
  return gsl_ran_binomial (r_global, 0.3, 55);
}

double
test_binomial_large_pdf (unsigned int n)
{
  return gsl_ran_binomial_pdf (n, 0.3, 55);
}

double
test_cauchy (void)
{
  return gsl_ran_cauchy (r_global, 2.0);
}

double
test_cauchy_pdf (double x)
{
  return gsl_ran_cauchy_pdf (x, 2.0);
}

double
test_chisq (void)
{
  return gsl_ran_chisq (r_global, 13.0);
}

double
test_chisq_pdf (double x)
{
  return gsl_ran_chisq_pdf (x, 13.0);
}

double
test_dir2d (void)
{
  double x=0, y=0, theta;
  gsl_ran_dir_2d (r_global, &x, &y);
  theta = atan2(x,y);
  return theta;
}

double
test_dir2d_pdf (double x)
{
  if (x > -M_PI && x <= M_PI)
    {
      return 1 / (2 * M_PI) ;
    }
  else
    {
      return 0 ;
    }
}

double
test_dir2d_trig_method (void)
{
  double x=0, y=0, theta;
  gsl_ran_dir_2d_trig_method (r_global, &x, &y);
  theta = atan2(x,y);
  return theta;
}

double
test_dir2d_trig_method_pdf (double x)
{
  if (x > -M_PI && x <= M_PI)
    {
      return 1 / (2 * M_PI) ;
    }
  else
    {
      return 0 ;
    }
}

double
test_dir3dxy (void)
{
  double x=0, y=0, z=0, theta;
  gsl_ran_dir_3d (r_global, &x, &y, &z);
  theta = atan2(x,y);
  return theta;
}

double
test_dir3dxy_pdf (double x)
{
  if (x > -M_PI && x <= M_PI)
    {
      return 1 / (2 * M_PI) ;
    }
  else
    {
      return 0 ;
    }
}

double
test_dir3dyz (void)
{
  double x=0, y=0, z=0, theta;
  gsl_ran_dir_3d (r_global, &x, &y, &z);
  theta = atan2(y,z);
  return theta;
}

double
test_dir3dyz_pdf (double x)
{
  if (x > -M_PI && x <= M_PI)
    {
      return 1 / (2 * M_PI) ;
    }
  else
    {
      return 0 ;
    }
}

double
test_dir3dzx (void)
{
  double x=0, y=0, z=0, theta;
  gsl_ran_dir_3d (r_global, &x, &y, &z);
  theta = atan2(z,x);
  return theta;
}

double
test_dir3dzx_pdf (double x)
{
  if (x > -M_PI && x <= M_PI)
    {
      return 1 / (2 * M_PI) ;
    }
  else
    {
      return 0 ;
    }
}

static gsl_ran_discrete_t *g1 = NULL;
static gsl_ran_discrete_t *g2 = NULL;

double
test_discrete1 (void)
{
    static double P[3]={0.59, 0.4, 0.01};
    if (g1==NULL) {
        g1 = gsl_ran_discrete_preproc(3,P);
    }
    return gsl_ran_discrete(r_global,g1);
}
double test_discrete1_pdf (unsigned int n)
{
    return gsl_ran_discrete_pdf((size_t)n,g1);
}

double
test_discrete2 (void)
{
    static double P[10]={ 1, 9, 3, 4, 5, 8, 6, 7, 2, 0 };
    if (g2==NULL) {
        g2 = gsl_ran_discrete_preproc(10,P);
    }
    return gsl_ran_discrete(r_global,g2);
}
double test_discrete2_pdf (unsigned int n)
{
    return gsl_ran_discrete_pdf((size_t)n,g2);
}

double
test_erlang (void)
{
  return gsl_ran_erlang (r_global, 3.0, 4.0);
}

double
test_erlang_pdf (double x)
{
  return gsl_ran_erlang_pdf (x, 3.0, 4.0);
}

double
test_exponential (void)
{
  return gsl_ran_exponential (r_global, 2.0);
}

double
test_exponential_pdf (double x)
{
  return gsl_ran_exponential_pdf (x, 2.0);
}

double
test_exppow0 (void)
{
  return gsl_ran_exppow (r_global, 3.7, 0.3);
}

double
test_exppow0_pdf (double x)
{
  return gsl_ran_exppow_pdf (x, 3.7, 0.3);
}

double
test_exppow1 (void)
{
  return gsl_ran_exppow (r_global, 3.7, 1.0);
}

double
test_exppow1_pdf (double x)
{
  return gsl_ran_exppow_pdf (x, 3.7, 1.0);
}

double
test_exppow1a (void)
{
  return gsl_ran_exppow (r_global, 3.7, 1.9);
}

double
test_exppow1a_pdf (double x)
{
  return gsl_ran_exppow_pdf (x, 3.7, 1.9);
}

double
test_exppow2 (void)
{
  return gsl_ran_exppow (r_global, 3.7, 2.0);
}

double
test_exppow2_pdf (double x)
{
  return gsl_ran_exppow_pdf (x, 3.7, 2.0);
}

double
test_exppow2a (void)
{
  return gsl_ran_exppow (r_global, 3.7, 7.5);
}

double
test_exppow2a_pdf (double x)
{
  return gsl_ran_exppow_pdf (x, 3.7, 7.5);
}

double
test_fdist (void)
{
  return gsl_ran_fdist (r_global, 3.0, 4.0);
}

double
test_fdist_pdf (double x)
{
  return gsl_ran_fdist_pdf (x, 3.0, 4.0);
}

double
test_flat (void)
{
  return gsl_ran_flat (r_global, 3.0, 4.0);
}

double
test_flat_pdf (double x)
{
  return gsl_ran_flat_pdf (x, 3.0, 4.0);
}

double
test_gamma (void)
{
  return gsl_ran_gamma (r_global, 2.5, 2.17);
}

double
test_gamma_pdf (double x)
{
  return gsl_ran_gamma_pdf (x, 2.5, 2.17);
}

double
test_gamma1 (void)
{
  return gsl_ran_gamma (r_global, 1.0, 2.17);
}

double
test_gamma1_pdf (double x)
{
  return gsl_ran_gamma_pdf (x, 1.0, 2.17);
}

double
test_gamma_int (void)
{
  return gsl_ran_gamma (r_global, 10.0, 2.17);
}

double
test_gamma_int_pdf (double x)
{
  return gsl_ran_gamma_pdf (x, 10.0, 2.17);
}

double
test_gamma_large (void)
{
  return gsl_ran_gamma (r_global, 20.0, 2.17);
}

double
test_gamma_large_pdf (double x)
{
  return gsl_ran_gamma_pdf (x, 20.0, 2.17);
}

double
test_gaussian (void)
{
  return gsl_ran_gaussian (r_global, 3.0);
}

double
test_gaussian_pdf (double x)
{
  return gsl_ran_gaussian_pdf (x, 3.0);
}

double
test_gaussian_ratio_method (void)
{
  return gsl_ran_gaussian_ratio_method (r_global, 3.0);
}

double
test_gaussian_ratio_method_pdf (double x)
{
  return gsl_ran_gaussian_pdf (x, 3.0);
}

double
test_gaussian_tail (void)
{
  return gsl_ran_gaussian_tail (r_global, 1.7, 0.25);
}

double
test_gaussian_tail_pdf (double x)
{
  return gsl_ran_gaussian_tail_pdf (x, 1.7, 0.25) ;
}

double
test_gaussian_tail1 (void)
{
  return gsl_ran_gaussian_tail (r_global, -1.7, 5.0);
}

double
test_gaussian_tail1_pdf (double x)
{
  return gsl_ran_gaussian_tail_pdf (x, -1.7, 5.0) ;
}

double
test_gaussian_tail2 (void)
{
  return gsl_ran_gaussian_tail (r_global, 0.1, 2.0);
}

double
test_gaussian_tail2_pdf (double x)
{
  return gsl_ran_gaussian_tail_pdf (x, 0.1, 2.0) ;
}

double
test_ugaussian (void)
{
  return gsl_ran_ugaussian (r_global);
}

double
test_ugaussian_pdf (double x)
{
  return gsl_ran_ugaussian_pdf (x);
}

double
test_ugaussian_ratio_method (void)
{
  return gsl_ran_ugaussian_ratio_method (r_global);
}

double
test_ugaussian_ratio_method_pdf (double x)
{
  return gsl_ran_ugaussian_pdf (x);
}

double
test_ugaussian_tail (void)
{
  return gsl_ran_ugaussian_tail (r_global, 3.0);
}

double
test_ugaussian_tail_pdf (double x)
{
  return gsl_ran_ugaussian_tail_pdf (x, 3.0) ;
}

double
test_bivariate_gaussian1 (void)
{
  double x = 0, y = 0;
  gsl_ran_bivariate_gaussian (r_global, 3.0, 2.0, 0.3, &x, &y);
  return x ;
}

double
test_bivariate_gaussian1_pdf (double x)
{
  return gsl_ran_gaussian_pdf (x, 3.0);
}

double
test_bivariate_gaussian2 (void)
{
  double x = 0, y = 0;
  gsl_ran_bivariate_gaussian (r_global, 3.0, 2.0, 0.3, &x, &y);
  return y ;
}

double
test_bivariate_gaussian2_pdf (double y)
{
  int i, n = 10 ;
  double sum = 0 ;
  double a = -10, b = 10, dx = (b - a)/n ;
  for (i = 0; i < n ; i++)
    {
      double x = a + i * dx ;
      sum += gsl_ran_bivariate_gaussian_pdf (x, y, 3.0, 2.0, 0.3) * dx ;
    }
  return sum ;
}

double
test_bivariate_gaussian3 (void)
{
  double x = 0, y = 0;
  gsl_ran_bivariate_gaussian (r_global, 3.0, 2.0, 0.3, &x, &y);
  return x + y ;
}

double
test_bivariate_gaussian3_pdf (double x)
{
  double sx = 3.0, sy = 2.0, r = 0.3;
  double su = (sx+r*sy) ;
  double sv = sy*sqrt(1-r*r) ;
  double sigma = sqrt(su*su + sv*sv) ;

  return gsl_ran_gaussian_pdf (x, sigma);
}

double
test_bivariate_gaussian4 (void)
{
  double x = 0, y = 0;
  gsl_ran_bivariate_gaussian (r_global, 3.0, 2.0, -0.5, &x, &y);
  return x + y ;
}

double
test_bivariate_gaussian4_pdf (double x)
{
  double sx = 3.0, sy = 2.0, r = -0.5;
  double su = (sx+r*sy) ;
  double sv = sy*sqrt(1-r*r) ;
  double sigma = sqrt(su*su + sv*sv) ;

  return gsl_ran_gaussian_pdf (x, sigma);
}

double
test_geometric (void)
{
  return gsl_ran_geometric (r_global, 0.5);
}

double
test_geometric_pdf (unsigned int n)
{
  return gsl_ran_geometric_pdf (n, 0.5);
}

double
test_geometric1 (void)
{
  return gsl_ran_geometric (r_global, 1.0);
}

double
test_geometric1_pdf (unsigned int n)
{
  return gsl_ran_geometric_pdf (n, 1.0);
}

double
test_hypergeometric1 (void)
{
  return gsl_ran_hypergeometric (r_global, 5, 7, 4);
}

double
test_hypergeometric1_pdf (unsigned int n)
{
  return gsl_ran_hypergeometric_pdf (n, 5, 7, 4);
}

double
test_hypergeometric2 (void)
{
  return gsl_ran_hypergeometric (r_global, 5, 7, 11);
}

double
test_hypergeometric2_pdf (unsigned int n)
{
  return gsl_ran_hypergeometric_pdf (n, 5, 7, 11);
}

double
test_hypergeometric3 (void)
{
  return gsl_ran_hypergeometric (r_global, 5, 7, 1);
}

double
test_hypergeometric3_pdf (unsigned int n)
{
  return gsl_ran_hypergeometric_pdf (n, 5, 7, 1);
}

double
test_hypergeometric4 (void)
{
  return gsl_ran_hypergeometric (r_global, 5, 7, 20);
}

double
test_hypergeometric4_pdf (unsigned int n)
{
  return gsl_ran_hypergeometric_pdf (n, 5, 7, 20);
}

double
test_hypergeometric5 (void)
{
  return gsl_ran_hypergeometric (r_global, 2, 7, 5);
}

double
test_hypergeometric5_pdf (unsigned int n)
{
  return gsl_ran_hypergeometric_pdf (n, 2, 7, 5);
}

double
test_hypergeometric6 (void)
{
  return gsl_ran_hypergeometric (r_global, 2, 10, 3);
}

double
test_hypergeometric6_pdf (unsigned int n)
{
  return gsl_ran_hypergeometric_pdf (n, 2, 10, 3);
}

double
test_gumbel1 (void)
{
  return gsl_ran_gumbel1 (r_global, 3.12, 4.56);
}

double
test_gumbel1_pdf (double x)
{
  return gsl_ran_gumbel1_pdf (x, 3.12, 4.56);
}

double
test_gumbel2 (void)
{
  return gsl_ran_gumbel2 (r_global, 3.12, 4.56);
}

double
test_gumbel2_pdf (double x)
{
  return gsl_ran_gumbel2_pdf (x, 3.12, 4.56);
}

double
test_landau (void)
{
  return gsl_ran_landau (r_global);
}

double
test_landau_pdf (double x)
{
  return gsl_ran_landau_pdf (x);
}

double
test_levy1 (void)
{
  return gsl_ran_levy (r_global, 5.0, 1.0);
}

double
test_levy1_pdf (double x)
{
  return gsl_ran_cauchy_pdf (x, 5.0);
}

double
test_levy2 (void)
{
  return gsl_ran_levy (r_global, 5.0, 2.0);
}

double
test_levy2_pdf (double x)
{
  return gsl_ran_gaussian_pdf (x, sqrt(2.0) * 5.0 );
}

double
test_levy1a (void)
{
  return gsl_ran_levy (r_global, 5.0, 1.01);
}

double
test_levy1a_pdf (double x)
{
  return gsl_ran_cauchy_pdf (x, 5.0);
}

double
test_levy2a (void)
{
  return gsl_ran_levy (r_global, 5.0, 1.99);
}

double
test_levy2a_pdf (double x)
{
  return gsl_ran_gaussian_pdf (x, sqrt(2.0) * 5.0 );
}

double
test_levy_skew1 (void)
{
  return gsl_ran_levy_skew (r_global, 5.0, 1.0, 0.0);
}

double
test_levy_skew1_pdf (double x)
{
  return gsl_ran_cauchy_pdf (x, 5.0);
}

double
test_levy_skew2 (void)
{
  return gsl_ran_levy_skew (r_global, 5.0, 2.0, 0.0);
}

double
test_levy_skew2_pdf (double x)
{
  return gsl_ran_gaussian_pdf (x, sqrt(2.0) * 5.0 );
}

double
test_levy_skew1a (void)
{
  return gsl_ran_levy_skew (r_global, 5.0, 1.01, 0.0);
}

double
test_levy_skew1a_pdf (double x)
{
  return gsl_ran_cauchy_pdf (x, 5.0);
}

double
test_levy_skew2a (void)
{
  return gsl_ran_levy_skew (r_global, 5.0, 1.99, 0.0);
}

double
test_levy_skew2a_pdf (double x)
{
  return gsl_ran_gaussian_pdf (x, sqrt(2.0) * 5.0 );
}

double
test_levy_skew1b (void)
{
  return gsl_ran_levy_skew (r_global, 5.0, 1.01, 0.001);
}

double
test_levy_skew1b_pdf (double x)
{
  return gsl_ran_cauchy_pdf (x, 5.0);
}

double
test_levy_skew2b (void)
{
  return gsl_ran_levy_skew (r_global, 5.0, 1.99, 0.001);
}

double
test_levy_skew2b_pdf (double x)
{
  return gsl_ran_gaussian_pdf (x, sqrt(2.0) * 5.0 );
}

double
test_logistic (void)
{
  return gsl_ran_logistic (r_global, 3.1);
}

double
test_logistic_pdf (double x)
{
  return gsl_ran_logistic_pdf (x, 3.1);
}

double
test_logarithmic (void)
{
  return gsl_ran_logarithmic (r_global, 0.4);
}

double
test_logarithmic_pdf (unsigned int n)
{
  return gsl_ran_logarithmic_pdf (n, 0.4);
}

double
test_lognormal (void)
{
  return gsl_ran_lognormal (r_global, 2.7, 1.3);
}

double
test_lognormal_pdf (double x)
{
  return gsl_ran_lognormal_pdf (x, 2.7, 1.3);
}

double
test_negative_binomial (void)
{
  return gsl_ran_negative_binomial (r_global, 0.3, 20.0);
}

double
test_negative_binomial_pdf (unsigned int n)
{
  return gsl_ran_negative_binomial_pdf (n, 0.3, 20.0);
}

double
test_pascal (void)
{
  return gsl_ran_pascal (r_global, 0.8, 3);
}

double
test_pascal_pdf (unsigned int n)
{
  return gsl_ran_pascal_pdf (n, 0.8, 3);
}

double
test_pareto (void)
{
  return gsl_ran_pareto (r_global, 1.9, 2.75);
}

double
test_pareto_pdf (double x)
{
  return gsl_ran_pareto_pdf (x, 1.9, 2.75);
}

double
test_rayleigh (void)
{
  return gsl_ran_rayleigh (r_global, 1.9);
}

double
test_rayleigh_pdf (double x)
{
  return gsl_ran_rayleigh_pdf (x, 1.9);
}

double
test_rayleigh_tail (void)
{
  return gsl_ran_rayleigh_tail (r_global, 2.7, 1.9);
}

double
test_rayleigh_tail_pdf (double x)
{
  return gsl_ran_rayleigh_tail_pdf (x, 2.7, 1.9);
}

double
test_poisson (void)
{
  return gsl_ran_poisson (r_global, 5.0);
}

double
test_poisson_pdf (unsigned int n)
{
  return gsl_ran_poisson_pdf (n, 5.0);
}

double
test_poisson_large (void)
{
  return gsl_ran_poisson (r_global, 30.0);
}

double
test_poisson_large_pdf (unsigned int n)
{
  return gsl_ran_poisson_pdf (n, 30.0);
}

double
test_tdist1 (void)
{
  return gsl_ran_tdist (r_global, 1.75);
}

double
test_tdist1_pdf (double x)
{
  return gsl_ran_tdist_pdf (x, 1.75);
}

double
test_tdist2 (void)
{
  return gsl_ran_tdist (r_global, 12.75);
}

double
test_tdist2_pdf (double x)
{
  return gsl_ran_tdist_pdf (x, 12.75);
}

double
test_laplace (void)
{
  return gsl_ran_laplace (r_global, 2.75);
}

double
test_laplace_pdf (double x)
{
  return gsl_ran_laplace_pdf (x, 2.75);
}

double
test_weibull (void)
{
  return gsl_ran_weibull (r_global, 3.14, 2.75);
}

double
test_weibull_pdf (double x)
{
  return gsl_ran_weibull_pdf (x, 3.14, 2.75);
}

double
test_weibull1 (void)
{
  return gsl_ran_weibull (r_global, 2.97, 1.0);
}

double
test_weibull1_pdf (double x)
{
  return gsl_ran_weibull_pdf (x, 2.97, 1.0);
}


text:
#pragma once
#include "../utils/utils.h"
#include <gsl/gsl>

namespace Halley {
	class Compression {
	public:
		static Bytes deflate(const Bytes& bytes);
		static Bytes deflate(gsl::span<const gsl::byte> bytes);
		static Bytes inflate(const Bytes& bytes);
		static Bytes inflate(gsl::span<const gsl::byte> bytes);
		static std::shared_ptr<const char> inflateToSharedPtr(gsl::span<const gsl::byte> bytes, size_t& outSize);
		static unsigned char* inflateRaw(gsl::span<const gsl::byte> bytes, size_t& outSize);
	};
}

ext: h

max_stars_repo_name: Healthire/halley

max_issues_repo_name: Healthire/halley

max_forks_repo_name: Healthire/halley


text:
/**
 * @file blas_mangle.h
 * @brief Macros to "de-mangle" BLAS routines and includes and allow generic calls to them in the code.
 * @author Moritz Kreutzer <moritz.kreutzer@fau.de>
 */
#ifndef GHOST_BLAS_MANGLE_H
#define GHOST_BLAS_MANGLE_H

#include <strings.h>

#ifdef GHOST_HAVE_MKL
#include <mkl_cblas.h>
#elif defined(GHOST_HAVE_GSL)
#include <gsl_cblas.h>
#else
#include <cblas.h>
#endif

#define blas_order(order) order==GHOST_DENSEMAT_COLMAJOR?CblasColMajor:CblasRowMajor
#define blas_trans(trans) trans[0]=='N'?CblasNoTrans:(trans[0]=='C'?CblasConjTrans:CblasTrans)

#define sgemm(order,transa,transb,m,n,k,alpha,a,lda,b,ldb,beta,c,ldc) \
    cblas_sgemm(blas_order(order),blas_trans(transa),blas_trans(transb),*m,*n,*k,*alpha,a,*lda,b,*ldb,*beta,c,*ldc)
#define dgemm(order,transa,transb,m,n,k,alpha,a,lda,b,ldb,beta,c,ldc) \
    cblas_dgemm(blas_order(order),blas_trans(transa),blas_trans(transb),*m,*n,*k,*alpha,a,*lda,b,*ldb,*beta,c,*ldc)
#define cgemm(order,transa,transb,m,n,k,alpha,a,lda,b,ldb,beta,c,ldc) \
    cblas_cgemm(blas_order(order),blas_trans(transa),blas_trans(transb),*m,*n,*k,alpha,a,*lda,b,*ldb,beta,c,*ldc)
#define zgemm(order,transa,transb,m,n,k,alpha,a,lda,b,ldb,beta,c,ldc) \
    cblas_zgemm(blas_order(order),blas_trans(transa),blas_trans(transb),*m,*n,*k,alpha,a,*lda,b,*ldb,beta,c,*ldc)

#endif

max_stars_repo_name: RRZE-HPC/GHOST

max_issues_repo_name: RRZE-HPC/GHOST

max_forks_repo_name: RRZE-HPC/GHOST


text:
#ifndef BFM_FRANCIS_H
#define BFM_FRANCIS_H

#include <cstdlib>
#include <string>
#include <cmath>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <fstream>
#include <complex>
#include <algorithm>
#include "UTSolve.h"
#include "Householder.h"
#include <util/time_cps.h>

//#define USE_LAPACK
//#define USE_EIGEN

#ifdef USE_LAPACK
#warning "Using LAPACK in Lanczos"
#define lapack_complex_float std::complex<float>
#define lapack_complex_double std::complex<double>
#include <lapacke.h>
#endif
#ifdef USE_EIGEN
#warning "Using EIGEN package in Lanczos"
#include <Eigen/Dense>
#endif

namespace BFM_Krylov{

  // template <class T> int SymmEigensystem(Matrix<T > &Ain, std::vector<T> &evals, std::vector<std::vector<T> > &evecs, double small);
  // template <class T> int Eigensystem(Matrix<T > &Ain, std::vector<T> &evals, std::vector<std::vector<T> > &evecs, double small);

  /**

     Find the eigenvalues of an upper hessenberg matrix using the Francis QR algorithm.

     H =
     x  x  x  x  x  x  x  x  x
     x  x  x  x  x  x  x  x  x
     0  x  x  x  x  x  x  x  x
     0  0  x  x  x  x  x  x  x
     0  0  0  x  x  x  x  x  x
     0  0  0  0  x  x  x  x  x
     0  0  0  0  0  x  x  x  x
     0  0  0  0  0  0  x  x  x
     0  0  0  0  0  0  0  x  x

     Factorization is P T P^H where T is upper triangular (mod cc blocks) and P is orthagonal/unitary.

  **/
  template <class T> int QReigensystem(Matrix<T> &Hin, std::vector<T> &evals, std::vector<std::vector<T> > &evecs, double small){
    Matrix<T > H(Hin.dim); H = Hin; ///I don't want to modify the input but matricies must be passed by reference
    for(int i=0;i<evecs.size();i++){evals[i] = 0;for(int j=0;j<evecs[i].size();j++){evecs[i][j] = 0;}}

    int N = H.dim;
    int M = N;
    T s,t,x=0,y=0,z=0;
    T u,d;
    T apd,amd,bc;
    std::vector<T > p(N,0);
    T nrm = H.Norm();		///Matrix Norm
    int n, m;
    int e = 0;
    int it = 0;
    int tot_it = 0;
    int l = 0;
    int r = 0;
    Matrix<T > P(N); P.Unity();
    std::vector<int> trows(N,0);

    /// Check if the matrix is really hessenberg, if not abort
    double sth = 0;
    for(int j=0;j<N;j++){for(int i=j+2;i<N;i++){sth = abs(H(i,j));
	if(sth > small){
	  std::cout << "Non hessenberg H = " << sth << " > " << small << std::endl;
	  exit(1);
	}
      }}

    /** Check for convergence

	x  x  x  x  x
	0  x  x  x  x
	0  0  x  x  x
	0  0  x  x  x
	0  0  0  0  x

	for this matrix l = 4
    **/
    do{
      do{

	l = H.Chop_subdiag(nrm,e,small);

	r = 0;		///May have converged on more than one eval

	///Single eval
	if(l == N-1){
	  evals[e] = H(l,l);
	  N--; e++; r++; it = 0;
	}
	///Double eval
	if(l == N-2){
	  trows[l+1] = 1;		///Needed for UTSolve
	  apd = H(l,l) + H(l+1,l+1);
	  amd = H(l,l) - H(l+1,l+1);
	  bc =  (T)4.0*H(l+1,l)*H(l,l+1);

	  evals[e] = (T)0.5*( apd + sqrt(amd*amd + bc) );
	  evals[e+1] = (T)0.5*( apd - sqrt(amd*amd + bc) );
	  N-=2; e+=2; r++; it = 0;
	}
      }while(r>0);

      if(N ==0){break;}

      std::vector<T > ck(3), v(3);
      for(int m = N-3; m >= l; m--){
	///Starting vector essentially random shift.
	if(it%10 == 0 && N >= 3 && it > 0){
	  s = (T)1.618033989*( abs( H(N-1,N-2) ) + abs( H(N-2,N-3) ) );
	  t = (T)0.618033989*( abs( H(N-1,N-2) ) + abs( H(N-2,N-3) ) );
	  x = H(m,m)*H(m,m) + H(m,m+1)*H(m+1,m) - s*H(m,m) + t;
	  y = H(m+1,m)*(H(m,m) + H(m+1,m+1) - s);
	  z = H(m+1,m)*H(m+2,m+1);
	}
	///Starting vector implicit Q theorem
	else{
	  s = (H(N-2,N-2) + H(N-1,N-1));
	  t = (H(N-2,N-2)*H(N-1,N-1) - H(N-2,N-1)*H(N-1,N-2));
	  x = H(m,m)*H(m,m) + H(m,m+1)*H(m+1,m) - s*H(m,m) + t;
	  y = H(m+1,m)*(H(m,m) + H(m+1,m+1) - s);
	  z = H(m+1,m)*H(m+2,m+1);
	}

	ck[0] = x; ck[1] = y; ck[2] = z;
	if(m == l) break;

	/** Some stupid thing from numerical recipies, seems to work**/
	u=abs(H(m,m-1))*(abs(y)+abs(z));
	d=abs(x)*(abs(H(m-1,m-1))+abs(H(m,m))+abs(H(m+1,m+1)));
	if ((T)abs(u+d) == (T)abs(d) ){l = m; break;}
	//if (u < small){l = m; break;}
      }

      if(it > 100000){
	std::cout << "QReigensystem: bugger it got stuck after 100000 iterations" << std::endl;
	std::cout << "got " << e << " evals " << l << " " << N << std::endl;
	exit(1);
      }

      normalize(ck);		///Normalization cancels in PHP anyway
      T beta;

      Householder_vector<T >(ck, 0, 2, v, beta);
      Householder_mult<T >(H,v,beta,0,l,l+2,0);
      Householder_mult<T >(H,v,beta,0,l,l+2,1);
      ///Accumulate eigenvector
      Householder_mult<T >(P,v,beta,0,l,l+2,1);

      int sw = 0;			///Are we on the last row?

      for(int k=l;k<N-2;k++){
	x = H(k+1,k);
	y = H(k+2,k);
	z = (T)0.0;
	if(k+3 <= N-1){z = H(k+3,k);}
	else{sw = 1; v[2] = (T)0.0;}

	ck[0] = x; ck[1] = y; ck[2] = z;

	normalize(ck);

	Householder_vector<T >(ck, 0, 2-sw, v, beta);
	Householder_mult<T >(H,v, beta,0,k+1,k+3-sw,0);
	Householder_mult<T >(H,v, beta,0,k+1,k+3-sw,1);
	///Accumulate eigenvector
	Householder_mult<T >(P,v, beta,0,k+1,k+3-sw,1);
      }

      it++;
      tot_it++;
    }while(N > 1);

    N = evals.size();

    ///Annoying - UT solves in reverse order;
    std::vector<T> tmp(N);for(int i=0;i<N;i++){tmp[i] = evals[N-i-1];} evals = tmp;
    UTeigenvectors(H, trows, evals, evecs);

    for(int i=0;i<evals.size();i++){evecs[i] = P*evecs[i]; normalize(evecs[i]);}

    return tot_it;
  }

  /**
     turn a matrix A =
     x  x  x  x  x
     x  x  x  x  x
     x  x  x  x  x
     x  x  x  x  x
     x  x  x  x  x
     into
     x  x  x  x  x
     x  x  x  x  x
     0  x  x  x  x
     0  0  x  x  x
     0  0  0  x  x
     with householder rotations
     Slow.
  */

  template <class T> void Hess(Matrix<T > &A, Matrix<T> &Q, int start){
    int N = A.dim;	//Matrix Size

    std::vector<T > p(N,0);
    for(int k=start;k<N-2;k++){
      std::vector<T > ck(N-k-1), v(N-k-1);
      for(int i=k+1;i<N;i++){ck[i-k-1] = A(i,k);}	///kth column

      normalize(ck);		///Normalization cancels in PHP anyway
      T beta;
      Householder_vector<T >(ck, 0, ck.size()-1, v, beta);	///Householder vector

      Householder_mult<T>(A,v,beta,start,k+1,N-1,0);	///A -> PA
      Householder_mult<T >(A,v,beta,start,k+1,N-1,1);  ///PA -> PAP^H
      ///Accumulate eigenvector
      Householder_mult<T >(Q,v,beta,start,k+1,N-1,1);  ///Q -> QP^H
    }
  }

  ///Tridiagonalize a matrix
  template <class T> void Tri(Matrix<T > &A, Matrix<T> &Q, int start){
    int N = A.dim;	//Matrix Size
    Hess(A,Q,start);
  }

  ///Tridiagonalize a matrix
  template <class T>
  void ForceTridiagonal(Matrix<T > &A){
    int N = A.dim;	//Matrix Size
    for(int l=0;l<N-2;l++){
      for(int k=l+2;k<N;k++){
	A(0,l,k);
	A(0,k,l);
      }
    }
  }

  /**
      Find the eigenvalues of an upper Hessenberg matrix using the Wilkinson QR algorithm.

      H =
      x  x  0  0  0  0
      x  x  x  0  0  0
      0  x  x  x  0  0
      0  0  x  x  x  0
      0  0  0  x  x  x
      0  0  0  0  x  x

      Factorization is P T P^H where T is upper triangular (mod cc blocks) and P is orthagonal/unitary.  **/

  //Rudy's original version with modifications passed down the ages.....

  template <class T>
  int basic_Wilkinson(Matrix<T> &Hin, std::vector<T> &evals, std::vector<std::vector<T> > &evecs, double small)
  {
    return basic_Wilkinson(Hin, evals, evecs, small, small);
  }

  template <class T>
  int basic_Wilkinson(Matrix<T> &Hin, std::vector<T> &evals, std::vector<std::vector<T> > &evecs, double small, double tol)
  {
    Matrix<T> H(Hin.dim);
    H = Hin; ///I don't want to modify the input but matricies must be passed by reference (CK: I don't think the author really understands why we usually pass by reference!)

    //Scale a matrix by its "norm"
    //double Hnorm = abs( Hin.LargestDiag() ); H =  H*(1.0/Hnorm);
    double Hnorm = abs(Hin.Norm());
    H = H * (1.0 / Hnorm);

    for(int i = 0; i < evecs.size(); ++i){
      evals[i] = 0;
      for(int j = 0; j < evecs[i].size(); ++j)
	evecs[i][j] = 0;
    }

    int N = H.dim;
    int M = N;
    T s, t, x = 0, y = 0, z = 0;
    T u, d;
    T apd, amd, bc;
    std::vector<T> p(N, 0);
    T nrm = H.Norm();		///Matrix Norm
    int n, m;
    int e = 0;
    int it = 0;
    int tot_it = 0;
    int l = 0;
    int r = 0;
    Matrix<T> P(N);
    P.Unity();
    std::vector<int> trows(N, 0);

    /// Check if the matrix is really symm tridiag
    double sth = 0;
    for(int j = 0; j < N; ++j){
	for(int i = j + 2; i < N; ++i){
	  if(abs(H(i, j)) > tol || abs(H(j, i)) > tol){
	    QDPIO::cout << "Non Tridiagonal H(" << i << ","<< j << ") = |" << Real( real( H(j,i) ) ) << "| > " << tol << std::endl;
	    QDPIO::cout << "Warning tridiagonalize and call again" << std::endl;
	    exit(1);
	  }
	}
    }
    do{
      do{
	//Jasper
	//Check if the subdiagonal term is small enough (<small)
	//if true then it is converged.
	//check start from H.dim - e - 1
	//How to deal with more than 2 are converged?
	//What if Chop_symm_subdiag return something int the middle?
	//--------------
	l = H.Chop_symm_subdiag(nrm, e, small);
	r = 0;		///May have converged on more than one eval

	//Jasper
	//In this case
	// x  x  0  0  0  0
	// x  x  x  0  0  0
	// 0  x  x  x  0  0
	// 0  0  x  x  x  0
	// 0  0  0  x  x  0
	// 0  0  0  0  0  x  <- l
	//--------------
	///Single eval
	if(l == N - 1){
	  evals[e] = H(l, l);
	  N--;
	  e++;
	  r++;
	  it = 0;
	}
	//Jasper
	// x  x  0  0  0  0
	// x  x  x  0  0  0
	// 0  x  x  x  0  0
	// 0  0  x  x  0  0
	// 0  0  0  0  x  x  <- l
	// 0  0  0  0  x  x
	//--------------
	///Double eval
	if(l == N - 2){
	  trows[l + 1] = 1;		///Needed for UTSolve
	  apd = H(l, l) + H(l + 1, l + 1);
	  amd = H(l, l) - H(l + 1, l + 1);
	  bc =  (T) 4.0 * H(l + 1, l) * H(l, l + 1);

	  evals[e] = (T) 0.5 * (apd + sqrt(amd * amd + bc));
	  evals[e + 1] = (T) 0.5 * (apd - sqrt(amd * amd + bc));
	  N -= 2;
	  e += 2;
	  r++;
	  it = 0;
	}
      }while(r > 0);

      //Jasper
      //Already converged
      //--------------
      if(N == 0)
	break;

      std::vector<T> ck(2), v(2);
      for(int m = N - 3; m >= l; m--){
	///Starting vector essentially random shift.
	if(it%10 == 0 && N >= 3 && it > 0){
	  t = abs(H(N - 1, N - 2)) + abs(H(N - 2, N - 3));
	  x = H(m, m) - t;
	  z = H(m + 1, m);
	}
	///Starting vector implicit Q theorem
	else{
	  d = (H(N - 2, N - 2) - H(N - 1, N - 1)) * (T) 0.5;
	  t = H(N - 1, N - 1) - H(N - 1, N - 2) * H(N - 1, N - 2) / (d + sign(d) * sqrt(d * d + H(N - 1, N - 2) * H(N - 1, N - 2)));
	  x = H(m, m) - t;
	  z = H(m + 1, m);
	}
	//Jasper
	//why it is here????
	//-----------------------
	if(m == l)
	  break;

	u = abs(H(m, m - 1)) * (abs(y) + abs(z));
	d = abs(x) * (abs(H(m - 1, m - 1)) + abs(H(m, m)) + abs(H(m + 1, m + 1)));
	if ((T)abs(u + d) == (T)abs(d)){
	  l = m;
	  break;
	}
      }
      //Jasper
      if(it > 1000000){
	QDPIO::cout << "Wilkinson: bugger it got stuck after 100000 iterations" << std::endl;
	QDPIO::cout << "got " << e << " evals " << l << " " << N << std::endl;
	exit(1);
      }

      T s, c;
      Givens_calc<T>(x, z, c, s);
      Givens_mult<T>(H, l, l + 1, c, -s, 0);
      Givens_mult<T>(H, l, l + 1, c,  s, 1);
      Givens_mult<T>(P, l, l + 1, c,  s, 1);

      for(int k = l; k < N - 2; ++k){
	x = H(k + 1,k);
	z = H(k + 2,k);

	Givens_calc<T>(x, z, c, s);
	Givens_mult<T>(H, k + 1, k + 2, c, -s, 0);
	Givens_mult<T>(H, k + 1, k + 2, c,  s, 1);
	Givens_mult<T>(P, k + 1, k + 2, c,  s, 1);
      }
      it++;
      tot_it++;
    }while(N > 1);

    N = evals.size();

    ///Annoying - UT solves in reverse order;
    std::vector<T> tmp(N);
    for(int i = 0; i < N; ++i)
      tmp[i] = evals[N-i-1];

    evals = tmp;

    UTeigenvectors(H, trows, evals, evecs);

    for(int i = 0; i < evals.size(); ++i){
      evecs[i] = P * evecs[i];
      normalize(evecs[i]);
      evals[i] = evals[i] * Hnorm;
    }

    return tot_it;
  }

  //An implementation using the Eigen package courtesy of Luchang
#ifdef USE_EIGEN
  template <class T>
  int eigen_Wilkinson(Matrix<T> &Ain, std::vector<T> &evals, std::vector<std::vector<T> > &evecs, double small){
    using namespace Eigen;
    const int size = Ain.dim;
    evals.resize(size);
    evecs.resize(size);
    for (int i = 0; i < size; i++) {
      evecs[i].resize(size);
    }
    MatrixXd A(size, size);
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
	A(i, j) = 0;
      }
    }
    for (int i = 0; i < size; i++) {
      for (int j = std::max(0, i-1); j <= std::min(size-1, i+1); j++) {
	A(i, j) = Ain(i, j);
      }
    }
    SelfAdjointEigenSolver<MatrixXd> eigensolver(A);
    for (int i = 0; i < size; i++) {
      evals[i] = eigensolver.eigenvalues()(i);
    }
    for (int i = 0; i < size; i++) {
      std::vector<T>& vec = evecs[i];
      for (int j = 0; j < size; j++) {
	vec[j] = eigensolver.eigenvectors()(j, i);
      }
    }
    return 0;
  }
#endif

 //An implementation using the LAPACK package courtesy of Luchang
#ifdef USE_LAPACK
  template <class T>
  int lapackcj_Wilkinson(Matrix<T> &AH, vector<T> &tevals, vector<vector<T> > &tevecs, double small) {
    double time = -cps::dclock();
    const int size = AH.dim;
    tevals.resize(size);
    tevecs.resize(size);
    int NN = tevals.size();
    for(int i=0;i<NN;i++) tevecs[i].resize(NN);
    double evals_tmp[NN];
    double evec_tmp[NN][NN];
    memset(evec_tmp[0],0,sizeof(double)*NN*NN);
    double AA[NN][NN];
    //    double ZZ[NN][NN];
    double DD[NN];
    double EE[NN];
    memset(AA, 0, sizeof(double)*NN*NN);
    //    memset(ZZ, 0, sizeof(double)*NN*NN);
    for (int i = 0; i< NN; i++)
      for (int j = i - 1; j <= i + 1; j++)
	if ( j < NN && j >= 0 ) {
	  AA[i][j] = AH(i,j);
	  if (i==j) DD[i] = AA[i][j];
	  if (i==j) evals_tmp[i] = AA[i][j];
	  if (j==(i-1)) EE[j] = AA[i][j];
	  //        if (i<20 && j<20)
	  QDPIO:: cout << "AA["<<i<<"]["<<j<<"]="<<AA[i][j]<<endl;
	}
    int evals_found;
    int lwork = ( (18*NN) > (1+4*NN+NN*NN)? (18*NN):(1+4*NN+NN*NN)) ;
    int liwork =  3+NN*10 ;
    int iwork[liwork];
    double work[lwork];
    int isuppz[2*NN];
    char jobz = 'V'; // calculate evals & evecs
    char range = 'I'; // calculate all evals
    //    char range = 'A'; // calculate all evals
    char uplo = 'U'; // refer to upper half of original matrix
    char compz = 'I'; // Compute eigenvectors of tridiagonal matrix
    int ifail[NN];
    int info;
    // dsyevx(&jobz, &range, &uplo, NN,
    //     AA, NN,
    //     0, 0, 0, 0, // these four are ignored if second parameteris 'A'
    //     0, // tolerance
    //     &evals_found, evals_tmp, evec_tmp, NN,
    //     work, lwork, iwork,
    //     ifail, info);
    int total = QMP_get_number_of_nodes();
    int node = QMP_get_node_number();
    int interval = (NN/total)+1;
    double vl = 0.0, vu = 0.0;
    int il = interval*node+1 , iu = interval*(node+1);
    if (iu > NN)  iu=NN;
    double tol = 0.0;
    if (0)
      {
	memset(evals_tmp,0,sizeof(double)*NN);
	if ( il <= NN){
	  printf("total=%d node=%d il=%d iu=%d\n",total,node,il,iu);
	  LAPACK_dsyevx(&jobz, &range, &uplo, &NN,
			(double*)AA, &NN,
			&vl, &vu, &il, &iu, // these four are ignored if second parameteris 'A'
			&tol, // tolerance
			&evals_found, evals_tmp, (double*)evec_tmp, &NN,
			work, &lwork, iwork,
			ifail, &info);
	  for (int i = iu-1; i>= il-1; i--){
	    printf("node=%d evals_found=%d evals_tmp[%d] = %g\n",node,evals_found, i - (il-1),evals_tmp[i - (il-1)]);
	    evals_tmp[i] = evals_tmp[i - (il-1)];
	    if (il>1) evals_tmp[i-(il-1)]=0.;
	    for (int j = 0; j< NN; j++){
	      evec_tmp[i][j] = evec_tmp[i - (il-1)][j];
	      if (il>1) evec_tmp[i-(il-1)][j]=0.;
	    }
	  }
	}
	{
	  QMP_sum_double_array(evals_tmp,NN);
	  QMP_sum_double_array((double *)evec_tmp,NN*NN);
	}
      } else
      if (1) {
	memset(evals_tmp,0,sizeof(double)*NN);
	if ( il <= NN){
	  printf("total=%d node=%d il=%d iu=%d\n",total,node,il,iu);
	  LAPACK_dstegr(&jobz, &range, &NN,
			(double*)DD, (double*)EE,
			&vl, &vu, &il, &iu, // these four are ignored if second parameteris 'A'
			&tol, // tolerance
			&evals_found, evals_tmp, (double*)evec_tmp, &NN,
			isuppz,
			work, &lwork, iwork, &liwork,
			&info);
	  for (int i = iu-1; i>= il-1; i--){
	    printf("node=%d evals_found=%d evals_tmp[%d] = %g\n",node,evals_found, i - (il-1),evals_tmp[i - (il-1)]);
	    evals_tmp[i] = evals_tmp[i - (il-1)];
	    if (il>1) evals_tmp[i-(il-1)]=0.;
	    for (int j = 0; j< NN; j++){
	      evec_tmp[i][j] = evec_tmp[i - (il-1)][j];
	      if (il>1) evec_tmp[i-(il-1)][j]=0.;
	    }
	  }
	}
	{
	  QMP_sum_double_array(evals_tmp,NN);
	  QMP_sum_double_array((double *)evec_tmp,NN*NN);
	}
      } else
	if(0) {
	  LAPACK_dsteqr(&compz, &NN,
			(double*)evals_tmp, (double*)EE,
			(double*)evec_tmp, &NN,
			work, &info);
	} else
	  {
	    LAPACK_dstedc(&compz, &NN,
			  (double*)evals_tmp, (double*)EE,
			  (double*)evec_tmp, &NN,
			  work, &lwork,
			  iwork, &liwork,
			  &info);
	  }
    for(int i=0;i<NN;i++){
      for(int j=0;j<NN;j++)
	tevecs[i][j]=evec_tmp[i][j];
      tevals[i]=evals_tmp[i];
    }

    cps::print_time("lapackcj_Wilkinson","run time",time+cps::dclock());
  }
#endif

  template <class T>
  int Wilkinson(Matrix<T> &Ain, std::vector<T> &evals, std::vector<std::vector<T> > &evecs, double small){
#ifdef USE_LAPACK
    return lapackcj_Wilkinson(Ain, evals, evecs, small);
#elif USE_EIGEN
    return eigen_Wilkinson(Ain, evals, evecs, small);
#else
    return basic_Wilkinson(Ain, evals, evecs, small);
#endif
  }

  ///Solve a symmetric eigensystem, not necessarily in tridiagonal form
  //Rudy's original implementation
  template <class T> int basic_SymmEigensystem(Matrix<T > &Ain, std::vector<T> &evals, std::vector<std::vector<T> > &evecs, double small){
    int N = Ain.dim;
    Matrix<T > A(N); A = Ain;
    Matrix<T > Q(N);Q.Unity();
    Tri(A,Q,0);
    int it = basic_Wilkinson<T>(A, evals, evecs, small);
    for(int k=0;k<N;k++){evecs[k] = Q*evecs[k];}
    return it;
  }

  //I guess Luchang wrote this....
#ifdef USE_EIGEN
  template <class T>
  int eigen_SymmEigensystem(Matrix<T> &Ain, std::vector<T> &evals, std::vector<std::vector<T> > &evecs, double small){
    using namespace Eigen;
    const int size = Ain.dim;
    evals.resize(size);
    evecs.resize(size);
    for (int i = 0; i < size; i++) {
      evecs[i].resize(size);
    }
    MatrixXd A(size, size);
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
	A(i, j) = Ain(i, j);
      }
    }
    SelfAdjointEigenSolver<MatrixXd> eigensolver(A);
    for (int i = 0; i < size; i++) {
      evals[i] = eigensolver.eigenvalues()(i);
    }
    for (int i = 0; i < size; i++) {
      std::vector<T>& vec = evecs[i];
      for (int j = 0; j < size; j++) {
	vec[j] = eigensolver.eigenvectors()(j, i);
      }
    }
    return 0;
  }
#endif

  template <class T>
  int SymmEigensystem(Matrix<T> &Ain, std::vector<T> &evals, std::vector<std::vector<T> > &evecs, double small){
#ifdef USE_EIGEN
    return eigen_SymmEigensystem(Ain, evals, evecs, small);
#else
    return basic_SymmEigensystem(Ain, evals, evecs, small);
#endif
  }

  ///Solve a general eigensystem, not necessarily in tridiagonal form
  template <class T> int Eigensystem(Matrix<T > &Ain, std::vector<T> &evals, std::vector<std::vector<T> > &evecs, double small){

    int N = Ain.dim;
    Matrix<T > A(N); A = Ain;
    Matrix<T > Q(N);Q.Unity();
    Hess(A,Q,0);
    int it = QReigensystem<T>(A, evals, evecs, small);
    for(int k=0;k<N;k++){evecs[k] = Q*evecs[k];}

    return it;
  }

}

#endif

max_stars_repo_name: RBC-UKQCD/CPS

max_issues_repo_name: RBC-UKQCD/CPS

max_forks_repo_name: RBC-UKQCD/CPS


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 std_9154454e_8e1d_4505_9034_5c7a71b77dd0_h
#define std_9154454e_8e1d_4505_9034_5c7a71b77dd0_h

#include <xutility>
#include <algorithm>
#include <memory>
#include <list>
#include <vector>
#include <deque>
#include <queue>
#include <map>
#include <set>
#include <stack>
#include <unordered_map>
#include <unordered_set>
#include <comip.h>

#include <gslib/type.h>
#include <gslib/string.h>

__gslib_begin__

class string;

/* used for base type delegate */
template<class _inner>
struct base_delegator
{
    _inner          _puppet;

    base_delegator() {}
    base_delegator(const _inner& i): _puppet(i) {}
    operator _inner() const { return _puppet; }
    _inner& operator=(const _inner& i) { return _puppet = i; }
};

template<class _inner>
struct base_ptr_delegator
{
    _inner*         _puppet;

    base_ptr_delegator() { _puppet = 0; }
    base_ptr_delegator(_inner* ptr) { _puppet = ptr; }
    operator _inner*() { return _puppet; }
    operator const _inner*() const { return _puppet; }
    _inner* operator=(_inner* ptr) { return _puppet = ptr; }
};

#if defined(_MSC_VER) && (_MSC_VER < 1910)
/* the same as the default std version */
inline size_t hash_bytes(const byte str[], int size)
{
    size_t val = 2166136261U;
    size_t first = 0;
    size_t last = (size_t)size;
    size_t stride = 1 + last / 10;
    for( ; first < last; first += stride)
        val = 16777619U * val ^ (size_t)str[first];
    return val;
}
#elif defined(_MSC_VER) && (_MSC_VER < 1916)
inline size_t hash_bytes(const byte str[], int size)
{ return std::_Hash_bytes(str, size); }
#else
inline size_t hash_bytes(const byte str[], int size)
{ return std::_Fnv1a_append_bytes(std::_FNV_offset_basis, str, size); }
#endif

#ifdef _UNICODE
inline size_t string_hash(const string& str) { return hash_bytes((const byte*)str.c_str(), str.size() * sizeof(wchar)); }
#else
inline size_t string_hash(const string& str) { return hash_bytes((const byte*)str.c_str(), str.size()); }
#endif

#if defined(_MSC_VER) && (_MSC_VER <= 1600)
class hasher
{
public:
    hasher() { _val = 2166136261u; }
    size_t add_bytes(const byte* str, int len)
    {
        size_t first = 0u;
        size_t last = (size_t)len;
        size_t stride = 1 + last / 10;
        for(; first < last; first += stride)
            _val = 16777619u * _val ^ (size_t)str[first];
        return _val;
    }
    size_t get_value() const { return _val; }

private:
    size_t          _val;
};
#elif defined(_MSC_VER) && (_MSC_VER < 1916)
class hasher:
    public std::_Fnv1a_hasher
{
public:
    size_t add_bytes(const byte* first, int len) { return _Add_bytes(first, first + len); }
    size_t get_value() const { return _Val; }
};
#else
class hasher
{
public:
    hasher() { _val = std::_FNV_offset_basis; }
    size_t add_bytes(const byte* first, int len) { return std::_Fnv1a_append_bytes(_val, first, len); }
    size_t get_value() const { return _val; }

private:
    size_t          _val;
};
#endif

/* switch to namespace tr1 and replace the hash for string */
__gslib_end__

namespace std {

template<>
class hash<gs::string>
#if defined(_MSC_VER) && (_MSC_VER < 1914)
    : public unary_function<gs::string, size_t>
#endif
{
public:
    size_t operator()(const gs::string& kval) const { return string_hash(kval); }
};

};
__gslib_begin__

/* wrappers for all kinds of vessels */
template<class _ty>
using list = std::list<_ty>;
template<class _ty>
using vector = std::vector<_ty>;
template<class _ty>
using queue = std::queue<_ty>;
template<class _ty>
using deque = std::deque<_ty>;
template<class _ty>
using stack = std::stack<_ty>;

template<class _kty, class _ty,
    class _pr = std::less<_kty>
    >
using map = std::map<_kty, _ty, _pr>;
template<class _ty, class _pr = std::less<_ty> >
using set = std::set<_ty, _pr>;

template<class _kty, class _ty,
    class _pr = std::less<_kty>
    >
using multimap = std::multimap<_kty, _ty, _pr>;
template<class _ty, class _pr = std::less<_ty> >
using multiset = std::multiset<_ty, _pr>;

template<class _kty, class _ty,
    class _hs = std::hash<_kty>,
    class _eq = std::equal_to<_kty>
    >
using unordered_map = std::unordered_map<_kty, _ty, _hs, _eq>;

template<class _kty,
    class _hs = std::hash<_kty>,
    class _equ = std::equal_to<_kty>
    >
using unordered_set = std::unordered_set<_kty, _hs, _equ>;

template<class _kty, class _ty,
    class _hs = std::hash<_kty>
    >
using unordered_multimap = std::unordered_multimap<_kty, _ty, _hs>;

template<class _kty, class _hs, class _equ>
using unordered_multiset = std::unordered_multiset<_kty, _hs, _equ>;

template<class _cls>
class com_ptr
{
public:
    typedef _cls* ptr_type;
    typedef com_ptr<_cls> myref;

public:
    com_ptr() { _ptr = nullptr; }
    ~com_ptr()
    {
        if(_ptr) {
            _ptr->Release();
            _ptr = nullptr;
        }
    }
    com_ptr(_cls* p)
    {
        _ptr = p;
        if(_ptr)
            _ptr->AddRef();
    }
    com_ptr(const myref& p)
    {
        _ptr = p.get();
        if(_ptr)
            _ptr->AddRef();
    }
    void attach(_cls* p)
    {
        if(_ptr)
            _ptr->Release();
        _ptr = p;
    }
    _cls* detach()
    {
        _cls* p = _ptr;
        _ptr = nullptr;
        return p;
    }
    void clear()
    {
        if(_ptr) {
            _ptr->Release();
            _ptr = nullptr;
        }
    }
    _cls* operator=(_cls* p)
    {
        if(p) p->AddRef();
        if(_ptr)
            _ptr->Release();
        return _ptr = p;
    }
    _cls* operator=(const myref& p)
    {
        if(p._ptr)
            p._ptr->AddRef();
        if(_ptr)
            _ptr->Release();
        return _ptr = p._ptr;
    }
    _cls** operator&()
    {
        assert(!_ptr);
        return &_ptr;
    }
    _cls* get() const { return _ptr; }
    _cls* operator->() const { return _ptr; }
    bool operator!() const { return !_ptr; }
    bool operator == (_cls* p) const { return _ptr == p; }
    bool operator != (_cls* p) const { return _ptr != p; }

private:
    _cls*           _ptr;
};

__gslib_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


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

/*
 * Author:  G. Jungman
 */
#ifndef __GSL_BLAS_H__
#define __GSL_BLAS_H__

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

#include <gsl/gsl_blas_types.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

/* ========================================================================
 * Level 1
 * ========================================================================
 */

GSL_EXPORT int gsl_blas_sdsdot (float alpha,
                                const gsl_vector_float * X,
                                const gsl_vector_float * Y,
                                float * result
                                );

GSL_EXPORT int gsl_blas_dsdot (const gsl_vector_float * X,
                               const gsl_vector_float * Y,
                               double * result
                               );

GSL_EXPORT int gsl_blas_sdot (const gsl_vector_float * X,
                              const gsl_vector_float * Y,
                              float * result
                              );

GSL_EXPORT int gsl_blas_ddot (const gsl_vector * X,
                              const gsl_vector * Y,
                              double * result
                              );

GSL_EXPORT int  gsl_blas_cdotu (const gsl_vector_complex_float * X,
                                const gsl_vector_complex_float * Y,
                                gsl_complex_float * dotu);

GSL_EXPORT int  gsl_blas_cdotc (const gsl_vector_complex_float * X,
                                const gsl_vector_complex_float * Y,
                                gsl_complex_float * dotc);

GSL_EXPORT int  gsl_blas_zdotu (const gsl_vector_complex * X,
                                const gsl_vector_complex * Y,
                                gsl_complex * dotu);

GSL_EXPORT int  gsl_blas_zdotc (const gsl_vector_complex * X,
                                const gsl_vector_complex * Y,
                                gsl_complex * dotc);

GSL_EXPORT float  gsl_blas_snrm2  (const gsl_vector_float * X);
GSL_EXPORT float  gsl_blas_sasum  (const gsl_vector_float * X);
GSL_EXPORT double gsl_blas_dnrm2  (const gsl_vector * X);
GSL_EXPORT double gsl_blas_dasum  (const gsl_vector * X);
GSL_EXPORT float  gsl_blas_scnrm2 (const gsl_vector_complex_float * X);
GSL_EXPORT float  gsl_blas_scasum (const gsl_vector_complex_float * X);
GSL_EXPORT double gsl_blas_dznrm2 (const gsl_vector_complex * X);
GSL_EXPORT double gsl_blas_dzasum (const gsl_vector_complex * X);

GSL_EXPORT CBLAS_INDEX_t gsl_blas_isamax (const gsl_vector_float * X);
GSL_EXPORT CBLAS_INDEX_t gsl_blas_idamax (const gsl_vector * X);
GSL_EXPORT CBLAS_INDEX_t gsl_blas_icamax (const gsl_vector_complex_float * X);
GSL_EXPORT CBLAS_INDEX_t gsl_blas_izamax (const gsl_vector_complex * X);

GSL_EXPORT int  gsl_blas_sswap (gsl_vector_float * X,
                                gsl_vector_float * Y);

GSL_EXPORT int  gsl_blas_scopy (const gsl_vector_float * X,
                                gsl_vector_float * Y);

GSL_EXPORT int  gsl_blas_saxpy (float alpha,
                                const gsl_vector_float * X,
                                gsl_vector_float * Y);

GSL_EXPORT int  gsl_blas_dswap (gsl_vector * X,
                                gsl_vector * Y);

GSL_EXPORT int  gsl_blas_dcopy (const gsl_vector * X,
                                gsl_vector * Y);

GSL_EXPORT int  gsl_blas_daxpy (double alpha,
                                const gsl_vector * X,
                                gsl_vector * Y);

GSL_EXPORT int  gsl_blas_cswap (gsl_vector_complex_float * X,
                                gsl_vector_complex_float * Y);

GSL_EXPORT int  gsl_blas_ccopy (const gsl_vector_complex_float * X,
                                gsl_vector_complex_float * Y);

GSL_EXPORT int  gsl_blas_caxpy (const gsl_complex_float alpha,
                                const gsl_vector_complex_float * X,
                                gsl_vector_complex_float * Y);

GSL_EXPORT int  gsl_blas_zswap (gsl_vector_complex * X,
                                gsl_vector_complex * Y);

GSL_EXPORT int  gsl_blas_zcopy (const gsl_vector_complex * X,
                                gsl_vector_complex * Y);

GSL_EXPORT int  gsl_blas_zaxpy (const gsl_complex alpha,
                                const gsl_vector_complex * X,
                                gsl_vector_complex * Y);

GSL_EXPORT int  gsl_blas_srotg (float a[], float b[], float c[], float s[]);

GSL_EXPORT int  gsl_blas_srotmg (float d1[], float d2[], float b1[], float b2, float P[]);

GSL_EXPORT int  gsl_blas_srot (gsl_vector_float * X,
                               gsl_vector_float * Y,
                               float c, float s);

GSL_EXPORT int  gsl_blas_srotm (gsl_vector_float * X,
                                gsl_vector_float * Y,
                                const float P[]);

GSL_EXPORT int  gsl_blas_drotg (double a[], double b[], double c[], double s[]);

GSL_EXPORT int  gsl_blas_drotmg (double d1[], double d2[], double b1[],
                                 double b2, double P[]);

GSL_EXPORT int  gsl_blas_drot (gsl_vector * X,
                               gsl_vector * Y,
                               const double c, const double s);

GSL_EXPORT int  gsl_blas_drotm (gsl_vector * X,
                                gsl_vector * Y,
                                const double P[]);

GSL_EXPORT void gsl_blas_sscal  (float  alpha, gsl_vector_float * X);
GSL_EXPORT void gsl_blas_dscal  (double alpha, gsl_vector * X);
GSL_EXPORT void gsl_blas_cscal  (const gsl_complex_float alpha, gsl_vector_complex_float * X);
GSL_EXPORT void gsl_blas_zscal  (const gsl_complex alpha, gsl_vector_complex * X);
GSL_EXPORT void gsl_blas_csscal (float  alpha, gsl_vector_complex_float * X);
GSL_EXPORT void gsl_blas_zdscal (double alpha, gsl_vector_complex * X);

/* ===========================================================================
 * Level 2
 * ===========================================================================
 */

/*
 * Routines with standard 4 prefixes (S, D, C, Z)
 */
GSL_EXPORT int  gsl_blas_sgemv (CBLAS_TRANSPOSE_t TransA,
                                float alpha,
                                const gsl_matrix_float * A,
                                const gsl_vector_float * X,
                                float beta,
                                gsl_vector_float * Y);

GSL_EXPORT int  gsl_blas_strmv (CBLAS_UPLO_t Uplo,
                                CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                                const gsl_matrix_float * A,
                                gsl_vector_float * X);

GSL_EXPORT int  gsl_blas_strsv (CBLAS_UPLO_t Uplo,
                                CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                                const gsl_matrix_float * A,
                                gsl_vector_float * X);

GSL_EXPORT int  gsl_blas_dgemv (CBLAS_TRANSPOSE_t TransA,
                                double alpha,
                                const gsl_matrix * A,
                                const gsl_vector * X,
                                double beta,
                                gsl_vector * Y);

GSL_EXPORT int  gsl_blas_dtrmv (CBLAS_UPLO_t Uplo,
                                CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                                const gsl_matrix * A,
                                gsl_vector * X);

GSL_EXPORT int  gsl_blas_dtrsv (CBLAS_UPLO_t Uplo,
                                CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                                const gsl_matrix * A,
                                gsl_vector * X);

GSL_EXPORT int  gsl_blas_cgemv (CBLAS_TRANSPOSE_t TransA,
                                const gsl_complex_float alpha,
                                const gsl_matrix_complex_float * A,
                                const gsl_vector_complex_float * X,
                                const gsl_complex_float beta,
                                gsl_vector_complex_float * Y);

GSL_EXPORT int  gsl_blas_ctrmv (CBLAS_UPLO_t Uplo,
                                CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                                const gsl_matrix_complex_float * A,
                                gsl_vector_complex_float * X);

GSL_EXPORT int  gsl_blas_ctrsv (CBLAS_UPLO_t Uplo,
                                CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                                const gsl_matrix_complex_float * A,
                                gsl_vector_complex_float * X);

GSL_EXPORT int  gsl_blas_zgemv (CBLAS_TRANSPOSE_t TransA,
                                const gsl_complex alpha,
                                const gsl_matrix_complex * A,
                                const gsl_vector_complex * X,
                                const gsl_complex beta,
                                gsl_vector_complex * Y);

GSL_EXPORT int  gsl_blas_ztrmv (CBLAS_UPLO_t Uplo,
                                CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                                const gsl_matrix_complex * A,
                                gsl_vector_complex * X);

GSL_EXPORT int  gsl_blas_ztrsv (CBLAS_UPLO_t Uplo,
                                CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                                const gsl_matrix_complex * A,
                                gsl_vector_complex *X);

/*
 * Routines with S and D prefixes only
 */
GSL_EXPORT int  gsl_blas_ssymv (CBLAS_UPLO_t Uplo,
                                float alpha,
                                const gsl_matrix_float * A,
                                const gsl_vector_float * X,
                                float beta,
                                gsl_vector_float * Y);

GSL_EXPORT int  gsl_blas_sger (float alpha,
                               const gsl_vector_float * X,
                               const gsl_vector_float * Y,
                               gsl_matrix_float * A);

GSL_EXPORT int  gsl_blas_ssyr (CBLAS_UPLO_t Uplo,
                               float alpha,
                               const gsl_vector_float * X,
                               gsl_matrix_float * A);

GSL_EXPORT int  gsl_blas_ssyr2 (CBLAS_UPLO_t Uplo,
                                float alpha,
                                const gsl_vector_float * X,
                                const gsl_vector_float * Y,
                                gsl_matrix_float * A);

GSL_EXPORT int  gsl_blas_dsymv (CBLAS_UPLO_t Uplo,
                                double alpha,
                                const gsl_matrix * A,
                                const gsl_vector * X,
                                double beta,
                                gsl_vector * Y);

GSL_EXPORT int  gsl_blas_dger (double alpha,
                               const gsl_vector * X,
                               const gsl_vector * Y,
                               gsl_matrix * A);

GSL_EXPORT int  gsl_blas_dsyr (CBLAS_UPLO_t Uplo,
                               double alpha,
                               const gsl_vector * X,
                               gsl_matrix * A);

GSL_EXPORT int  gsl_blas_dsyr2 (CBLAS_UPLO_t Uplo,
                                double alpha,
                                const gsl_vector * X,
                                const gsl_vector * Y,
                                gsl_matrix * A);

/*
 * Routines with C and Z prefixes only
 */

GSL_EXPORT int  gsl_blas_chemv (CBLAS_UPLO_t Uplo,
                                const gsl_complex_float alpha,
                                const gsl_matrix_complex_float * A,
                                const gsl_vector_complex_float * X,
                                const gsl_complex_float beta,
                                gsl_vector_complex_float * Y);

GSL_EXPORT int  gsl_blas_cgeru (const gsl_complex_float alpha,
                                const gsl_vector_complex_float * X,
                                const gsl_vector_complex_float * Y,
                                gsl_matrix_complex_float * A);

GSL_EXPORT int  gsl_blas_cgerc (const gsl_complex_float alpha,
                                const gsl_vector_complex_float * X,
                                const gsl_vector_complex_float * Y,
                                gsl_matrix_complex_float * A);

GSL_EXPORT int  gsl_blas_cher (CBLAS_UPLO_t Uplo,
                              float alpha,
                              const gsl_vector_complex_float * X,
                              gsl_matrix_complex_float * A);

GSL_EXPORT int  gsl_blas_cher2 (CBLAS_UPLO_t Uplo,
                                const gsl_complex_float alpha,
                                const gsl_vector_complex_float * X,
                                const gsl_vector_complex_float * Y,
                                gsl_matrix_complex_float * A);

GSL_EXPORT int  gsl_blas_zhemv (CBLAS_UPLO_t Uplo,
                                const gsl_complex alpha,
                                const gsl_matrix_complex * A,
                                const gsl_vector_complex * X,
                                const gsl_complex beta,
                                gsl_vector_complex * Y);

GSL_EXPORT int  gsl_blas_zgeru (const gsl_complex alpha,
                                const gsl_vector_complex * X,
                                const gsl_vector_complex * Y,
                                gsl_matrix_complex * A);

GSL_EXPORT int  gsl_blas_zgerc (const gsl_complex alpha,
                                const gsl_vector_complex * X,
                                const gsl_vector_complex * Y,
                                gsl_matrix_complex * A);

GSL_EXPORT int  gsl_blas_zher (CBLAS_UPLO_t Uplo,
                              double alpha,
                              const gsl_vector_complex * X,
                              gsl_matrix_complex * A);

GSL_EXPORT int  gsl_blas_zher2 (CBLAS_UPLO_t Uplo,
                                const gsl_complex alpha,
                                const gsl_vector_complex * X,
                                const gsl_vector_complex * Y,
                                gsl_matrix_complex * A);

/*
 * ===========================================================================
 * Prototypes for level 3 BLAS
 * ===========================================================================
 */

/*
 * Routines with standard 4 prefixes (S, D, C, Z)
 */
GSL_EXPORT int  gsl_blas_sgemm (CBLAS_TRANSPOSE_t TransA,
                                CBLAS_TRANSPOSE_t TransB,
                                float alpha,
                                const gsl_matrix_float * A,
                                const gsl_matrix_float * B,
                                float beta,
                                gsl_matrix_float * C);

GSL_EXPORT int  gsl_blas_ssymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                                float alpha,
                                const gsl_matrix_float * A,
                                const gsl_matrix_float * B,
                                float beta,
                                gsl_matrix_float * C);

GSL_EXPORT int  gsl_blas_ssyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
                                float alpha,
                                const gsl_matrix_float * A,
                                float beta,
                                gsl_matrix_float * C);

GSL_EXPORT int  gsl_blas_ssyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
                                 float alpha,
                                 const gsl_matrix_float * A,
                                 const gsl_matrix_float * B,
                                 float beta,
                                 gsl_matrix_float * C);

GSL_EXPORT int  gsl_blas_strmm (CBLAS_SIDE_t Side,
                                CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                                CBLAS_DIAG_t Diag,
                                float alpha,
                                const gsl_matrix_float * A,
                                gsl_matrix_float * B);

GSL_EXPORT int  gsl_blas_strsm (CBLAS_SIDE_t Side,
                                CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                                CBLAS_DIAG_t Diag,
                                float alpha,
                                const gsl_matrix_float * A,
                                gsl_matrix_float * B);

GSL_EXPORT int  gsl_blas_dgemm (CBLAS_TRANSPOSE_t TransA,
                                CBLAS_TRANSPOSE_t TransB,
                                double alpha,
                                const gsl_matrix * A,
                                const gsl_matrix * B,
                                double beta,
                                gsl_matrix * C);

GSL_EXPORT int  gsl_blas_dsymm (CBLAS_SIDE_t Side,
                                CBLAS_UPLO_t Uplo,
                                double alpha,
                                const gsl_matrix * A,
                                const gsl_matrix * B,
                                double beta,
                                gsl_matrix * C);

GSL_EXPORT int  gsl_blas_dsyrk (CBLAS_UPLO_t Uplo,
                                CBLAS_TRANSPOSE_t Trans,
                                double alpha,
                                const gsl_matrix * A,
                                double beta,
                                gsl_matrix * C);

GSL_EXPORT int  gsl_blas_dsyr2k (CBLAS_UPLO_t Uplo,
                                 CBLAS_TRANSPOSE_t Trans,
                                 double alpha,
                                 const  gsl_matrix * A,
                                 const  gsl_matrix * B,
                                 double beta,
                                 gsl_matrix * C);

GSL_EXPORT int  gsl_blas_dtrmm (CBLAS_SIDE_t Side,
                                CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                                CBLAS_DIAG_t Diag,
                                double alpha,
                                const gsl_matrix * A,
                                gsl_matrix * B);

GSL_EXPORT int  gsl_blas_dtrsm (CBLAS_SIDE_t Side,
                                CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                                CBLAS_DIAG_t Diag,
                                double alpha,
                                const gsl_matrix * A,
                                gsl_matrix * B);

GSL_EXPORT int  gsl_blas_cgemm (CBLAS_TRANSPOSE_t TransA,
                                CBLAS_TRANSPOSE_t TransB,
                                const gsl_complex_float alpha,
                                const gsl_matrix_complex_float * A,
                                const gsl_matrix_complex_float * B,
                                const gsl_complex_float beta,
                                gsl_matrix_complex_float * C);

GSL_EXPORT int  gsl_blas_csymm (CBLAS_SIDE_t Side,
                                CBLAS_UPLO_t Uplo,
                                const gsl_complex_float alpha,
                                const gsl_matrix_complex_float * A,
                                const gsl_matrix_complex_float * B,
                                const gsl_complex_float beta,
                                gsl_matrix_complex_float * C);

GSL_EXPORT int  gsl_blas_csyrk (CBLAS_UPLO_t Uplo,
                                CBLAS_TRANSPOSE_t Trans,
                                const gsl_complex_float alpha,
                                const gsl_matrix_complex_float * A,
                                const gsl_complex_float beta,
                                gsl_matrix_complex_float * C);

GSL_EXPORT int  gsl_blas_csyr2k (CBLAS_UPLO_t Uplo,
                                 CBLAS_TRANSPOSE_t Trans,
                                 const gsl_complex_float alpha,
                                 const gsl_matrix_complex_float * A,
                                 const gsl_matrix_complex_float * B,
                                 const gsl_complex_float beta,
                                 gsl_matrix_complex_float * C);

GSL_EXPORT int  gsl_blas_ctrmm (CBLAS_SIDE_t Side,
                                CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                                CBLAS_DIAG_t Diag,
                                const gsl_complex_float alpha,
                                const gsl_matrix_complex_float * A,
                                gsl_matrix_complex_float * B);

GSL_EXPORT int  gsl_blas_ctrsm (CBLAS_SIDE_t Side,
                                CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                                CBLAS_DIAG_t Diag,
                                const gsl_complex_float alpha,
                                const gsl_matrix_complex_float * A,
                                gsl_matrix_complex_float * B);

GSL_EXPORT int  gsl_blas_zgemm (CBLAS_TRANSPOSE_t TransA,
                                CBLAS_TRANSPOSE_t TransB,
                                const gsl_complex alpha,
                                const gsl_matrix_complex * A,
                                const gsl_matrix_complex * B,
                                const gsl_complex beta,
                                gsl_matrix_complex * C);

GSL_EXPORT int  gsl_blas_zsymm (CBLAS_SIDE_t Side,
                                CBLAS_UPLO_t Uplo,
                                const gsl_complex alpha,
                                const gsl_matrix_complex * A,
                                const gsl_matrix_complex * B,
                                const gsl_complex beta,
                                gsl_matrix_complex * C);

GSL_EXPORT int  gsl_blas_zsyrk (CBLAS_UPLO_t Uplo,
                                CBLAS_TRANSPOSE_t Trans,
                                const gsl_complex alpha,
                                const gsl_matrix_complex * A,
                                const gsl_complex beta,
                                gsl_matrix_complex * C);

GSL_EXPORT int  gsl_blas_zsyr2k (CBLAS_UPLO_t Uplo,
                                 CBLAS_TRANSPOSE_t Trans,
                                 const gsl_complex alpha,
                                 const gsl_matrix_complex * A,
                                 const gsl_matrix_complex * B,
                                 const gsl_complex beta,
                                 gsl_matrix_complex *C);

GSL_EXPORT int  gsl_blas_ztrmm (CBLAS_SIDE_t Side,
                                CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                                CBLAS_DIAG_t Diag,
                                const gsl_complex alpha,
                                const gsl_matrix_complex * A,
                                gsl_matrix_complex * B);

GSL_EXPORT int  gsl_blas_ztrsm (CBLAS_SIDE_t Side,
                                CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                                CBLAS_DIAG_t Diag,
                                const gsl_complex alpha,
                                const gsl_matrix_complex * A,
                                gsl_matrix_complex * B);

/*
 * Routines with prefixes C and Z only
 */
GSL_EXPORT int  gsl_blas_chemm (CBLAS_SIDE_t Side,
                                CBLAS_UPLO_t Uplo,
                                const gsl_complex_float alpha,
                                const gsl_matrix_complex_float * A,
                                const gsl_matrix_complex_float * B,
                                const gsl_complex_float beta,
                                gsl_matrix_complex_float * C);

GSL_EXPORT int  gsl_blas_cherk (CBLAS_UPLO_t Uplo,
                                CBLAS_TRANSPOSE_t Trans,
                                float alpha,
                                const gsl_matrix_complex_float * A,
                                float beta,
                                gsl_matrix_complex_float * C);

GSL_EXPORT int  gsl_blas_cher2k (CBLAS_UPLO_t Uplo,
                                 CBLAS_TRANSPOSE_t Trans,
                                 const gsl_complex_float alpha,
                                 const gsl_matrix_complex_float * A,
                                 const gsl_matrix_complex_float * B,
                                 float beta,
                                 gsl_matrix_complex_float * C);

GSL_EXPORT int  gsl_blas_zhemm (CBLAS_SIDE_t Side,
                                CBLAS_UPLO_t Uplo,
                                const gsl_complex alpha,
                                const gsl_matrix_complex * A,
                                const gsl_matrix_complex * B,
                                const gsl_complex beta,
                                gsl_matrix_complex * C);

GSL_EXPORT int  gsl_blas_zherk (CBLAS_UPLO_t Uplo,
                                CBLAS_TRANSPOSE_t Trans,
                                double alpha,
                                const gsl_matrix_complex * A,
                                double beta,
                                gsl_matrix_complex * C);

GSL_EXPORT int  gsl_blas_zher2k (CBLAS_UPLO_t Uplo,
                                 CBLAS_TRANSPOSE_t Trans,
                                 const gsl_complex alpha,
                                 const gsl_matrix_complex * A,
                                 const gsl_matrix_complex * B,
                                 double beta,
                                 gsl_matrix_complex * C);

__END_DECLS

#endif /* __GSL_BLAS_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws