text:
/* To execute:
 * ./driver <n> # for some positive integer n
 * */
#include "q_incs.h"
#include <lapacke.h>
#include "matrix_helpers.h"
#include "positive_solver.h"

/* testing and benchmarking for AxEqualsB solver.
 * generates a random matrix A, and tests the positive semidefinite solver
 * by running iton AtA and tests the full solver by running it on A.
 *
 * the current method of generating random matrices isn't perfect, but it's
 * better than it seems because solving Ax = b is independent of any scaling
 * factor, so the size of the matrix is irrelevant.
 */

static int _bench_drop_n = 0;
static int _bench_iters = 1;
static bool _verbose = false;

static void
_print_input(
            double **A,
            double *x,
            double *b,
            int n
            )
{
  for ( int i = 0; i < n; i++ ) {
    fprintf(stderr, "[ ");
    for ( int j = 0; j < n; j++ ) {
      if ( j == 0 ) {
        fprintf(stderr, " %2d ", (int)index_symm_matrix(A, i, j));
      }
      else {
        fprintf(stderr, ", %2d ", (int)index_symm_matrix(A, i, j));
      }
    }
    fprintf(stderr, "] ");
    // print x
    fprintf(stderr, " [ ");
    fprintf(stderr, " %2d ", (int)x[i]);
    fprintf(stderr, "] ");
    // print b
    fprintf(stderr, "  = [ ");
    fprintf(stderr, " %3d ", (int)b[i]);
    fprintf(stderr, "] ");

    fprintf(stderr, "\n");
  }
}

/* assembly code to read the TSC */
static inline uint64_t
_RDTSC(void)
{
  unsigned int hi, lo;
  __asm__ volatile("rdtsc" : "=a" (lo), "=d" (hi));
  return ((uint64_t)hi << 32) | lo;
}

static void
_print_result(double *x_expected, double *x_returned,
              double *b_expected, double *b_returned,
              char *name, int n, double runtime, bool checker_success)
{
  bool error = false;
  fprintf(stderr, "CHECKING %s RESULTS\n", name);
  if (_verbose) {
    fprintf(stderr, "x(expect), x(%s), b(expect), b(%s)\n", name, name);
  }
  for (int i=0; i < n; i++) {
    if (_verbose) {
      fprintf(stderr, " %lf, %lf, %lf, %lf",
              x_expected[i], x_returned[i], b_expected[i], b_returned[i]);
    }

    if ( abs(b_returned[i] - b_expected[i]) < 0.001 ) {
      if (_verbose) {
        fprintf(stderr, " %s_MATCH\n", name);
      }
    }
    else {
      error = true;
      if (_verbose) {
        fprintf(stderr, " %s_ERROR\n", name);
      }
    }
  }
  fprintf(stderr, "%s %s in %.3e cycles, checker said %s.\n\n",
          name,
          error ? "FAILED" : "SUCCEEDED",
          runtime,
          checker_success ? "SUCCEEDED" : "FAILED");
}

static int
_run_our_tests(
    double **A,
    double *x_expected,
    double *b,
    double **A_posdef,
    double *b_posdef,
    int n
    )
{
  int status = 0;

  double **A_posdef_copy = NULL;
  double **A_posdef_full = NULL;
  double *b_posdef_copy = NULL;
  double *b_returned = NULL;
  double *x_returned = NULL;
  uint64_t begin, total;
  double avg_cycles;
  bool checker_success;

  status = alloc_symm_matrix(&A_posdef_copy, n); cBYE(status);
  status = alloc_matrix(&A_posdef_full, n); cBYE(status);
  b_posdef_copy = malloc(n * sizeof(double)); return_if_malloc_failed(b_posdef_copy);
  b_returned = malloc(n * sizeof(double)); return_if_malloc_failed(b_returned);
  x_returned = malloc(n * sizeof(double)); return_if_malloc_failed(x_returned);

  for (int i = 0; i < n; i++) {
    b_posdef_copy[i] = b_posdef[i];
    for(int j = 0; j < n; j++) {
      A_posdef_full[i][j] = index_symm_matrix(A_posdef, i, j);
      if (j < n - i) {
        A_posdef_copy[i][j] = A_posdef[i][j];
      }
    }
  }

  // time computations for each solver
  total = 0;
  for (int i = 0; i < _bench_iters; i++) {
    begin = _RDTSC();
    status = posdef_positive_solver(A_posdef, x_returned, b_posdef, n); cBYE(status);
    if (i >= _bench_drop_n) {
      total += _RDTSC() - begin;
    }
  }
  avg_cycles = (double)total / (_bench_iters - _bench_drop_n);
  multiply_symm_matrix_vector(A_posdef, x_returned, n, b_returned);
  checker_success = symm_positive_solver_check(A_posdef, x_returned, b_posdef, n, 0.0);
  _print_result(x_expected, x_returned, b_posdef, b_returned, "POSDEF_SLOW", n, avg_cycles, checker_success);

  total = 0;
  for (int i = 0; i < _bench_iters; i++) {
    for (int i = 0; i < n; i++) {
      b_posdef_copy[i] = b_posdef[i];
      for(int j = 0; j < n - i; j++) {
        A_posdef_copy[i][j] = A_posdef[i][j];
      }
    }
    begin = _RDTSC();
    status = posdef_positive_solver_fast(A_posdef_copy, x_returned, b_posdef_copy, n); cBYE(status); cBYE(status);
    if (i >= _bench_drop_n) {
      total += _RDTSC() - begin;
    }
  }
  avg_cycles = (double)total / (_bench_iters - _bench_drop_n);
  multiply_symm_matrix_vector(A_posdef, x_returned, n, b_returned);
  checker_success = symm_positive_solver_check(A_posdef, x_returned, b_posdef, n, 0.0);
  _print_result(x_expected, x_returned, b_posdef, b_returned, "POSDEF_FAST", n, avg_cycles, checker_success);

  total = 0;
  for (int i = 0; i < _bench_iters; i++) {
    for (int i = 0; i < n; i++) {
      for(int j = 0; j < n; j++) {
        A_posdef_full[i][j] = index_symm_matrix(A_posdef, i, j);
      }
    }
    begin = _RDTSC();
    status = full_posdef_positive_solver(A_posdef_full, x_returned, b_posdef, n); cBYE(status);
    if (i >= _bench_drop_n) {
      total += _RDTSC() - begin;
    }
  }
  avg_cycles = (double)total / (_bench_iters - _bench_drop_n);
  multiply_symm_matrix_vector(A_posdef, x_returned, n, b_returned);
  checker_success = symm_positive_solver_check(A_posdef, x_returned, b_posdef, n, 0.0);
  _print_result(x_expected, x_returned, b_posdef, b_returned, "FULL_POSDEF_SLOW", n, avg_cycles, checker_success);

  total = 0;
  for (int i = 0; i < _bench_iters; i++) {
    for (int i = 0; i < n; i++) {
      b_posdef_copy[i] = b_posdef[i];
      for(int j = 0; j < n; j++) {
        A_posdef_full[i][j] = index_symm_matrix(A_posdef, i, j);
      }
    }
    begin = _RDTSC();
    status = full_posdef_positive_solver_fast(A_posdef_full, x_returned, b_posdef_copy, n); cBYE(status);
    if (i >= _bench_drop_n) {
      total += _RDTSC() - begin;
    }
  }
  avg_cycles = (double)total / (_bench_iters - _bench_drop_n);
  multiply_symm_matrix_vector(A_posdef, x_returned, n, b_returned);
  checker_success = symm_positive_solver_check(A_posdef, x_returned, b_posdef, n, 0.0);
  _print_result(x_expected, x_returned, b_posdef, b_returned, "FULL_POSDEF_FAST", n, avg_cycles, checker_success);

  total = 0;
  for (int i = 0; i < _bench_iters; i++) {
    begin = _RDTSC();
    status = positive_solver(A, x_returned, b, n); cBYE(status);
    if (i >= _bench_drop_n) {
      total += _RDTSC() - begin;
    }
  }
  avg_cycles = (double)total / (_bench_iters - _bench_drop_n);
  multiply_matrix_vector(A, x_returned, n, b_returned);
  checker_success = full_positive_solver_check(A, x_returned, b, n, 0.0);
  _print_result(x_expected, x_returned, b, b_returned, "FULL", n, avg_cycles, checker_success);

BYE:
  free_matrix(A_posdef_copy, n);
  free_matrix(A_posdef_full, n);
  free_if_non_null(b_posdef_copy);
  free_if_non_null(b_returned);
  free_if_non_null(x_returned);

  return status;
}

static int
_run_lapack_tests(
    double **A,
    double *x_expected,
    double *b,
    double **A_posdef,
    double *b_posdef,
    int n
    )
{
  int status = 0;
  lapack_int N = n;
  lapack_int NRHS = 1;
  lapack_int LDA = N;
  lapack_int LDB = N;

  double *A_unrolled = NULL;
  double *A_posdef_unrolled = NULL;
  double *b_copy = NULL;
  double *b_posdef_copy = NULL;
  double *b_returned = NULL;
  lapack_int *ipiv = NULL;

  uint64_t begin, total;
  double avg_cycles;
  bool checker_success;

  A_unrolled = malloc(n * n * sizeof(double)); return_if_malloc_failed(A_unrolled);
  A_posdef_unrolled = malloc(n * n * sizeof(double)); return_if_malloc_failed(A_posdef_unrolled);
  b_copy = malloc(n * sizeof(double)); return_if_malloc_failed(b_copy);
  b_posdef_copy = malloc(n * sizeof(double)); return_if_malloc_failed(b_posdef_copy);
  b_returned = malloc(n * sizeof(double)); return_if_malloc_failed(b_returned);
  ipiv = malloc(n * sizeof(lapack_int)); return_if_malloc_failed(ipiv);

  total = 0;
  for (int i = 0; i < _bench_iters; i++) {
    for (int i = 0; i < n; i++) {
      b_copy[i] = b[i];
      for (int j = 0; j < n; j++) {
        A_unrolled[n * i + j] = A[i][j];
      }
    }
    begin = _RDTSC();
    LAPACKE_dgesv(LAPACK_COL_MAJOR, N, NRHS, A_unrolled, LDA, ipiv, b_copy, LDB);
    if (i >= _bench_drop_n) {
      total += _RDTSC() - begin;
    }
  }
  avg_cycles = (double)total / (_bench_iters - _bench_drop_n);
  multiply_matrix_vector(A, b_copy, n, b_returned);
  checker_success = full_positive_solver_check(A, b_copy, b, n, 0.0);
  _print_result(x_expected, b_copy, b, b_returned, "LAPACK_FULL", n, avg_cycles, checker_success);

  total = 0;
  for (int i = 0; i < _bench_iters; i++) {
    for (int i = 0; i < n; i++) {
      b_posdef_copy[i] = b_posdef[i];
      for (int j = 0; j <= i; j++) {
        A_posdef_unrolled[n * i + j] = A_posdef[j][i - j];
      }
      for (int j = i + 1; j < n; j++) {
        A_posdef_unrolled[n * i + j] = 0;
      }
    }
    begin = _RDTSC();
    LAPACKE_dposv(LAPACK_COL_MAJOR, 'U', N, NRHS, A_posdef_unrolled, LDA, b_posdef_copy, LDB);
    if (i >= _bench_drop_n) {
      total += _RDTSC() - begin;
    }
  }
  avg_cycles = (double)total / (_bench_iters - _bench_drop_n);
  multiply_symm_matrix_vector(A_posdef, b_posdef_copy, n, b_returned);
  checker_success = symm_positive_solver_check(A_posdef, b_posdef_copy, b_posdef, n, 0.0);
  _print_result(x_expected, b_posdef_copy, b_posdef, b_returned, "LAPACK_POSDEF", n, avg_cycles, checker_success);

BYE:
  free_if_non_null(A_unrolled);
  free_if_non_null(A_posdef_unrolled);
  free_if_non_null(b_copy);
  free_if_non_null(b_posdef_copy);
  free_if_non_null(b_returned);
  free_if_non_null(ipiv);

  return status;
}

int
main(
    int argc,
    char **argv
    )
{
  int status = 0;
  int n = 0;                 // dimension of matrices
  double **A = NULL;         // randomly generated n * n matrix
  double **A_posdef = NULL;  // randomly generated positive definite matrix (right now just A transpose * A)
  double *x_expected = NULL; // randomly generated solution
  double *b_posdef = NULL;   // A_posdef * x_expected
  double *b_full = NULL;     // A * x_expected

  srand48(_RDTSC());

  bool print_usage = false, bench = false;
  switch ( argc ) {
  case 3 :
    if (strcmp("-vb", argv[2]) == 0 || strcmp("-bv", argv[2]) == 0) {
      _verbose = true;
      bench = true;
    } else if (strcmp("-v", argv[2]) == 0) {
      _verbose = true;
    } else if (strcmp("-b", argv[2]) == 0) {
      bench = true;
    } else {
      print_usage = true;
    } // fall through
  case 2 :
    n = atoi(argv[1]);
    break;
  default :
    print_usage = true;
    break;
  }
  if (print_usage || n <= 0) {
    printf("Usage: ./test_driver <n> [-v]\n");
    printf("where n is a positive integer and -v provides verbose output.\n");
    go_BYE(-1);
  }
  if (bench == true) {
    _bench_iters = 15;
    _bench_drop_n = 3;
  }

  status = alloc_matrix(&A, n); cBYE(status);
  status = alloc_symm_matrix(&A_posdef, n); cBYE(status);

  for ( int i = 0; i < n; i++ )  {
    for ( int j = 0; j < n; j++ ) {
      A[i][j] = (drand48() - 0.5) * 100;
    }
  }
  transpose_and_multiply(A, A_posdef, n);
  b_posdef = (double *) malloc(n * sizeof(double)); return_if_malloc_failed(b_posdef);
  b_full = (double *) malloc(n * sizeof(double)); return_if_malloc_failed(b_full);
  x_expected = (double *) malloc(n * sizeof(double)); return_if_malloc_failed(x_expected);

  // Initialize x and b
  for ( int i = 0; i < n; i++ )  {
    x_expected[i] = (lrand48() % 16) - 16/2;
  }
  multiply_symm_matrix_vector(A_posdef, x_expected, n, b_posdef);
  multiply_matrix_vector(A, x_expected, n, b_full);

  if (_verbose) {
    _print_input(A_posdef, x_expected, b_posdef, n);
  }

  _run_our_tests(A, x_expected, b_full, A_posdef, b_posdef, n);

  _run_lapack_tests(A, x_expected, b_full, A_posdef, b_posdef, n);

BYE:
  free_matrix(A, n);
  free_matrix(A_posdef, n);
  free_if_non_null(x_expected);
  free_if_non_null(b_posdef);
  free_if_non_null(b_full);

  return status;
}

max_stars_repo_name: subramon/qlu

max_issues_repo_name: subramon/qlu

max_forks_repo_name: subramon/qlu


text:
#pragma once

#include <bitset>
#include <type_traits>

#include <cuda/define_specifiers.hpp>

#include <gsl-lite/gsl-lite.hpp>

namespace thrustshift {

namespace detail {

// to be compiler independent std::bitset is used here
template <typename I>
int count_leading_zeros_cpu(I i) {

	std::bitset<sizeof(I) * 8> bs(i);
	static_assert(bs.size() > 0);
	int lz = 0;
	for (int bi = gsl_lite::narrow<int>(sizeof(I) * 8) - 1; bi >= 0; --bi) {
		if (bs[bi]) {
			break;
		}
		++lz;
	}
	return lz;
}

template <typename I>
int count_leading_zeros_gpu(I i);

template<> CUDA_FD int count_leading_zeros_gpu<int>(int i) {
	return __clz(i);
}

template<> CUDA_FD int count_leading_zeros_gpu<long long int>(long long int i) {
	return __clzll(i);
}

} // namespace detail

template <typename I>
CUDA_FHD int count_leading_zeros(I i) {

#ifdef __CUDA_ARCH__
	return detail::count_leading_zeros_gpu(i);
#else
	return detail::count_leading_zeros_cpu(i);
#endif

}

} // namespace thrustshift

ext: h

max_stars_repo_name: pauleonix/thrustshift

max_issues_repo_name: pauleonix/thrustshift

max_forks_repo_name: pauleonix/thrustshift


text:
#pragma once

#include <stdint.h>

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

#include "nn_math.h"
#include "nn_utils.h"

typedef struct {
    size_t size;
    uintmax_t *layers;
    gsl_vector **biases;
    gsl_matrix **weights;
} nn_network;

nn_network * nn_network_create(const size_t n, const uintmax_t *layers);
void nn_network_destroy(nn_network *network);
void nn_network_set_biases(nn_network *network, gsl_vector **biases);
void nn_network_get_biases(nn_network *network, gsl_vector **biases);
void nn_network_set_weights(nn_network *network, gsl_matrix **weights);
void nn_network_get_weights(nn_network *network, gsl_matrix **weights);
gsl_vector * nn_network_ff(nn_network *network, gsl_vector *x);

max_stars_repo_name: mygulamali/neural-net

max_issues_repo_name: mygulamali/neural-net

max_forks_repo_name: mygulamali/neural-net


text:
/*
** t-tests
**
** G.Lohmann, Jan 2017
*/

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

#include <gsl/gsl_cdf.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#define TINY 1.0e-10
#define ABS(x) ((x) > 0 ? (x) : -(x))
#define SQR(x) ((x)*(x))

/*
** convert t to z values
*/
double t2z(double t,double df)
{
  double p=0,z=0;
  double a,b,x;
  extern double gsl_sf_beta_inc(double,double,double);

  /* t to p */
  x = df/(df+(t*t));
  if (x <= 0 || x > 1) return 0;

  a = 0.5*df;
  b = 0.5;
  p = 0.5 * gsl_sf_beta_inc(a,b,x);

  /* p to z */
  z = gsl_cdf_ugaussian_Qinv(p);
  return z;
}

void avevar(double *data,int n,double *a,double *v)
{
  int j;
  double ave,var,nx,s,u;

  nx = (double)n;
  ave = 0;
  for (j=0; j<n; j++) ave += data[j];
  ave /= nx;

  var = u = 0;
  for (j=0; j<n; j++) {
    s = data[j]-ave;
    u   += s;
    var += s*s;
  }
  var=(var-u*u/nx)/(nx-1);

  *v = var;
  *a = ave;
}

/* welch test, unequal variance */
double welchtest(double *data1,double *data2,int n1,int n2)
{
  double tiny=TINY;
  double ave1,ave2,var1,var2;
  double z=0,t=0,df=0;
  double nx1 = (double)n1;
  double nx2 = (double)n2;

  avevar(data1,n1,&ave1,&var1);
  if (var1 < tiny) return 0;
  avevar(data2,n2,&ave2,&var2);
  if (var2 < tiny) return 0;
  t = (ave1 - ave2)/sqrt(var1/nx1 + var2/nx2);
  df = SQR(var1/nx1+var2/nx2)/(SQR(var1/nx1)/(nx1-1)+SQR(var2/nx2)/(nx2-1));
  z  = t2z((double) t,(double) df);
  if (t < 0) z = -z;
  return z;
}

/* paired twosample test */
double paired_ttest(double *data1,double *data2,int n)
{
  int j;
  double ave,var,tiny=TINY;
  for (j=0; j<n; j++) {
    data1[j] -= data2[j];
  }
  avevar(data1,n,&ave,&var);
  if (var < tiny) return 0;
  double nx = (double)n;
  double df = nx-1.0;
  double t = sqrt(nx) * ave/sqrt(var);
  double z = t2z(t,df);
  if (t < 0) z = -z;
  return z;
}

/* paired twosample test */
double xtest2(double *data1,double *data2,int n)
{
  int j;
  double nx,ave1,ave2,var1,var2,sd,u,df,cov;
  double tiny=TINY;
  double t=0,z=0;

  avevar(data1,n,&ave1,&var1);
  avevar(data2,n,&ave2,&var2);
  if (var1 < tiny || var2 < tiny) return 0;

  nx = (double)n;
  df = nx-1;
  cov = 0;
  for (j=0; j<n; j++)
    cov += (data1[j]-ave1)*(data2[j]-ave2);
  cov /= df;

  t = 0;
  u = (var1+var2-2.0*cov);
  if (u < tiny) return 0;
  sd = sqrt(u/nx);
  if (sd < tiny) return 0;
  t = (double)(ave1-ave2)/sd;
  z = t2z(t,df);
  if (t < 0) z = -z;
  return z;
}

/* twosample t-test, pooled variance  */
double ttest2(double *data1,double *data2,int n1,int n2)
{
  int j;
  double nx1,nx2,ave1,ave2,sd,s1,s2,df;
  double tiny=TINY;
  double t=0,z=0;

  nx1 = (double)n1;
  nx2 = (double)n2;

  ave1 = ave2 = 0;
  for (j=0; j<n1; j++) ave1 += data1[j];
  for (j=0; j<n2; j++) ave2 += data2[j];
  if (ABS(ave1) < tiny) return 0;
  if (ABS(ave2) < tiny) return 0;
  ave1 /= nx1;
  ave2 /= nx2;

  s1 = 0;
  for (j=0; j<n1; j++) s1 += SQR(data1[j]-ave1);

  s2 = 0;
  for (j=0; j<n2; j++) s2 += SQR(data2[j]-ave2);

  df = nx1 + nx2 - 2.0;
  sd = sqrt((1.0/nx1 + 1.0/nx2) * (s1+s2)/df);
  if (sd < tiny) return 0;

  t = (double)(ave1-ave2)/sd;
  z = t2z(t,df);
  if (t < 0) z = -z;
  return z;
}

/* onesample t-test */
double ttest1(double *data,int n)
{
  double ave=0,var=0,tiny=TINY;
  avevar(data,n,&ave,&var);
  if (var < tiny) return 0;
  double nx = (double)n;
  double df = nx-1.0;
  double t  = sqrt(nx) * ave/sqrt(var);
  double z  = t2z(t,df);
  if (t < 0) z = -z;
  return z;
}

ext: c

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>

#include "tests.h"

void
test_swap (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int N = 1;
   float X[] = { 0.539f };
   int incX = 1;
   float Y[] = { -0.262f };
   int incY = -1;
   float expected1[] = { -0.262f };
   float expected2[] = { 0.539f };
   cblas_sswap(N, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], expected1[i], flteps, "sswap(case 88)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], expected2[i], flteps, "sswap(case 89)");
     }
   };
  };

  {
   int N = 1;
   double X[] = { 0.906 };
   int incX = 1;
   double Y[] = { 0.373 };
   int incY = -1;
   double expected1[] = { 0.373 };
   double expected2[] = { 0.906 };
   cblas_dswap(N, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], expected1[i], dbleps, "dswap(case 90)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], expected2[i], dbleps, "dswap(case 91)");
     }
   };
  };

  {
   int N = 1;
   float X[] = { -0.316f, -0.529f };
   int incX = 1;
   float Y[] = { -0.313f, 0.363f };
   int incY = -1;
   float expected1[] = { -0.313f, 0.363f };
   float expected2[] = { -0.316f, -0.529f };
   cblas_cswap(N, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[2*i], expected1[2*i], flteps, "cswap(case 92) real");
       gsl_test_rel(X[2*i+1], expected1[2*i+1], flteps, "cswap(case 92) imag");
     };
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[2*i], expected2[2*i], flteps, "cswap(case 93) real");
       gsl_test_rel(Y[2*i+1], expected2[2*i+1], flteps, "cswap(case 93) imag");
     };
   };
  };

  {
   int N = 1;
   double X[] = { 0.512, -0.89 };
   int incX = 1;
   double Y[] = { -0.225, -0.511 };
   int incY = -1;
   double expected1[] = { -0.225, -0.511 };
   double expected2[] = { 0.512, -0.89 };
   cblas_zswap(N, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[2*i], expected1[2*i], dbleps, "zswap(case 94) real");
       gsl_test_rel(X[2*i+1], expected1[2*i+1], dbleps, "zswap(case 94) imag");
     };
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[2*i], expected2[2*i], dbleps, "zswap(case 95) real");
       gsl_test_rel(Y[2*i+1], expected2[2*i+1], dbleps, "zswap(case 95) imag");
     };
   };
  };

  {
   int N = 1;
   float X[] = { 0.336f };
   int incX = -1;
   float Y[] = { -0.431f };
   int incY = 1;
   float expected1[] = { -0.431f };
   float expected2[] = { 0.336f };
   cblas_sswap(N, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], expected1[i], flteps, "sswap(case 96)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], expected2[i], flteps, "sswap(case 97)");
     }
   };
  };

  {
   int N = 1;
   double X[] = { 0.764 };
   int incX = -1;
   double Y[] = { -0.293 };
   int incY = 1;
   double expected1[] = { -0.293 };
   double expected2[] = { 0.764 };
   cblas_dswap(N, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], expected1[i], dbleps, "dswap(case 98)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], expected2[i], dbleps, "dswap(case 99)");
     }
   };
  };

  {
   int N = 1;
   float X[] = { -0.239f, 0.361f };
   int incX = -1;
   float Y[] = { 0.149f, 0.347f };
   int incY = 1;
   float expected1[] = { 0.149f, 0.347f };
   float expected2[] = { -0.239f, 0.361f };
   cblas_cswap(N, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[2*i], expected1[2*i], flteps, "cswap(case 100) real");
       gsl_test_rel(X[2*i+1], expected1[2*i+1], flteps, "cswap(case 100) imag");
     };
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[2*i], expected2[2*i], flteps, "cswap(case 101) real");
       gsl_test_rel(Y[2*i+1], expected2[2*i+1], flteps, "cswap(case 101) imag");
     };
   };
  };

  {
   int N = 1;
   double X[] = { -0.171, -0.936 };
   int incX = -1;
   double Y[] = { 0.495, -0.835 };
   int incY = 1;
   double expected1[] = { 0.495, -0.835 };
   double expected2[] = { -0.171, -0.936 };
   cblas_zswap(N, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[2*i], expected1[2*i], dbleps, "zswap(case 102) real");
       gsl_test_rel(X[2*i+1], expected1[2*i+1], dbleps, "zswap(case 102) imag");
     };
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[2*i], expected2[2*i], dbleps, "zswap(case 103) real");
       gsl_test_rel(Y[2*i+1], expected2[2*i+1], dbleps, "zswap(case 103) imag");
     };
   };
  };

  {
   int N = 1;
   float X[] = { -0.405f };
   int incX = -1;
   float Y[] = { -0.213f };
   int incY = -1;
   float expected1[] = { -0.213f };
   float expected2[] = { -0.405f };
   cblas_sswap(N, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], expected1[i], flteps, "sswap(case 104)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], expected2[i], flteps, "sswap(case 105)");
     }
   };
  };

  {
   int N = 1;
   double X[] = { -0.761 };
   int incX = -1;
   double Y[] = { -0.585 };
   int incY = -1;
   double expected1[] = { -0.585 };
   double expected2[] = { -0.761 };
   cblas_dswap(N, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[i], expected1[i], dbleps, "dswap(case 106)");
     }
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], expected2[i], dbleps, "dswap(case 107)");
     }
   };
  };

  {
   int N = 1;
   float X[] = { 0.853f, 0.146f };
   int incX = -1;
   float Y[] = { 0.009f, -0.178f };
   int incY = -1;
   float expected1[] = { 0.009f, -0.178f };
   float expected2[] = { 0.853f, 0.146f };
   cblas_cswap(N, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[2*i], expected1[2*i], flteps, "cswap(case 108) real");
       gsl_test_rel(X[2*i+1], expected1[2*i+1], flteps, "cswap(case 108) imag");
     };
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[2*i], expected2[2*i], flteps, "cswap(case 109) real");
       gsl_test_rel(Y[2*i+1], expected2[2*i+1], flteps, "cswap(case 109) imag");
     };
   };
  };

  {
   int N = 1;
   double X[] = { -0.228, 0.386 };
   int incX = -1;
   double Y[] = { 0.988, -0.084 };
   int incY = -1;
   double expected1[] = { 0.988, -0.084 };
   double expected2[] = { -0.228, 0.386 };
   cblas_zswap(N, X, incX, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(X[2*i], expected1[2*i], dbleps, "zswap(case 110) real");
       gsl_test_rel(X[2*i+1], expected1[2*i+1], dbleps, "zswap(case 110) imag");
     };
   };
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[2*i], expected2[2*i], dbleps, "zswap(case 111) real");
       gsl_test_rel(Y[2*i+1], expected2[2*i+1], dbleps, "zswap(case 111) imag");
     };
   };
  };

}

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
//
// Created by robin on 2018/9/1.
// Copyright (c) 2018 Robin. All rights reserved.
//

#pragma once

#include <ostream>
#include <algorithm>
#include <vector>
#include <gsl/gsl>
#include "arithmetic.h"

namespace satz::gf2::v2 {

/// @brief Polynomial over GF(2).
class Polynomial {
public:
  using int_type = int;
  using container_type = std::vector<int_type>;

  Polynomial () = default;

  template <typename T>
  Polynomial (std::initializer_list<T> degrees)
      : Polynomial(degrees.begin(), degrees.end()) { }

  template <typename InputIt>
  Polynomial (InputIt first, InputIt last) : deg_(first, last) {
    Expects(std::all_of(deg_.begin(),
                        deg_.end(),
                        [] (const auto& x) { return x >= 0; }));

    if (!std::is_sorted(deg_.begin(), deg_.end())) {
      std::sort(deg_.begin(), deg_.end());
    }
    deg_.erase(std::unique(deg_.begin(), deg_.end()), deg_.end());
  }

  static Polynomial from_ulong (unsigned long);
  static Polynomial from_bytes (gsl::span<uint8_t>);
  static Polynomial from_bytes (gsl::span<uint8_t>, int_type degree);
  static Polynomial make_random (int_type degree);
  static Polynomial make_irreducible (int_type degree);

  /**
   * @brief Return a bit representation of the coefficients.
   *
   *    When the degree of this polynomial is -1, return an empty vector.
   *    In other cases, return a bit representation with just enough bytes,
   *    i.e., the most significant byte should be non-zero.
   */
  [[nodiscard]] std::vector<uint8_t> to_bytes () const;

  explicit operator bool () const;

  friend bool operator == (const Polynomial& lhs, const Polynomial& rhs);
  friend bool operator != (const Polynomial& lhs, const Polynomial& rhs);
  friend bool operator < (const Polynomial& lhs, const Polynomial& rhs);
  friend bool operator > (const Polynomial& lhs, const Polynomial& rhs);
  friend bool operator <= (const Polynomial& lhs, const Polynomial& rhs);
  friend bool operator >= (const Polynomial& lhs, const Polynomial& rhs);
  Polynomial& operator ^= (const Polynomial& rhs);
  Polynomial& operator |= (const Polynomial& rhs);
  Polynomial& operator &= (const Polynomial& rhs);
  Polynomial& operator += (const Polynomial& rhs);
  Polynomial& operator -= (const Polynomial& rhs);
  Polynomial& operator *= (const Polynomial& rhs);
  Polynomial& operator %= (const Polynomial& rhs);
  Polynomial& operator <<= (int_type n);
  Polynomial& operator >>= (int_type n);

  friend Polynomial operator ^ (const Polynomial& lhs, const Polynomial& rhs);
  friend Polynomial operator | (const Polynomial& lhs, const Polynomial& rhs);
  friend Polynomial operator & (const Polynomial& lhs, const Polynomial& rhs);
  friend Polynomial operator + (const Polynomial& lhs, const Polynomial& rhs);
  friend Polynomial operator - (const Polynomial& lhs, const Polynomial& rhs);
  friend Polynomial operator * (const Polynomial& lhs, const Polynomial& rhs);
  friend Polynomial operator % (const Polynomial& lhs, const Polynomial& rhs);
  Polynomial operator << (int_type n) const;
  Polynomial operator >> (int_type n) const;

  friend std::ostream& operator << (std::ostream& os, const Polynomial& obj);

  /**
   * @brief Indicates whether it is a zero polynomial.
   */
  [[nodiscard]] bool empty () const;

  /**
   * @brief Returns the degree of the polynomial.
   */
  [[nodiscard]] int_type degree () const;

  /**
   * @brief Returns the number of non-zero coefficients.
   */
  [[nodiscard]] int_type nnz () const;

  /**
   * @brief Given an interger n, return whether this polynomial contains an
   *    x^n term.
   * @param n
   */
  [[nodiscard]] bool contains (int_type n) const;

protected:
  explicit Polynomial (container_type&& v);

private:
  container_type deg_; // `deg_` for degrees
};

}

namespace satz::gf2 {
using namespace v2;
}

namespace satz {

template <>
const gf2::v2::Polynomial zero<gf2::v2::Polynomial> = gf2::v2::Polynomial{};

template <>
const gf2::v2::Polynomial one<gf2::v2::Polynomial> = gf2::v2::Polynomial{0};

template <>
const gf2::v2::Polynomial X<gf2::v2::Polynomial> = gf2::v2::Polynomial{1};

}

max_stars_repo_name: lie-yan/rabin-fingerprint

max_issues_repo_name: lie-yan/rabin-fingerprint

max_forks_repo_name: lie-yan/rabin-fingerprint


text:
#ifndef INCLUDED_TELEPATH_BLAS_2_GEMV_H
#define INCLUDED_TELEPATH_BLAS_2_GEMV_H

#include <cblas.h>
#include <complex>

namespace telepath{

    template< typename T >
    void gemv( CBLAS_ORDER order, CBLAS_TRANSPOSE trans, std::size_t m,
        std::size_t n, T alpha, const T* a, std::size_t lda, const T* x,
        std::size_t incx, T beta, T* y, std::size_t incy );
    template< typename T >
    void gemv( CBLAS_ORDER order, CBLAS_TRANSPOSE trans, std::size_t m,
        std::size_t n, const T* alpha, const T* a, std::size_t lda, const T* x,
        std::size_t incx, const T* beta, T* y, std::size_t incy );

    template<>
    void gemv<float>( CBLAS_ORDER order, CBLAS_TRANSPOSE trans, std::size_t m,
        std::size_t n, float alpha, const float* a, std::size_t lda,
        const float* x, std::size_t incx, float beta, float* y,
        std::size_t incy ){
        return cblas_sgemv( order, trans, m, n, alpha, a, lda, x, incx, beta,
            y, incy );
    }

    template<>
    void gemv<double>( CBLAS_ORDER order, CBLAS_TRANSPOSE trans,
        std::size_t m, std::size_t n, double alpha, const double* a,
        std::size_t lda, const double* x, std::size_t incx, double beta,
        double* y, std::size_t incy ){
        return cblas_dgemv( order, trans, m, n, alpha, a, lda, x, incx, beta,
            y, incy );
    }

    template<>
    void gemv<std::complex<float> >( CBLAS_ORDER order, CBLAS_TRANSPOSE trans,
        std::size_t m, std::size_t n, const std::complex<float>* alpha,
        const std::complex<float>* a, std::size_t lda,
        const std::complex<float>* x, std::size_t incx,
        const std::complex<float>* beta, std::complex<float>* y,
        std::size_t incy ){
        return cblas_cgemv( order, trans, m, n, (float*) alpha, (float*) a, lda,
            (float*) x, incx, (float*) beta, (float*) y, incy );
    }

    template<>
    void gemv<std::complex<double> >( CBLAS_ORDER order,
        CBLAS_TRANSPOSE trans, std::size_t m, std::size_t n,
        const std::complex<double>* alpha, const std::complex<double>* a,
        std::size_t lda, const std::complex<double>* x, std::size_t incx,
        const std::complex<double>* beta, std::complex<double>* y,
        std::size_t incy ){
        return cblas_zgemv( order, trans, m, n, (double*) alpha, (double*) a,
            lda, (double*) x, incx, (double*) beta, (double*) y, incy );
    }

} //namespace telepath

#endif

max_stars_repo_name: tbepler/telepath

max_issues_repo_name: tbepler/telepath

max_forks_repo_name: tbepler/telepath