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

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_machine.h>
#include <gsl/gsl_histogram2d.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

#define M 107
#define N 239
#define M1 17
#define N1 23
#define MR 10
#define NR 5

void
test2d (void)
{
  double xr[MR + 1] =
    { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 };

  double yr[NR + 1] = { 90.0, 91.0, 92.0, 93.0, 94.0, 95.0 };

  gsl_histogram2d *h, *h1, *g, *hr;
  size_t i, j, k;

  gsl_ieee_env_setup ();

  h = gsl_histogram2d_calloc (M, N);
  h1 = gsl_histogram2d_calloc (M, N);
  g = gsl_histogram2d_calloc (M, N);

  gsl_test (h->xrange == 0,
            "gsl_histogram2d_calloc returns valid xrange pointer");
  gsl_test (h->yrange == 0,
            "gsl_histogram2d_calloc returns valid yrange pointer");
  gsl_test (h->bin == 0, "gsl_histogram2d_calloc returns valid bin pointer");

  gsl_test (h->nx != M, "gsl_histogram2d_calloc returns valid nx");
  gsl_test (h->ny != N, "gsl_histogram2d_calloc returns valid ny");

  hr = gsl_histogram2d_calloc_range (MR, NR, xr, yr);

  gsl_test (hr->xrange == 0,
            "gsl_histogram2d_calloc_range returns valid xrange pointer");
  gsl_test (hr->yrange == 0,
            "gsl_histogram2d_calloc_range returns valid yrange pointer");
  gsl_test (hr->bin == 0,
            "gsl_histogram2d_calloc_range returns valid bin pointer");

  gsl_test (hr->nx != MR, "gsl_histogram2d_calloc_range returns valid nx");
  gsl_test (hr->ny != NR, "gsl_histogram2d_calloc_range returns valid ny");

  {
    int status = 0;
    for (i = 0; i <= MR; i++)
      {
        if (hr->xrange[i] != xr[i])
          {
            status = 1;
          }
      };

    gsl_test (status,
              "gsl_histogram2d_calloc_range creates xrange");
  }

  {
    int status = 0;
    for (i = 0; i <= NR; i++)
      {
        if (hr->yrange[i] != yr[i])
          {
            status = 1;
          }
      };

    gsl_test (status,
              "gsl_histogram2d_calloc_range creates yrange");
  }

  for (i = 0; i <= MR; i++)
    {
      hr->xrange[i] = 0.0;
    }

  for (i = 0; i <= NR; i++)
    {
      hr->yrange[i] = 0.0;
    }

  {
    int status = gsl_histogram2d_set_ranges (hr, xr, MR + 1, yr, NR + 1);

    for (i = 0; i <= MR; i++)
      {
        if (hr->xrange[i] != xr[i])
          {
            status = 1;
          }
      };

    gsl_test (status, "gsl_histogram2d_set_ranges sets xrange");
  }

  {
    int status = 0;
    for (i = 0; i <= NR; i++)
      {
        if (hr->yrange[i] != yr[i])
          {
            status = 1;
          }
      };

    gsl_test (status, "gsl_histogram2d_set_ranges sets yrange");
  }

  k = 0;
  for (i = 0; i < M; i++)
    {
      for (j = 0; j < N; j++)
        {
          k++;
          gsl_histogram2d_accumulate (h, (double) i, (double) j, (double) k);
        };
    }

  {
    int status = 0;
    k = 0;
    for (i = 0; i < M; i++)
      {
        for (j = 0; j < N; j++)
          {
            k++;
            if (h->bin[i * N + j] != (double) k)
              {
                status = 1;
              }
          }
      }
    gsl_test (status,
              "gsl_histogram2d_accumulate writes into array");
  }

  {
    int status = 0;
    k = 0;
    for (i = 0; i < M; i++)
      {
        for (j = 0; j < N; j++)
          {
            k++;
            if (gsl_histogram2d_get (h, i, j) != (double) k)
              status = 1;
          };
      }
    gsl_test (status, "gsl_histogram2d_get reads from array");
  }

  for (i = 0; i <= M; i++)
    {
      h1->xrange[i] = 100.0 + i;
    }

  for (i = 0; i <= N; i++)
    {
      h1->yrange[i] = 900.0 + i * i;
    }

  gsl_histogram2d_memcpy (h1, h);

  {
    int status = 0;
    for (i = 0; i <= M; i++)
      {
        if (h1->xrange[i] != h->xrange[i])
          status = 1;
      };
    gsl_test (status, "gsl_histogram2d_memcpy copies bin xranges");
  }

  {
    int status = 0;
    for (i = 0; i <= N; i++)
      {
        if (h1->yrange[i] != h->yrange[i])
          status = 1;
      };
    gsl_test (status, "gsl_histogram2d_memcpy copies bin yranges");
  }

  {
    int status = 0;
    for (i = 0; i < M; i++)
      {
        for (j = 0; j < N; j++)
          {
            if (gsl_histogram2d_get (h1, i, j) !=
                gsl_histogram2d_get (h, i, j))
              status = 1;
          }
      }
    gsl_test (status, "gsl_histogram2d_memcpy copies bin values");
  }

  gsl_histogram2d_free (h1);

  h1 = gsl_histogram2d_clone (h);

  {
    int status = 0;
    for (i = 0; i <= M; i++)
      {
        if (h1->xrange[i] != h->xrange[i])
          status = 1;
      };
    gsl_test (status, "gsl_histogram2d_clone copies bin xranges");
  }

  {
    int status = 0;
    for (i = 0; i <= N; i++)
      {
        if (h1->yrange[i] != h->yrange[i])
          status = 1;
      };
    gsl_test (status, "gsl_histogram2d_clone copies bin yranges");
  }

  {
    int status = 0;
    for (i = 0; i < M; i++)
      {
        for (j = 0; j < N; j++)
          {
            if (gsl_histogram2d_get (h1, i, j) !=
                gsl_histogram2d_get (h, i, j))
              status = 1;
          }
      }
    gsl_test (status, "gsl_histogram2d_clone copies bin values");
  }

  gsl_histogram2d_reset (h);

  {
    int status = 0;

    for (i = 0; i < M * N; i++)
      {
        if (h->bin[i] != 0)
          status = 1;
      }
    gsl_test (status, "gsl_histogram2d_reset zeros array");
  }

  gsl_histogram2d_free (h);
  h = gsl_histogram2d_calloc (M1, N1);

  {

    int status = 0;

    for (i = 0; i < M1; i++)
      {
        for (j = 0; j < N1; j++)
          {
            gsl_histogram2d_increment (h, (double) i, (double) j);

            for (k = 0; k <= i * N1 + j; k++)
              {
                if (h->bin[k] != 1)
                  {
                    status = 1;
                  }
              }

            for (k = i * N1 + j + 1; k < M1 * N1; k++)
              {
                if (h->bin[k] != 0)
                  {
                    status = 1;
                  }
              }
          }
      }
    gsl_test (status, "gsl_histogram2d_increment increases bin value");
  }

  gsl_histogram2d_free (h);
  h = gsl_histogram2d_calloc (M, N);

  {
    int status = 0;
    for (i = 0; i < M; i++)
      {
        double x0 = 0, x1 = 0;
        gsl_histogram2d_get_xrange (h, i, &x0, &x1);

        if (x0 != i || x1 != i + 1)
          {
            status = 1;
          }
      }
    gsl_test (status,
              "gsl_histogram2d_get_xlowerlimit and xupperlimit");
  }

  {
    int status = 0;
    for (i = 0; i < N; i++)
      {
        double y0 = 0, y1 = 0;
        gsl_histogram2d_get_yrange (h, i, &y0, &y1);

        if (y0 != i || y1 != i + 1)
          {
            status = 1;
          }
      }
    gsl_test (status,
              "gsl_histogram2d_get_ylowerlimit and yupperlimit");
  }

  {
    int status = 0;
    if (gsl_histogram2d_xmax (h) != M)
      status = 1;
    gsl_test (status, "gsl_histogram2d_xmax");
  }

  {
    int status = 0;
    if (gsl_histogram2d_xmin (h) != 0)
      status = 1;
    gsl_test (status, "gsl_histogram2d_xmin");
  }

  {
    int status = 0;
    if (gsl_histogram2d_nx (h) != M)
      status = 1;
    gsl_test (status, "gsl_histogram2d_nx");
  }

  {
    int status = 0;
    if (gsl_histogram2d_ymax (h) != N)
      status = 1;
    gsl_test (status, "gsl_histogram2d_ymax");
  }

  {
    int status = 0;
    if (gsl_histogram2d_ymin (h) != 0)
      status = 1;
    gsl_test (status, "gsl_histogram2d_ymin");
  }

  {
    int status = 0;
    if (gsl_histogram2d_ny (h) != N)
      status = 1;
    gsl_test (status, "gsl_histogram2d_ny");
  }

  h->bin[3 * N + 2] = 123456.0;
  h->bin[4 * N + 3] = -654321;

  {
    double max = gsl_histogram2d_max_val (h);
    gsl_test (max != 123456.0, "gsl_histogram2d_max_val finds maximum value");
  }

  {
    double min = gsl_histogram2d_min_val (h);
    gsl_test (min != -654321.0,
              "gsl_histogram2d_min_val finds minimum value");
  }

  {
    size_t imax, jmax;
    gsl_histogram2d_max_bin (h, &imax, &jmax);
    gsl_test (imax != 3
              || jmax != 2,
              "gsl_histogram2d_max_bin finds maximum value bin");
  }

  {
    size_t imin, jmin;
    gsl_histogram2d_min_bin (h, &imin, &jmin);
    gsl_test (imin != 4
              || jmin != 3, "gsl_histogram2d_min_bin find minimum value bin");
  }

  for (i = 0; i < M * N; i++)
    {
      h->bin[i] = i + 27;
      g->bin[i] = (i + 27) * (i + 1);
    }

  {
    double sum = gsl_histogram2d_sum (h);
    gsl_test (sum != N * M * 27 + ((N * M - 1) * N * M) / 2,
              "gsl_histogram2d_sum sums all bin values");
  }

  {
    /* first test... */
    const double xpos = 0.6;
    const double ypos = 0.85;
    double xmean;
    double ymean;
    size_t xbin;
    size_t ybin;
    gsl_histogram2d *h3 = gsl_histogram2d_alloc (M, N);
    gsl_histogram2d_set_ranges_uniform (h3, 0, 1, 0, 1);
    gsl_histogram2d_increment (h3, xpos, ypos);
    gsl_histogram2d_find (h3, xpos, ypos, &xbin, &ybin);
    xmean = gsl_histogram2d_xmean (h3);
    ymean = gsl_histogram2d_ymean (h3);

    {
      double expected_xmean = (h3->xrange[xbin] + h3->xrange[xbin + 1]) / 2.0;
      double expected_ymean = (h3->yrange[ybin] + h3->yrange[ybin + 1]) / 2.0;
      gsl_test_abs (xmean, expected_xmean, 100.0 * GSL_DBL_EPSILON,
                    "gsl_histogram2d_xmean");
      gsl_test_abs (ymean, expected_ymean, 100.0 * GSL_DBL_EPSILON,
                    "gsl_histogram2d_ymean");
    };
    gsl_histogram2d_free (h3);
  }

  {
    /* test it with bivariate normal distribution */
    const double xmean = 0.7;
    const double ymean = 0.7;
    const double xsigma = 0.1;
    const double ysigma = 0.1;
    const double correl = 0.5;
    const double norm =
      10.0 / M_PI / xsigma / ysigma / sqrt (1.0 - correl * correl);
    size_t xbin;
    size_t ybin;
    gsl_histogram2d *h3 = gsl_histogram2d_alloc (M, N);
    gsl_histogram2d_set_ranges_uniform (h3, 0, 1, 0, 1);
    /* initialize with 2d gauss pdf in two directions */
    for (xbin = 0; xbin < M; xbin++)
      {
        double xi =
          ((h3->xrange[xbin] + h3->xrange[xbin + 1]) / 2.0 - xmean) / xsigma;
        for (ybin = 0; ybin < N; ybin++)
          {
            double yi =
              ((h3->yrange[ybin] + h3->yrange[ybin + 1]) / 2.0 -
               ymean) / ysigma;
            double prob =
              norm * exp (-(xi * xi - 2.0 * correl * xi * yi + yi * yi) /
                          2.0 / (1 - correl * correl));
            h3->bin[xbin * N + ybin] = prob;
          }
      }
    {
      double xs = gsl_histogram2d_xsigma (h3);
      double ys = gsl_histogram2d_ysigma (h3);
      /* evaluate results and compare with parameters */

      gsl_test_abs (gsl_histogram2d_xmean (h3), xmean, 2.0/M,
                    "gsl_histogram2d_xmean histogram mean(x)");
      gsl_test_abs (gsl_histogram2d_ymean (h3), ymean, 2.0/N,
                    "gsl_histogram2d_ymean histogram mean(y)");
      gsl_test_abs (xs, xsigma, 2.0/M,
                    "gsl_histogram2d_xsigma histogram stdev(x)");
      gsl_test_abs (ys, ysigma, 2.0/N,
                    "gsl_histogram2d_ysigma histogram stdev(y)");
      gsl_test_abs (gsl_histogram2d_cov (h3) / xs / ys, correl,
                    2.0/((M < N) ? M : N),
                    "gsl_histogram2d_cov histogram covariance");
    }
    gsl_histogram2d_free (h3);
  }

  gsl_histogram2d_memcpy (h1, g);
  gsl_histogram2d_add (h1, h);

  {
    int status = 0;
    for (i = 0; i < M * N; i++)
      {
        if (h1->bin[i] != g->bin[i] + h->bin[i])
          status = 1;
      }
    gsl_test (status, "gsl_histogram2d_add histogram addition");
  }

  gsl_histogram2d_memcpy (h1, g);
  gsl_histogram2d_sub (h1, h);

  {
    int status = 0;
    for (i = 0; i < M * N; i++)
      {
        if (h1->bin[i] != g->bin[i] - h->bin[i])
          status = 1;
      }
    gsl_test (status, "gsl_histogram2d_sub histogram subtraction");
  }

  gsl_histogram2d_memcpy (h1, g);
  gsl_histogram2d_mul (h1, h);

  {
    int status = 0;
    for (i = 0; i < M * N; i++)
      {
        if (h1->bin[i] != g->bin[i] * h->bin[i])
          status = 1;
      }
    gsl_test (status, "gsl_histogram2d_mul histogram multiplication");
  }

  gsl_histogram2d_memcpy (h1, g);
  gsl_histogram2d_div (h1, h);

  {
    int status = 0;
    for (i = 0; i < M * N; i++)
      {
        if (h1->bin[i] != g->bin[i] / h->bin[i])
          status = 1;
      }
    gsl_test (status, "gsl_histogram2d_div histogram division");
  }

  gsl_histogram2d_memcpy (h1, g);
  gsl_histogram2d_scale (h1, 0.5);

  {
    int status = 0;
    for (i = 0; i < M * N; i++)
      {
        if (h1->bin[i] != 0.5 * g->bin[i])
          status = 1;
      }
    gsl_test (status, "gsl_histogram2d_scale histogram scaling");
  }

  gsl_histogram2d_memcpy (h1, g);
  gsl_histogram2d_shift (h1, 0.25);

  {
    int status = 0;
    for (i = 0; i < M * N; i++)
      {
        if (h1->bin[i] != 0.25 + g->bin[i])
          status = 1;
      }
    gsl_test (status, "gsl_histogram2d_shift histogram shift");
  }

  gsl_histogram2d_free (h);     /* free whatever is in h */

  h = gsl_histogram2d_calloc_uniform (M1, N1, 0.0, 5.0, 0.0, 5.0);

  gsl_test (h->xrange == 0,
            "gsl_histogram2d_calloc_uniform returns valid range pointer");
  gsl_test (h->yrange == 0,
            "gsl_histogram2d_calloc_uniform returns valid range pointer");
  gsl_test (h->bin == 0,
            "gsl_histogram2d_calloc_uniform returns valid bin pointer");
  gsl_test (h->nx != M1, "gsl_histogram2d_calloc_uniform returns valid nx");
  gsl_test (h->ny != N1, "gsl_histogram2d_calloc_uniform returns valid ny");

  gsl_histogram2d_accumulate (h, 0.0, 3.01, 1.0);
  gsl_histogram2d_accumulate (h, 0.1, 2.01, 2.0);
  gsl_histogram2d_accumulate (h, 0.2, 1.01, 3.0);
  gsl_histogram2d_accumulate (h, 0.3, 0.01, 4.0);

  {
    size_t i1, i2, i3, i4;
    size_t j1, j2, j3, j4;
    double expected;
    int status;
    status = gsl_histogram2d_find (h, 0.0, 3.01, &i1, &j1);
    status = gsl_histogram2d_find (h, 0.1, 2.01, &i2, &j2);
    status = gsl_histogram2d_find (h, 0.2, 1.01, &i3, &j3);
    status = gsl_histogram2d_find (h, 0.3, 0.01, &i4, &j4);

    for (i = 0; i < M1; i++)
      {
        for (j = 0; j < N1; j++)
          {
            if (i == i1 && j == j1)
              {
                expected = 1.0;
              }
            else if (i == i2 && j == j2)
              {
                expected = 2.0;
              }
            else if (i == i3 && j == j3)
              {
                expected = 3.0;
              }
            else if (i == i4 && j == j4)
              {
                expected = 4.0;
              }
            else
              {
                expected = 0.0;
              }

            if (h->bin[i * N1 + j] != expected)
              {
                status = 1;
              }
          }
      }
    gsl_test (status, "gsl_histogram2d_find returns index");
  }

  {
    FILE *f = fopen ("test.txt", "w");
    gsl_histogram2d_fprintf (f, h, "%.19e", "%.19e");
    fclose (f);
  }

  {
    FILE *f = fopen ("test.txt", "r");
    gsl_histogram2d *hh = gsl_histogram2d_calloc (M1, N1);
    int status = 0;

    gsl_histogram2d_fscanf (f, hh);

    for (i = 0; i <= M1; i++)
      {
        if (h->xrange[i] != hh->xrange[i])
          {
            printf ("xrange[%d] : %g orig vs %g\n",
                    (int) i, h->xrange[i], hh->xrange[i]);
            status = 1;
          }
      }

    for (j = 0; j <= N1; j++)
      {
        if (h->yrange[j] != hh->yrange[j])
          {
            printf ("yrange[%d] : %g orig vs %g\n",
                    (int) j, h->yrange[j], hh->yrange[j]);
            status = 1;
          }
      }

    for (i = 0; i < M1 * N1; i++)
      {
        if (h->bin[i] != hh->bin[i])
          {
            printf ("bin[%d] : %g orig vs %g\n",
                    (int) i, h->bin[i], hh->bin[i]);
            status = 1;
          }
      }

    gsl_test (status, "gsl_histogram2d_fprintf and fscanf");

    gsl_histogram2d_free (hh);
    fclose (f);
  }

  {
    FILE *f = fopen ("test.dat", "wb");
    gsl_histogram2d_fwrite (f, h);
    fclose (f);
  }

  {
    FILE *f = fopen ("test.dat", "rb");
    gsl_histogram2d *hh = gsl_histogram2d_calloc (M1, N1);
    int status = 0;

    gsl_histogram2d_fread (f, hh);

    for (i = 0; i <= M1; i++)
      {
        if (h->xrange[i] != hh->xrange[i])
          {
            printf ("xrange[%d] : %g orig vs %g\n",
                    (int) i, h->xrange[i], hh->xrange[i]);
            status = 1;
          }
      }

    for (j = 0; j <= N1; j++)
      {
        if (h->yrange[j] != hh->yrange[j])
          {
            printf ("yrange[%d] : %g orig vs %g\n",
                    (int) j, h->yrange[j], hh->yrange[j]);
            status = 1;
          }
      }

    for (i = 0; i < M1 * N1; i++)
      {
        if (h->bin[i] != hh->bin[i])
          {
            printf ("bin[%d] : %g orig vs %g\n",
                    (int) i, h->bin[i], hh->bin[i]);
            status = 1;
          }
      }

    gsl_test (status, "gsl_histogram2d_fwrite and fread");

    gsl_histogram2d_free (hh);
    fclose (f);
  }

  gsl_histogram2d_free (h);
  gsl_histogram2d_free (h1);
  gsl_histogram2d_free (g);
  gsl_histogram2d_free (hr);
}

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
#pragma once

#ifndef __CONSTANTS_H__
#define __CONSTANTS_H__

/***************************************************
 *************** Auto-Generated File ***************
 ***************************************************/

#include <cstring>
#include <string>

#include <gsl/gsl_math.h>

#include "../../Utils/BinarySearch.h"

constexpr const int numConstants = 4;
constexpr const char* constants[numConstants] = {
    "e", "gamma", "pi", "vphi"
};
constexpr const double constantValues[numConstants] = {
    M_E, M_EULER, M_PI, 1.618033988749895
};
const std::string shortConstants[numConstants] = {
    "e", "γ", "π", "ϕ"
};
constexpr const char * constantLongValues[numConstants] = {
    "2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932",
	"0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495146314472498",
	"3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328",
	"1.61803398874989484820458683436563811772030917980576286213544862270526046281890244970720720418939113748475408807538"
};
constexpr const int longestConstantName = 5;
constexpr const int shortestConstantName = 1;

/*
Returns index of the constant in the constants array.

Uses binary search under the hood to search for the index.

Parameters
----------
name: The name of the constant

Returns
-------
The index or -1 if the provided name is not a constant.
*/
CONSTEXPR_BINARY_SEARCH(getConstantIndex, constants, numConstants)

/*
Returns the value of the constant at the provided index.

Parameters
----------
index: the index of the constant name

Returns
-------
If the index is valid, it will return the number of arguments.
If the index is not valid, the function will return 0.
Note: a return value of -1 means that it accepts a variable number of parameters.
*/
static double getConstantValue(const std::string& name){
    int index = getConstantIndex(name.c_str());
    if (index != -1){
        return constantValues[index];
    }
    return GSL_NAN;
}
constexpr double getConstantValue(int index){
    if (index >= 0 && index < numConstants){
        return constantValues[index];
    }
    return GSL_NAN;
}

#endif // __CONSTANTS_H__

ext: h

max_stars_repo_name: antoniojkim/CalcPlusPlus

max_issues_repo_name: antoniojkim/CalcPlusPlus

max_forks_repo_name: antoniojkim/CalcPlusPlus


text:
#pragma once

#include "rev/gl/Context.h"
#include "rev/gl/Resource.h"
#include "rev/gl/Uniform.h"

#include <array>
#include <gsl/gsl_assert>
#include <string>

namespace rev {

template <GLenum shaderType>
GLuint createShader()
{
    return glCreateShader(shaderType);
}

namespace detail {
    template <void (*propertyGetter)(GLuint, GLenum, GLint*),
        void (*logGetter)(GLuint, GLsizei, GLsizei*, GLchar*)>
    std::string extractLog(GLuint objectId)
    {
        GLint logLength;
        propertyGetter(objectId, GL_INFO_LOG_LENGTH, &logLength);

        std::string log;
        log.resize(logLength);

        GLint fetchedLength;
        logGetter(objectId, logLength, &fetchedLength, log.data());

        log.resize(fetchedLength);

        return log;
    }
} // namespace detail

template <GLenum shaderType>
class Shader : public Resource<createShader<shaderType>, gl::deleteShader> {
public:
    void setSource(const std::string_view& source)
    {
        setSource(std::array<std::string_view, 1>{ source });
    }

    template <size_t arrayLength>
    void setSource(const std::array<std::string_view, arrayLength>& source)
    {
        std::array<const char*, arrayLength> pointers;
        std::array<GLint, arrayLength> lengths;
        for (size_t i = 0; i < arrayLength; i++) {
            pointers[i] = source[i].data();
            lengths[i] = static_cast<GLint>(source[i].size());
            Expects(source[i].size() <= std::numeric_limits<GLint>::max());
        }
        glShaderSource(this->getId(), arrayLength, pointers.data(), lengths.data());
    }

    void compile() { glCompileShader(this->getId()); }

    bool getCompileStatus()
    {
        GLint status;
        glGetShaderiv(this->getId(), GL_COMPILE_STATUS, &status);
        return (status == GL_TRUE);
    }

    std::string getCompileLog()
    {
        return detail::extractLog<gl::getShaderiv, gl::getShaderInfoLog>(this->getId());
    }
};

using VertexShader = Shader<GL_VERTEX_SHADER>;
using FragmentShader = Shader<GL_FRAGMENT_SHADER>;

class ProgramResource : public Resource<gl::createProgram, gl::deleteProgram> {
public:
    template <GLenum shaderType>
    void attachShader(const Shader<shaderType>& shader)
    {
        glAttachShader(getId(), shader.getId());
    }

    void link() { glLinkProgram(getId()); }

    bool getLinkStatus()
    {
        GLint status;
        glGetProgramiv(getId(), GL_LINK_STATUS, &status);
        return (status == GL_TRUE);
    }

    std::string getLinkLog()
    {
        return detail::extractLog<gl::getProgramiv, gl::getProgramInfoLog>(getId());
    }

    template <typename VertexSourceType, typename FragmentSourceType>
    void buildWithSource(
        const VertexSourceType& vertexSource, const FragmentSourceType& fragmentSource)
    {
        VertexShader vShader;
        vShader.setSource(vertexSource);
        vShader.compile();
        if (!vShader.getCompileStatus()) {
            throw vShader.getCompileLog();
        }

        FragmentShader fShader;
        fShader.setSource(fragmentSource);
        fShader.compile();
        if (!fShader.getCompileStatus()) {
            throw fShader.getCompileLog();
        }

        attachShader(vShader);
        attachShader(fShader);
        link();
        if (!getLinkStatus()) {
            throw getLinkLog();
        }
    }

    template <typename VariableType>
    Uniform<VariableType> getUniform(const char* name)
    {
        GLint location = glGetUniformLocation(getId(), name);

        return Uniform<VariableType>(location);
    }
};

using ProgramContext = ResourceContext<ProgramResource, gl::useProgram>;

} // namespace rev

max_stars_repo_name: eyebrowsoffire/rev

max_issues_repo_name: eyebrowsoffire/rev

max_forks_repo_name: eyebrowsoffire/rev


text:
#include <iostream>
#include <math.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_wavelet2d.h>

#include <cgi/mosaic/Mosaic.hpp>

void normalizeBySum(double* data, int n)
{
  int sum = 0;
  for(int i=0; i<n; ++i)
    sum += data[i];

  for(int i=0; i<n; ++i)
    data[i] /= sum;
}

int
main (int argc, char **argv)
{

  const int n = 1024;

  if(argc < 2)
    {
      std::cout << "ERROR" << std::endl;
      return 1;
    }

  cgi::Mosaic m;
  if(cgi::NoError != m.open(argv[1]))
    {
      std::cout << "Failed to open" << std::endl;
      return 2;
    }

  m.setTileSize(cgi::core::Size2d<int>(n,n));

  double* data = new double [n * n];
  if(cgi::NoError != m.getTile<double>(0, 0, data))
    {
      std::cout << "getTile failed" << std::endl;
      return 3;
    }

  gsl_wavelet *w;
  gsl_wavelet_workspace *work;

  w = gsl_wavelet_alloc (gsl_wavelet_daubechies, 4);
  work = gsl_wavelet_workspace_alloc (n);

  if(GSL_SUCCESS != gsl_wavelet2d_transform_forward (w, data, n, n, n, work))
    {
      std::cout << "Transform failed" << std::endl;
      return 4;
    }

  /////////////////////////////////////////
  cgi::Mosaic moDwt;
  moDwt.create("dwt.tif", "GTiff", cgi::core::Size2d<int>(n,n), 1, cgi::Depth64F);
  moDwt.setTileSize(cgi::core::Size2d<int>(n,n));
  moDwt.putTile<double>(data, 0, 0);
  moDwt.close();

  /////////////////////////////////////////
  // normalization
  normalizeBySum(data, n);

  /////////////////////////////////////////
  // do the inverse and reconstruct the image

  if(GSL_SUCCESS != gsl_wavelet2d_transform_inverse (w, data, n, n, n, work))
    {
      std::cout << "Transform failed" << std::endl;
      return 5;
    }

  /////////////////////////////////////////
  // convert the reconstructed image to unsigned short
  unsigned short* recData = new unsigned short [n * n];
  int i;
  for(i=0; i<(n*n); ++i)
    {
      recData[i] = static_cast<unsigned short>(std::max(std::min(data[i], 65535.0), 0.0));
    }

  cgi::Mosaic moRec;
  moRec.create("rec.tif", "GTiff", cgi::core::Size2d<int>(n,n), 1, cgi::Depth16U);
  moRec.setTileSize(cgi::core::Size2d<int>(n,n));
  moRec.putTile<unsigned short>(recData, 0, 0);
  moRec.close();

  /////////////////////////////////////////
  gsl_wavelet_free (w);
  gsl_wavelet_workspace_free (work);

  delete [] recData;
  delete [] data;

  return 0;
}

ext: c

max_stars_repo_name: klaricmn/snippets

max_issues_repo_name: klaricmn/snippets

max_forks_repo_name: klaricmn/snippets


text:
#ifndef _CALCLAMBDA_CHISQ_DIST
#define _CALCLAMBDA_CHISQ_DIST

#include <gsl/gsl_blas.h>

#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION

#define MIN_EIGENVAL 1E-6
#define SIGINT_DEFAULT 0.001

struct chisq_dist {
    int mode;
    double *covmat;
    double *c;
    double *slope;
    double *pivotmag;
    double *refmag;
    double *refmagerr;
    double *magerr;
    double *color;
    double *col_err_ratio;
    double *lupcorr;
    int ncalc;
    int ngal;
    int nz;
    int ncol;

    double sigint;
    //    long do_chisq;
    //long nophotoerr;
};

int chisq_dist(int mode, int do_chisq, int nophotoerr, int ncalc, int ncol, double *covmat, double *c, double *slope, double *pivotmag, double *refmag, double *refmagerr, double *magerr, double *color, double *lupcorr, double *dist, double sigint);

int check_and_fix_covmat(gsl_matrix *covmat);

#endif

max_stars_repo_name: erykoff/redmapper

max_issues_repo_name: erykoff/redmapper

max_forks_repo_name: erykoff/redmapper


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

#define A(m,n) BLKADDR( A, double, m, n )

/***************************************************************************//**
 *
 * @ingroup CORE_double
 *
 * CORE_dgeqp3_pivot finds next pivot, pvt, based on maximum column norm.
 * It applies the swap to the matrices A, F, and vectors jpvt, norms1, norms2.
 * If info != 0, it returns immediately, doing no work.
 *
 *******************************************************************************
 *
 * @param[in,out] A
 *         On entry, descriptor for m by n matrix A.
 *         On exit, column k of jj-th block column is swapped with column pvt.
 *
 * @param[in,out] F
 *         On entry, n by nb matrix F.
 *         On exit, row k is swapped with row pvt - jj*nb.
 *         Currently, F is stored column-wise, not tile-wise.
 *
 * @param[in] ldf
 *         Leading dimension of F. ldf >= max(1,A.n).
 *
 * @param[in] jj
 *         Index of current block column, 0 <= jj < A.nt.
 *
 * @param[in] k
 *         Index of current column within block column, 0 <= k < A.nb.
 *
 * @param[in,out] jpvt
 *         Permutation vector, dimension n.
 *         On exit, swaps entries jpvt[k+jj*nb] and jpvt[pvt].
 *
 * @param[in,out] norms1
 *         On entry, vector of partial column norms, dimension n.
 *         On exit, sets norms1[pvt] = norms1[k+jj*nb].
 *
 * @param[in,out] norms2
 *         On entry, vector of original column norms, dimension n.
 *         On exit, sets norms2[pvt] = norms2[k+jj*nb].
 *
 * @param[in] info
 *         Error code from dgeqp3_update; zero if no error.
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dgeqp3_pivot = PCORE_dgeqp3_pivot
#define CORE_dgeqp3_pivot PCORE_dgeqp3_pivot
#endif
void CORE_dgeqp3_pivot( PLASMA_desc A, double *F, int ldf,
                        int jj, int k, int *jpvt,
                        double *norms1, double *norms2, int *info )
{
    double *Aij, *Aip;
    int pvt, ii, pp, p, mb, lda, tmp, jk;

    /* since pivot depends on all the norm updates, check info here
     * to detect errors from dgeqp3_updates. */
    if ( *info != 0 )
        return;

    /* jk and pvt are indices in global matrix */
    jk = jj*A.nb + k;
    pvt = jk + cblas_idamax( A.n - jk, &norms1[jk], 1 );
    if ( pvt != jk ) {
        tmp       = jpvt[jk];
        jpvt[jk]  = jpvt[pvt];
        jpvt[pvt] = tmp;

        norms1[pvt] = norms1[jk];  /* don't need to save norms[pvt] */
        norms2[pvt] = norms2[jk];

        cblas_dswap( A.nb, &F[k], ldf, &F[pvt-jj*A.nb], ldf );  /* rows k and pvt */

        pp = pvt / A.nb;  /* tile containing pvt  */
        p  = pvt % A.nb;  /* index within pp tile */
        for( ii = 0; ii < A.mt; ++ii ) {
            mb  = min( A.mb, A.m - ii*A.mb );
            lda = BLKLDD( A, ii );
            Aij = A(ii,jj);
            Aip = A(ii,pp);
            cblas_dswap( mb, &Aij[k*lda], 1, &Aip[p*lda], 1 );  /* cols k and pvt */
        }
    }
}


text:
	#pragma once

#include <DirectXMath.h>
#include <d3d11.h>
#include <gsl\gsl>

namespace Library
{
	class Mesh;

	template <typename T>
	struct VertexDeclaration
	{
		static constexpr uint32_t VertexSize() { return gsl::narrow_cast<uint32_t>(sizeof(T)); }
		static constexpr uint32_t VertexBufferByteWidth(size_t vertexCount) { return gsl::narrow_cast<uint32_t>(sizeof(T) * vertexCount); }
		static void CreateVertexBuffer(gsl::not_null<ID3D11Device*> device, const gsl::span<const T>& vertices, gsl::not_null<ID3D11Buffer**> vertexBuffer);
	};

	class VertexPosition : public VertexDeclaration<VertexPosition>
	{
	private:
		inline static const D3D11_INPUT_ELEMENT_DESC _InputElements[]
		{
			{ "POSITION", 0, DXGI_FORMAT_R32G32B32A32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
		};

	public:
		VertexPosition() = default;

		VertexPosition(const DirectX::XMFLOAT4& position) :
			Position(position) { }

		DirectX::XMFLOAT4 Position;

		inline static const gsl::span<const D3D11_INPUT_ELEMENT_DESC> InputElements { _InputElements };

		static void CreateVertexBuffer(gsl::not_null<ID3D11Device*> device, const Library::Mesh& mesh, gsl::not_null<ID3D11Buffer**> vertexBuffer);
		static void CreateVertexBuffer(gsl::not_null<ID3D11Device*> device, const gsl::span<const VertexPosition>& vertices, gsl::not_null<ID3D11Buffer**> vertexBuffer)
		{
			VertexDeclaration::CreateVertexBuffer(device, vertices, vertexBuffer);
		}
	};

	class VertexPositionColor : public VertexDeclaration<VertexPositionColor>
	{
	private:
		inline static const D3D11_INPUT_ELEMENT_DESC _InputElements[]
		{
			{ "POSITION", 0, DXGI_FORMAT_R32G32B32A32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
			{ "COLOR", 0, DXGI_FORMAT_R32G32B32A32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
		};

	public:
		VertexPositionColor() = default;

		VertexPositionColor(const DirectX::XMFLOAT4& position, const DirectX::XMFLOAT4& color) :
			Position(position), Color(color) { }

		DirectX::XMFLOAT4 Position;
		DirectX::XMFLOAT4 Color;

		inline static const gsl::span<const D3D11_INPUT_ELEMENT_DESC> InputElements{ _InputElements };

		static void CreateVertexBuffer(gsl::not_null<ID3D11Device*> device, const Library::Mesh& mesh, gsl::not_null<ID3D11Buffer**> vertexBuffer);
		static void CreateVertexBuffer(gsl::not_null<ID3D11Device*> device, const gsl::span<const VertexPositionColor>& vertices, gsl::not_null<ID3D11Buffer**> vertexBuffer)
		{
			VertexDeclaration::CreateVertexBuffer(device, vertices, vertexBuffer);
		}
	};

	class VertexPositionTexture : public VertexDeclaration<VertexPositionTexture>
	{
	private:
		inline static const D3D11_INPUT_ELEMENT_DESC _InputElements[]
		{
			{ "POSITION", 0, DXGI_FORMAT_R32G32B32A32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
			{ "TEXCOORD", 0, DXGI_FORMAT_R32G32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
		};
	public:
		VertexPositionTexture() = default;

		VertexPositionTexture(const DirectX::XMFLOAT4& position, const DirectX::XMFLOAT2& textureCoordinates) :
			Position(position), TextureCoordinates(textureCoordinates) { }

		DirectX::XMFLOAT4 Position;
		DirectX::XMFLOAT2 TextureCoordinates;

		inline static const gsl::span<const D3D11_INPUT_ELEMENT_DESC> InputElements{ _InputElements };

		static void CreateVertexBuffer(gsl::not_null<ID3D11Device*> device, const Library::Mesh& mesh, gsl::not_null<ID3D11Buffer**> vertexBuffer);
		static void CreateVertexBuffer(gsl::not_null<ID3D11Device*> device, const gsl::span<const VertexPositionTexture>& vertices, gsl::not_null<ID3D11Buffer**> vertexBuffer)
		{
			VertexDeclaration::CreateVertexBuffer(device, vertices, vertexBuffer);
		}
	};

	class VertexPositionSize : public VertexDeclaration<VertexPositionSize>
	{
	private:
		inline static const D3D11_INPUT_ELEMENT_DESC _InputElements[]
		{
			{ "POSITION", 0, DXGI_FORMAT_R32G32B32A32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
			{ "SIZE", 0, DXGI_FORMAT_R32G32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
		};

	public:
		VertexPositionSize() = default;

		VertexPositionSize(const DirectX::XMFLOAT4& position, const DirectX::XMFLOAT2& size) :
			Position(position), Size(size) { }

		DirectX::XMFLOAT4 Position;
		DirectX::XMFLOAT2 Size;

		inline static const gsl::span<const D3D11_INPUT_ELEMENT_DESC> InputElements{ _InputElements };

		static void CreateVertexBuffer(gsl::not_null<ID3D11Device*> device, const gsl::span<const VertexPositionSize>& vertices, gsl::not_null<ID3D11Buffer**> vertexBuffer)
		{
			VertexDeclaration::CreateVertexBuffer(device, vertices, vertexBuffer);
		}
	};

	class VertexPositionNormal : public VertexDeclaration<VertexPositionNormal>
	{
	private:
		inline static const D3D11_INPUT_ELEMENT_DESC _InputElements[]
		{
			{ "POSITION", 0, DXGI_FORMAT_R32G32B32A32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
			{ "NORMAL", 0, DXGI_FORMAT_R32G32B32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
		};

	public:
		VertexPositionNormal() = default;

		VertexPositionNormal(const DirectX::XMFLOAT4& position, const DirectX::XMFLOAT3& normal) :
			Position(position), Normal(normal) { }

		DirectX::XMFLOAT4 Position;
		DirectX::XMFLOAT3 Normal;

		inline static const gsl::span<const D3D11_INPUT_ELEMENT_DESC> InputElements{ _InputElements };

		static void CreateVertexBuffer(gsl::not_null<ID3D11Device*> device, const Library::Mesh& mesh, gsl::not_null<ID3D11Buffer**> vertexBuffer);
		static void CreateVertexBuffer(gsl::not_null<ID3D11Device*> device, const gsl::span<const VertexPositionNormal>& vertices, gsl::not_null<ID3D11Buffer**> vertexBuffer)
		{
			VertexDeclaration::CreateVertexBuffer(device, vertices, vertexBuffer);
		}
	};

	class VertexPositionTextureNormal : public VertexDeclaration<VertexPositionTextureNormal>
	{
	private:
		inline static const D3D11_INPUT_ELEMENT_DESC _InputElements[]
		{
			{ "POSITION", 0, DXGI_FORMAT_R32G32B32A32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
			{ "TEXCOORD", 0, DXGI_FORMAT_R32G32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
			{ "NORMAL", 0, DXGI_FORMAT_R32G32B32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
		};

	public:
		VertexPositionTextureNormal() = default;

		VertexPositionTextureNormal(const DirectX::XMFLOAT4& position, const DirectX::XMFLOAT2& textureCoordinates, const DirectX::XMFLOAT3& normal) :
			Position(position), TextureCoordinates(textureCoordinates), Normal(normal) { }

		DirectX::XMFLOAT4 Position;
		DirectX::XMFLOAT2 TextureCoordinates;
		DirectX::XMFLOAT3 Normal;

		inline static const gsl::span<const D3D11_INPUT_ELEMENT_DESC> InputElements{ _InputElements };

		static void CreateVertexBuffer(gsl::not_null<ID3D11Device*> device, const Library::Mesh& mesh, gsl::not_null<ID3D11Buffer**> vertexBuffer);
		static void CreateVertexBuffer(gsl::not_null<ID3D11Device*> device, const gsl::span<const VertexPositionTextureNormal>& vertices, gsl::not_null<ID3D11Buffer**> vertexBuffer)
		{
			VertexDeclaration::CreateVertexBuffer(device, vertices, vertexBuffer);
		}
	};

	class VertexPositionTextureNormalTangent : public VertexDeclaration<VertexPositionTextureNormalTangent>
	{
	private:
		inline static const D3D11_INPUT_ELEMENT_DESC _InputElements[]
		{
			{ "POSITION", 0, DXGI_FORMAT_R32G32B32A32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
			{ "TEXCOORD", 0, DXGI_FORMAT_R32G32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
			{ "NORMAL", 0, DXGI_FORMAT_R32G32B32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
			{ "TANGENT", 0, DXGI_FORMAT_R32G32B32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
		};

	public:
		VertexPositionTextureNormalTangent() = default;

		VertexPositionTextureNormalTangent(const DirectX::XMFLOAT4& position, const DirectX::XMFLOAT2& textureCoordinates, const DirectX::XMFLOAT3& normal, const DirectX::XMFLOAT3& tangent) :
			Position(position), TextureCoordinates(textureCoordinates), Normal(normal), Tangent(tangent) { }

		DirectX::XMFLOAT4 Position;
		DirectX::XMFLOAT2 TextureCoordinates;
		DirectX::XMFLOAT3 Normal;
		DirectX::XMFLOAT3 Tangent;

		inline static const gsl::span<const D3D11_INPUT_ELEMENT_DESC> InputElements{ _InputElements };

		static void CreateVertexBuffer(gsl::not_null<ID3D11Device*> device, const Library::Mesh& mesh, gsl::not_null<ID3D11Buffer**> vertexBuffer);
		static void CreateVertexBuffer(gsl::not_null<ID3D11Device*> device, const gsl::span<const VertexPositionTextureNormalTangent>& vertices, gsl::not_null<ID3D11Buffer**> vertexBuffer)
		{
			VertexDeclaration::CreateVertexBuffer(device, vertices, vertexBuffer);
		}
	};

	class VertexSkinnedPositionTextureNormal : public VertexDeclaration<VertexSkinnedPositionTextureNormal>
	{
	private:
		inline static const D3D11_INPUT_ELEMENT_DESC _InputElements[]
		{
			{ "POSITION", 0, DXGI_FORMAT_R32G32B32A32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
			{ "TEXCOORD", 0, DXGI_FORMAT_R32G32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
			{ "NORMAL", 0, DXGI_FORMAT_R32G32B32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
			{ "BONEINDICES", 0, DXGI_FORMAT_R32G32B32A32_UINT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 },
			{ "BONEWEIGHTS", 0, DXGI_FORMAT_R32G32B32A32_FLOAT, 0, D3D11_APPEND_ALIGNED_ELEMENT, D3D11_INPUT_PER_VERTEX_DATA, 0 }
		};

	public:
		VertexSkinnedPositionTextureNormal() = default;

		VertexSkinnedPositionTextureNormal(const DirectX::XMFLOAT4& position, const DirectX::XMFLOAT2& textureCoordinates, const DirectX::XMFLOAT3& normal, const DirectX::XMUINT4& boneIndices, const DirectX::XMFLOAT4& boneWeights) :
			Position(position), TextureCoordinates(textureCoordinates), Normal(normal), BoneIndices(boneIndices), BoneWeights(boneWeights) { }

		DirectX::XMFLOAT4 Position;
		DirectX::XMFLOAT2 TextureCoordinates;
		DirectX::XMFLOAT3 Normal;
		DirectX::XMUINT4 BoneIndices;
		DirectX::XMFLOAT4 BoneWeights;

		inline static const gsl::span<const D3D11_INPUT_ELEMENT_DESC> InputElements{ _InputElements };

		static void CreateVertexBuffer(gsl::not_null<ID3D11Device*> device, const Library::Mesh& mesh, gsl::not_null<ID3D11Buffer**> vertexBuffer);
		static void CreateVertexBuffer(gsl::not_null<ID3D11Device*> device, const gsl::span<const VertexSkinnedPositionTextureNormal>& vertices, gsl::not_null<ID3D11Buffer**> vertexBuffer)
		{
			VertexDeclaration::CreateVertexBuffer(device, vertices, vertexBuffer);
		}
	};
}

#include "VertexDeclarations.inl"


text:
/* -*- linux-c -*- */
/* fewbody_int.c

   Copyright (C) 2002-2004 John M. Fregeau

   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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include "fewbody.h"

/* allocate memory for ks_params */
void fb_malloc_ks_params(fb_ks_params_t *ks_params)
{
	ks_params->m = fb_malloc_vector(ks_params->nstar);
	ks_params->M = fb_malloc_vector(ks_params->kstar);
	ks_params->amat = fb_malloc_matrix(ks_params->nstar, ks_params->kstar);
	ks_params->Tmat = fb_malloc_matrix(ks_params->kstar, ks_params->kstar);
}

/* initialize ks_params; assumes ks_params is already malloc()ed */
void fb_init_ks_params(fb_ks_params_t *ks_params, fb_hier_t hier)
{
	int i, j, k;
	double *y;

	/* exit if hier is not consistent with ks_params */
	if (ks_params->nstar != hier.nobj) {
		fprintf(stderr, "fb_init_ks_params(): ks_params->nstar != hier.nobj: ks_params->nstar=%d  hier.nobj=%d\n", \
			ks_params->nstar, hier.nobj);
		exit(1);
	}

	/* first set the mass matrix */
	for (i=0; i<hier.nobj; i++) {
		ks_params->m[i] = hier.obj[i]->m;
	}

	/* calculate the M_k */
	k = -1;
	for (i=0; i<hier.nobj-1; i++) {
		for (j=i+1; j<hier.nobj; j++) {
			k++;
			ks_params->M[k] = ks_params->m[i] * ks_params->m[j];
		}
	}

	/* calculate and set the a and T matrices */
	fb_calc_amat(ks_params->amat, ks_params->nstar, ks_params->kstar);
	fb_calc_Tmat(ks_params->amat, ks_params->m, ks_params->Tmat, ks_params->nstar, ks_params->kstar);

	/* set Einit */
	y = fb_malloc_vector(8*ks_params->kstar+1);
	fb_euclidean_to_ks(hier.obj, y, ks_params->nstar, ks_params->kstar);
	ks_params->Einit = fb_ks_Einit(y, *ks_params);

	fb_free_vector(y);
}

/* free memory for ks_params */
void fb_free_ks_params(fb_ks_params_t ks_params)
{
	fb_free_vector(ks_params.m);
	fb_free_vector(ks_params.M);
	fb_free_matrix(ks_params.amat);
	fb_free_matrix(ks_params.Tmat);
}

/* allocate memory for nonks_params */
void fb_malloc_nonks_params(fb_nonks_params_t *nonks_params)
{
	nonks_params->m = fb_malloc_vector(nonks_params->nstar);
}

/* initialize nonks_params; assumes nonks_params is already malloc()ed */
void fb_init_nonks_params(fb_nonks_params_t *nonks_params, fb_hier_t hier)
{
	int i;

	/* exit if hier is not consistent with nonks_params */
	if (nonks_params->nstar != hier.nobj) {
		fprintf(stderr, "fb_init_nonks_params(): nonks_params->nstar != hier.nobj: nonks_params->nstar=%d  hier.nobj=%d\n", \
			nonks_params->nstar, hier.nobj);
		exit(1);
	}

	/* set the mass vector */
	for (i=0; i<hier.nobj; i++) {
		nonks_params->m[i] = hier.obj[i]->m;
	}
}

/* free memory for ks_params */
void fb_free_nonks_params(fb_nonks_params_t nonks_params)
{
	fb_free_vector(nonks_params.m);
}

max_stars_repo_name: gnodvi/cosmos

max_issues_repo_name: gnodvi/cosmos

max_forks_repo_name: gnodvi/cosmos


text:
#pragma once

#include <Runtime/Runtime.h>
#include <napi/napi.h>
#include <memory>
#include <gsl/gsl>

namespace babylon
{
    class ScriptHost final
    {
    public:
        explicit ScriptHost(RuntimeImpl&);
        ScriptHost(const ScriptHost&) = delete;
        ~ScriptHost();

        void RunScript(gsl::czstring<> script, gsl::czstring<> url);

        Napi::Env& Env();

    private:
        class Impl;
        std::unique_ptr<Impl> m_impl;
    };
}

max_stars_repo_name: TrevorDev/BabylonNative

max_issues_repo_name: TrevorDev/BabylonNative

max_forks_repo_name: TrevorDev/BabylonNative


text:
#ifndef math_utility_c__H
#define math_utility_c__H
/********************************* TRICK HEADER *******************************
PURPOSE:
      (Math functions C Version)
LIBRARY DEPENDENCY:
      ((../src/math_utility_c.c)
        (../../cad/src/global_constants.c))
*******************************************************************************/

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

#ifdef __cplusplus
extern "C" {
#endif

extern const double ___PI;
extern const double __RAD;
extern const double __EPS;

int pol_from_cart_C(const gsl_vector *cart, gsl_vector *pol);
int build_psivg_thtvg_TM_C(const double psivg, const double thtvg, gsl_matrix *AMAT);
int build_psi_tht_phi_TM_C(const double psi, const double tht, const double phi, gsl_matrix *TM);
int Matrix2Quaternion_C(gsl_matrix *Matrix_in, gsl_vector *Quaternion_out);
int Quaternion2Matrix_C(const gsl_vector *Quaternion_in, gsl_matrix *Matrix_out);
int Quaternion_cross_C(const gsl_vector *Quaternion_in1, const gsl_vector *Quaternion_in2, gsl_vector *Quaternion_out);
int Quaternion2Euler_C(const gsl_vector *Quaternion_in, double *Roll, double *Pitch, double *Yaw);
int Euler2Quaternion_C(const double Roll, const double Pitch, const double Yaw, gsl_vector *Quaternion_out);
int QuaternionMultiply_C(const gsl_vector *Quaternion_in1, const gsl_vector *Quaternion_in2, gsl_vector *Quaternion_out);
int Quaternion_rotation_C(const gsl_vector *Vector_in, const gsl_vector *Quaternion_in, gsl_vector *Vector_out);
int build_321_rotation_matrix(const gsl_vector *angle, gsl_matrix *mat_out);
int sign(const double variable);
int euler_angle(gsl_matrix *TBD, gsl_vector *euler_ang);

gsl_matrix* moore_penrose_pinv(gsl_matrix *A, const double rcond);
gsl_matrix *skew_sym_C(const gsl_vector *vec);
gsl_vector *Quaternion_conjugate_C(const gsl_vector *Quaternion_in);
gsl_vector *QuaternionInverse_C(const gsl_vector *Quaternion_in);
#ifdef __cplusplus
}
#endif
#endif


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

#ifndef SIM_MATH_UTIL_H_
#define SIM_MATH_UTIL_H_

#include <stdint.h>

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

#include <functional>
#include <limits>
#include <random>
#include <string>

#include "common/c_math/util.h"

// Generates random numbers using a generator seeded by a string hash. This
// facilitates binding of generators to simulator models.
//
// A static seed offset can be specified using SetSeedOffset, supporting
// deterministic perturbation of the generators used by all models.
class NamedRandomNumberGenerator {
 public:
  explicit NamedRandomNumberGenerator(const std::string &name)
      : generator_(static_cast<uint32_t>(std::hash<std::string>{}(name)) +
                   seed_offset_),
        normal_(0.0, 1.0),
        uniform_real_(0.0, 1.0),
        uniform_int32_(0, std::numeric_limits<int32_t>::max()) {}

  // Returns a random number drawn from a normal distribution with mean=0.0
  // and stddev=1.0.
  double GetNormal() { return normal_(generator_); }

  // Returns a random number drawn from a uniform distribution on [0.0, 1.0).
  double GetUniformReal() { return uniform_real_(generator_); }

  // Returns a random integer drawn from a uniform distribution on
  // [0, std::numeric_limits<int32_t>::max()].
  int32_t GetUniformInt32() { return uniform_int32_(generator_); }

  static void SetSeedOffset(uint32_t offset) { seed_offset_ = offset; }

 private:
  static uint32_t seed_offset_;

  // 32-bit Mersenne Twister generator.
  std::mt19937 generator_;

  std::normal_distribution<double> normal_;
  std::uniform_real_distribution<double> uniform_real_;
  std::uniform_int_distribution<int32_t> uniform_int32_;
};

// The hypergeometric function, {_2}F_1(a, b; c; z)
double HyperGeom(double a, double b, double c, double z);
uint32_t Factorial(uint32_t n);

// The rising factorial, denoted by the Pochhammer symbol (q)_n, is
// defined by:
//
//   (q)_n = 1                     if n == 0
//           q*(q+1)*...*(q+n-1)   if n > 0
//
// Warning: The Pochhammer symbol can denote either rising or falling
// factorial depending on context.
double RisingFactorial(double q, int32_t n);

// Computes the intersection point of two circles in a coordinate
// system centered at the first circle and with the positive x-axis
// along the line between the centers of the first and second circles.
// For both the non-degenerate intersecting and single point
// intersecting cases, this returns the intersection point with
// non-negative y.  For the identical center case, this returns the
// point (r1, 0) where r1 is the radius of the first circle.  For the
// remaining non-intersecting case, this returns the point:
//
//   ((d^2 + r1^2 - r2^2) / (2 * d), 0)
//
// For a derivation of the formulas used here, see:
// http://mathworld.wolfram.com/Circle-CircleIntersection.html
//
// Args:
//   d: Distance, always positive, between circle centers.
//   r1: Radius of first circle.
//   r2: Radius of second circle.
//   x: Output pointer to x coordinate of intersection (i.e. the
//       distance along line between first and second circle).
//   y: Output pointer to y coordinate of intersection (i.e. the
//       perpendicular distance, always positive, from line connecting
//       circle centers).
//
// Returns:
//   True if the circles intersect.
bool IntersectCircles(double d, double r1, double r2, double *x, double *y);

// Finds the fractional index of a value x_i in a GSL vector x
// assuming x is monotonically increasing.
double gsl_interp_index(const gsl_vector *x, double x_i);

// Performs a linear interpolation along a pair of GSL vectors.  The
// function to interpolate on is described by the input GSL vector x
// and the output GSL vector y.  Interp1 linearly interpolates to find
// the value of y for input x_i.  If opt is kInterpOptionDefault,
// extrapolate linearly beyond the end of the vectors.  If opt is
// kInterpOptionSaturate, saturate the output value at the
// ends of the output (y) vector.
//
// Warning: This currently assumes a monotonically increasing x.
double gsl_interp1(const gsl_vector *x, const gsl_vector *y, double xi,
                   InterpOption opt);

// Linear interpolation along a pair of vectors and a matrix.
//
// The function to interpolate on is described by the input vectors x
// and y and the output matrix z, with corresponding entries for its
// value at each (x,y) input pair.  Interp2 bilinearly interpolates to
// find the value of z for inputs x_i and y_i.  If opt is
// kInterpOptionDefault, extrapolate linearly beyond the end of the
// vectors.  If opt is kInterpOptionSaturate, saturate the
// output value at the edges of the output (z) matrix.
//
// Warning: This currently assumes a monotonically increasing x, y.
double gsl_interp2(const gsl_vector *x, const gsl_vector *y,
                   const gsl_matrix *z, double x_i, double y_i,
                   InterpOption opt);
double gsl_interp2_scaled(double s, double t, const gsl_matrix *z,
                          InterpOption opt);
void gsl_interp2_indices(double s, double t, int32_t size1, int32_t size2,
                         InterpOption opt, int32_t *s0, int32_t *t0, double *ds,
                         double *dt);
double gsl_interp2_scaled_helper(const gsl_matrix *z, int32_t s0, int32_t t0,
                                 double ds, double dt);

// Linear interpolation along a set of three vectors and a
// three-dimensional matrix.
//
// The function to interpolate on is described by the input vectors x,
// y, and z and the output matrix mat, with corresponding entries for
// its value at each (x,y,z) input pair.  Interp3 bilinearly
// interpolates to find the value of mat for inputs x_i, y_i, and z_i.
// If opt is kInterpOptionDefault, extrapolate linearly beyond the
// end of the vectors.  If opt is kInterpOptionSaturate,
// saturate the output value at the edges of the output (mat) matrix.
//
// Warning: This currently assumes a monotonically increasing x, y, z.
double gsl_interp3(const gsl_vector *x, const gsl_vector *y,
                   const gsl_vector *z, const gsl_vector *mat, double x_i,
                   double y_i, double z_i, InterpOption opt);
double gsl_interp3_scaled(double u, double v, double w, const gsl_vector *z,
                          int32_t size1, int32_t size2, int32_t size3,
                          InterpOption opt);
void gsl_interp3_indices(double u, double v, double w, int32_t size1,
                         int32_t size2, int32_t size3, InterpOption opt,
                         int32_t *indices, double *du, double *dv, double *dw);
double gsl_interp3_scaled_helper(const gsl_vector *z, int32_t indices[],
                                 double du, double dv, double dw);

#endif  // SIM_MATH_UTIL_H_


text:
/**
 * @file localization.h
 * @brief source localization
 * @author Ulrich Klee
 */

#ifndef LOCALIZATION_H
#define LOCALIZATION_H

#include <stdio.h>
#include <assert.h>

#include <gsl/gsl_block.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_fit.h>

#include "common/jexception.h"

gsl_vector* getSrpPhat(double delta_f, gsl_matrix_complex *mFramePerChannel, gsl_vector *searchRangeX, gsl_vector *searchRangeY, gsl_matrix *arrgeom, int zPos);

void calcDelays(int x, int y, int z, const gsl_matrix* mpos, gsl_vector* delays);

void calcDelaysOfLinearMicrophoneArray(float azimuth, const gsl_matrix* mpos, gsl_vector* delays);

void calcDelaysOfCircularMicrophoneArray(float azimuth, float elevation, const gsl_matrix* mpos, gsl_vector* delays);

gsl_vector* getDelays(double delta_f, gsl_matrix_complex* mFramePerChannel, gsl_vector* searchRange);

double getAzimuth(double delta_f, gsl_matrix_complex *mFramePerChannel, gsl_matrix *arrgeom, gsl_vector *delays);

double getPlaneWaveSrp(double delta_f, gsl_matrix_complex *mFramePerChannel, gsl_vector *searchRangeY, gsl_matrix *arrgeom);

//void halfComplexPack(double* tgt, const gsl_vector_complex* src);

void getGCCRaw(const gsl_matrix_complex* spectralSample, double sampleRate, gsl_vector* gcc);

const gsl_vector* getGCC(gsl_matrix_complex *spectralSample, double sampleRate);
const gsl_vector* getWindowedGCC(gsl_matrix_complex *spectralSample, double sampleRate, double minDelay, double maxDelay);

const gsl_vector* getWindowedGCCratio(gsl_matrix_complex *spectralSample, double sampleRate, double minDelay, double maxDelay);

const gsl_vector* getWindowedGCCdirect(gsl_matrix_complex *spectralSample, double sampleRate, double minDelay, double maxDelay);

const gsl_vector* getWindowedGCCabs(gsl_matrix_complex *spectralSample, double sampleRate, double minDelay, double maxDelay);

const gsl_vector* getDynWindowedGCC(gsl_matrix_complex *spectralSample, double sampleRate, double minDelay, double maxDelay, double wMinDelay, double wMaxDelay, double threshold);

double getInterpolation(gsl_matrix *crossResult, int delayPos);

gsl_vector* get3DPosition(gsl_vector *yCoord, gsl_vector *azimuth1, gsl_vector *azimuth2, double xPos, double zPos);

gsl_vector* get3DPosition_T_shape(gsl_matrix *arrgeom1, int arrayNr1, gsl_matrix *arrgeom2, int arrayNr2,  gsl_matrix *arrgeom3, double azimuth1, double azimuth2, double azimuth3);

gsl_vector* getGCC_old(gsl_matrix_complex *spectralSample, double delta_f, gsl_vector *delays);

gsl_matrix* getLowerTriangMatrix(gsl_matrix* fullMatrix);

gsl_vector* getXi(gsl_matrix* D1_2);

class PhaseTransform {
public:
  PhaseTransform(unsigned sz);
  ~PhaseTransform();
private:
  gsl_vector_complex*					_crossSpectrum;
  double*						_crossCorrelation;
};

class NoisePowerSpectrum
{
 public:
  NoisePowerSpectrum(double alpha = 0.95);
  ~NoisePowerSpectrum() {
    if (powerSpectrum != NULL)
      gsl_vector_free(powerSpectrum);
  }
  void setAlpha(double alpha) {
    this->alpha = alpha;
    alpha1 = 1 - alpha;
  }
  double getAlpha() { return alpha; }
  void add(const gsl_vector_complex *noiseSpectrum, double timestamp);
  const gsl_vector *get() {
    return powerSpectrum;
  }
 private:
  gsl_vector *powerSpectrum;
  double alpha, alpha1;
  double timestamp;
};

class NoiseCrossSpectrum
{
 public:
  NoiseCrossSpectrum(double alpha = 0.95);
  ~NoiseCrossSpectrum() {
    if (crossSpectrum != NULL)
      gsl_vector_complex_free(crossSpectrum);
  }
  void setAlpha(double alpha) {
    this->alpha = alpha;
    alpha1 = 1 - alpha;
  }
  double getAlpha() { return alpha; }
  const gsl_vector_complex *get() {
    return crossSpectrum;
  }
  void add(const gsl_vector_complex *noiseSpectrum1, const gsl_vector_complex *noiseSpectrum2);
 private:
  gsl_vector_complex *crossSpectrum;
  double alpha, alpha1;
};

class GCC
{
 public:
  GCC(double sampleRate = 44100.0 , unsigned fftLen = 2048, unsigned nChan = 16, unsigned pairs = 6, double alpha = 0.95, double beta = 0.5, double q = 0.3, bool interpolate = true, bool noisereduction = true);
  ~GCC();
  void calculate(const gsl_vector_complex *spectralSample1, unsigned chan1, const gsl_vector_complex *spectralSample2, unsigned chan2, unsigned pair, double timestamp, bool sad = false, bool smooth = true);
  const gsl_vector* findMaximum(double minDelay = -HUGE_VAL, double maxDelay = HUGE_VAL);
  double getPeakDelay() { return delay; }
  double getPeakCorr() { return maxCorr; }
  double getRatio() { return ratio; }
  const gsl_vector* getNoisePowerSpectrum(unsigned chan) { return noisePowerSpectrum[chan].get(); }
  const gsl_vector_complex* getNoiseCrossSpectrum(unsigned pair) { return noiseCrossSpectrum[pair].get(); }
  const gsl_vector_complex* getCrossSpectrum() { return crossSpectrum; }
  const gsl_vector* getCrossCorrelation() { return crossCorrelation; }
  void setAlpha(double alpha) {
    for(unsigned i=0; i < pairs; i++)
      noiseCrossSpectrum[i].setAlpha(alpha);
    for(unsigned i=0; i < nChan; i++)
      noisePowerSpectrum[i].setAlpha(alpha);
  }
  double getAlpha() { return noiseCrossSpectrum[0].getAlpha(); }
 protected:
  virtual gsl_complex calcCrossSpectrumValue(gsl_complex x1, gsl_complex x2, const gsl_vector_complex* Gn1n2, const gsl_vector* N1, const gsl_vector* N2, unsigned i) { throw jdimension_error("Not implemented!!!\n"); };

  //auxiliary variables
  gsl_complex x1, x2, Gx1x2, Gs1s2, G;
  double W1, W2, X1, X12, X2, X22;

  double sampleRate;
  unsigned fftLen;
  unsigned fftLen2;
  unsigned len;
  unsigned nChan;
  unsigned pairs;
  unsigned delayPos;
  double q, q1, q2, beta, beta1;
  double ratio;
  double delay;
  double maxCorr;
  double maxCorr2;
  gsl_vector *crossCorrelation;
  gsl_matrix *interpolValues;
  gsl_vector* retValues;
  gsl_vector *powerSpectrum1, *powerSpectrum2;
  gsl_vector_complex *cSpectrum;
  gsl_vector_complex *crossSpectrum;
  NoisePowerSpectrum* noisePowerSpectrum;
  NoiseCrossSpectrum* noiseCrossSpectrum;
  bool noisereduction;
  bool sad;
  bool interpolate;
};

class GCCRaw : public GCC
{
 public:
  GCCRaw(double sampleRate = 44100.0 , unsigned fftLen = 2048, unsigned nChan = 16, unsigned pairs = 6, double alpha = 0.95, double beta = 0.5, double q = 0.3, bool interpolate = true, bool noisereduction = true) : GCC(sampleRate, fftLen, nChan, pairs, alpha, beta, q, interpolate, noisereduction) {};
 protected:
  gsl_complex calcCrossSpectrumValue(gsl_complex x1, gsl_complex x2, const gsl_vector_complex* Gn1n2, const gsl_vector* N1, const gsl_vector* N2, unsigned i);
};

class GCCGnnSub : public GCC
{
 public:
  GCCGnnSub(double sampleRate = 44100.0 , unsigned fftLen = 2048, unsigned nChan = 16, unsigned pairs = 6, double alpha = 0.95, double beta = 0.5, double q = 0.3, bool interpolate = true, bool noisereduction = true) : GCC(sampleRate, fftLen, nChan, pairs, alpha, beta, q, interpolate, noisereduction) {};
 protected:
  gsl_complex calcCrossSpectrumValue(gsl_complex x1, gsl_complex x2, const gsl_vector_complex* Gn1n2, const gsl_vector* N1, const gsl_vector* N2, unsigned i);
};

class GCCPhat : public GCC
{
 public:
  GCCPhat(double sampleRate = 44100.0 , unsigned fftLen = 2048, unsigned nChan = 16, unsigned pairs = 6, double alpha = 0.95, double beta = 0.5, double q = 0.3, bool interpolate = true, bool noisereduction = true) : GCC(sampleRate, fftLen, nChan, pairs, alpha, beta, q, interpolate, noisereduction) {};
 protected:
  gsl_complex calcCrossSpectrumValue(gsl_complex x1, gsl_complex x2, const gsl_vector_complex* Gn1n2, const gsl_vector* N1, const gsl_vector* N2, unsigned i);
};

class GCCGnnSubPhat : public GCC
{
 public:
  GCCGnnSubPhat(double sampleRate = 44100.0 , unsigned fftLen = 2048, unsigned nChan = 16, unsigned pairs = 6, double alpha = 0.95, double beta = 0.5, double q = 0.3, bool interpolate = true, bool noisereduction = true) : GCC(sampleRate, fftLen, nChan, pairs, alpha, beta, q, interpolate, noisereduction) {};
 protected:
  gsl_complex calcCrossSpectrumValue(gsl_complex x1, gsl_complex x2, const gsl_vector_complex* Gn1n2, const gsl_vector* N1, const gsl_vector* N2, unsigned i);
};

class GCCMLRRaw : public GCC
{
 public:
  GCCMLRRaw(double sampleRate = 44100.0 , unsigned fftLen = 2048, unsigned nChan = 16, unsigned pairs = 6, double alpha = 0.95, double beta = 0.5, double q = 0.3, bool interpolate = true, bool noisereduction = true) : GCC(sampleRate, fftLen, nChan, pairs, alpha, beta, q, interpolate, noisereduction) {};
 protected:
  gsl_complex calcCrossSpectrumValue(gsl_complex x1, gsl_complex x2, const gsl_vector_complex* Gn1n2, const gsl_vector* N1, const gsl_vector* N2, unsigned i);
};

class GCCMLRGnnSub : public GCC
{
 public:
  GCCMLRGnnSub(double sampleRate = 44100.0 , unsigned fftLen = 2048, unsigned nChan = 16, unsigned pairs = 6, double alpha = 0.95, double beta = 0.5, double q = 0.3, bool interpolate = true, bool noisereduction = true) : GCC(sampleRate, fftLen, nChan, pairs, alpha, beta, q, interpolate, noisereduction) {};
 protected:
  gsl_complex calcCrossSpectrumValue(gsl_complex x1, gsl_complex x2, const gsl_vector_complex* Gn1n2, const gsl_vector* N1, const gsl_vector* N2, unsigned i);
};

#endif


text:
/**
 *
 * @file core_stsmqr_sytra1.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Jakub Kurzak
 * @author Azzam Haidar
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:44:48 2014
 *
 **/
#include <lapacke.h>
#include "common.h"
#undef COMPLEX
#define REAL

/***************************************************************************//**
 *
 * @ingroup CORE_float
 *
 *  CORE_stsmqr_sytra1: see CORE_stsmqr
 *
 *  This kernel applies a left transformation on | A1'|
 *                                               | A2 |
 *
 * Needs therefore to make the explicit transpose of A1 before
 * and after the application of the block of reflectors
 * Can be further optimized by changing accordingly the underneath
 * kernel ztsrfb!
 *
 *******************************************************************************
 *
 * @param[in] side
 *         @arg PlasmaLeft  : apply Q or Q**T from the Left;
 *         @arg PlasmaRight : apply Q or Q**T from the Right.
 *
 * @param[in] trans
 *         @arg PlasmaNoTrans   :  No transpose, apply Q;
 *         @arg PlasmaTrans :  ConjTranspose, apply Q**T.
 *
 * @param[in] m1
 *         The number of rows of the tile A1. M1 >= 0.
 *
 * @param[in] n1
 *         The number of columns of the tile A1. N1 >= 0.
 *
 * @param[in] m2
 *         The number of rows of the tile A2. M2 >= 0.
 *         M2 = M1 if side == PlasmaRight.
 *
 * @param[in] n2
 *         The number of columns of the tile A2. N2 >= 0.
 *         N2 = N1 if side == PlasmaLeft.
 *
 * @param[in] k
 *         The number of elementary reflectors whose product defines
 *         the matrix Q.
 *
 * @param[in] ib
 *         The inner-blocking size.  IB >= 0.
 *
 * @param[in,out] A1
 *         On entry, the M1-by-N1 tile A1.
 *         On exit, A1 is overwritten by the application of Q.
 *
 * @param[in] lda1
 *         The leading dimension of the array A1. LDA1 >= max(1,M1).
 *
 * @param[in,out] A2
 *         On entry, the M2-by-N2 tile A2.
 *         On exit, A2 is overwritten by the application of Q.
 *
 * @param[in] lda2
 *         The leading dimension of the tile A2. LDA2 >= max(1,M2).
 *
 * @param[in] V
 *         The i-th row must contain the vector which defines the
 *         elementary reflector H(i), for i = 1,2,...,k, as returned by
 *         CORE_STSQRT in the first k columns of its array argument V.
 *
 * @param[in] ldv
 *         The leading dimension of the array V. LDV >= max(1,K).
 *
 * @param[in] T
 *         The IB-by-N1 triangular factor T of the block reflector.
 *         T is upper triangular by block (economic storage);
 *         The rest of the array is not referenced.
 *
 * @param[in] ldt
 *         The leading dimension of the array T. LDT >= IB.
 *
 * @param[out] WORK
 *         Workspace array of size
 *             LDWORK-by-N1 if side == PlasmaLeft
 *             LDWORK-by-IB if side == PlasmaRight
 *
 * @param[in] ldwork
 *         The leading dimension of the array WORK.
 *             LDWORK >= max(1,IB) if side == PlasmaLeft
 *             LDWORK >= max(1,M1) if side == PlasmaRight
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_stsmqr_sytra1 = PCORE_stsmqr_sytra1
#define CORE_stsmqr_sytra1 PCORE_stsmqr_sytra1
#define CORE_stsmqr PCORE_stsmqr
int  CORE_stsmqr(PLASMA_enum side, PLASMA_enum trans,
                 int M1, int N1, int M2, int N2, int K, int IB,
                 float *A1, int LDA1,
                 float *A2, int LDA2,
                 const float *V, int LDV,
                 const float *T, int LDT,
                 float *WORK, int LDWORK);
#endif
int CORE_stsmqr_sytra1( PLASMA_enum side, PLASMA_enum trans,
                        int m1, int n1, int m2, int n2,
                        int k, int ib,
                        float *A1, int lda1,
                        float *A2, int lda2,
                        const float *V, int ldv,
                        const float *T, int ldt,
                        float *WORK, int ldwork)
{
    int i, j;

    if ( (m1 != n1) ) {
        coreblas_error(3, "Illegal value of M1, N1");
        return -3;
    }

    /* in-place transposition of A1 */
    for (j = 0; j < n1; j++){
        A1[j + j*lda1] = (A1[j + j*lda1]);

        for (i = j+1; i < m1; i++){
            *WORK = *(A1 + i + j*lda1);
            *(A1 + i + j*lda1) = (*(A1 + j + i*lda1));
            *(A1 + j + i*lda1) = (*WORK);
        }
    }

    CORE_stsmqr(side, trans, m1, n1, m2, n2, k, ib, A1, lda1, A2, lda2, V, ldv, T, ldt, WORK, ldwork);

    /* in-place transposition of A1 */
    for (j = 0; j < n1; j++){
        A1[j + j*lda1] = (A1[j + j*lda1]);

        for (i = j+1; i < m1; i++){
            *WORK = *(A1 + i + j*lda1);
            *(A1 + i + j*lda1) = (*(A1 + j + i*lda1));
            *(A1 + j + i*lda1) = (*WORK);
        }
    }

    return PLASMA_SUCCESS;
}


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

#include "tests.h"

void
test_gbmv (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int order = 101;
   int trans = 111;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   float alpha = -1.0f;
   float beta = -1.0f;
   float A[] = { 0.423f, -0.143f, -0.182f, -0.076f, -0.855f, 0.599f, 0.389f, -0.473f, 0.493f, -0.902f, -0.889f, -0.256f, 0.112f, 0.128f, -0.277f, -0.777f };
   float X[] = { 0.488f, 0.029f, -0.633f, 0.84f };
   int incX = -1;
   float Y[] = { 0.874f, 0.322f, -0.477f };
   int incY = -1;
   float y_expected[] = { -0.101941f, 0.764086f, 0.481914f };
   cblas_sgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "sgbmv(case 794)");
     }
   };
  };

  {
   int order = 102;
   int trans = 111;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   float alpha = -1.0f;
   float beta = -1.0f;
   float A[] = { 0.423f, -0.143f, -0.182f, -0.076f, -0.855f, 0.599f, 0.389f, -0.473f, 0.493f, -0.902f, -0.889f, -0.256f, 0.112f, 0.128f, -0.277f, -0.777f };
   float X[] = { 0.488f, 0.029f, -0.633f, 0.84f };
   int incX = -1;
   float Y[] = { 0.874f, 0.322f, -0.477f };
   int incY = -1;
   float y_expected[] = { -0.656261f, 0.19575f, 0.055905f };
   cblas_sgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "sgbmv(case 795)");
     }
   };
  };

  {
   int order = 101;
   int trans = 112;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   float alpha = 0.0f;
   float beta = 0.1f;
   float A[] = { -0.066f, -0.153f, -0.619f, 0.174f, 0.777f, 0.543f, 0.614f, -0.446f, -0.138f, -0.767f, 0.725f, 0.222f, 0.165f, -0.063f, -0.047f, 0.267f };
   float X[] = { -0.096f, -0.007f, -0.657f };
   int incX = -1;
   float Y[] = { -0.88f, 0.102f, -0.278f, 0.403f };
   int incY = -1;
   float y_expected[] = { -0.088f, 0.0102f, -0.0278f, 0.0403f };
   cblas_sgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "sgbmv(case 796)");
     }
   };
  };

  {
   int order = 102;
   int trans = 112;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   float alpha = 0.0f;
   float beta = 0.1f;
   float A[] = { -0.066f, -0.153f, -0.619f, 0.174f, 0.777f, 0.543f, 0.614f, -0.446f, -0.138f, -0.767f, 0.725f, 0.222f, 0.165f, -0.063f, -0.047f, 0.267f };
   float X[] = { -0.096f, -0.007f, -0.657f };
   int incX = -1;
   float Y[] = { -0.88f, 0.102f, -0.278f, 0.403f };
   int incY = -1;
   float y_expected[] = { -0.088f, 0.0102f, -0.0278f, 0.0403f };
   cblas_sgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "sgbmv(case 797)");
     }
   };
  };

  {
   int order = 101;
   int trans = 111;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   double alpha = 0.1;
   double beta = 0;
   double A[] = { -0.688, 0.29, 0.442, -0.001, 0.313, -0.073, 0.991, -0.654, -0.12, 0.416, 0.571, 0.932, -0.179, -0.724, 0.492, -0.965 };
   double X[] = { 0.187, -0.338, -0.976, -0.052 };
   int incX = -1;
   double Y[] = { -0.101, 0.8, 0.026 };
   int incY = -1;
   double y_expected[] = { 0.0083289, -0.0279986, -0.0446472 };
   cblas_dgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dgbmv(case 798)");
     }
   };
  };

  {
   int order = 102;
   int trans = 111;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   double alpha = 0.1;
   double beta = 0;
   double A[] = { -0.688, 0.29, 0.442, -0.001, 0.313, -0.073, 0.991, -0.654, -0.12, 0.416, 0.571, 0.932, -0.179, -0.724, 0.492, -0.965 };
   double X[] = { 0.187, -0.338, -0.976, -0.052 };
   int incX = -1;
   double Y[] = { -0.101, 0.8, 0.026 };
   int incY = -1;
   double y_expected[] = { -0.1141297, 0.0088824, -0.0320568 };
   cblas_dgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dgbmv(case 799)");
     }
   };
  };

  {
   int order = 101;
   int trans = 112;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   double alpha = -0.3;
   double beta = -0.3;
   double A[] = { 0.746, 0.262, -0.449, -0.954, -0.093, 0.108, -0.496, 0.927, 0.177, 0.729, -0.92, -0.469, 0.87, -0.877, -0.308, -0.806 };
   double X[] = { 0.662, -0.887, 0.261 };
   int incX = -1;
   double Y[] = { 0.771, 0.637, -0.177, -0.018 };
   int incY = -1;
   double y_expected[] = { -0.048588, -0.467865, 0.0818433, -0.0398619 };
   cblas_dgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dgbmv(case 800)");
     }
   };
  };

  {
   int order = 102;
   int trans = 112;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   double alpha = -0.3;
   double beta = -0.3;
   double A[] = { 0.746, 0.262, -0.449, -0.954, -0.093, 0.108, -0.496, 0.927, 0.177, 0.729, -0.92, -0.469, 0.87, -0.877, -0.308, -0.806 };
   double X[] = { 0.662, -0.887, 0.261 };
   int incX = -1;
   double Y[] = { 0.771, 0.637, -0.177, -0.018 };
   int incY = -1;
   double y_expected[] = { -0.404082, -0.2887797, 0.1876263, -0.1345935 };
   cblas_dgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dgbmv(case 801)");
     }
   };
  };

  {
   int order = 101;
   int trans = 111;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   float alpha[2] = {0.0f, 1.0f};
   float beta[2] = {0.0f, 0.1f};
   float A[] = { -0.107f, 0.926f, -0.246f, -0.555f, -0.301f, 0.276f, 0.471f, -0.084f, -0.754f, 0.082f, -0.952f, -0.394f, 0.659f, 0.054f, 0.795f, 0.923f, 0.232f, -0.788f, 0.478f, 0.775f, -0.118f, 0.691f, -0.933f, 0.809f, 0.164f, -0.263f, -0.923f, -0.88f, 0.819f, -0.521f, -0.045f, 0.034f };
   float X[] = { 0.407f, 0.895f, 0.301f, 0.769f, -0.269f, -0.465f, 0.455f, -0.628f };
   int incX = -1;
   float Y[] = { -0.116f, -0.744f, -0.936f, -0.064f, -0.232f, -0.665f };
   int incY = -1;
   float y_expected[] = { -0.806176f, -1.559f, -1.57611f, -0.155463f, 0.098816f, -0.274361f };
   cblas_cgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], flteps, "cgbmv(case 802) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], flteps, "cgbmv(case 802) imag");
     };
   };
  };

  {
   int order = 102;
   int trans = 111;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   float alpha[2] = {0.0f, 1.0f};
   float beta[2] = {0.0f, 0.1f};
   float A[] = { -0.107f, 0.926f, -0.246f, -0.555f, -0.301f, 0.276f, 0.471f, -0.084f, -0.754f, 0.082f, -0.952f, -0.394f, 0.659f, 0.054f, 0.795f, 0.923f, 0.232f, -0.788f, 0.478f, 0.775f, -0.118f, 0.691f, -0.933f, 0.809f, 0.164f, -0.263f, -0.923f, -0.88f, 0.819f, -0.521f, -0.045f, 0.034f };
   float X[] = { 0.407f, 0.895f, 0.301f, 0.769f, -0.269f, -0.465f, 0.455f, -0.628f };
   int incX = -1;
   float Y[] = { -0.116f, -0.744f, -0.936f, -0.064f, -0.232f, -0.665f };
   int incY = -1;
   float y_expected[] = { -0.245235f, -0.313725f, -0.798094f, 0.691455f, -0.164015f, -0.242714f };
   cblas_cgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], flteps, "cgbmv(case 803) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], flteps, "cgbmv(case 803) imag");
     };
   };
  };

  {
   int order = 101;
   int trans = 112;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   float alpha[2] = {-1.0f, 0.0f};
   float beta[2] = {0.0f, 0.1f};
   float A[] = { -0.258f, 0.838f, -0.106f, -0.066f, 0.395f, 0.982f, -0.546f, 0.565f, 0.14f, -0.18f, 0.165f, -0.186f, 0.499f, -0.038f, -0.305f, -0.653f, -0.811f, -0.466f, -0.674f, -0.013f, -0.552f, -0.807f, -0.536f, 0.864f, -0.027f, -0.606f, 0.459f, 0.564f, -0.968f, 0.717f, -0.312f, -0.485f };
   float X[] = { -0.399f, 0.459f, 0.398f, 0.358f, -0.161f, -0.359f };
   int incX = -1;
   float Y[] = { 0.572f, 0.293f, -0.813f, -0.096f, -0.611f, -0.717f, 0.736f, 0.259f };
   int incY = -1;
   float y_expected[] = { -0.619961f, -0.011425f, -0.477499f, 0.059361f, -0.886984f, 0.44008f, -0.139432f, 0.04644f };
   cblas_cgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], flteps, "cgbmv(case 804) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], flteps, "cgbmv(case 804) imag");
     };
   };
  };

  {
   int order = 102;
   int trans = 112;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   float alpha[2] = {-1.0f, 0.0f};
   float beta[2] = {0.0f, 0.1f};
   float A[] = { -0.258f, 0.838f, -0.106f, -0.066f, 0.395f, 0.982f, -0.546f, 0.565f, 0.14f, -0.18f, 0.165f, -0.186f, 0.499f, -0.038f, -0.305f, -0.653f, -0.811f, -0.466f, -0.674f, -0.013f, -0.552f, -0.807f, -0.536f, 0.864f, -0.027f, -0.606f, 0.459f, 0.564f, -0.968f, 0.717f, -0.312f, -0.485f };
   float X[] = { -0.399f, 0.459f, 0.398f, 0.358f, -0.161f, -0.359f };
   int incX = -1;
   float Y[] = { 0.572f, 0.293f, -0.813f, -0.096f, -0.611f, -0.717f, 0.736f, 0.259f };
   int incY = -1;
   float y_expected[] = { -0.318227f, -0.172201f, -0.109343f, 0.698685f, 0.208261f, -0.269065f, 0.175074f, -0.507326f };
   cblas_cgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], flteps, "cgbmv(case 805) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], flteps, "cgbmv(case 805) imag");
     };
   };
  };

  {
   int order = 101;
   int trans = 113;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   float alpha[2] = {-1.0f, 0.0f};
   float beta[2] = {0.0f, 0.1f};
   float A[] = { -0.804f, 0.232f, -0.448f, -0.558f, -0.078f, -0.056f, -0.345f, -0.379f, 0.369f, -0.662f, -0.169f, -0.391f, -0.215f, 0.467f, 0.374f, 0.889f, -0.698f, 0.734f, 0.377f, -0.955f, 0.498f, 0.151f, -0.725f, -0.728f, -0.655f, -0.581f, 0.389f, 0.949f, -0.553f, -0.434f, 0.237f, 0.641f };
   float X[] = { -0.262f, -0.823f, -0.357f, -0.994f, -0.347f, -0.375f };
   int incX = -1;
   float Y[] = { -0.683f, -0.87f, -0.708f, 0.071f, 0.575f, -0.575f, 0.845f, 0.032f };
   int incY = -1;
   float y_expected[] = { 0.341749f, 0.301992f, -0.306848f, 0.109252f, -0.018347f, -0.747479f, -0.894201f, 0.713246f };
   cblas_cgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], flteps, "cgbmv(case 806) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], flteps, "cgbmv(case 806) imag");
     };
   };
  };

  {
   int order = 102;
   int trans = 113;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   float alpha[2] = {-1.0f, 0.0f};
   float beta[2] = {0.0f, 0.1f};
   float A[] = { -0.804f, 0.232f, -0.448f, -0.558f, -0.078f, -0.056f, -0.345f, -0.379f, 0.369f, -0.662f, -0.169f, -0.391f, -0.215f, 0.467f, 0.374f, 0.889f, -0.698f, 0.734f, 0.377f, -0.955f, 0.498f, 0.151f, -0.725f, -0.728f, -0.655f, -0.581f, 0.389f, 0.949f, -0.553f, -0.434f, 0.237f, 0.641f };
   float X[] = { -0.262f, -0.823f, -0.357f, -0.994f, -0.347f, -0.375f };
   int incX = -1;
   float Y[] = { -0.683f, -0.87f, -0.708f, 0.071f, 0.575f, -0.575f, 0.845f, 0.032f };
   int incY = -1;
   float y_expected[] = { -0.562773f, -0.455143f, -0.213881f, -0.466169f, -0.183683f, 0.097891f, -0.451416f, 0.052586f };
   cblas_cgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], flteps, "cgbmv(case 807) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], flteps, "cgbmv(case 807) imag");
     };
   };
  };

  {
   int order = 101;
   int trans = 111;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   double alpha[2] = {0, 0.1};
   double beta[2] = {1, 0};
   double A[] = { -0.919, -0.002, 0.105, -0.338, -0.358, -0.715, -0.157, 0.307, 0.334, 0.121, 0.366, 0.029, -0.006, -0.662, -0.314, 0.061, -0.322, -0.865, -0.586, 0.556, 0.507, 0.581, 0.855, -0.09, 0.836, -0.788, -0.209, -0.694, -0.695, 0.11, -0.234, 0.17 };
   double X[] = { 0.356, -0.76, -0.96, 0.437, -0.849, 0.397, -0.382, -0.826 };
   int incX = -1;
   double Y[] = { 0.288, -0.832, 0.889, 0.576, -0.809, 0.4 };
   int incY = -1;
   double y_expected[] = { 0.3241775, -0.6761577, 0.8458527, 0.5705165, -0.8597295, 0.4268499 };
   cblas_zgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], dbleps, "zgbmv(case 808) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], dbleps, "zgbmv(case 808) imag");
     };
   };
  };

  {
   int order = 102;
   int trans = 111;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   double alpha[2] = {0, 0.1};
   double beta[2] = {1, 0};
   double A[] = { -0.919, -0.002, 0.105, -0.338, -0.358, -0.715, -0.157, 0.307, 0.334, 0.121, 0.366, 0.029, -0.006, -0.662, -0.314, 0.061, -0.322, -0.865, -0.586, 0.556, 0.507, 0.581, 0.855, -0.09, 0.836, -0.788, -0.209, -0.694, -0.695, 0.11, -0.234, 0.17 };
   double X[] = { 0.356, -0.76, -0.96, 0.437, -0.849, 0.397, -0.382, -0.826 };
   int incX = -1;
   double Y[] = { 0.288, -0.832, 0.889, 0.576, -0.809, 0.4 };
   int incY = -1;
   double y_expected[] = { 0.4026074, -0.8033768, 0.7510795, 0.5671044, -0.8162255, 0.3349099 };
   cblas_zgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 3; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], dbleps, "zgbmv(case 809) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], dbleps, "zgbmv(case 809) imag");
     };
   };
  };

  {
   int order = 101;
   int trans = 112;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   double alpha[2] = {1, 0};
   double beta[2] = {1, 0};
   double A[] = { 0.511, -0.707, -0.906, 0.345, -0.524, -0.933, 0.154, -0.529, -0.651, -0.851, 0.104, 0.532, -0.297, 0.477, 0.511, 0.469, -0.888, -0.789, 0.656, 0.288, -0.749, 0.961, 0.571, 0.539, 0.465, 0.647, 0.653, -0.994, -0.515, 0.297, 0.35, -0.707 };
   double X[] = { -0.991, 0.658, -0.909, -0.99, -0.517, -0.071 };
   int incX = -1;
   double Y[] = { 0.451, 0.351, -0.113, -0.62, 0.983, 0.511, 0.142, -0.186 };
   int incY = -1;
   double y_expected[] = { 0.560921, -1.094193, -0.210397, -0.613323, 3.018979, 0.641612, 0.384166, 1.11801 };
   cblas_zgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], dbleps, "zgbmv(case 810) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], dbleps, "zgbmv(case 810) imag");
     };
   };
  };

  {
   int order = 102;
   int trans = 112;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   double alpha[2] = {1, 0};
   double beta[2] = {1, 0};
   double A[] = { 0.511, -0.707, -0.906, 0.345, -0.524, -0.933, 0.154, -0.529, -0.651, -0.851, 0.104, 0.532, -0.297, 0.477, 0.511, 0.469, -0.888, -0.789, 0.656, 0.288, -0.749, 0.961, 0.571, 0.539, 0.465, 0.647, 0.653, -0.994, -0.515, 0.297, 0.35, -0.707 };
   double X[] = { -0.991, 0.658, -0.909, -0.99, -0.517, -0.071 };
   int incX = -1;
   double Y[] = { 0.451, 0.351, -0.113, -0.62, 0.983, 0.511, 0.142, -0.186 };
   int incY = -1;
   double y_expected[] = { -0.435541, 0.015793, -0.926518, 1.122561, 1.671751, -0.257493, 0.187543, 1.066818 };
   cblas_zgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], dbleps, "zgbmv(case 811) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], dbleps, "zgbmv(case 811) imag");
     };
   };
  };

  {
   int order = 101;
   int trans = 113;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   double alpha[2] = {0, 0.1};
   double beta[2] = {-0.3, 0.1};
   double A[] = { 0.534, 0.67, -0.621, 0.143, -0.794, 0.073, 0.414, -0.9, 0.155, -0.368, 0.122, -0.583, 0.03, 0.646, -0.768, -0.892, -0.741, -0.397, 0.626, 0.004, -0.515, 0.355, 0.196, -0.989, -0.982, 0.985, 0.445, 0.63, -0.849, -0.528, 0.146, -0.319 };
   double X[] = { -0.199, -0.259, 0.386, -0.131, -0.867, 0.888 };
   int incX = -1;
   double Y[] = { 0.106, 0.874, 0.962, 0.636, -0.759, 0.415, -0.053, 0.315 };
   int incY = -1;
   double y_expected[] = { -0.139603, -0.250546, -0.3107376, -0.1144656, 0.2181809, -0.0877031, 0.0149724, -0.0224571 };
   cblas_zgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], dbleps, "zgbmv(case 812) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], dbleps, "zgbmv(case 812) imag");
     };
   };
  };

  {
   int order = 102;
   int trans = 113;
   int M = 3;
   int N = 4;
   int KL = 1;
   int KU = 1;
   int lda = 4;
   double alpha[2] = {0, 0.1};
   double beta[2] = {-0.3, 0.1};
   double A[] = { 0.534, 0.67, -0.621, 0.143, -0.794, 0.073, 0.414, -0.9, 0.155, -0.368, 0.122, -0.583, 0.03, 0.646, -0.768, -0.892, -0.741, -0.397, 0.626, 0.004, -0.515, 0.355, 0.196, -0.989, -0.982, 0.985, 0.445, 0.63, -0.849, -0.528, 0.146, -0.319 };
   double X[] = { -0.199, -0.259, 0.386, -0.131, -0.867, 0.888 };
   int incX = -1;
   double Y[] = { 0.106, 0.874, 0.962, 0.636, -0.759, 0.415, -0.053, 0.315 };
   int incY = -1;
   double y_expected[] = { -0.1642353, -0.2575697, -0.3610975, -0.1305629, 0.1713576, -0.2514988, 0.0195631, -0.0648656 };
   cblas_zgbmv(order, trans, M, N, KU, KL, alpha, A,                                  lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(Y[2*i], y_expected[2*i], dbleps, "zgbmv(case 813) real");
       gsl_test_rel(Y[2*i+1], y_expected[2*i+1], dbleps, "zgbmv(case 813) imag");
     };
   };
  };

}


text:
#include <math.h>
#include <gsl/gsl_ntuple.h>
#include <gsl/gsl_histogram.h>

struct data
{
  double x;
  double y;
  double z;
};

int sel_func (void *ntuple_data, void *params);
double val_func (void *ntuple_data, void *params);

int
main (void)
{
  struct data ntuple_row;

  gsl_ntuple *ntuple
    = gsl_ntuple_open ("test.dat", &ntuple_row,
                       sizeof (ntuple_row));
  double lower = 1.5;

  gsl_ntuple_select_fn S;
  gsl_ntuple_value_fn V;

  gsl_histogram *h = gsl_histogram_alloc (100);
  gsl_histogram_set_ranges_uniform(h, 0.0, 10.0);

  S.function = &sel_func;
  S.params = &lower;

  V.function = &val_func;
  V.params = 0;

  gsl_ntuple_project (h, ntuple, &V, &S);
  gsl_histogram_fprintf (stdout, h, "%f", "%f");
  gsl_histogram_free (h);
  gsl_ntuple_close (ntuple);

  return 0;
}

int
sel_func (void *ntuple_data, void *params)
{
  struct data * data = (struct data *) ntuple_data;
  double x, y, z, E2, scale;
  scale = *(double *) params;

  x = data->x;
  y = data->y;
  z = data->z;

  E2 = x * x + y * y + z * z;

  return E2 > scale;
}

double
val_func (void *ntuple_data, void *params)
{
  struct data * data = (struct data *) ntuple_data;
  double x, y, z;

  x = data->x;
  y = data->y;
  z = data->z;

  return x * x + y * y + z * z;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/* cheb/deriv.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 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_chebyshev.h>

int gsl_cheb_calc_deriv(gsl_cheb_series * deriv, const gsl_cheb_series * f)
{
  const size_t n = f->order + 1;
  const double con = 2.0 / (f->b - f->a);
  size_t i;

  if(deriv->order != f->order)
    {
      GSL_ERROR ("order of chebyshev series must be equal", GSL_ENOMEM);
    }

  /* set the other parameters in the chebyshev struct */

  deriv->a = f->a;
  deriv->b = f->b;

  /* error in derivative is n^2 c_n */
  /* deriv->err = n * n * f->c[n-1];*/

  /* FIXME:  should probably set deriv->f[] as well */

  deriv->c[n-1] = 0.0;

  if(n > 1) {
    deriv->c[n-2] = 2.0 *(n-1.0) * f->c[n-1];

    for(i = n; i>=3; i--)
      deriv->c[i-3] = deriv->c[i-1] + 2.0 *(i-2.0) * f->c[i-2];

    for(i = 0  ; i<n ; i++)
      deriv->c[i] *= con;
  }

  return GSL_SUCCESS;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/**
 * @file cblas.h
 * @brief Header file to handle Intel MKL, OpenBLAS, or Netlib includes.
 *
 * `npypacke` can be linked with a reference CBLAS implementation, OpenBLAS, or
 * Intel MKL, but Intel MKL uses a different header file, `mkl.h`, while
 * OpenBLAS uses the CBLAS header `cblas.h`. This file includes the correct
 * header file given an appropriate preprocessor macro is defined, and defines
 * some types appropriately so that the user can simply write in terms of the
 * Intel MKL interface, whether or not Intel MKL is actually linked.
 */

#ifndef NPY_LPK_CBLAS_H
#define NPY_LPK_CBLAS_H
// if linking with Netlib CBLAS
#if defined(CBLAS_INCLUDE)
#include <cblas.h>
// define MKL_INT used in extension modules as in mkl.h
#ifndef MKL_INT
#define MKL_INT int
#endif /* MKL_INT */
// else if linking with OpenBLAS CBLAS
#elif defined(OPENBLAS_INCLUDE)
#include <cblas.h>
// OpenBLAS has blasint typedef, so define MKL_INT using blasint
#ifndef MKL_INT
#define MKL_INT blasint
#endif /* MKL_INT */
// else if linking to Intel MKL
#elif defined(MKL_INCLUDE)
#include <mkl.h>
// else error, no LAPACKE includes specified
#else
// silence error squiggles in VS Code (__INTELLISENSE__ always defined)
#ifndef __INTELLISENSE__
#error "no CBLAS includes specified. try -D(CBLAS|OPENBLAS|MKL)_INCLUDE"
#endif /* __INTELLISENSE__ */
#endif

// silence error squiggles in VS Code. use mkl.h since it also define MKL types
#ifdef __INTELLISENSE__
#include <mkl.h>
#endif /* __INTELLISENSE__ */

#endif /* NPY_LPK_CBLAS_H */

max_issues_repo_name: phetdam/numpy-lapacke-demo

max_forks_repo_name: phetdam/numpy-lapacke-demo