text:
/*
 * testnonsymm.c
 * Patrick Alken
 *
 * Compile: gcc -g -O2 -Wall -o testnonsymm testnonsymm.c -lm -lgsl -lcblas -latlas
 *
 * Usage: testnonsymm [options]
 *
 * -i             : incremental matrices
 * -b             : balance the matrices
 * -z             : compute Schur vectors and test them
 * -n size        : size of matrices
 * -l lower-bound : lower bound for matrix elements
 * -u upper-bound : upper bound for matrix elements
 * -c num         : number of matrices to solve
 */

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <errno.h>
#include <math.h>
#include <getopt.h>
#include <sys/times.h>

#include <gsl/gsl_eigen.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_test.h>

typedef struct
{
  gsl_eigen_nonsymm_workspace *nonsymm_p;
  gsl_eigen_nonsymmv_workspace *nonsymmv_p;
  gsl_matrix *A;
  gsl_matrix *Av;
  gsl_vector_complex *eval;
  gsl_vector_complex *evalv;
  gsl_matrix_complex *evec;

  int compute_z;
  gsl_matrix *Z;
  gsl_matrix *Zv;

  size_t n_evals;
} nonsymm_workspace;

nonsymm_workspace *nonsymm_alloc(size_t n, int compute_z, int do_balance);
void nonsymm_free(nonsymm_workspace *w);
int nonsymm_proc(nonsymm_workspace *w);

/*
 * Global variables
 */
unsigned long count = 0;

/*
 * Prototypes
 */

void make_random_matrix(gsl_matrix *m, gsl_rng *r, int lower,
                        int upper);
void make_start_matrix(gsl_matrix *m, int lower);
int inc_matrix (gsl_matrix *m, int lower, int upper);
void output_matrix(gsl_matrix *m);
void print_octave(gsl_matrix *m, const char *str);
void print_matrix(gsl_matrix *m, const char *str);
void print_hess(gsl_matrix *H, const char *str);
void print_vector(gsl_vector_complex *eval, const char *str);
int cmp(double a, double b);
int compare(const void *a, const void *b);
void sort_complex_vector(gsl_vector_complex *v);
int test_evals(gsl_vector_complex *obs, gsl_vector_complex *expected,
               gsl_matrix *A, const char *obsname, const char *expname);
int test_Z(gsl_matrix *A, gsl_matrix *Z, gsl_matrix *T);
int test_eigenvectors(gsl_matrix *A, gsl_vector_complex *eval,
                      gsl_matrix_complex *evec);
void my_error_handler(const char *reason, const char *file, int line,
                      int err);

nonsymm_workspace *
nonsymm_alloc(size_t n, int compute_z, int do_balance)

{
  nonsymm_workspace *w;

  w = (nonsymm_workspace *) malloc(sizeof(nonsymm_workspace));

  memset(w, '\0', sizeof(nonsymm_workspace));

  w->nonsymm_p = gsl_eigen_nonsymm_alloc(n);
  w->nonsymmv_p = gsl_eigen_nonsymmv_alloc(n);

  w->A = gsl_matrix_alloc(n, n);
  w->Av = gsl_matrix_alloc(n, n);

  if (compute_z)
    {
      w->Z = gsl_matrix_alloc(n, n);
      w->Zv = gsl_matrix_alloc(n, n);
      w->compute_z = 1;
      gsl_eigen_nonsymm_params(1, do_balance, w->nonsymm_p);
    }
  else
    gsl_eigen_nonsymm_params(0, do_balance, w->nonsymm_p);

  w->eval = gsl_vector_complex_alloc(n);
  w->evalv = gsl_vector_complex_alloc(n);

  w->evec = gsl_matrix_complex_alloc(n, n);

  return (w);
} /* nonsymm_alloc() */

void
nonsymm_free(nonsymm_workspace *w)

{
  if (w->nonsymm_p)
    gsl_eigen_nonsymm_free(w->nonsymm_p);

  if (w->nonsymmv_p)
    gsl_eigen_nonsymmv_free(w->nonsymmv_p);

  if (w->A)
    gsl_matrix_free(w->A);

  if (w->Av)
    gsl_matrix_free(w->Av);

  if (w->Z)
    gsl_matrix_free(w->Z);

  if (w->Zv)
    gsl_matrix_free(w->Zv);

  if (w->eval)
    gsl_vector_complex_free(w->eval);

  if (w->evalv)
    gsl_vector_complex_free(w->evalv);

  if (w->evec)
    gsl_matrix_complex_free(w->evec);

  free(w);
} /* nonsymm_free() */

int
nonsymm_proc(nonsymm_workspace *w)

{
  int s1, s2, s;

  if (w->compute_z)
    {
      s1 = gsl_eigen_nonsymm_Z(w->A, w->eval, w->Z, w->nonsymm_p);
      s2 = gsl_eigen_nonsymmv_Z(w->Av, w->evalv, w->evec, w->Zv, w->nonsymmv_p);
    }
  else
    {
      s1 = gsl_eigen_nonsymm(w->A, w->eval, w->nonsymm_p);
      s2 = gsl_eigen_nonsymmv(w->Av, w->evalv, w->evec, w->nonsymmv_p);
    }

  w->n_evals = w->nonsymm_p->n_evals;

  s = 0;
  if (s1)
    s = s1;
  else if (s2)
    s = s2;

  return s;
} /* nonsymm_proc() */

/**********************************************
 * General routines
 **********************************************/

void
make_random_matrix(gsl_matrix *m, gsl_rng *r, int lower, int upper)

{
  size_t i, j;
  size_t N = m->size1;

  for (i = 0; i < N; ++i)
  {
    for (j = 0; j < N; ++j)
    {
      gsl_matrix_set(m,
                     i,
                     j,
                     gsl_rng_uniform(r) * (upper - lower) + lower);
    }
  }
} /* make_random_matrix() */

void
make_start_matrix(gsl_matrix *m, int lower)

{
  size_t i, j;
  size_t N = m->size1;

  for (i = 0; i < N; ++i)
    for (j = 0; j < N; ++j)
      gsl_matrix_set(m, i, j, (double)lower);
} /* make_start_matrix() */

int
inc_matrix (gsl_matrix *m, int lower, int upper)
{
  size_t i = 0;
  size_t N = m->size1 * m->size2;
  int carry = 1;

  for (i = 0; carry > 0 && i < N; i++)
    {
      double v = m->data[i] + carry;
      carry = (v > upper) ? 1 : 0;
      m->data[i] = (v > upper) ? lower : v;
    }

  return carry;
} /* inc_matrix() */

void
output_matrix(gsl_matrix *m)
{
  size_t i, j;
  size_t N = m->size1;
  size_t M = m->size2;

  for (i = 0; i < N; ++i)
  {
    for (j = 0; j < M; ++j)
      {
        /*printf("%10.9f%s",*/
        printf("%10.18e%s",
               gsl_matrix_get(m, i, j),
               (j < M - 1) ? "," : ";\n");
      }
  }
}

void
print_octave(gsl_matrix *m, const char *str)
{
  FILE *fp;
  size_t i, j;
  const size_t N = m->size1;
  const size_t M = m->size2;

  fp = fopen(str, "w");

  if (!fp)
    return;

  fprintf(fp, "# Created by Octave 2.1.73, Tue Aug 01 15:00:27 2006 MDT <blah@blah>\n");
  fprintf(fp, "# name: %s\n", str);
  fprintf(fp, "# type: matrix\n");
  fprintf(fp, "# rows: %u\n", N);
  fprintf(fp, "# columns: %u\n", N);

  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < M; ++j)
        {
          fprintf(fp,
                  "%10.9f%s",
                  gsl_matrix_get(m, i, j),
                  (j < M - 1) ? " " : "\n");
        }
    }

  fclose(fp);
}

void
print_matrix(gsl_matrix *m, const char *str)

{
  size_t i, j;
  size_t N = m->size1;
  size_t M = m->size2;
  gsl_matrix_view v;
  size_t rows, cols;
  size_t r, c;
  char buf[100];

  /*print_octave(m, str);
  return;*/

  /*rows = GSL_MIN(15, N);*/
  rows = N;
  cols = GSL_MIN(15, N);
  /*cols = N;*/

  for (i = 0; i < N; i += rows)
  {
    for (j = 0; j < M; j += cols)
    {
      r = GSL_MIN(rows, N - i);
      c = GSL_MIN(cols, N - j);

      v = gsl_matrix_submatrix(m,
                               i,
                               j,
                               r,
                               c);

      sprintf(buf, "%s(%u:%u,%u:%u)",
              str,
              i + 1,
              i + r,
              j + 1,
              j + c);

      printf("%s = [\n", buf);

      output_matrix(&v.matrix);

      printf("]\n");
    }
  }
} /* print_matrix() */

void
print_vector(gsl_vector_complex *eval, const char *str)

{
  size_t N = eval->size;
  size_t i;
  gsl_complex z;

  printf("%s = [\n", str);

  for (i = 0; i < N; ++i)
    {
      z = gsl_vector_complex_get(eval, i);
      printf("%.18e %.18e;\n", GSL_REAL(z), GSL_IMAG(z));
    }

  printf("\n]\n");
} /* print_vector() */

int
cmp(double a, double b)

{
  return ((a > b) ? 1 : ((a < b) ? -1 : 0));
} /* cmp() */

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

{
  const double *x = a;
  const double *y = b;
  int r1 = cmp(y[0], x[0]);
  int r2 = cmp(y[1], x[1]);

  if (fabs(x[0] - y[0]) < 1.0e-8)
    {
      /* real parts are very close to each other */
      return r2;
    }
  else
    {
      return r1 ? r1 : r2;
    }
} /* compare() */

void
sort_complex_vector(gsl_vector_complex *v)

{
  qsort(v->data, v->size, 2 * sizeof(double), &compare);
} /* sort_complex_vector() */

int
test_evals(gsl_vector_complex *obs, gsl_vector_complex *expected,
           gsl_matrix *A, const char *obsname, const char *expname)

{
  size_t N = expected->size;
  size_t i, k;
  double max, max_abserr, max_relerr;

  max = 0.0;
  max_abserr = 0.0;
  max_relerr = 0.0;
  k = 0;

  for (i = 0; i < N; ++i)
    {
      gsl_complex z = gsl_vector_complex_get(expected, i);
      max = GSL_MAX_DBL(max, gsl_complex_abs(z));
    }

  for (i = 0; i < N; ++i)
    {
      gsl_complex z_obs = gsl_vector_complex_get(obs, i);
      gsl_complex z_exp = gsl_vector_complex_get(expected, i);

      double x_obs = GSL_REAL(z_obs);
      double y_obs = GSL_IMAG(z_obs);
      double x_exp = GSL_REAL(z_exp);
      double y_exp = GSL_IMAG(z_exp);

      double abserr_x = fabs(x_obs - x_exp);
      double abserr_y = fabs(y_obs - y_exp);
      double noise = max * GSL_DBL_EPSILON * N * N;

      max_abserr = GSL_MAX_DBL(max_abserr, abserr_x + abserr_y);

      if (abserr_x < noise && abserr_y < noise)
        continue;

      if (abserr_x > 1.0e-6 || abserr_y > 1.0e-6)
        ++k;
    }

    if (k)
      {
        printf("==== CASE %lu ===========================\n\n", count);

        print_matrix(A, "A");

        printf("=== eval - %s ===\n", expname);
        print_vector(expected, expname);

        printf("=== eval - %s ===\n", obsname);
        print_vector(obs, obsname);

        printf("max abserr = %g  max relerr = %g\n", max_abserr, max_relerr);

        printf("=========================================\n\n");
      }

    return k;
} /* test_evals() */

/* test if A = ZTZ^t (or AZ = ZT) */

int
test_Z(gsl_matrix *A, gsl_matrix *Z, gsl_matrix *T)

{
  size_t N = A->size1;
  size_t i, j, k;
  double rhs, lhs;
  double abserr;
  gsl_matrix *T1, *T2;

  T1 = gsl_matrix_alloc(N, N);
  T2 = gsl_matrix_alloc(N, N);

  /* zero the lower triangle of T */
  if (N > 3)
    {
      for (i = 2; i < N; ++i)
        {
          for (j = 0; j < (i - 1); ++j)
            {
              gsl_matrix_set(T, i, j, 0.0);
            }
        }
    }
  else if (N == 3)
    {
      gsl_matrix_set(T, 2, 0, 0.0);
    }

  /* compute T1 = A Z */
  gsl_blas_dgemm(CblasNoTrans,
                 CblasNoTrans,
                 1.0,
                 A,
                 Z,
                 0.0,
                 T1);

  /* compute T2 = Z T */
  gsl_blas_dgemm(CblasNoTrans,
                 CblasNoTrans,
                 1.0,
                 Z,
                 T,
                 0.0,
                 T2);

  k = 0;
  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          lhs = gsl_matrix_get(T1, i, j);
          rhs = gsl_matrix_get(T2, i, j);

          abserr = fabs(lhs - rhs);

          if (abserr > 1.0e-6)
            ++k;
        }
    }

  if (k)
    {
      printf("==== CASE %lu ===========================\n\n", count);

      print_matrix(A, "A");

      printf("=== Schur Form matrix ===\n");
      print_matrix(T, "T");

      printf("=== Similarity matrix ===\n");
      print_matrix(Z, "Z");

      printf("=== A Z ===\n");
      print_matrix(T1, "A Z");

      printf("=== Z T ===\n");
      print_matrix(T2, "Z T");

      printf("=== A Z - Z T ===\n");
      gsl_matrix_sub(T1, T2);
      print_matrix(T1, "A Z - Z T");

      printf("=========================================\n\n");
    }

  gsl_matrix_free(T1);
  gsl_matrix_free(T2);

  return k;
} /* test_Z() */

int
test_eigenvectors(gsl_matrix *A, gsl_vector_complex *eval,
                  gsl_matrix_complex *evec)
{
  const size_t N = A->size1;
  size_t i, j;
  int k, s;
  gsl_matrix_complex *m;
  gsl_vector_complex *x, *y;
  gsl_complex z_one, z_zero;

  m = gsl_matrix_complex_alloc(N, N);
  y = gsl_vector_complex_alloc(N);
  x = gsl_vector_complex_alloc(N);

  /* m <- A */
  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          gsl_complex z;

          GSL_SET_COMPLEX(&z, gsl_matrix_get(A, i, j), 0.0);
          gsl_matrix_complex_set(m, i, j, z);
        }
    }

  GSL_SET_COMPLEX(&z_one, 1.0, 0.0);
  GSL_SET_COMPLEX(&z_zero, 0.0, 0.0);

  s = 0;

  /* check eigenvalues */
  for (i = 0; i < N; ++i)
    {
      gsl_vector_complex_view vi = gsl_matrix_complex_column(evec, i);
      gsl_complex ei = gsl_vector_complex_get(eval, i);
      double norm = gsl_blas_dznrm2(&vi.vector);

      /* check that eigenvector is normalized */
      gsl_test_rel(norm, 1.0, N * GSL_DBL_EPSILON, "case %u, normalized",
                   count);

      /* compute x = lambda * v */
      gsl_vector_complex_memcpy(x, &vi.vector);
      gsl_blas_zscal(ei, x);

      /* compute y = A v */
      gsl_blas_zgemv(CblasNoTrans, z_one, m, &vi.vector, z_zero, y);

      k = 0;

      /* now test if y = x */
      for (j = 0; j < N; ++j)
        {
          gsl_complex z;
          double lhs_r, lhs_i;
          double rhs_r, rhs_i;

          z = gsl_vector_complex_get(y, j);
          lhs_r = GSL_REAL(z);
          lhs_i = GSL_IMAG(z);

          z = gsl_vector_complex_get(x, j);
          rhs_r = GSL_REAL(z);
          rhs_i = GSL_IMAG(z);

          if (fabs(lhs_r - rhs_r) > 1e8 * GSL_DBL_EPSILON)
            ++k;
          if (fabs(lhs_i - rhs_i) > 1e8 * GSL_DBL_EPSILON)
            ++k;
        }

      if (k)
        {
          s++;

          printf("==== CASE %lu ===========================\n\n", count);

          print_matrix(A, "A");

          printf("eval = %f + %fi\n", GSL_REAL(ei), GSL_IMAG(ei));

          print_vector(&vi.vector, "v");

          print_vector(y, "Av");

          print_vector(x, "ev_v");

          printf("=========================================\n\n");
        }
    }

  gsl_matrix_complex_free(m);
  gsl_vector_complex_free(y);
  gsl_vector_complex_free(x);

  return s;
} /* test_eigenvectors() */

void
my_error_handler(const char *reason, const char *file, int line,
                 int err)

{
  printf("[caught: %s:%d: errno=%d %s]\n", file, line, err, reason);
} /* my_error_handler() */

int
main(int argc, char *argv[])

{
  nonsymm_workspace *nonsymm_workspace_p;
  size_t N;
  int c;
  gsl_matrix *A;
  gsl_rng *r;
  int incremental; /* incremental/random matrices */
  int lower;  /* lower bound */
  int upper;  /* upper bound */
  unsigned int nmat;   /* number of matrices to solve */
  int status;
  int compute_z;
  int do_balance;

  gsl_ieee_env_setup();
  gsl_rng_env_setup();

  /*gsl_set_error_handler(&my_error_handler);*/

  N = 30;
  incremental = 0;
  compute_z = 0;
  do_balance = 0;
  lower = -10;
  upper = 10;
  nmat = 0;

  while ((c = getopt(argc, argv, "izbc:n:l:u:t:")) != (-1))
    {
      switch (c)
        {
          case 'i':
            incremental = 1;
            break;

          case 'n':
            N = strtol(optarg, NULL, 0);
            break;

          case 'l':
            lower = strtol(optarg, NULL, 0);
            break;

          case 'u':
            upper = strtol(optarg, NULL, 0);
            break;

          case 'c':
            nmat = strtoul(optarg, NULL, 0);
            break;

          case 'b':
            do_balance = 1;
            break;

          case 'z':
            compute_z = 1;
            break;

          case '?':
          default:
            printf("usage: %s [-i] [-b] [-z] [-n size] [-l lower-bound] [-u upper-bound] [-c num]\n", argv[0]);
            exit(1);
            break;
        } /* switch (c) */
    }

  A = gsl_matrix_alloc(N, N);
  nonsymm_workspace_p = nonsymm_alloc(N, compute_z, do_balance);

  if (!incremental)
    r = gsl_rng_alloc(gsl_rng_default);
  else
  {
    r = 0;
    make_start_matrix(A, lower);
  }

  fprintf(stderr, "testing N = %d", N);
  if (incremental)
    fprintf(stderr, " incrementally");
  else
    fprintf(stderr, " randomly");

  fprintf(stderr, " on element range [%d, %d]", lower, upper);

  if (compute_z)
    fprintf(stderr, ", with Schur vectors");

  if (do_balance)
    fprintf(stderr, ", with balancing");

  fprintf(stderr, "\n");

  while (1)
    {
      if (nmat && (count >= nmat))
        break;

      ++count;

      if (!incremental)
        make_random_matrix(A, r, lower, upper);
      else
        {
          status = inc_matrix(A, lower, upper);
          if (status)
            break; /* all done */
        }

      /* make copies of the matrix */
      gsl_matrix_memcpy(nonsymm_workspace_p->A, A);
      gsl_matrix_memcpy(nonsymm_workspace_p->Av, A);

      status = nonsymm_proc(nonsymm_workspace_p);

      if (status)
        {
          printf("=========== CASE %lu ============\n", count);
          printf("Failed to converge: found %u eigenvalues\n",
                 nonsymm_workspace_p->n_evals);
          print_matrix(A, "A");
        }

      gsl_eigen_nonsymmv_sort(nonsymm_workspace_p->evalv,
                              nonsymm_workspace_p->evec,
                              GSL_EIGEN_SORT_ABS_ASC);

      status = test_eigenvectors(A,
                                 nonsymm_workspace_p->evalv,
                                 nonsymm_workspace_p->evec);

      sort_complex_vector(nonsymm_workspace_p->eval);
      sort_complex_vector(nonsymm_workspace_p->evalv);

      status = test_evals(nonsymm_workspace_p->eval,
                          nonsymm_workspace_p->evalv,
                          A,
                          "nonsymm",
                          "nonsymmv");

      if (compute_z)
        {
          test_Z(A, nonsymm_workspace_p->Z, nonsymm_workspace_p->A);
          test_Z(A, nonsymm_workspace_p->Zv, nonsymm_workspace_p->Av);
        }
    }

  gsl_matrix_free(A);
  nonsymm_free(nonsymm_workspace_p);

  if (r)
    gsl_rng_free(r);

  return 0;
} /* main() */

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#ifndef BASIC_DATA_STRUCTURE_H
#define BASIC_DATA_STRUCTURE_H

#include <vector>
#include <array>
#include <string>
#include <fstream>
#include <iostream>
#include <petsc.h>
#include <petscksp.h>

#include "petscsys.h"
#include "petscmat.h"
using namespace std;
/////////////////////////////////

// const int degree = 3;
// const int dim = 3;
const int phydim = 3;
// const int bzpt_num = (dim == 3) ? 64 : 16;
// const int dof_all = 3;

class Vertex3D
{
public:
	double coor[3];
	int label; //ID for inlet and outlet
	Vertex3D();
};

class Element3D
{
public:
	int degree;
	int order;
	int nbf;
	int type;	//0 for interior and 1 for boundary, for visualization purpose
	int bzflag; //0 for spline element, 1 for Bezier element
	int label;

	vector<int> IEN;
	vector<array<double, 3>> pts; //tmp
	vector<array<double, 64>> cmat;
	vector<int> IDBC;
	double velocity[3];

	Element3D(int p = 3);
	void BezierPolyn(double u, vector<double> &Nu, vector<double> &dNdu) const;
	void Basis(double u, double v, double w, vector<double> &Nt, vector<array<double, 3>> &dNdt) const;
	void Para2Phys(double u, double v, double w, double pt[3]) const;
};

class Vertex2D
{
public:
	double coor[3];
	int label; //ID for inlet and outlet
	Vertex2D();
};

class Element2D
{
public:
	int degree;
	int order;
	int nbf;
	int type;	//0 for interior and 1 for boundary, for visualization purpose
	int bzflag; //0 for spline element, 1 for Bezier element
	int label;

	vector<int> IEN;
	vector<array<double, phydim>> pts; //tmp
	vector<array<double, 16>> cmat;
	vector<int> IDBC;
	double velocity[2];

	Element2D(int p = 3);
	void BezierPolyn(double u, vector<double> &Nu, vector<double> &dNdu) const;
	void Basis(double u, double v, vector<double> &Nt, vector<array<double, 2>> &dNdt) const;
	void Para2Phys(double u, double v, double pt[phydim]) const;
};

//mesh format conversion

void Raw2Vtk_hex(string fn);

void Rawn2Vtk_hex(string fn);

#endif

ext: h

max_stars_repo_name: CMU-CBML/NeuronTransportOptimization

max_issues_repo_name: CMU-CBML/NeuronTransportOptimization

max_forks_repo_name: CMU-CBML/NeuronTransportOptimization


text:
/*  Copyright 2017 International Business Machines Corporation

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.  */

#pragma once

#include "kernelpp/config.h"

#include <gsl.h>
#include <memory>

#if defined(kernelpp_WITH_CUDA)
#   define checkCudaErrors(err)  __checkCudaErrors (err, __FILE__, __LINE__)
#   define checkCudaLastError()  __checkLastCudaError (__FILE__, __LINE__)
#endif

namespace kernelpp
{
    /*  Initialize the cuda runtime, and ensure at least
     *  one device is available.
     */
    bool init_cudart();

    /*  `device_ptr<T>` is a smart pointer which owns and manages one
     *  or more objects of type `T` in contiguous memory on a CUDA device,
     *  and disposes of these objects when it goes out of scope.
     */
    template<typename T> class device_ptr final
    {
        struct cuda_deleter;
        std::unique_ptr<T, cuda_deleter> m_ptr;
        size_t m_size;

      public:
        device_ptr();
        device_ptr(size_t n);
        device_ptr(device_ptr<T>&& other);
        device_ptr(gsl::span<T> span);

        gsl::span<T> span();
        const gsl::span<T> span() const;

        T* get();

        void copy_from(const gsl::span<T> from);
        void copy_to(gsl::span<T> to) const;
    };
}

#if defined(kernelpp_WITH_CUDA) && defined(__CUDACC__)
#   include "cuda_util-inl.cuh"
#endif

max_stars_repo_name: rayglover-ibm/kernelpp

max_issues_repo_name: rayglover-ibm/kernelpp

max_forks_repo_name: rayglover-ibm/kernelpp


text:
#pragma once

#include <gsl/span>
#include <rapidjson/fwd.h>
#include <spdlog/fwd.h>

#include <memory>

namespace CesiumGltf {
struct Model;
}

namespace Cesium3DTilesSelection {

/**
 * @brief Parses the provided B3DM batch table and adds an equivalent
 * EXT_feature_metadata extension to the provided glTF.
 *
 * @param pLogger
 * @param gltf
 * @param featureTable
 * @param batchTableJson
 * @param batchTableBinaryData
 */
void upgradeBatchTableToFeatureMetadata(
    const std::shared_ptr<spdlog::logger>& pLogger,
    CesiumGltf::Model& gltf,
    const rapidjson::Document& featureTable,
    const rapidjson::Document& batchTableJson,
    const gsl::span<const std::byte>& batchTableBinaryData);

} // namespace Cesium3DTilesSelection

max_stars_repo_name: yieryi/cesium-native

max_issues_repo_name: yieryi/cesium-native

max_forks_repo_name: yieryi/cesium-native


text:
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include "mathfunctions.h"
#include "globals.h"
#include "constitutive_equations.h"

/*****************************************************************************/
/***************************************************************************/
/************ Compute the spatial discretization of the equation ***********/
/***************************************************************************/

extern const double g;
extern void Compute_InnerProducts(double IP1[], double IP2[], double IP3[], double IP4[], double IP5[], double IP6[], double IP7[], double IP8[], double IP9[], double IP10[]);
extern void RoeFlux( double *Fhat, double A_L, double A_R, double Q_L, double Q_R, double b, double m1val, double m2val, double g);
extern void LF(double *Fhat, double A_L, double A_R, double Q_L, double Q_R, double b, double g);

extern void* xcalloc(int items, int size);
extern void GetLGLWeights(int N, double *w);

void computeL()
{
	double* Fhat1L = xcalloc(NumNodes, sizeof(double));
	double* Fhat2L = xcalloc(NumNodes, sizeof(double));
	double* Fhat1R = xcalloc(NumNodes, sizeof(double));
	double* Fhat2R = xcalloc(NumNodes, sizeof(double));

#ifdef WDON
/**** Compute the mass over an element to use later for ensuring ********/
/**** positive mass with the flux in a wetting and drying treatment *************/
	double mass[NumEl];
	double LGLWeight[Np];
	GetLGLWeights(P, LGLWeight);

	for (int i = 0; i < NumEl; ++i)
	{
		double avgArea = 0;
		int begNode = i*Np + 1;
		for (int j =0; j < Np; j++)
		{
			avgArea += LGLWeight[j]*A[begNode+1];
		}
		mass[i] = avgArea;
	}
#endif

/************** Compute the numerical flux at the faces *****************/
	for (int i=0; i < NumEdges; ++i)
	{
		double A_L, Q_L, A_R, Q_R, m1val, m2val;

		int leftNode = i*Np;
		int rightNode = leftNode+1;

		m1val = m1[i];
		m2val = m2[i];

		double tmpF[2];

		#ifdef WDON
		// Check to see if the elements separated by this boundary are both dry
		if (i > 0 && i < NumNodes-1)
		{
			if (WD[i-1] == 0 && WD[i] == 0)
			{
				// Reflection flux for the left element
				A_L = A[leftNode];
				A_R = A_L;
				Q_L = Q[leftNode];
				Q_R = -Q_L;
				//LF(tmpF, A_L, A_R, Q_L, Q_R,b[i], m1val, m2val, 0);
				RoeFlux(tmpF, A_L, A_R, Q_L, Q_R,b[i], m1val, m2val,0);
				Fhat1L[i] = tmpF[0];
				Fhat2L[i] = tmpF[1];

				// Reflection flux for the right element
				A_R = A[rightNode];
				A_L = A_R;
				Q_R = Q[rightNode];
				Q_L = -Q_R;
				//LF(tmpF, A_L, A_R, Q_L, Q_R, b[i], m1val, m2val, 0);
				RoeFlux(tmpF, A_L, A_R, Q_L, Q_R, b[i], m1val, m2val, 0);
				Fhat1R[i] = tmpF[0];
				Fhat2R[i] = tmpF[1];

				if (isnan(tmpF[0]) || isnan(tmpF[1]))
				{
					printf("both elements dry flux not a number, edge %d \n",i);
					exit(EXIT_FAILURE);
				}

			}

		}

		// if the elements are not both dry
		if ((i == 0) || (i == NumEl) || WD[i-1] == 1 || WD[i] == 1)
		{
			A_L = A[leftNode];
			A_R = A[rightNode];

			Q_L = Q[leftNode];
			Q_R = Q[rightNode];

			RoeFlux(tmpF, A_L, A_R, Q_L,Q_R,b[i], m1val, m2val, g);
			//LF(tmpF, A_L, A_R, Q_L,Q_R,b[i], m1val, m2val, g);

			if (isnan(tmpF[0]) || isnan(tmpF[1]))
			{
				printf(" both elements wet flux not a number, edge %d \n",i);
				exit(EXIT_FAILURE);
			}

			if (i==0 || WD[i-1] == 1)
			{
				Fhat1L[i] = tmpF[0];
				Fhat2L[i] = tmpF[1];

			}
			else
			{
				double newtmpF[2];
				//LF(newtmpF, A_L, A_R, Q_L, Q_R, b[i], m1val, m2val, 0);
				RoeFlux(newtmpF, A_L, A_R, Q_L, Q_R, b[i], m1val, m2val, 0);
				Fhat1L[i] = newtmpF[0];
				Fhat2L[i] = newtmpF[1];
				if (isnan(newtmpF[0]) || isnan(newtmpF[1]))
				{
					printf("left element dry flux not a number, edge %d \n",i);
					exit(EXIT_FAILURE);
				}

			}
			if (i == NumEl || WD[i]==1)
			{
				Fhat1R[i] = tmpF[0];
				Fhat2R[i] = tmpF[1];
			}
			else
			{
				//LF(tmpF, A_L,A_R,Q_L,Q_R,b[i], m1val, m2val, 0);
				RoeFlux(tmpF, A_L,A_R,Q_L,Q_R,b[i], m1val, m2val, 0);
				Fhat1R[i] = tmpF[0];
				Fhat2R[i] = tmpF[1];
				if (isnan(tmpF[0]) || isnan(tmpF[1]))
				{
					printf("right element dry flux not a number, edge %d \n",i);
					exit(EXIT_FAILURE);
				}

			}
		}

		#else
		A_L = A[leftNode];
		A_R = A[rightNode];

		Q_L = Q[leftNode];
		Q_R = Q[rightNode];

		RoeFlux(tmpF, A_L, A_R, Q_L,Q_R,b[i],m1val,m2val,g);
		//LF(tmpF, A_L, A_R, Q_L,Q_R,b[i],m1val,m2val,g);

		if (isnan(tmpF[0]) || isnan(tmpF[1]))
		{
			printf("flux not a number, edge %d \n",i);
			printf("A_L = %lf A_R = %lf Q_L = %lf Q_R = %lf, b = %lf, m1 = %lf, m2 = %lf \n", A_L, A_R, Q_L, Q_R, b[i], m1val, m2val);
			exit(EXIT_FAILURE);
		}

		Fhat1L[i] = tmpF[0];
		Fhat2L[i] = tmpF[1];

		Fhat1R[i] = tmpF[0];
		Fhat2R[i] = tmpF[1];

		#endif

		/**************************************************************************************************/
		double u_L = Q_L/A_L;
		double u_R = Q_R/A_R;
		double c_L = sqrt(g*A_L/b[i]);
		double c_R = sqrt(g*A_R/b[i]);

		if (i==0)
			max_lambda = fmax((fabs(u_L)+c_L), (fabs(u_R)+c_R));
		// Compute maximum eigenvalue for the next time step
		double current_max =fmax((fabs(u_L) + c_L), (fabs(u_R) + c_R));
		max_lambda = max((max_lambda),(current_max));
	}

	printf("max_lambda = %lf\n", max_lambda);
	#ifdef WDON
	for (int i =1; i < NumNodes-1; ++i)
	{
		int leftNode = i*Np;
		int rightNode = leftNode+1;
		double maxBetaOverAlpha = 2;
		// Check to see if this flux might possibly result in negative mass
		// If the mass will be negative on the left side
	      	if (Fhat1L[i]*maxBetaOverAlpha*dt > mass[i-1])
		{
	    double A_L = A[leftNode];
			double A_R = A_L;
			double Q_L = Q[leftNode];
	    double Q_R = -Q_L;
	    double tmpF[2];
			double localG = g;

			if (WD[i] == 0)
				localG = 0;

			//LF(tmpF, A_L, A_R, Q_L, Q_R, b[i],m1val,m2val,localG);
			RoeFlux(tmpF, A_L, A_R, Q_L, Q_R, b[i],m1val,m2val,localG);
	      		Fhat1L[i] = tmpF[0];
	      		Fhat2L[i] = tmpF[1];
			printf("element %d will be dry\n",i);
	    	}

		// If the mass will be negative on the right side
		if (-Fhat1R[i]*maxBetaOverAlpha*dt > mass[i])
		{
			double A_R = A[rightNode];
			double A_L = A_R;
			double Q_R = Q[rightNode];
			double Q_L = -Q_R;
			double tmpF[2];
			double localG = g;
			if (WD[i] == 0)
				localG = 0;
			//LF(tmpF, A_L, A_R, Q_L, Q_R, b[i],m1val,m2val,localG);
			RoeFlux(tmpF, A_L, A_R, Q_L, Q_R, b[i],m1val,m2val,localG);
			Fhat1R[i] = tmpF[0];
			Fhat2R[i] = tmpF[1];
			printf("element %d will be dry \n",i);

		}
	}
	#endif

	// Compute the right hand side
//	double IP1[NumEl], IP2[NumEl], IP3[NumEl], IP4[NumEl], IP5[NumEl], IP6[NumEl], IP7[NumEl], IP8[NumEl], IP9[NumEl], IP10[NumEl];
//	Compute_InnerProducts(IP1, IP2, IP3, IP4, IP5, IP6, IP7, IP8, IP9, IP10);

	int k;
	for (k=0; k < NumEl; ++k)
	{

		double h = dh[k];
		gsl_vector *F1 = gsl_vector_alloc(Np);
		gsl_vector *F2 = gsl_vector_alloc(Np);
		gsl_vector *ST21 = gsl_vector_alloc(Np);
		gsl_vector *ST22 = gsl_vector_alloc(Np);
		gsl_vector *ST23 = gsl_vector_alloc(Np);

		int begNode = k*Np;
		for (int i = 0; i < Np; i++)
		{
			double Aval = A[begNode+i+1];
			double Qval = Q[begNode+i+1];
			double bval = NodalB[begNode+i];
			double S0 = dz[begNode+i];
			double m1val = Nodalm1[begNode+i];
			double m2val = Nodalm2[begNode+i];
			double dm1val = dm1[begNode+i];
			double dm2val = dm2[begNode+i];
			double hval = getH(Aval, bval, m1val, m2val);
			double I1val = getI1(Aval, bval, m1val, m2val);
			double I2val = getI2(Aval, bval, db[k*Np+i], m1val, dm1val, m2val, dm2val);
			double Sfval = getS_f(Aval, Qval, bval, m1val, m2val, NodalnFriction[begNode+i]);
			double beta = 1.0;
			double localG;
			#ifdef WDON
			if(WD[k]==1)
				localG = g;
			else
				localG = 0;
			#else
				localG = g;
			#endif

			gsl_vector_set(F1, i, Qval);
			gsl_vector_set(F2, i, beta*Qval*Qval/Aval + g*I1val);
			gsl_vector_set(ST21, i, g*I2val);
			//gsl_vector_set(ST22, i, g*hval*bval*S0);
			gsl_vector_set(ST22, i, g*Aval*S0);
			gsl_vector_set(ST23, i, g*Aval*Sfval);

		}

		gsl_vector *localFhat1 = gsl_vector_alloc(2);
		gsl_vector_set(localFhat1, 0, Fhat1R[k]);
		gsl_vector_set(localFhat1, 1, Fhat1L[k+1]);

		gsl_vector *localFhat2 = gsl_vector_alloc(2);
		gsl_vector_set(localFhat2, 0, Fhat2R[k]);
		gsl_vector_set(localFhat2, 1, Fhat2L[k+1]);

		// cacluate the volume integral of the flux
		gsl_vector *localRHS1 = gsl_vector_calloc(Np);
		gsl_vector *localRHS2 = gsl_vector_calloc(Np);
		gsl_blas_dgemv(CblasNoTrans, 2.0/h, VolMat, F1, 1.0, localRHS1);
		gsl_blas_dgemv(CblasNoTrans, 2.0/h, VolMat, F2, 1.0, localRHS2);

		// calculate the surface integral of the flux
		gsl_vector *SurfPart1 = gsl_vector_calloc(Np);
		gsl_vector *SurfPart2 = gsl_vector_calloc(Np);
		gsl_blas_dgemv(CblasNoTrans, 2.0/h, LIFT, localFhat1, 1.0, SurfPart1);
		gsl_blas_dgemv(CblasNoTrans, 2.0/h, LIFT, localFhat2, 1.0, SurfPart2);

		// calculate the RHS
		gsl_vector_add(localRHS1, SurfPart1);
		gsl_vector_add(localRHS2, SurfPart2);
		gsl_vector_add(localRHS2, ST21);
		gsl_vector_add(localRHS2, ST22);
		gsl_vector_sub(localRHS2, ST23);

		gsl_vector_free(F1);
		gsl_vector_free(F2);
		gsl_vector_free(localFhat1);
		gsl_vector_free(localFhat2);
		gsl_vector_free(SurfPart1);
		gsl_vector_free(SurfPart2);
		gsl_vector_free(ST21);
		gsl_vector_free(ST22);
		gsl_vector_free(ST23);

		for(int i = 0; i < Np; i++)
		{
			RHSA[begNode+i+1] = gsl_vector_get(localRHS1, i);
			RHSQ[begNode+i+1] = gsl_vector_get(localRHS2, i);
		}

		gsl_vector_free(localRHS1);
		gsl_vector_free(localRHS2);

	/*
		double F_el1[2] = {Fhat1R[k], Fhat1L[k+1]};
		double F_el2[2] = {Fhat2R[k], Fhat2L[k+1]};

		double invMrow1[2] = {4/h, -2/h};
		double invMrow2[2] = {-2/h, 4/h};
		//double invM[2][2] = {{4/(b-a), 2/(a-b)},{2/(a-b), 4/(b-a)}};

		double S1 = IP1[k] + F_el1[0];
		double S2 = IP2[k] - F_el1[1];
		double S3 = IP3[k] + F_el2[0]  + IP4[k] + IP5[k] - IP6[k];
		double S4 = IP7[k] - F_el2[1] + IP8[k] + IP9[k] - IP10[k];

		// Compute RHSA = invM*[S1;S2] and RHSQ = invM*[S3;S4]
		RHSA[2*k + 1] = invMrow1[0]*S1 + invMrow1[1]*S2;
		RHSA[2*k + 2] = invMrow2[0]*S1 + invMrow2[1]*S2;
		RHSQ[2*k + 1] = invMrow1[0]*S3 + invMrow1[1]*S4;
		RHSQ[2*k + 2] = invMrow2[0]*S3 + invMrow2[1]*S4;
*/
	}

	// The value at the ghost nodes will be the same as the values at the other side of the face

	RHSA[0] = RHSA[1];
	RHSA[NumNodes - 1] = RHSA[NumNodes - 2];

	RHSQ[0] = RHSQ[1];
	RHSQ[NumNodes - 1] = RHSQ[NumNodes - 2];

	free(Fhat1L);
	free(Fhat1R);
	free(Fhat2L);
	free(Fhat2R);

}

ext: c

max_stars_repo_name: evalseth/DG-RAIN

max_issues_repo_name: evalseth/DG-RAIN

max_forks_repo_name: evalseth/DG-RAIN


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

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

#undef COMPLEX
#define REAL

static int check_solution(PLASMA_enum side, PLASMA_enum uplo, PLASMA_enum trans, PLASMA_enum diag,
                          int M, int N, double alpha,
                          double *A, int LDA,
                          double *Bref, double *Bplasma, int LDB);

int testing_dtrsm(int argc, char **argv)
{
    /* Check for number of arguments*/
    if ( argc != 5 ) {
        USAGE("TRSM", "alpha M N LDA LDB",
              "   - alpha  : alpha coefficient\n"
              "   - M      : number of rows of matrices B\n"
              "   - N      : number of columns of matrices B\n"
              "   - LDA    : leading dimension of matrix A\n"
              "   - LDB    : leading dimension of matrix B\n");
        return -1;
    }

    double alpha = (double) atol(argv[0]);
    int M     = atoi(argv[1]);
    int N     = atoi(argv[2]);
    int LDA   = atoi(argv[3]);
    int LDB   = atoi(argv[4]);

    double eps;
    int info_solution;
    int s, u, t, d, i;
    int LDAxM = LDA*max(M,N);
    int LDBxN = LDB*max(M,N);

    double *A      = (double *)malloc(LDAxM*sizeof(double));
	#pragma omp register([LDAxM]A)
    double *B      = (double *)malloc(LDBxN*sizeof(double));
	#pragma omp register([LDBxN]B)
    double *Binit  = (double *)malloc(LDBxN*sizeof(double));
	#pragma omp register([LDBxN]Binit)
    double *Bfinal = (double *)malloc(LDBxN*sizeof(double));
	#pragma omp register([LDBxN]Bfinal)

    /* Check if unable to allocate memory */
    if ( (!A) || (!B) || (!Binit) || (!Bfinal)){
        printf("Out of Memory \n ");
        return -2;
    }

    eps = LAPACKE_dlamch_work('e');

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

    /*----------------------------------------------------------
     *  TESTING DTRSM
     */

    /* Initialize A, B, C */
    LAPACKE_dlarnv_work(IONE, ISEED, LDAxM, A);
    LAPACKE_dlarnv_work(IONE, ISEED, LDBxN, B);
    for(i=0;i<max(M,N);i++)
      A[LDA*i+i] = A[LDA*i+i] + 2.0;

    for (s=0; s<2; s++) {
        for (u=0; u<2; u++) {
#ifdef COMPLEX
            for (t=0; t<3; t++) {
#else
            for (t=0; t<2; t++) {
#endif
                for (d=0; d<2; d++) {

                    memcpy(Binit,  B, LDBxN*sizeof(double));
                    memcpy(Bfinal, B, LDBxN*sizeof(double));

                    /* PLASMA DTRSM */
                    PLASMA_dtrsm(side[s], uplo[u], trans[t], diag[d],
                                 M, N, alpha, A, LDA, Bfinal, LDB);

                    /* Check the solution */
                    info_solution = check_solution(side[s], uplo[u], trans[t], diag[d],
                                                   M, N, alpha, A, LDA, Binit, Bfinal, LDB);

                    printf("***************************************************\n");
                    if (info_solution == 0) {
                        printf(" ---- TESTING DTRSM (%s, %s, %s, %s) ...... PASSED !\n",
                               sidestr[s], uplostr[u], transstr[t], diagstr[d]);
                    }
                    else {
                        printf(" ---- TESTING DTRSM (%s, %s, %s, %s) ... FAILED !\n",
                               sidestr[s], uplostr[u], transstr[t], diagstr[d]);
                    }
                    printf("***************************************************\n");
                }
            }
        }
    }

    free(A); free(B);
    free(Binit); free(Bfinal);

    return 0;
}

/*--------------------------------------------------------------
 * Check the solution
 */
static int check_solution(PLASMA_enum side, PLASMA_enum uplo, PLASMA_enum trans, PLASMA_enum diag,
                          int M, int N, double alpha,
                          double *A, int LDA,
                          double *Bref, double *Bplasma, int LDB)
{
    int info_solution;
    double Anorm, Binitnorm, Bplasmanorm, Blapacknorm, Rnorm, result;
    double eps;
    double mzone = (double)-1.0;

    double *work = (double *)malloc(max(M, N)* sizeof(double));
    int Am, An;

    if (side == PlasmaLeft) {
        Am = M; An = M;
    } else {
        Am = N; An = N;
    }

    Anorm       = LAPACKE_dlantr_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), lapack_const(uplo), lapack_const(diag),
                                Am, An, A, LDA, work);
    Binitnorm   = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, Bref,    LDB, work);
    Bplasmanorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, Bplasma, LDB, work);

    cblas_dtrsm(CblasColMajor, (CBLAS_SIDE)side, (CBLAS_UPLO)uplo, (CBLAS_TRANSPOSE)trans,
                (CBLAS_DIAG)diag, M, N, (alpha), A, LDA, Bref, LDB);

    Blapacknorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, Bref, LDB, work);

    cblas_daxpy(LDB * N, (mzone), Bplasma, 1, Bref, 1);

    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, Bref, LDB, work);

    eps = LAPACKE_dlamch_work('e');

    printf("Rnorm %e, Anorm %e, Binitnorm %e, Bplasmanorm %e, Blapacknorm %e\n",
           Rnorm, Anorm, Binitnorm, Bplasmanorm, Blapacknorm);

    result = Rnorm / ((Anorm + Blapacknorm) * max(M,N) * eps);

    printf("============\n");
    printf("Checking the norm of the difference against reference DTRSM \n");
    printf("-- ||Cplasma - Clapack||_oo/((||A||_oo+||B||_oo).N.eps) = %e \n", result);

    if ( isinf(Blapacknorm) || isinf(Bplasmanorm) || isnan(result) || isinf(result) || (result > 10.0) ) {
        printf("-- The solution is suspicious ! \n");
        info_solution = 1;
    }
    else {
        printf("-- The solution is CORRECT ! \n");
        info_solution= 0 ;
    }
    free(work);

    return info_solution;
}


text:
#pragma once
#include <vector>
#include <chrono>
#include <complex>
#include <boost/circular_buffer.hpp>
#include <gsl/gsl_odeiv.h>
#include "Tract.h"

class SpeechSynthesizer
{
public:
	double h_min;

private:
	// vocal chord
	const double m1;	// mass 1 of vocal cords
	const double m2;	// mass 2 of vocal cords
	const double d1;	// thickness of m1
	const double d2;	// thickness of m2
	const double lg;	// effective length of vocal chords
	const double es;	// nonlinear spring coeff.
    const double eh;
    const double ks1;	// linear spring coeff.
    const double ks2;
    const double kh1;
    const double kh2;
    const double kc;	// coupling spring coeff.
    const double om1;
	double Ag0;
	double p_s;
	double x1, x2, v1, v2;
	boost::circular_buffer<double> ug_buf;	// front : current (latest) value
	double dydt[5];
	// vocal tract
    const double toVelum;
    const double toSinus;
    const double toConstr;
    const double R_sin;
    const double L_sin;
    const double C_sin;
	Tract oral, nasal;
	double *z_in_curr, *h_out_curr;
	double *z_in_prev, *h_out_prev;
	// other params
	const double r_N;
	const double r_vib;
	// simulation params
	const double dt, df;
	const int N;
	const double calcSpan;
	double lastCalcSpan;
	double lastCalcTime;
	gsl_odeiv_system system;
	gsl_odeiv_step *step;
	double t;

public:
	SpeechSynthesizer(int sampling, int freqResolution, std::vector<double> &tractArea, double dl);
	~SpeechSynthesizer(void);
	double Step(void);
	double GetTime(void);

private:
	static int Differentiate(double t, const double y[], double dydt[], void *params);
	void CalcTract(void);
	void CalcResponse(double res[], const std::complex<double> spec[]);
	void GetTractImpedance(double res[], double t);
	void GetOutputImpedance(double res[], double t);
};

max_stars_repo_name: kinoh/VoiceSynthesis

max_issues_repo_name: kinoh/VoiceSynthesis

max_forks_repo_name: kinoh/VoiceSynthesis


text:
#ifndef matops_h
#define matops_h

#include <ceed.h>
#include <petsc.h>
#include <petscdmplex.h>

#include "structs.h"

PetscErrorCode MatGetDiag(Mat A, Vec D);
PetscErrorCode ApplyLocal_Ceed(Vec X, Vec Y, UserO user);
PetscErrorCode MatMult_Ceed(Mat A, Vec X, Vec Y);
PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx);
PetscErrorCode MatMult_Prolong(Mat A, Vec X, Vec Y);
PetscErrorCode MatMult_Restrict(Mat A, Vec X, Vec Y);
PetscErrorCode ComputeErrorMax(UserO user, CeedOperator op_error,
                               Vec X, CeedVector target, PetscReal *max_error);

#endif // matops_h

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "compearth.h"
#ifdef COMPEARTH_USE_MKL
#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wreserved-id-macro"
#pragma clang diagnostic ignored "-Wstrict-prototypes"
#endif
#include <mkl_cblas.h>
#ifdef __clang__
#pragma clang diagnostic pop
#endif
#else
#include <cblas.h>
#endif

/*!
 * @brief Computes the angle between moment tensors.  This is from
 *        Equation 15 of Tape and Tape 2016 - A confidence parameter for
 *        seismic moment tensors
 *
 * @param[in] n       Number of moment tensors.
 * @param[in] M1in    First moment tensor packed {m11, m22, m33, m12, m13, m23}.
 *                    This is an array of dimension [6 x n] with leading
 *                    dimension 6.
 * @param[in] M2in    Second moment tensor packed {m11, m22, m33, m12, m13, m23}.
 *                    This is an array of dimension [6 x n] with leading
 *                    dimension 6.
 *
 * @param[out] theta  Angle between moment tensors M1 and M2 in radians.
 *                    This is an array of dimension [n].
 *
 * @result 0 indicates success.
 *
 * @author Ben Baker
 *
 * @copyright MIT
 *
 */
int compearth_angleMT(const int n,
                      const double *__restrict__ M1in,
                      const double *__restrict__ M2in,
                      double *__restrict__ theta)
{
    double M1[9*CE_CHUNKSIZE] __attribute__((aligned(64)));
    double M2[9*CE_CHUNKSIZE] __attribute__((aligned(64)));
    double M1_mag[CE_CHUNKSIZE] __attribute__((aligned(64)));
    double M2_mag[CE_CHUNKSIZE] __attribute__((aligned(64)));
    double xden[CE_CHUNKSIZE] __attribute__((aligned(64)));
    double arg[CE_CHUNKSIZE] __attribute__((aligned(64)));
    double xnum[CE_CHUNKSIZE] __attribute__((aligned(64)));
    int i, ierr, imt, nmtLoc;
    ierr = 0;
    for (imt=0; imt<n; imt=imt+CE_CHUNKSIZE)
    {
        nmtLoc = MIN(CE_CHUNKSIZE, n - imt);
        // Convert to 3x3 matrices
        compearth_Mvec2Mmat(nmtLoc, &M1in[6*imt], 1, M1);
        compearth_Mvec2Mmat(nmtLoc, &M2in[6*imt], 1, M2);
        // Compute norms
        ierr = compearth_normMat(nmtLoc, M1, CE_TWO_NORM, 2.0, M1_mag);
        ierr = compearth_normMat(nmtLoc, M2, CE_TWO_NORM, 2.0, M2_mag);
        // Compute the denomimators
        for (i=0; i<nmtLoc; i++)
        {
            xden[i] = M1_mag[i]*M2_mag[i];
        }
        // Compute the numerators
        for (i=0; i<nmtLoc; i++)
        {
            xnum[i] = cblas_ddot(9, &M1[9*i], 1, &M2[9*i], 1);
        }
        // Compute argument and do all checks
        for (i=0; i<nmtLoc; i++)
        {
            // Can't divide by zero - force this to work b/c M1 or M2 is
            // equivalently 0.
            if (xden[i] == 0.0)
            {
                xnum[i] = 0.0;
                xden[i] = 1.0;
                fprintf(stderr, "%s: Denominator is 0\n", __func__);
            }
            arg[i] = xnum[i]/xden[i];
            // Correct for numerical errors
            if (arg[i] <-1.0)
            {
                fprintf(stdout, "%s: Warning arg is %f setting to -1\n",
                         __func__, arg[i]);
                arg[i] =-1.0;
            }
            if (arg[i] > 1.0)
            {
                fprintf(stdout, "%s: Warning arg is %f setting to +1\n",
                        __func__, arg[i]);
                arg[i] = 1.0;
            }
        }
        // Finally compute the angle
        for (i=0; i<nmtLoc; i++)
        {
            theta[imt+i] = acos(arg[i]);
        }
    }
    return ierr;
}

max_stars_repo_name: OUCyf/mtbeach

max_issues_repo_name: carltape/mtbeach

max_forks_repo_name: carltape/mtbeach


text:
//CELL (~soma) stage: IIR filtering of each row or col of X according to dim.
//The IIR filters are specified by an Nx(Q+1) or (Q+1)xN matrix A,
//where Q is the IIR filter order (Q=0 means only a0; Q=1 means a0 and a1; etc.).

//The calling program must ensure that the sizes are correct, the filter is stable, etc.

//I just started this... finish later!! (Or skip.)

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

#ifdef __cplusplus
namespace ov {
extern "C" {
#endif

int iir_s (float *Y, const float *X, const float *A, const size_t N, const size_t T, const int Q, const char iscolmajor, const size_t dim);
int iir_d (double *Y, const double *X, const double *A, const size_t N, const size_t T, const int Q, const char iscolmajor, const size_t dim);
int iir_c (float *Y, const float *X, const float *A, const size_t N, const size_t T, const int Q, const char iscolmajor, const size_t dim);
int iir_z (double *Y, const double *X, const double *A, const size_t N, const size_t T, const int Q, const char iscolmajor, const size_t dim);

int iir_inplace_s (float *X, const float *A, const size_t N, const size_t T, const int Q, const char iscolmajor, const size_t dim);
int iir_inplace_d (double *X, const double *A, const size_t N, const size_t T, const int Q, const char iscolmajor, const size_t dim);
int iir_inplace_c (float *X, const float *A, const size_t N, const size_t T, const int Q, const char iscolmajor, const size_t dim);
int iir_inplace_z (double *X, const double *A, const size_t N, const size_t T, const int Q, const char iscolmajor, const size_t dim);

int iir_inplace_s (float *X, const float *A, const size_t N, const size_t T, const int Q, const char iscolmajor, const size_t dim)
{
    const size_t M = Q - 1;
    int n, t;

    //Checks
    if (N<1) { fprintf(stderr,"error in iir_s: N (num neurons) must be positive\n"); return 1; }
    if (T<1) { fprintf(stderr,"error in iir_s: T (num time points) must be positive\n"); return 1; }
    if (Q<0) { fprintf(stderr,"error in iir_s: Q (filter order) must be nonnegative\n"); return 1; }

    if (N==1u)
    {
        if (A[0]!=1.0f) { cblas_sscal((int)T,1.0f/A[0],A,1); cblas_sscal((int)(N*T),1.0f/A[0],X,1); }
        for (size_t t=1; t<M; ++t) { X[t] -= cblas_sdot(t,&A[M-t],1,&X[0],1); }
        for (size_t t=M; t<T; ++t) { X[t] -= cblas_sdot(M,&A[0],1,&X[t-M],1); }
    }
    else if (dim==0u)
    {
        if (N<3)
        {
            if (iscolmajor)
            {
                for (size_t n=0u; n<N; ++n)
                {
                    if (A[n*Q1]!=1.0f) { cblas_sscal((int)T,1.0f/A[n*Q1],&A[n*Q1],1); }
                }
                for (size_t n=1; n<M; ++n)
                {
                    X[n] -= cblas_sdot(n,&A[M-n],1,&X[0],1);
                    X[n+R] -= cblas_sdot(n,&A[M-n],1,&X[R],1);
                }
                for (size_t n=M; n<R; ++n)
                {
                    X[n] -= cblas_sdot(M,&A[0],1,&X[n-M],1);
                    X[n+R] -= cblas_sdot(M,&A[0],1,&X[n-M+R],1);
                }
            }
            else
            {
                for (size_t n=1; n<M; ++n)
                {
                    X[n*C] -= cblas_sdot(n,&A[M-n],1,&X[0],(int)C);
                    X[1+n*C] -= cblas_sdot(n,&A[M-n],1,&X[1],(int)C);
                }
                for (size_t n=M; n<R; ++n)
                {
                    X[n*C] -= cblas_sdot(M,&A[0],1,&X[(n-M)*C],(int)C);
                    X[1+n*C] -= cblas_sdot(M,&A[0],1,&X[1+(n-M)*C],(int)C);
                }
            }
        }
        else
        {
            if (iscolmajor)
            {
                for (size_t n=1; n<M; ++n) { cblas_sgemv(CblasColMajor,CblasTrans,n,(int)C,-1.0f,&X[0],(int)R,&A[M-n],1,1.0f,&X[n],(int)R); }
                for (size_t n=M; n<R; ++n) { cblas_sgemv(CblasColMajor,CblasTrans,M,(int)C,-1.0f,&X[n-M],(int)R,&A[0],1,1.0f,&X[n],(int)R); }
            }
            else
            {
                for (size_t n=1; n<M; ++n) { cblas_sgemv(CblasRowMajor,CblasTrans,n,(int)C,-1.0f,&X[0],(int)C,&A[M-n],1,1.0f,&X[n*C],1); }
                for (size_t n=M; n<R; ++n) { cblas_sgemv(CblasRowMajor,CblasTrans,M,(int)C,-1.0f,&X[(n-M)*C],(int)C,&A[0],1,1.0f,&X[n*C],1); }
            }
        }
    }
    else if (dim==1u)
    {
        if (R==1)
        {
            for (size_t n=1; n<M; ++n) { X[n] -= cblas_sdot(n,&A[M-n],1,&X[0],1); }
            for (size_t n=M; n<C; ++n) { X[n] -= cblas_sdot(M,&A[0],1,&X[n-M],1); }
        }
        else if (R==2)
        {
            if (iscolmajor)
            {
                for (size_t n=1; n<M; ++n)
                {
                    X[n*R] -= cblas_sdot(n,&A[M-n],1,&X[0],(int)R);
                    X[1+n*R] -= cblas_sdot(n,&A[M-n],1,&X[1],(int)R);
                }
                for (size_t n=M; n<C; ++n)
                {
                    X[n*R] -= cblas_sdot(M,&A[0],1,&X[(n-M)*R],(int)R);
                    X[1+n*R] -= cblas_sdot(M,&A[0],1,&X[1+(n-M)*R],(int)R);
                }
            }
            else
            {
                for (size_t n=1; n<M; ++n)
                {
                    X[n] -= cblas_sdot(n,&A[M-n],1,&X[0],1);
                    X[n+C] -= cblas_sdot(n,&A[M-n],1,&X[C],1);
                }
                for (size_t n=M; n<C; ++n)
                {
                    X[n] -= cblas_sdot(M,&A[0],1,&X[n-M],1);
                    X[n+C] -= cblas_sdot(M,&A[0],1,&X[n-M+C],1);
                }
            }
        }
        else
        {
            if (iscolmajor)
            {
                for (size_t n=1; n<M; ++n) { cblas_sgemv(CblasColMajor,CblasNoTrans,(int)R,n,-1.0f,&X[0],(int)R,&A[M-n],1,1.0f,&X[n*R],1); }
                for (size_t n=M; n<C; ++n) { cblas_sgemv(CblasColMajor,CblasNoTrans,(int)R,M,-1.0f,&X[(n-M)*R],(int)R,&A[0],1,1.0f,&X[n*R],1); }
            }
            else
            {
                for (size_t n=1; n<M; ++n) { cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)R,n,-1.0f,&X[0],(int)C,&A[M-n],1,1.0f,&X[n],(int)C); }
                for (size_t n=M; n<C; ++n) { cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)R,M,-1.0f,&X[n-M],(int)C,&A[0],1,1.0f,&X[n],(int)C); }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in iir_s: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int iir_d (double *X, const char iscolmajor, const size_t R, const size_t C, const double *A, const size_t N, const size_t dim)
{
    const size_t M = N - 1;
    int n;

    //Checks
    if (R<1) { fprintf(stderr,"error in iir_d: R (nrows X) must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in iir_d: C (ncols X) must be positive\n"); return 1; }
    if (N<1) { fprintf(stderr,"error in iir_d: N (filter order) must be positive\n"); return 1; }

    if (dim==0u)
    {
        if (C==1)
        {
            for (size_t n=1; n<M; ++n) { X[n] -= cblas_ddot(n,&A[M-n],1,&X[0],1); }
            for (size_t n=M; n<R; ++n) { X[n] -= cblas_ddot(M,&A[0],1,&X[n-M],1); }
        }
        else if (C==2)
        {
            if (iscolmajor)
            {
                for (size_t n=1; n<M; ++n)
                {
                    X[n] -= cblas_ddot(n,&A[M-n],1,&X[0],1);
                    X[n+R] -= cblas_ddot(n,&A[M-n],1,&X[R],1);
                }
                for (size_t n=M; n<R; ++n)
                {
                    X[n] -= cblas_ddot(M,&A[0],1,&X[n-M],1);
                    X[n+R] -= cblas_ddot(M,&A[0],1,&X[n-M+R],1);
                }
            }
            else
            {
                for (size_t n=1; n<M; ++n)
                {
                    X[n*C] -= cblas_ddot(n,&A[M-n],1,&X[0],(int)C);
                    X[1+n*C] -= cblas_ddot(n,&A[M-n],1,&X[1],(int)C);
                }
                for (size_t n=M; n<R; ++n)
                {
                    X[n*C] -= cblas_ddot(M,&A[0],1,&X[(n-M)*C],(int)C);
                    X[1+n*C] -= cblas_ddot(M,&A[0],1,&X[1+(n-M)*C],(int)C);
                }
            }
        }
        else
        {
            if (iscolmajor)
            {
                for (size_t n=1; n<M; ++n) { cblas_dgemv(CblasColMajor,CblasTrans,n,(int)C,-1.0,&X[0],(int)R,&A[M-n],1,1.0,&X[n],(int)R); }
                for (size_t n=M; n<R; ++n) { cblas_dgemv(CblasColMajor,CblasTrans,M,(int)C,-1.0,&X[n-M],(int)R,&A[0],1,1.0,&X[n],(int)R); }
            }
            else
            {
                for (size_t n=1; n<M; ++n) { cblas_dgemv(CblasRowMajor,CblasTrans,n,(int)C,-1.0,&X[0],(int)C,&A[M-n],1,1.0,&X[n*C],1); }
                for (size_t n=M; n<R; ++n) { cblas_dgemv(CblasRowMajor,CblasTrans,M,(int)C,-1.0,&X[(n-M)*C],(int)C,&A[0],1,1.0,&X[n*C],1); }
            }
        }
    }
    else if (dim==1u)
    {
        if (R==1)
        {
            for (size_t n=1; n<M; ++n) { X[n] -= cblas_ddot(n,&A[M-n],1,&X[0],1); }
            for (size_t n=M; n<C; ++n) { X[n] -= cblas_ddot(M,&A[0],1,&X[n-M],1); }
        }
        else if (R==2)
        {
            if (iscolmajor)
            {
                for (size_t n=1; n<M; ++n)
                {
                    X[n*R] -= cblas_ddot(n,&A[M-n],1,&X[0],(int)R);
                    X[1+n*R] -= cblas_ddot(n,&A[M-n],1,&X[1],(int)R);
                }
                for (size_t n=M; n<C; ++n)
                {
                    X[n*R] -= cblas_ddot(M,&A[0],1,&X[(n-M)*R],(int)R);
                    X[1+n*R] -= cblas_ddot(M,&A[0],1,&X[1+(n-M)*R],(int)R);
                }
            }
            else
            {
                for (size_t n=1; n<M; ++n)
                {
                    X[n] -= cblas_ddot(n,&A[M-n],1,&X[0],1);
                    X[n+C] -= cblas_ddot(n,&A[M-n],1,&X[C],1);
                }
                for (size_t n=M; n<C; ++n)
                {
                    X[n] -= cblas_ddot(M,&A[0],1,&X[n-M],1);
                    X[n+C] -= cblas_ddot(M,&A[0],1,&X[n-M+C],1);
                }
            }
        }
        else
        {
            if (iscolmajor)
            {
                for (size_t n=1; n<M; ++n) { cblas_dgemv(CblasColMajor,CblasNoTrans,(int)R,n,-1.0,&X[0],(int)R,&A[M-n],1,1.0,&X[n*R],1); }
                for (size_t n=M; n<C; ++n) { cblas_dgemv(CblasColMajor,CblasNoTrans,(int)R,M,-1.0,&X[(n-M)*R],(int)R,&A[0],1,1.0,&X[n*R],1); }
            }
            else
            {
                for (size_t n=1; n<M; ++n) { cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)R,n,-1.0,&X[0],(int)C,&A[M-n],1,1.0,&X[n],(int)C); }
                for (size_t n=M; n<C; ++n) { cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)R,M,-1.0,&X[n-M],(int)C,&A[0],1,1.0,&X[n],(int)C); }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in iir_d: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int iir_c (float *X, const char iscolmajor, const size_t R, const size_t C, const float *A, const size_t N, const size_t dim)
{
    const float a[2] = {-1.0f,0.0f}, b[2] = {1.0f,0.0f};
    const size_t M = N - 1;
    int n;

    //Checks
    if (R<1) { fprintf(stderr,"error in iir_c: R (nrows X) must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in iir_c: C (ncols X) must be positive\n"); return 1; }
    if (N<1) { fprintf(stderr,"error in iir_c: N (filter order) must be positive\n"); return 1; }

    if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=1; n<M; ++n) { cblas_cgemv(CblasColMajor,CblasTrans,n,(int)C,&a[0],&X[0],(int)R,&A[2*(M-n)],1,&b[0],&X[2*n],(int)R); }
            for (size_t n=M; n<R; ++n) { cblas_cgemv(CblasColMajor,CblasTrans,M,(int)C,&a[0],&X[2*(n-M)],(int)R,&A[0],1,&b[0],&X[2*n],(int)R); }
        }
        else
        {
            for (size_t n=1; n<M; ++n) { cblas_cgemv(CblasRowMajor,CblasTrans,n,(int)C,&a[0],&X[0],(int)C,&A[2*(M-n)],1,&b[0],&X[2*n*C],1); }
            for (size_t n=M; n<R; ++n) { cblas_cgemv(CblasRowMajor,CblasTrans,M,(int)C,&a[0],&X[2*(n-M)*C],(int)C,&A[0],1,&b[0],&X[2*n*C],1); }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=1; n<M; ++n) { cblas_cgemv(CblasColMajor,CblasNoTrans,(int)R,n,&a[0],&X[0],(int)R,&A[2*(M-n)],1,&b[0],&X[2*n*R],1); }
            for (size_t n=M; n<C; ++n) { cblas_cgemv(CblasColMajor,CblasNoTrans,(int)R,M,&a[0],&X[2*(n-M)*R],(int)R,&A[0],1,&b[0],&X[2*n*R],1); }
        }
        else
        {
            for (size_t n=1; n<M; ++n) { cblas_cgemv(CblasRowMajor,CblasNoTrans,(int)R,n,&a[0],&X[0],(int)C,&A[2*(M-n)],1,&b[0],&X[2*n],(int)C); }
            for (size_t n=M; n<C; ++n) { cblas_cgemv(CblasRowMajor,CblasNoTrans,(int)R,M,&a[0],&X[2*(n-M)],(int)C,&A[0],1,&b[0],&X[2*n],(int)C); }
        }
    }
    else
    {
        fprintf(stderr,"error in iir_c: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int iir_z (double *X, const char iscolmajor, const size_t R, const size_t C, const double *A, const size_t N, const size_t dim)
{
    const double a[2] = {-1.0,0.0}, b[2] = {1.0,0.0};
    const size_t M = N - 1;
    int n;

    //Checks
    if (R<1) { fprintf(stderr,"error in iir_z: R (nrows X) must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in iir_z: C (ncols X) must be positive\n"); return 1; }
    if (N<1) { fprintf(stderr,"error in iir_z: N (filter order) must be positive\n"); return 1; }

    if (dim==0u)
    {
        if (iscolmajor)
        {
            for (size_t n=1; n<M; ++n) { cblas_zgemv(CblasColMajor,CblasTrans,n,(int)C,&a[0],&X[0],(int)R,&A[2*(M-n)],1,&b[0],&X[2*n],(int)R); }
            for (size_t n=M; n<R; ++n) { cblas_zgemv(CblasColMajor,CblasTrans,M,(int)C,&a[0],&X[2*(n-M)],(int)R,&A[0],1,&b[0],&X[2*n],(int)R); }
        }
        else
        {
            for (size_t n=1; n<M; ++n) { cblas_zgemv(CblasRowMajor,CblasTrans,n,(int)C,&a[0],&X[0],(int)C,&A[2*(M-n)],1,&b[0],&X[2*n*C],1); }
            for (size_t n=M; n<R; ++n) { cblas_zgemv(CblasRowMajor,CblasTrans,M,(int)C,&a[0],&X[2*(n-M)*C],(int)C,&A[0],1,&b[0],&X[2*n*C],1); }
        }
    }
    else if (dim==1u)
    {
        if (iscolmajor)
        {
            for (size_t n=1; n<M; ++n) { cblas_zgemv(CblasColMajor,CblasNoTrans,(int)R,n,&a[0],&X[0],(int)R,&A[2*(M-n)],1,&b[0],&X[2*n*R],1); }
            for (size_t n=M; n<C; ++n) { cblas_zgemv(CblasColMajor,CblasNoTrans,(int)R,M,&a[0],&X[2*(n-M)*R],(int)R,&A[0],1,&b[0],&X[2*n*R],1); }
        }
        else
        {
            for (size_t n=1; n<M; ++n) { cblas_zgemv(CblasRowMajor,CblasNoTrans,(int)R,n,&a[0],&X[0],(int)C,&A[2*(M-n)],1,&b[0],&X[2*n],(int)C); }
            for (size_t n=M; n<C; ++n) { cblas_zgemv(CblasRowMajor,CblasNoTrans,(int)R,M,&a[0],&X[2*(n-M)],(int)C,&A[0],1,&b[0],&X[2*n],(int)C); }
        }
    }
    else
    {
        fprintf(stderr,"error in iir_z: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif


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

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_sf_trig.h>

#include "error.h"

#include "chebyshev.h"
#include "cheb_eval.c"

/* sinh(x) series
 * double-precision for |x| < 1.0
 */
inline
static
int
sinh_series(const double x, double * result)
{
  const double y = x*x;
  const double c0 = 1.0/6.0;
  const double c1 = 1.0/120.0;
  const double c2 = 1.0/5040.0;
  const double c3 = 1.0/362880.0;
  const double c4 = 1.0/39916800.0;
  const double c5 = 1.0/6227020800.0;
  const double c6 = 1.0/1307674368000.0;
  const double c7 = 1.0/355687428096000.0;
  *result = x*(1.0 + y*(c0+y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*(c6+y*c7))))))));
  return GSL_SUCCESS;
}

/* cosh(x)-1 series
 * double-precision for |x| < 1.0
 */
inline
static
int
cosh_m1_series(const double x, double * result)
{
  const double y = x*x;
  const double c0 = 0.5;
  const double c1 = 1.0/24.0;
  const double c2 = 1.0/720.0;
  const double c3 = 1.0/40320.0;
  const double c4 = 1.0/3628800.0;
  const double c5 = 1.0/479001600.0;
  const double c6 = 1.0/87178291200.0;
  const double c7 = 1.0/20922789888000.0;
  const double c8 = 1.0/6402373705728000.0;
  *result = y*(c0+y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*(c6+y*(c7+y*c8))))))));
  return GSL_SUCCESS;
}

/* Chebyshev expansion for f(t) = sinc((t+1)/2), -1 < t < 1
 */
static double sinc_data[17] = {
  1.133648177811747875422,
 -0.532677564732557348781,
 -0.068293048346633177859,
  0.033403684226353715020,
  0.001485679893925747818,
 -0.000734421305768455295,
 -0.000016837282388837229,
  0.000008359950146618018,
  0.000000117382095601192,
 -0.000000058413665922724,
 -0.000000000554763755743,
  0.000000000276434190426,
  0.000000000001895374892,
 -0.000000000000945237101,
 -0.000000000000004900690,
  0.000000000000002445383,
  0.000000000000000009925
};
static cheb_series sinc_cs = {
  sinc_data,
  16,
  -1, 1,
  10
};

/* Chebyshev expansion for f(t) = g((t+1)Pi/8), -1<t<1
 * g(x) = (sin(x)/x - 1)/(x*x)
 */
static double sin_data[12] = {
  -0.3295190160663511504173,
   0.0025374284671667991990,
   0.0006261928782647355874,
  -4.6495547521854042157541e-06,
  -5.6917531549379706526677e-07,
   3.7283335140973803627866e-09,
   3.0267376484747473727186e-10,
  -1.7400875016436622322022e-12,
  -1.0554678305790849834462e-13,
   5.3701981409132410797062e-16,
   2.5984137983099020336115e-17,
  -1.1821555255364833468288e-19
};
static cheb_series sin_cs = {
  sin_data,
  11,
  -1, 1,
  11
};

/* Chebyshev expansion for f(t) = g((t+1)Pi/8), -1<t<1
 * g(x) = (2(cos(x) - 1)/(x^2) + 1) / x^2
 */
static double cos_data[11] = {
  0.165391825637921473505668118136,
 -0.00084852883845000173671196530195,
 -0.000210086507222940730213625768083,
  1.16582269619760204299639757584e-6,
  1.43319375856259870334412701165e-7,
 -7.4770883429007141617951330184e-10,
 -6.0969994944584252706997438007e-11,
  2.90748249201909353949854872638e-13,
  1.77126739876261435667156490461e-14,
 -7.6896421502815579078577263149e-17,
 -3.7363121133079412079201377318e-18
};
static cheb_series cos_cs = {
  cos_data,
  10,
  -1, 1,
  10
};

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

/* I would have prefered just using the library sin() function.
 * But after some experimentation I decided that there was
 * no good way to understand the error; library sin() is just a black box.
 * So we have to roll our own.
 */
int
gsl_sf_sin_e(double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  {
    const double P1 = 7.85398125648498535156e-1;
    const double P2 = 3.77489470793079817668e-8;
    const double P3 = 2.69515142907905952645e-15;

    const double sgn_x = GSL_SIGN(x);
    const double abs_x = fabs(x);

    if(abs_x < GSL_ROOT4_DBL_EPSILON) {
      const double x2 = x*x;
      result->val = x * (1.0 - x2/6.0);
      result->err = fabs(x*x2*x2 / 100.0);
      return GSL_SUCCESS;
    }
    else {
      double sgn_result = sgn_x;
      double y = floor(abs_x/(0.25*M_PI));
      int octant = y - ldexp(floor(ldexp(y,-3)),3);
      int stat_cs;
      double z;

      if(GSL_IS_ODD(octant)) {
        octant += 1;
        octant &= 07;
        y += 1.0;
      }

      if(octant > 3) {
        octant -= 4;
        sgn_result = -sgn_result;
      }

      z = ((abs_x - y * P1) - y * P2) - y * P3;

      if(octant == 0) {
        gsl_sf_result sin_cs_result;
        const double t = 8.0*fabs(z)/M_PI - 1.0;
        stat_cs = cheb_eval_e(&sin_cs, t, &sin_cs_result);
        result->val = z * (1.0 + z*z * sin_cs_result.val);
      }
      else { /* octant == 2 */
        gsl_sf_result cos_cs_result;
        const double t = 8.0*fabs(z)/M_PI - 1.0;
        stat_cs = cheb_eval_e(&cos_cs, t, &cos_cs_result);
        result->val = 1.0 - 0.5*z*z * (1.0 - z*z * cos_cs_result.val);
      }

      result->val *= sgn_result;

      if(abs_x > 1.0/GSL_DBL_EPSILON) {
        result->err = fabs(result->val);
      }
      else if(abs_x > 100.0/GSL_SQRT_DBL_EPSILON) {
        result->err = 2.0 * abs_x * GSL_DBL_EPSILON * fabs(result->val);
      }
      else if(abs_x > 0.1/GSL_SQRT_DBL_EPSILON) {
        result->err = 2.0 * GSL_SQRT_DBL_EPSILON * fabs(result->val);
      }
      else {
        result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      }

      return stat_cs;
    }
  }
}

int
gsl_sf_cos_e(double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  {
    const double P1 = 7.85398125648498535156e-1;
    const double P2 = 3.77489470793079817668e-8;
    const double P3 = 2.69515142907905952645e-15;

    const double abs_x = fabs(x);

    if(abs_x < GSL_ROOT4_DBL_EPSILON) {
      const double x2 = x*x;
      result->val = 1.0 - 0.5*x2;
      result->err = fabs(x2*x2/12.0);
      return GSL_SUCCESS;
    }
    else {
      double sgn_result = 1.0;
      double y = floor(abs_x/(0.25*M_PI));
      int octant = y - ldexp(floor(ldexp(y,-3)),3);
      int stat_cs;
      double z;

      if(GSL_IS_ODD(octant)) {
        octant += 1;
        octant &= 07;
        y += 1.0;
      }

      if(octant > 3) {
        octant -= 4;
        sgn_result = -sgn_result;
      }

      if(octant > 1) {
        sgn_result = -sgn_result;
      }

      z = ((abs_x - y * P1) - y * P2) - y * P3;

      if(octant == 0) {
        gsl_sf_result cos_cs_result;
        const double t = 8.0*fabs(z)/M_PI - 1.0;
        stat_cs = cheb_eval_e(&cos_cs, t, &cos_cs_result);
        result->val = 1.0 - 0.5*z*z * (1.0 - z*z * cos_cs_result.val);
      }
      else { /* octant == 2 */
        gsl_sf_result sin_cs_result;
        const double t = 8.0*fabs(z)/M_PI - 1.0;
        stat_cs = cheb_eval_e(&sin_cs, t, &sin_cs_result);
        result->val = z * (1.0 + z*z * sin_cs_result.val);
      }

      result->val *= sgn_result;

      if(abs_x > 1.0/GSL_DBL_EPSILON) {
        result->err = fabs(result->val);
      }
      else if(abs_x > 100.0/GSL_SQRT_DBL_EPSILON) {
        result->err = 2.0 * abs_x * GSL_DBL_EPSILON * fabs(result->val);
      }
      else if(abs_x > 0.1/GSL_SQRT_DBL_EPSILON) {
        result->err = 2.0 * GSL_SQRT_DBL_EPSILON * fabs(result->val);
      }
      else {
        result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      }

      return stat_cs;
    }
  }
}

int
gsl_sf_hypot_e(const double x, const double y, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x == 0.0 && y == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    const double a = fabs(x);
    const double b = fabs(y);
    const double min = GSL_MIN_DBL(a,b);
    const double max = GSL_MAX_DBL(a,b);
    const double rat = min/max;
    const double root_term = sqrt(1.0 + rat*rat);

    if(max < GSL_DBL_MAX/root_term) {
      result->val = max * root_term;
      result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      return GSL_SUCCESS;
    }
    else {
      OVERFLOW_ERROR(result);
    }
  }
}

int
gsl_sf_complex_sin_e(const double zr, const double zi,
                        gsl_sf_result * szr, gsl_sf_result * szi)
{
  /* CHECK_POINTER(szr) */
  /* CHECK_POINTER(szi) */

  if(fabs(zi) < 1.0) {
    double ch_m1, sh;
    sinh_series(zi, &sh);
    cosh_m1_series(zi, &ch_m1);
    szr->val = sin(zr)*(ch_m1 + 1.0);
    szi->val = cos(zr)*sh;
    szr->err = 2.0 * GSL_DBL_EPSILON * fabs(szr->val);
    szi->err = 2.0 * GSL_DBL_EPSILON * fabs(szi->val);
    return GSL_SUCCESS;
  }
  else if(fabs(zi) < GSL_LOG_DBL_MAX) {
    double ex = exp(zi);
    double ch = 0.5*(ex+1.0/ex);
    double sh = 0.5*(ex-1.0/ex);
    szr->val = sin(zr)*ch;
    szi->val = cos(zr)*sh;
    szr->err = 2.0 * GSL_DBL_EPSILON * fabs(szr->val);
    szi->err = 2.0 * GSL_DBL_EPSILON * fabs(szi->val);
    return GSL_SUCCESS;
  }
  else {
    OVERFLOW_ERROR_2(szr, szi);
  }
}

int
gsl_sf_complex_cos_e(const double zr, const double zi,
                        gsl_sf_result * czr, gsl_sf_result * czi)
{
  /* CHECK_POINTER(czr) */
  /* CHECK_POINTER(czi) */

  if(fabs(zi) < 1.0) {
    double ch_m1, sh;
    sinh_series(zi, &sh);
    cosh_m1_series(zi, &ch_m1);
    czr->val =  cos(zr)*(ch_m1 + 1.0);
    czi->val = -sin(zr)*sh;
    czr->err = 2.0 * GSL_DBL_EPSILON * fabs(czr->val);
    czi->err = 2.0 * GSL_DBL_EPSILON * fabs(czi->val);
    return GSL_SUCCESS;
  }
  else if(fabs(zi) < GSL_LOG_DBL_MAX) {
    double ex = exp(zi);
    double ch = 0.5*(ex+1.0/ex);
    double sh = 0.5*(ex-1.0/ex);
    czr->val =  cos(zr)*ch;
    czi->val = -sin(zr)*sh;
    czr->err = 2.0 * GSL_DBL_EPSILON * fabs(czr->val);
    czi->err = 2.0 * GSL_DBL_EPSILON * fabs(czi->val);
    return GSL_SUCCESS;
  }
  else {
    OVERFLOW_ERROR_2(czr,czi);
  }
}

int
gsl_sf_complex_logsin_e(const double zr, const double zi,
                           gsl_sf_result * lszr, gsl_sf_result * lszi)
{
  /* CHECK_POINTER(lszr) */
  /* CHECK_POINTER(lszi) */

  if(zi > 60.0) {
    lszr->val = -M_LN2 + zi;
    lszi->val =  0.5*M_PI - zr;
    lszr->err = 2.0 * GSL_DBL_EPSILON * fabs(lszr->val);
    lszi->err = 2.0 * GSL_DBL_EPSILON * fabs(lszi->val);
  }
  else if(zi < -60.0) {
    lszr->val = -M_LN2 - zi;
    lszi->val = -0.5*M_PI + zr;
    lszr->err = 2.0 * GSL_DBL_EPSILON * fabs(lszr->val);
    lszi->err = 2.0 * GSL_DBL_EPSILON * fabs(lszi->val);
  }
  else {
    gsl_sf_result sin_r, sin_i;
    int status;
    gsl_sf_complex_sin_e(zr, zi, &sin_r, &sin_i); /* ok by construction */
    status = gsl_sf_complex_log_e(sin_r.val, sin_i.val, lszr, lszi);
    if(status == GSL_EDOM) {
      DOMAIN_ERROR_2(lszr, lszi);
    }
  }
  return gsl_sf_angle_restrict_symm_e(&(lszi->val));
}

int
gsl_sf_lnsinh_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(fabs(x) < 1.0) {
    double eps;
    sinh_series(x, &eps);
    result->val = log(eps);
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < -0.5*GSL_LOG_DBL_EPSILON) {
    result->val = x + log(0.5*(1.0 - exp(-2.0*x)));
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    result->val = -M_LN2 + x;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int gsl_sf_lncosh_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(fabs(x) < 1.0) {
    double eps;
    cosh_m1_series(x, &eps);
    return gsl_sf_log_1plusx_e(eps, result);
  }
  else if(x < -0.5*GSL_LOG_DBL_EPSILON) {
    result->val = x + log(0.5*(1.0 + exp(-2.0*x)));
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    result->val = -M_LN2 + x;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

/*
inline int gsl_sf_sincos_e(const double theta, double * s, double * c)
{
  double tan_half = tan(0.5 * theta);
  double den = 1. + tan_half*tan_half;
  double cos_theta = (1.0 - tan_half*tan_half) / den;
  double sin_theta = 2.0 * tan_half / den;
}
*/

int
gsl_sf_polar_to_rect(const double r, const double theta,
                          gsl_sf_result * x, gsl_sf_result * y)
{
  double t   = theta;
  int status = gsl_sf_angle_restrict_symm_e(&t);
  double c = cos(t);
  double s = sin(t);
  x->val = r * cos(t);
  y->val = r * sin(t);
  x->err  = r * fabs(s * GSL_DBL_EPSILON * t);
  x->err += 2.0 * GSL_DBL_EPSILON * fabs(x->val);
  y->err  = r * fabs(c * GSL_DBL_EPSILON * t);
  y->err += 2.0 * GSL_DBL_EPSILON * fabs(y->val);
  return status;
}

int
gsl_sf_rect_to_polar(const double x, const double y,
                          gsl_sf_result * r, gsl_sf_result * theta)
{
  int stat_h = gsl_sf_hypot_e(x, y, r);
  if(r->val > 0.0) {
    theta->val = atan2(y, x);
    theta->err = 2.0 * GSL_DBL_EPSILON * fabs(theta->val);
    return stat_h;
  }
  else {
    DOMAIN_ERROR(theta);
  }
}

int gsl_sf_angle_restrict_symm_err_e(const double theta, gsl_sf_result * result)
{
  /* synthetic extended precision constants */
  const double P1 = 4 * 7.8539812564849853515625e-01;
  const double P2 = 4 * 3.7748947079307981766760e-08;
  const double P3 = 4 * 2.6951514290790594840552e-15;
  const double TwoPi = 2*(P1 + P2 + P3);

  const double y = GSL_SIGN(theta) * 2 * floor(fabs(theta)/TwoPi);
  double r = ((theta - y*P1) - y*P2) - y*P3;

  if(r >  M_PI) { r = (((r-2*P1)-2*P2)-2*P3); }  /* r-TwoPi */
  else if (r < -M_PI) r = (((r+2*P1)+2*P2)+2*P3); /* r+TwoPi */

  result->val = r;

  if(fabs(theta) > 0.0625/GSL_DBL_EPSILON) {
    result->val = GSL_NAN;
    result->err = GSL_NAN;
    GSL_ERROR ("error", GSL_ELOSS);
  }
  else if(fabs(theta) > 0.0625/GSL_SQRT_DBL_EPSILON) {
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val - theta);
    return GSL_SUCCESS;
  }
  else {
    double delta = fabs(result->val - theta);
    result->err = 2.0 * GSL_DBL_EPSILON * ((delta < M_PI) ? delta : M_PI);
    return GSL_SUCCESS;
  }
}

int gsl_sf_angle_restrict_pos_err_e(const double theta, gsl_sf_result * result)
{
  /* synthetic extended precision constants */
  const double P1 = 4 * 7.85398125648498535156e-01;
  const double P2 = 4 * 3.77489470793079817668e-08;
  const double P3 = 4 * 2.69515142907905952645e-15;
  const double TwoPi = 2*(P1 + P2 + P3);

  const double y = 2*floor(theta/TwoPi);

  double r = ((theta - y*P1) - y*P2) - y*P3;

  if(r > TwoPi) {r = (((r-2*P1)-2*P2)-2*P3); }  /* r-TwoPi */
  else if (r < 0) { /* may happen due to FP rounding */
    r = (((r+2*P1)+2*P2)+2*P3); /* r+TwoPi */
  }

  result->val = r;

  if(fabs(theta) > 0.0625/GSL_DBL_EPSILON) {
    result->val = GSL_NAN;
    result->err = fabs(result->val);
    GSL_ERROR ("error", GSL_ELOSS);
  }
  else if(fabs(theta) > 0.0625/GSL_SQRT_DBL_EPSILON) {
    result->err = GSL_DBL_EPSILON * fabs(result->val - theta);
    return GSL_SUCCESS;
  }
  else {
    double delta = fabs(result->val - theta);
    result->err = 2.0 * GSL_DBL_EPSILON * ((delta < M_PI) ? delta : M_PI);
    return GSL_SUCCESS;
  }
}

int gsl_sf_angle_restrict_symm_e(double * theta)
{
  gsl_sf_result r;
  int stat = gsl_sf_angle_restrict_symm_err_e(*theta, &r);
  *theta = r.val;
  return stat;
}

int gsl_sf_angle_restrict_pos_e(double * theta)
{
  gsl_sf_result r;
  int stat = gsl_sf_angle_restrict_pos_err_e(*theta, &r);
  *theta = r.val;
  return stat;
}

int gsl_sf_sin_err_e(const double x, const double dx, gsl_sf_result * result)
{
  int stat_s = gsl_sf_sin_e(x, result);
  result->err += fabs(cos(x) * dx);
  result->err += GSL_DBL_EPSILON * fabs(result->val);
  return stat_s;
}

int gsl_sf_cos_err_e(const double x, const double dx, gsl_sf_result * result)
{
  int stat_c = gsl_sf_cos_e(x, result);
  result->err += fabs(sin(x) * dx);
  result->err += GSL_DBL_EPSILON * fabs(result->val);
  return stat_c;
}

#if 0
int
gsl_sf_sin_pi_x_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(-100.0 < x && x < 100.0) {
    result->val = sin(M_PI * x) / (M_PI * x);
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    const double N = floor(x + 0.5);
    const double f = x - N;

    if(N < INT_MAX && N > INT_MIN) {
      /* Make it an integer if we can. Saves another
       * call to floor().
       */
      const int intN    = (int)N;
      const double sign = ( GSL_IS_ODD(intN) ? -1.0 : 1.0 );
      result->val = sign * sin(M_PI * f);
      result->err = GSL_DBL_EPSILON * fabs(result->val);
    }
    else if(N > 2.0/GSL_DBL_EPSILON || N < -2.0/GSL_DBL_EPSILON) {
      /* All integer-valued floating point numbers
       * bigger than 2/eps=2^53 are actually even.
       */
      result->val = 0.0;
      result->err = 0.0;
    }
    else {
      const double resN = N - 2.0*floor(0.5*N); /* 0 for even N, 1 for odd N */
      const double sign = ( fabs(resN) > 0.5 ? -1.0 : 1.0 );
      result->val = sign * sin(M_PI*f);
      result->err = GSL_DBL_EPSILON * fabs(result->val);
    }

    return GSL_SUCCESS;
  }
}
#endif

int gsl_sf_sinc_e(double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  {
    const double ax = fabs(x);

    if(ax < 0.8) {
      /* Do not go to the limit of the fit since
       * there is a zero there and the Chebyshev
       * accuracy will go to zero.
       */
      return cheb_eval_e(&sinc_cs, 2.0*ax-1.0, result);
    }
    else if(ax < 100.0) {
      /* Small arguments are no problem.
       * We trust the library sin() to
       * roughly machine precision.
       */
      result->val = sin(M_PI * ax)/(M_PI * ax);
      result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      return GSL_SUCCESS;
    }
    else {
      /* Large arguments must be handled separately.
       */
      const double r = M_PI*ax;
      gsl_sf_result s;
      int stat_s = gsl_sf_sin_e(r, &s);
      result->val = s.val/r;
      result->err = s.err/r + 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      return stat_s;
    }
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_sin(const double x)
{
  EVAL_RESULT(gsl_sf_sin_e(x, &result));
}

double gsl_sf_cos(const double x)
{
  EVAL_RESULT(gsl_sf_cos_e(x, &result));
}

double gsl_sf_hypot(const double x, const double y)
{
  EVAL_RESULT(gsl_sf_hypot_e(x, y, &result));
}

double gsl_sf_lnsinh(const double x)
{
  EVAL_RESULT(gsl_sf_lnsinh_e(x, &result));
}

double gsl_sf_lncosh(const double x)
{
  EVAL_RESULT(gsl_sf_lncosh_e(x, &result));
}

double gsl_sf_angle_restrict_symm(const double theta)
{
  double result = theta;
  EVAL_DOUBLE(gsl_sf_angle_restrict_symm_e(&result));
}

double gsl_sf_angle_restrict_pos(const double theta)
{
  double result = theta;
  EVAL_DOUBLE(gsl_sf_angle_restrict_pos_e(&result));
}

#if 0
double gsl_sf_sin_pi_x(const double x)
{
  EVAL_RESULT(gsl_sf_sin_pi_x_e(x, &result));
}
#endif

double gsl_sf_sinc(const double x)
{
  EVAL_RESULT(gsl_sf_sinc_e(x, &result));
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#ifndef I3SPRNGRANDOMSERVICE_H
#define I3SPRNGRANDOMSERVICE_H

#include "phys-services/I3RandomService.h"

#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_test.h>

#include <string>

/**
 * copyright  (C) 2004
 * the icecube collaboration
 * $Id: I3SPRNGRandomService.h 127790 2015-01-14 04:03:06Z olivas $
 *
 * @brief SPRNG Implementation of the I3RandomService interface.
 * This implementation uses a combination of SPRNG and GSL to generate
 * statistically independent streams of pseudo-random number distributions.
 * See gsl-sprng.h for more details.
 *
 * NB : It's important that you use the same seed for different jobs.  Set
 * nstreams to the number of jobs and use a different streamnum for each job.
 * Otherwise you'll get correlations between the RNG streams.  I know this is
 * counterintuitive, but this is how SPRNG works.
 *
 * The code for this class is based on John Pretz's implementation of
 * I3GSLRandomService.
 *
 * @version $Revision: 127790 $
 * @date $Date: 2015-01-13 21:03:06 -0700 (Tue, 13 Jan 2015) $
 * @author juancarlos
 *
 * @todo Add ability to save state of rng after run is complete
 *  SPRNG has the functions:
 *
 *      int pack_sprng(char *bytes); // returns size of bytes
 *      void unpack_sprng(char bytes[MAX_PACKED_LENGTH]);
 *
 *  which can be used to save and retrieve the state of an rng
 */

class I3SPRNGRandomService : public I3RandomService{
 public:
  /**
   * default constructor
   */
  I3SPRNGRandomService();

  /**
   * constructors
   */
  I3SPRNGRandomService(int seed, int nstreams, int streamnum,
		       std::string instatefile=std::string(),
		       std::string outstatefile=std::string());

  /**
   * destructor
   */
  virtual ~I3SPRNGRandomService();

  /**
   * Binomial distribution
   */
  virtual int Binomial(int ntot, double prob);

  // As with John Pretz's GSL implementation, I have left this out for now.
  /* virtual double BreitWigner(double mean = 0, double gamma = 1)=0; */

  /**
   * Exponential distribution
   */
  virtual double Exp(double tau);

  /**
   * Uniform int distribution with range [0,imax)
   */
  virtual unsigned int Integer(unsigned int imax);

  /**
   * int Poisson distribution
   */
  virtual int Poisson(double mean);

  /**
   * double Poisson distribution
   */
  virtual double PoissonD(double mean);

  /**
   * double uniform distribution with range (0,x1)
   */
  virtual double Uniform(double x1 = 1);

  /**
   * double  uniform distribution with range (x1,x2)
   */
  virtual double Uniform(double x1, double x2);

  /**
   * double Gaussian distribution given mean and StdD
   */
  virtual double Gaus(double mean, double stddev);

  virtual I3FrameObjectPtr GetState() const;
  virtual void RestoreState(I3FrameObjectConstPtr state);

 private:

  // private copy constructors and assignment
  I3SPRNGRandomService(const I3SPRNGRandomService& );
  I3SPRNGRandomService operator=(const I3SPRNGRandomService& );

  gsl_rng* rng_;
  std::string instatefile_;
  std::string outstatefile_;

  int seed_;
  int streamnum_;
  int nstreams_;

  SET_LOGGER("I3SPRNGRandomService");

};

I3_POINTER_TYPEDEFS(I3SPRNGRandomService);

#endif // I3SPRNGRANDOMSERVICE_H


text:
/*
 *  Player - One Hell of a Robot Server
 *  Copyright (C) 2000  Brian Gerkey   &  Kasper Stoy
 *                      gerkey@usc.edu    kaspers@robotics.usc.edu
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This library is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */
/**************************************************************************
 * Desc: Useful pdf functions
 * Author: Andrew Howard
 * Date: 10 Dec 2002
 * CVS: $Id: pf_pdf.h 6345 2008-04-17 01:36:39Z gerkey $
 *************************************************************************/

#ifndef PF_PDF_H
#define PF_PDF_H

#include "nav2d_localizer/pf_vector.h"

//#include <gsl/gsl_rng.h>
//#include <gsl/gsl_randist.h>

#ifdef __cplusplus
extern "C" {
#endif

/**************************************************************************
 * Gaussian
 *************************************************************************/

// Gaussian PDF info
typedef struct
{
  // Mean, covariance and inverse covariance
  pf_vector_t x;
  pf_matrix_t cx;
  //pf_matrix_t cxi;
  double cxdet;

  // Decomposed covariance matrix (rotation * diagonal)
  pf_matrix_t cr;
  pf_vector_t cd;

  // A random number generator
  //gsl_rng *rng;

} pf_pdf_gaussian_t;

// Create a gaussian pdf
pf_pdf_gaussian_t *pf_pdf_gaussian_alloc(pf_vector_t x, pf_matrix_t cx);

// Destroy the pdf
void pf_pdf_gaussian_free(pf_pdf_gaussian_t *pdf);

// Compute the value of the pdf at some point [z].
//double pf_pdf_gaussian_value(pf_pdf_gaussian_t *pdf, pf_vector_t z);

// Draw randomly from a zero-mean Gaussian distribution, with standard
// deviation sigma.
// We use the polar form of the Box-Muller transformation, explained here:
//   http://www.taygeta.com/random/gaussian.html
double pf_ran_gaussian(double sigma);

// Generate a sample from the the pdf.
pf_vector_t pf_pdf_gaussian_sample(pf_pdf_gaussian_t *pdf);

#if 0

/**************************************************************************
 * Discrete
 *************************************************************************/

// Discrete PDF info
typedef struct
{
  // The list of discrete probs
  int prob_count;
  double *probs;

  // A random number generator
  gsl_rng *rng;

  // The discrete prob generator
  gsl_ran_discrete_t *ran;

} pf_pdf_discrete_t;

// Create a discrete pdf
pf_pdf_discrete_t *pf_pdf_discrete_alloc(int count, double *probs);

// Destroy the pdf
void pf_pdf_discrete_free(pf_pdf_discrete_t *pdf);

// Compute the value of the probability of some element [i]
double pf_pdf_discrete_value(pf_pdf_discrete_t *pdf, int i);

// Generate a sample from the the pdf.
int pf_pdf_discrete_sample(pf_pdf_discrete_t *pdf);
#endif

#ifdef __cplusplus
}
#endif

#endif

max_issues_repo_name: pplankton/MRSLAM

max_forks_repo_name: pplankton/MRSLAM


text:
#ifndef PMVS3_OPTIM_H
#define PMVS3_OPTIM_H

#include <vector>
#include "patch.h"
#include <gsl/gsl_multimin.h>

namespace PMVS3 {

class CfindMatch;

class Coptim {
 public:
  Coptim(CfindMatch& findMatch);

  void init(void);

  //-----------------------------------------------------------------
  // Image manipulation
  //-----------------------------------------------------------------
  void collectImages(const int index, std::vector<int>& indexes) const;
  void addImages(Patch::Cpatch& patch) const;
  void removeImagesEdge(Patch::Cpatch& patch) const;

  float getUnit(const int index, const Vec4f& coord) const;

  void computeUnits(const Patch::Cpatch& patch,
                    std::vector<int>& indexes,
                    std::vector<float>& fineness,
                    std::vector<Vec4f>& rays) const;
  void computeUnits(const Patch::Cpatch& patch,
                    std::vector<float>& fineness) const;

  //-----------------------------------------------------------------
  // Optimization
  //-----------------------------------------------------------------

  int preProcess(Patch::Cpatch& patch, const int id, const int seed);
  void refinePatch(Patch::Cpatch& patch, const int id, const int time);

  void refinePatchBFGS(Patch::Cpatch& patch, const int id, const int time);
  void refinePatchBFGS(Patch::Cpatch& patch, const int id, const int time,
                       const int ncc);
  bool refinePatchBFGS2(Patch::Cpatch& patch, const int id, const int time,
                        const int ncc); // LM version
  void refineDepthBFGS(Patch::Cpatch& patch, const int id, const int time,
                       const int ncc);

  int postProcess(Patch::Cpatch& patch, const int id, const int seed);

  void setRefImage(Patch::Cpatch& patch, const int id);

  int check(Patch::Cpatch& patch);

  std::vector<int> m_status;

 protected:
  void filterImagesByAngle(Patch::Cpatch& patch);

  void sortImages(Patch::Cpatch& patch) const;
  void constraintImages(Patch::Cpatch& patch, const float nccThreshold,
                        const int id);
  void setRefConstraintImages(Patch::Cpatch& patch, const float nccThreshold,
                              const int id);

  void setINCCs(const Patch::Cpatch& patch,
                std::vector<float> & nccs,
                const std::vector<int>& indexes,
                const int id, const int robust);

  void setINCCs(const Patch::Cpatch& patch,
                std::vector<std::vector<float> >& nccs,
                const std::vector<int>& indexes,
                const int id, const int robust);

  int grabTex(const Vec4f& coord, const Vec4f& pxaxis, const Vec4f& pyaxis,
              const Vec4f& pzaxis, const int index, const int size,
              std::vector<float>& tex) const;

  int grabSafe(const int index, const int size, const Vec3f& center,
               const Vec3f& dx, const Vec3f& dy, const int level) const;

  double computeSSD(const Vec4f& coord, const Vec4f& normal,
                    const std::vector<int>& indexes, const int id);

  double computeSSD(const Vec4f& coord, const Vec4f& normal,
                    const std::vector<int>& indexes, const Vec4f& pxaxis,
                    const Vec4f& pyaxis, const int id);
  /*
  double computeINCC(const Vec4f& coord, const Vec4f& normal,
                     const std::vector<int>& indexes, const int id,
                     const int robust);
  */
  double computeINCC(const Vec4f& coord, const Vec4f& normal,
                     const std::vector<int>& indexes, const Vec4f& pxaxis,
                     const Vec4f& pyaxis, const int id,
                     const int robust);

 public:
  static void normalize(std::vector<float>& tex);
  static void normalize(std::vector<std::vector<float> >& texs, const int size);

  float dot(const std::vector<float>& tex0, const std::vector<float>& tex1) const;
  float ssd(const std::vector<float>& tex0, const std::vector<float>& tex1) const;
 protected:
  static void lfunc(double* p, double* hx, int m, int n, void* adata);
  void func(int m, int n, double* x, double* fvec, int* iflag, void* arg);

  //BFGS
  static double my_f(const gsl_vector *v, void *params);
  static void my_f_lm(const double *par, int m_dat, const void *data, double *fvec, int *info); // LM version
  static void my_df(const gsl_vector *v, void *params,
                    gsl_vector *df);
  static void my_fdf(const gsl_vector *x, void *params,
                     double *f, gsl_vector *df);
  // for derivative computation
  static double my_f0(double x, void* params);
  static double my_f1(double x, void* params);
  static double my_f2(double x, void* params);
  //----------------------------------------------------------------------
  // For ssd
  static double my_f_ssd(const gsl_vector *v, void *params);
  static void my_df_ssd(const gsl_vector *v, void *params,
                        gsl_vector *df);
  static void my_fdf_ssd(const gsl_vector *x, void *params,
                         double *f, gsl_vector *df);
  // for derivative computation
  static double my_f_ssd0(double x, void* params);
  static double my_f_ssd1(double x, void* params);
  static double my_f_ssd2(double x, void* params);
  //----------------------------------------------------------------------
  // For debugging depth
  static double my_f_depth(const gsl_vector *v, void *params);
  static void my_df_depth(const gsl_vector *v, void *params,
                    gsl_vector *df);
  static void my_fdf_depth(const gsl_vector *x, void *params,
                     double *f, gsl_vector *df);
  // for derivative computation
  static double my_f0_depth(double x, void* params);

  void encode(const Vec4f& coord,
              double* const vect, const int id) const;
  void encode(const Vec4f& coord, const Vec4f& normal,
              double* const vect, const int id) const;
  void decode(Vec4f& coord, Vec4f& normal,
	      const double* const vect, const int id) const;
  void decode(Vec4f& coord, const double* const vect, const int id) const;

 public:
  void setWeightsT(const Patch::Cpatch& patch, const int id);

  double computeINCC(const Vec4f& coord, const Vec4f& normal,
                     const std::vector<int>& indexes, const int id,
                     const int robust);
  void getPAxes(const int index, const Vec4f& coord, const Vec4f& normal,
                Vec4f& pxaxis, Vec4f& pyaxis) const;

  static inline float robustincc(const float rhs) {
    return rhs / (1 + 3 * rhs);
  }

  static inline float unrobustincc(const float rhs) {
    return rhs / (1 - 3 * rhs);
  }

 protected:

  void setAxesScales(void);

  static Coptim* m_one;
  CfindMatch& m_fm;

  //-----------------------------------------------------------------
  // Axes
  std::vector<Vec3f> m_xaxes;
  std::vector<Vec3f> m_yaxes;
  std::vector<Vec3f> m_zaxes;
  // Scales
  std::vector<float> m_ipscales;

  //-----------------------------------------------------------------
  // For threads
  std::vector<float> m_vect0T;
  std::vector<Vec4f> m_centersT;
  std::vector<Vec4f> m_raysT;
  std::vector<std::vector<int> > m_indexesT;
  std::vector<float> m_dscalesT;
  std::vector<float> m_ascalesT;

  // stores current parameters for derivative computation
  std::vector<Vec3f> m_paramsT;

  // Grabbed texture
  std::vector<std::vector<std::vector<float> > > m_texsT; // last is 7x7x3 patch
  // weights for refineDepthOrientationWeighed
  std::vector<std::vector<float> > m_weightsT;
  // Working array for levmar
  std::vector<std::vector<double> > m_worksT;

};
};

#endif // PMVS3_OPTIM_H

max_stars_repo_name: mattjr/structured