text:
#pragma once
#include <typeinfo>
#include <gsl/gsl_assert>
#include <algorithm>
#include "halleystring.h"
#include <halley/support/exception.h>

namespace Halley
{

	template <typename T>
	struct EnumNames {
	};

	struct UserConverter
	{
		template<typename T, typename std::enable_if<std::is_enum<T>::value, int>::type = 0>
		static String toString(const T& v)
		{
			return EnumNames<T>()()[int(v)];
		}

		template<typename T, typename std::enable_if<!std::is_enum<T>::value, int>::type = 0>
		static String toString(const T& v)
		{
			return v.toString();
		}

		template<typename T, typename std::enable_if<std::is_enum<T>::value, int>::type = 0>
		static T fromString(const String& str)
		{
			EnumNames<T> n;
			auto names = n();
			auto res = std::find_if(std::begin(names), std::end(names), [&](const char* v) { return str == v; });
			if (res == std::end(names)) {
				throw Exception("String \"" + str + "\" does not exist in enum \"" + typeid(T).name() + "\".", HalleyExceptions::Utils);
			}
			return T(res - std::begin(names));
		}

		template<typename T, typename std::enable_if<!std::is_enum<T>::value, int>::type = 0>
		static T fromString(const String& str)
		{
			return T(str);
		}
	};

	template <typename T>
	struct ToStringConverter {
		String operator()(const T& s) const
		{
			return UserConverter::toString(s);
		}
	};

	template <typename T>
	struct FromStringConverter {
		T operator()(const String& s) const
		{
			return UserConverter::fromString<T>(s);
		}
	};

	template<size_t N>
	struct ToStringConverter<char[N]>
	{
		String operator()(const char s[N]) const
		{
			return String(s);
		}
	};

	template<>
	struct ToStringConverter<const char*>
	{
		String operator()(const char* s) const
		{
			return String(s);
		}
	};

	template<>
	struct ToStringConverter<char*>
	{
		String operator()(const char* s) const
		{
			return String(s);
		}
	};

	template<size_t N>
	struct ToStringConverter<wchar_t[N]>
	{
		String operator()(const wchar_t s[N]) const
		{
			return String(s);
		}
	};

	template<>
	struct ToStringConverter<const wchar_t*>
	{
		String operator()(const wchar_t* s) const
		{
			return String(s);
		}
	};

	template<>
	struct ToStringConverter<wchar_t*>
	{
		String operator()(const wchar_t* s) const
		{
			return String(s);
		}
	};

	template<>
	struct ToStringConverter<std::string>
	{
		String operator()(const std::string& s) const
		{
			return String(s);
		}
	};

	template<>
	struct ToStringConverter<bool>
	{
		String operator()(bool s) const
		{
			return String(s ? "true" : "false");
		}
	};

	template<>
	struct FromStringConverter<bool>
	{
		bool operator()(const String& s) const
		{
			return s == "true";
		}
	};

	template <typename T, typename std::enable_if<std::is_floating_point<T>::value, int>::type = 0>
	String toString(T src, int precisionDigits = -1, char decimalSeparator = '.')
	{
		Expects(precisionDigits >= -1 && precisionDigits <= 20);
		std::stringstream str;
		if (precisionDigits != -1) {
			str << std::fixed << std::setprecision(precisionDigits);
		}
		str << src;

		String result;
		if (precisionDigits == -1) {
			result = String::prettyFloat(str.str());
		} else {
			result = str.str();
		}

		if (decimalSeparator != '.') {
			result = result.replaceAll(String("."), String(decimalSeparator));
		}

		return result;
	}

	template <typename T, typename std::enable_if<std::is_integral<T>::value, int>::type = 0>
	String toString(T value, int base = 10, int width = 1)
	{
		Expects(base == 10 || base == 16 || base == 8);
		std::stringstream ss;
		if (base == 16) {
			ss.setf(std::ios::hex, std::ios::basefield);
		} else if (base == 8) {
			ss.setf(std::ios::oct, std::ios::basefield);
		}
		if (width > 1) {
			ss << std::setfill('0') << std::setw(width);
		}
		ss << value;
		return ss.str();
	}

	template <typename T, typename std::enable_if<!std::is_integral<T>::value && !std::is_floating_point<T>::value, int>::type = 0>
	String toString(const T& value)
	{
		return ToStringConverter<typename std::remove_cv<T>::type>()(value);
	}

	template <typename T>
	T fromString(const String& value)
	{
		return FromStringConverter<T>()(value);
	}

	template <typename T, typename std::enable_if<!std::is_same<T, String>::value, int>::type = 0>
	String operator+ (const String& lhp, const T& rhp)
	{
		return lhp + toString<typename std::remove_cv<T>::type>(rhp);
	}

	template <typename T, typename std::enable_if<!std::is_same<T, String>::value, int>::type = 0>
	String operator+ (const T& lhp, const String& rhp)
	{
		return toString<typename std::remove_cv<T>::type>(lhp) + rhp;
	}

}

max_stars_repo_name: sunhay/halley

max_issues_repo_name: sunhay/halley

max_forks_repo_name: sunhay/halley


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

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sum.h>

int
gsl_sum_levin_utrunc_accel (const double *array,
                            const size_t array_size,
                            gsl_sum_levin_utrunc_workspace * w,
                            double *sum_accel, double *abserr_trunc)
{
  return gsl_sum_levin_utrunc_minmax (array, array_size,
                                      0, array_size - 1,
                                      w, sum_accel, abserr_trunc);
}

int
gsl_sum_levin_utrunc_minmax (const double *array,
                             const size_t array_size,
                             const size_t min_terms,
                             const size_t max_terms,
                             gsl_sum_levin_utrunc_workspace * w,
                             double *sum_accel, double *abserr_trunc)
{
  if (array_size == 0)
    {
      *sum_accel = 0.0;
      *abserr_trunc = 0.0;
      w->sum_plain = 0.0;
      w->terms_used = 0;
      return GSL_SUCCESS;
    }
  else if (array_size == 1)
    {
      *sum_accel = array[0];
      *abserr_trunc = GSL_POSINF;
      w->sum_plain = array[0];
      w->terms_used = 1;
      return GSL_SUCCESS;
    }
  else
    {
      const double SMALL = 0.01;
      const size_t nmax = GSL_MAX (max_terms, array_size) - 1;
      double trunc_n = 0.0, trunc_nm1 = 0.0;
      double actual_trunc_n = 0.0, actual_trunc_nm1 = 0.0;
      double result_n = 0.0, result_nm1 = 0.0;
      size_t n;
      int better = 0;
      int before = 0;
      int converging = 0;
      double least_trunc = GSL_DBL_MAX;
      double result_least_trunc;

      /* Calculate specified minimum number of terms. No convergence
         tests are made, and no truncation information is stored. */

      for (n = 0; n < min_terms; n++)
        {
          const double t = array[n];

          result_nm1 = result_n;
          gsl_sum_levin_utrunc_step (t, n, w, &result_n);
        }

      /* Assume the result after the minimum calculation is the best. */

      result_least_trunc = result_n;

      /* Calculate up to maximum number of terms. Check truncation
         condition. */

      for (; n <= nmax; n++)
        {
          const double t = array[n];

          result_nm1 = result_n;
          gsl_sum_levin_utrunc_step (t, n, w, &result_n);

          /* Compute the truncation error directly */

          actual_trunc_nm1 = actual_trunc_n;
          actual_trunc_n = fabs (result_n - result_nm1);

          /* Average results to make a more reliable estimate of the
             real truncation error */

          trunc_nm1 = trunc_n;
          trunc_n = 0.5 * (actual_trunc_n + actual_trunc_nm1);

          /* Determine if we are in the convergence region. */

          better = (trunc_n < trunc_nm1 || trunc_n < SMALL * fabs (result_n));
          converging = converging || (better && before);
          before = better;

          if (converging)
            {
              if (trunc_n < least_trunc)
                {
                  /* Found a low truncation point in the convergence
                     region. Save it. */

                  least_trunc = trunc_n;
                  result_least_trunc = result_n;
                }

              if (fabs (trunc_n / result_n) < 10.0 * GSL_MACH_EPS)
                break;
            }
        }

      if (converging)
        {
          /* Stopped in the convergence region. Return result and
             error estimate. */

          *sum_accel = result_least_trunc;
          *abserr_trunc = least_trunc;
          w->terms_used = n;
          return GSL_SUCCESS;
        }
      else
        {
          /* Never reached the convergence region. Use the last
             calculated values. */

          *sum_accel = result_n;
          *abserr_trunc = trunc_n;
          w->terms_used = n;
          return GSL_SUCCESS;
        }
    }
}

int
gsl_sum_levin_utrunc_step (const double term,
                           const size_t n,
                           gsl_sum_levin_utrunc_workspace * w, double *sum_accel)
{
  if (term == 0.0)
    {
      /* This is actually harmless when treated in this way. A term
         which is exactly zero is simply ignored; the state is not
         changed. We return GSL_EZERODIV as an indicator that this
         occured. */

      return GSL_EZERODIV;
    }
  else if (n == 0)
    {
      *sum_accel = term;
      w->sum_plain = term;
      w->q_den[0] = 1.0 / term;
      w->q_num[0] = 1.0;
      return GSL_SUCCESS;
    }
  else
    {
      double factor = 1.0;
      double ratio = (double) n / (n + 1.0);
      int j;

      w->sum_plain += term;
      w->q_den[n] = 1.0 / (term * (n + 1.0) * (n + 1.0));
      w->q_num[n] = w->sum_plain * w->q_den[n];

      for (j = n - 1; j >= 0; j--)
        {
          double c = factor * (j + 1) / (n + 1);
          factor *= ratio;
          w->q_den[j] = w->q_den[j + 1] - c * w->q_den[j];
          w->q_num[j] = w->q_num[j + 1] - c * w->q_num[j];
        }

      *sum_accel = w->q_num[0] / w->q_den[0];
      return GSL_SUCCESS;
    }
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/*
GENETIC - A simple genetic algorithm.

Copyright 2014, Javier Burguete Tolosa.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this
	list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice,
	this list of conditions and the following disclaimer in the documentation
	and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``AS IS'' AND ANY EXPRESS
OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT SHALL Javier Burguete Tolosa OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/**
 * \file entity.c
 * \brief Source file to define the entity functions.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2014 Javier Burguete Tolosa. All rights reserved.
 */
#define _GNU_SOURCE
#include <gsl/gsl_rng.h>
#include <glib.h>
#include "entity.h"

/**
 * Function to create an entity.
 */
void
entity_new (Entity * entity,    ///< Entity struct.
            unsigned int genome_nbytes,
            ///< Number of bytes of the entity genome.
            unsigned int id)    ///< Identifier number.
{
  entity->id = id;
  // Aligning in 4 bytes
  entity->nbytes = ((genome_nbytes + 3) / 4) * 4;
  entity->genome = (char *) g_slice_alloc (entity->nbytes);
}

/**
 * Function to init randomly the genome of an entity.
 */
void
entity_init (Entity * entity,   ///< Entity struct.
             gsl_rng * rng)     ///< GSL random numbers generator.
{
  unsigned int i;
  for (i = 0; i < entity->nbytes; ++i)
    entity->genome[i] = (char) gsl_rng_uniform_int (rng, 256);
}

/**
 * Function to free the memory used by an entity.
 */
void
entity_free (Entity * entity)   ///< Entity struct.
{
  g_slice_free1 (entity->nbytes, entity->genome);
}

max_stars_repo_name: jburguete/genetic

max_issues_repo_name: jburguete/genetic

max_forks_repo_name: jburguete/genetic


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

#include <config.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_ieee_utils.h>

void
test_bspline (gsl_bspline_workspace * bw)
{
  gsl_vector *B;
  size_t i, j;
  size_t n = 100;
  size_t ncoeffs = gsl_bspline_ncoeffs (bw);
  size_t order = gsl_bspline_order (bw);
  size_t nbreak = gsl_bspline_nbreak (bw);
  double a = gsl_bspline_breakpoint (0, bw);
  double b = gsl_bspline_breakpoint (nbreak - 1, bw);

  B = gsl_vector_alloc (ncoeffs);

  for (i = 0; i < n; i++)
    {
      double xi = a + (b - a) * (i / (n - 1.0));
      double sum = 0;
      gsl_bspline_eval (xi, B, bw);

      for (j = 0; j < ncoeffs; j++)
        {
          double Bj = gsl_vector_get (B, j);
          int s = (Bj < 0 || Bj > 1);
          gsl_test (s,
                    "basis-spline coefficient %u is in range [0,1] for x=%g",
                    j, xi);
          sum += Bj;
        }

      gsl_test_rel (sum, 1.0, order * GSL_DBL_EPSILON,
                    "basis-spline order %u is normalized for x=%g", order,
                    xi);
    }

  gsl_vector_free (B);
}

int
main (int argc, char **argv)
{
  int status = 0;
  size_t order, breakpoints, i;

  gsl_ieee_env_setup ();

  argc = 0;                     /* prevent warnings about unused parameters */
  argv = 0;

  for (order = 1; order < 10; order++)
    {
      for (breakpoints = 2; breakpoints < 100; breakpoints++)
        {
          double a = -1.23 * order, b = 45.6 * order;
          gsl_bspline_workspace *bw = gsl_bspline_alloc (order, breakpoints);
          gsl_bspline_knots_uniform (a, b, bw);
          test_bspline (bw);
          gsl_bspline_free (bw);
        }
    }

  for (order = 1; order < 10; order++)
    {
      for (breakpoints = 2; breakpoints < 100; breakpoints++)
        {
          double a = -1.23 * order, b = 45.6 * order;
          gsl_bspline_workspace *bw = gsl_bspline_alloc (order, breakpoints);
          gsl_vector *k = gsl_vector_alloc (breakpoints);
          for (i = 0; i < breakpoints; i++)
            {
              double f, x;
              f = sqrt (i / (breakpoints - 1.0));
              x = (1 - f) * a + f * b;
              gsl_vector_set (k, i, x);
            };
          gsl_bspline_knots (k, bw);
          test_bspline (bw);
          gsl_vector_free (k);
          gsl_bspline_free (bw);
        }
    }

  exit (gsl_test_summary ());
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


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

#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>

#include "ccl.h"

ccl_a_finder *ccl_a_finder_new(int na, double *a_arr)
{
  if(na<=0)
    return NULL;

  ccl_a_finder *finda=malloc(sizeof(ccl_a_finder));
  if(finda == NULL)
    return NULL;

  finda->ia_last=0;
  finda->amin = a_arr[0];
  finda->amax = a_arr[na-1];
  finda->na = na;
  finda->a_arr = malloc(na*sizeof(double));
  if(finda->a_arr==NULL) {
    free(finda);
    return NULL;
  }

  memcpy(finda->a_arr, a_arr, na*sizeof(double));

  return finda;
}

void ccl_a_finder_free(ccl_a_finder *finda)
{
  if(finda!=NULL) {
    if(finda->na>0)
      free(finda->a_arr);
    free(finda);
  }
}

int ccl_find_a_index(ccl_a_finder *finda, double a)
{
  int ia_0;

  if(a>=finda->amax)
    ia_0=finda->na-1;
  else if(a<=finda->amin)
    ia_0=0;
  else {
    int gotit=0;
    ia_0=finda->ia_last;
    while(!gotit) {
      if(ia_0==0) {
        if(a<finda->a_arr[1])
          gotit=1;
        else
          ia_0++;
      }
      else if(ia_0==finda->na-1) {
        if(a>=finda->a_arr[ia_0-1])
          gotit=1;
        ia_0--;
      }
      else {
        if(a<finda->a_arr[ia_0])
          ia_0--;
        else {
          if(a>=finda->a_arr[ia_0+1])
            ia_0++;
          else
            gotit=1;
        }
      }
    }
  }

  finda->ia_last = ia_0;
  return ia_0;
}

ccl_f3d_t *ccl_f3d_t_copy(ccl_f3d_t *f3d_o, int *status)
{
  int ia, s2dstatus=0;
  ccl_f3d_t *f3d = malloc(sizeof(ccl_f3d_t));
  if (f3d == NULL)
    *status = CCL_ERROR_MEMORY;

  if(*status==0) {
    f3d->lkmin = f3d_o->lkmin;
    f3d->lkmax = f3d_o->lkmax;
    f3d->na = f3d_o->na;
    f3d->is_product = f3d_o->is_product;
    f3d->extrap_linear_growth = f3d_o->extrap_linear_growth;
    f3d->extrap_order_lok = f3d_o->extrap_order_lok;
    f3d->extrap_order_hik = f3d_o->extrap_order_hik;
    f3d->is_log = f3d_o->is_log;
    f3d->growth_factor_0 = f3d_o->growth_factor_0;
    f3d->growth_exponent = f3d_o->growth_exponent;

    f3d->a_arr = malloc(f3d->na*sizeof(double));
    if(f3d->a_arr == NULL)
      *status = CCL_ERROR_MEMORY;
  }

  if(*status==0) {
    memcpy(f3d->a_arr, f3d_o->a_arr, f3d->na*sizeof(double));

    if(f3d_o->fka_1 != NULL)
      f3d->fka_1 = ccl_f2d_t_copy(f3d_o->fka_1, status);
    else
      f3d->fka_1 = NULL;
  }

  if(*status==0) {
    if(f3d_o->fka_2 != NULL)
      f3d->fka_2 = ccl_f2d_t_copy(f3d_o->fka_2, status);
    else
      f3d->fka_2 = NULL;
  }

  if(*status==0) {
    if(f3d_o->tkka != NULL) {
      f3d->tkka = malloc(f3d->na*sizeof(gsl_spline2d));
      if (f3d->tkka == NULL)
        *status = CCL_ERROR_MEMORY;

      if(*status == 0) {
        s2dstatus = 0;
        for(ia=0; ia<f3d->na; ia++) {
          f3d->tkka[ia] = gsl_spline2d_alloc(gsl_interp2d_bicubic,
                                             f3d_o->tkka[ia]->interp_object.xsize,
                                             f3d_o->tkka[ia]->interp_object.ysize);
          if(f3d->tkka[ia] == NULL)
            *status = CCL_ERROR_MEMORY;

          if(*status==0) {
            s2dstatus |= gsl_spline2d_init(f3d->tkka[ia],
                                           f3d_o->tkka[ia]->xarr,
                                           f3d_o->tkka[ia]->yarr,
                                           f3d_o->tkka[ia]->zarr,
                                           f3d_o->tkka[ia]->interp_object.xsize,
                                           f3d_o->tkka[ia]->interp_object.ysize);
            if(s2dstatus)
              *status = CCL_ERROR_SPLINE;
          }
        }
      }
    }
    else
      f3d->tkka = NULL;
  }

  return f3d;
}

ccl_f3d_t *ccl_f3d_t_new(int na,double *a_arr,
                         int nk,double *lk_arr,
                         double *tkka_arr,
                         double *fka1_arr,
                         double *fka2_arr,
                         int is_product,
                         int extrap_order_lok,
                         int extrap_order_hik,
                         ccl_f2d_extrap_growth_t extrap_linear_growth,
                         int is_tkka_log,
                         double growth_factor_0,
                         int growth_exponent,
                         ccl_f2d_interp_t interp_type,
                         int *status) {
  int ia, s2dstatus;
  ccl_f3d_t *f3d = malloc(sizeof(ccl_f3d_t));
  if (f3d == NULL)
    *status = CCL_ERROR_MEMORY;

  if (*status == 0) {
    is_product = is_product || (tkka_arr == NULL);
    f3d->is_product = is_product;
    f3d->extrap_order_lok = extrap_order_lok;
    f3d->extrap_order_hik = extrap_order_hik;
    f3d->extrap_linear_growth = extrap_linear_growth;
    f3d->is_log = is_tkka_log;
    f3d->growth_factor_0 = growth_factor_0;
    f3d->growth_exponent = growth_exponent;
    f3d->fka_1 = NULL;
    f3d->fka_2 = NULL;
    f3d->tkka = NULL;

    f3d->lkmin = lk_arr[0];
    f3d->lkmax = lk_arr[nk-1];
    f3d->na = na;
    f3d->a_arr = malloc(na*sizeof(double));
    if(f3d->a_arr == NULL)
      *status = CCL_ERROR_MEMORY;
  }

  if (*status == 0)
    memcpy(f3d->a_arr, a_arr, na*sizeof(double));

  if ((extrap_order_lok > 1) || (extrap_order_lok < 0) ||
      (extrap_order_hik > 1) || (extrap_order_hik < 0))
    *status = CCL_ERROR_INCONSISTENT;

  if ((extrap_linear_growth != ccl_f2d_cclgrowth) &&
      (extrap_linear_growth != ccl_f2d_constantgrowth) &&
      (extrap_linear_growth != ccl_f2d_no_extrapol))
    *status = CCL_ERROR_INCONSISTENT;

  if (f3d->is_product) {
    if ((fka1_arr == NULL) || (fka2_arr == NULL))
      *status = CCL_ERROR_INCONSISTENT;
  }
  else {
    if (tkka_arr == NULL)
      *status = CCL_ERROR_INCONSISTENT;
  }

  if(*status == 0) {
    if (f3d->is_product) {
      f3d->fka_1 = ccl_f2d_t_new(na, a_arr,
                                 nk, lk_arr,
                                 fka1_arr, NULL, NULL, 0,
                                 extrap_order_lok, extrap_order_hik,
                                 extrap_linear_growth,
                                 is_tkka_log,
                                 growth_factor_0, growth_exponent/2,
                                 interp_type, status);
      f3d->fka_2 = ccl_f2d_t_new(na, a_arr,
                                 nk, lk_arr,
                                 fka2_arr, NULL, NULL, 0,
                                 extrap_order_lok, extrap_order_hik,
                                 extrap_linear_growth,
                                 is_tkka_log,
                                 growth_factor_0, growth_exponent/2,
                                 interp_type, status);
    }
    else {
      switch(interp_type) {
      case(ccl_f2d_3):
        f3d->tkka = malloc(na*sizeof(gsl_spline2d));
        if (f3d->tkka == NULL)
          *status = CCL_ERROR_MEMORY;
        if(*status == 0) {
          for(ia=0; ia<na; ia++) {
            double *tkk = &(tkka_arr[ia*nk*nk]);
            f3d->tkka[ia] = gsl_spline2d_alloc(gsl_interp2d_bicubic, nk, nk);
            if (f3d->tkka[ia] == NULL) {
              *status = CCL_ERROR_MEMORY;
              break;
            }
            s2dstatus = gsl_spline2d_init(f3d->tkka[ia],
                                          lk_arr, lk_arr, tkk, nk, nk);
            if (s2dstatus) {
              *status = CCL_ERROR_SPLINE;
              break;
            }
          }
        }
        break;
      default:
        f3d->tkka = NULL;
      }
    }
  }

  return f3d;
}

double ccl_f3d_t_eval(ccl_f3d_t *f3d,double lk1,double lk2,double a,ccl_a_finder *finda,
                      void *cosmo, int *status) {
  double tkka_post;

  double a_ev = a;
  int is_hiz = a < f3d->a_arr[0];
  int is_loz = a > f3d->a_arr[f3d->na-1];
  if (is_loz) { // Are we above the interpolation range in a?
    if (f3d->extrap_linear_growth == ccl_f2d_no_extrapol) {
      *status=CCL_ERROR_SPLINE_EV;
      return NAN;
    }
    a_ev = f3d->a_arr[f3d->na-1];
  }
  else if (is_hiz) { // Are we below the interpolation range in a?
    if (f3d->extrap_linear_growth == ccl_f2d_no_extrapol) {
      *status=CCL_ERROR_SPLINE_EV;
      return NAN;
    }
    a_ev = f3d->a_arr[0];
  }

  if (f3d->is_product) {
    double fka1 = ccl_f2d_t_eval(f3d->fka_1, lk1, a_ev, cosmo, status);
    double fka2 = ccl_f2d_t_eval(f3d->fka_2, lk2, a_ev, cosmo, status);
    if(*status != 0)
      return NAN;
    tkka_post = fka1*fka2;
  }
  else {
    double lk1_ev = lk1;
    int is_hik1 = lk1 > f3d->lkmax;
    int is_lok1 = lk1 < f3d->lkmin;
    int extrap_k1 = (is_hik1 & (f3d->extrap_order_hik > 0)) ||  (is_lok1 & (f3d->extrap_order_lok > 0));
    if (is_hik1) // Are we above the interpolation range in k?
      lk1_ev = f3d->lkmax;
    else if (is_lok1) // Are we below the interpolation range in k?
      lk1_ev = f3d->lkmin;

    double lk2_ev = lk2;
    int is_hik2 = lk2 > f3d->lkmax;
    int is_lok2 = lk2 < f3d->lkmin;
    int extrap_k2 = (is_hik2 & (f3d->extrap_order_hik > 0)) ||  (is_lok2 & (f3d->extrap_order_lok > 0));
    if (is_hik2) // Are we above the interpolation range in k?
      lk2_ev = f3d->lkmax;
    else if (is_lok2) // Are we below the interpolation range in k?
      lk2_ev = f3d->lkmin;

    int ia = ccl_find_a_index(finda, a_ev);

    if(*status == 0) {
      int spstatus = 0;
      double tkka, dtkka1, dtkka2;
      spstatus |= gsl_spline2d_eval_e(f3d->tkka[ia], lk1_ev, lk2_ev,
                                      NULL, NULL, &tkka);
      if(extrap_k1)
        spstatus |= gsl_spline2d_eval_deriv_x_e(f3d->tkka[ia], lk1_ev, lk2_ev,
                                                NULL, NULL, &dtkka1);
      if(extrap_k2)
        spstatus |= gsl_spline2d_eval_deriv_y_e(f3d->tkka[ia], lk1_ev, lk2_ev,
                                                NULL, NULL, &dtkka2);
      if(ia < f3d->na-1) {
        double tkka_p1, dtkka1_p1, dtkka2_p1;
        double h = (a_ev-f3d->a_arr[ia])/(f3d->a_arr[ia+1]-f3d->a_arr[ia]);

        spstatus |= gsl_spline2d_eval_e(f3d->tkka[ia+1], lk1_ev, lk2_ev,
                                        NULL, NULL, &tkka_p1);
        if(!spstatus)
          tkka = tkka*(1-h) + tkka_p1*h;

        if(extrap_k1) {
          spstatus |= gsl_spline2d_eval_deriv_x_e(f3d->tkka[ia+1], lk1_ev, lk2_ev,
                                                  NULL, NULL, &dtkka1_p1);
          if(!spstatus)
            dtkka1 = dtkka1*(1-h) + dtkka1_p1*h;
        }

        if(extrap_k2) {
          spstatus |= gsl_spline2d_eval_deriv_y_e(f3d->tkka[ia+1], lk1_ev, lk2_ev,
                                                  NULL, NULL, &dtkka2_p1);
          if(!spstatus)
            dtkka2 = dtkka2*(1-h) + dtkka2_p1*h;
        }
      }
      if(spstatus) {
        *status = CCL_ERROR_SPLINE_EV;
        return NAN;
      }

      // Extrapolate if needed
      if(extrap_k1)
        tkka += dtkka1 * (lk1 - lk1_ev);
      if(extrap_k2)
        tkka += dtkka2 * (lk2 - lk2_ev);
      tkka_post = tkka;
    }

    // Exponentiate if needed
    if (f3d->is_log)
      tkka_post = exp(tkka_post);
  }

  // Extrapolate in a if needed
  if (is_hiz) {
    double gz;
    if (f3d->extrap_linear_growth == ccl_f2d_cclgrowth) { // Use CCL's growth function
      ccl_cosmology *csm = (ccl_cosmology *)cosmo;
      if (!csm->computed_growth) {
        *status = CCL_ERROR_GROWTH_INIT;
        ccl_cosmology_set_status_message(
          csm,
          "ccl_f3d.c: ccl_f3d_t_eval(): growth factor splines have not been precomputed!");
        return NAN;
      }
      gz = (
        ccl_growth_factor(csm, a, status) /
        ccl_growth_factor(csm, a_ev, status));
    }
    else // Use constant growth factor
      gz = f3d->growth_factor_0;

    tkka_post *= pow(gz, f3d->growth_exponent);
  }

  return tkka_post;
}

void ccl_f3d_t_free(ccl_f3d_t *f3d)
{
  if(f3d != NULL) {
    if(f3d->fka_1 != NULL)
      ccl_f2d_t_free(f3d->fka_1);
    if(f3d->fka_2 != NULL)
      ccl_f2d_t_free(f3d->fka_2);
    if(f3d->tkka != NULL) {
      int ia;
      for(ia=0; ia<f3d->na; ia++)
        gsl_spline2d_free(f3d->tkka[ia]);
      free(f3d->tkka);
    }
    if(f3d->na > 0)
      free(f3d->a_arr);
    free(f3d);
  }
}

ccl_a_finder *ccl_a_finder_new_from_f3d(ccl_f3d_t *f3d)
{
  return ccl_a_finder_new(f3d->na, f3d->a_arr);
}

max_stars_repo_name: Jappenn/CCL

max_issues_repo_name: Jappenn/CCL

max_forks_repo_name: Jappenn/CCL


text:
//This program takes a double array and returns the roots of the polynomial defined by them. The roots in question can be complex, and the complex data type will need usage in modelica.

#include <stdio.h>
#include <gsl/gsl_poly.h>
void poly_roots(double* a, const size_t size,double *z,size_t size_z);

void poly_roots(double* a, const size_t size,double *z,size_t size_z)
{
  /* coefficients of P(x) =  -1 + x^5  */

  gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc ((size+1));

  gsl_poly_complex_solve (a, size+1, w, z);

  gsl_poly_complex_workspace_free (w);
	for (int i = 0; i < size; i++)
    {
      printf ("z%d = %+.18f %+.18f\n",
              i, z[2*i], z[2*i+1]);
    }

}

//The main function has been provided solely for the sake of testing the function. It will be removed when the problem set is formally passed on.
/*
int main (void)
{
 	double a[6] = { -1, 0, 0, 0, 0, 1 };
	int size_a = 5;
	double *z;
	poly_roots(a, size_a,z, 10);

	return 0;
}

*/

max_stars_repo_name: RITUait/OpenModelica

max_issues_repo_name: RITUait/OpenModelica

max_forks_repo_name: RITUait/OpenModelica


text:
#pragma once

#if SEAL_COMPILER == SEAL_COMPILER_MSVC

// Require Visual Studio 2015 or newer and C++14 support
#if (_MSC_VER < 1900) || (_MSVC_LANG < 201402L)
#error "Microsoft Visual Studio 2015 (14.0) or newer required; C++14 support required"
#endif

// Can work with Visual Studio 2015 (C++14) with some limitations
#if _MSVC_LANG == 1900
#undef SEAL_USE_SHARED_MUTEX_FOR_RW_LOCK
#endif

// Deal with Debug mode in Visual Studio
#ifdef _DEBUG
#define SEAL_DEBUG
#endif

// Cannot use std::shared_mutex when compiling with /clr
#ifdef _M_CEE
#undef SEAL_USE_SHARED_MUTEX_FOR_RW_LOCK
#endif

// Try to check presence of additional headers using __has_include
#ifdef __has_include

// Check for MSGSL
#if __has_include(<gsl/gsl>)
#include <gsl/gsl>
#define SEAL_USE_MSGSL
#else
#undef SEAL_USE_MSGSL
#endif //__has_include(<gsl/gsl>)

#endif

// X64
#ifdef _M_X64

// Use compiler intrinsics for better performance
#define SEAL_USE_INTRIN

#ifdef SEAL_USE_INTRIN
#include <intrin.h>

#pragma intrinsic(_addcarry_u64)
#define SEAL_ADD_CARRY_UINT64(operand1, operand2, carry, result) _addcarry_u64(     \
    carry,                                                                          \
    static_cast<unsigned long long>(operand1),                                      \
    static_cast<unsigned long long>(operand2),                                      \
    reinterpret_cast<unsigned long long*>(result))

#pragma intrinsic(_subborrow_u64)
#define SEAL_SUB_BORROW_UINT64(operand1, operand2, borrow, result) _subborrow_u64(  \
    borrow,                                                                         \
    static_cast<unsigned long long>(operand1),                                      \
    static_cast<unsigned long long>(operand2),                                      \
    reinterpret_cast<unsigned long long*>(result))

#pragma intrinsic(_BitScanReverse64)
#define SEAL_MSB_INDEX_UINT64(result, value) _BitScanReverse64(result, value)

#pragma intrinsic(_umul128)
#define SEAL_MULTIPLY_UINT64(operand1, operand2, result128) {                       \
    result128[0] = _umul128(                                                        \
        static_cast<unsigned long long>(operand1),                                  \
        static_cast<unsigned long long>(operand2),                                  \
        reinterpret_cast<unsigned long long*>(result128 + 1));                      \
}
#define SEAL_MULTIPLY_UINT64_HW64(operand1, operand2, hw64) {                       \
    _umul128(                                                                       \
        static_cast<unsigned long long>(operand1),                                  \
        static_cast<unsigned long long>(operand2),                                  \
        reinterpret_cast<unsigned long long*>(hw64));                               \
}

#endif
#else
#undef SEAL_USE_INTRIN

#endif //_M_X64

#endif

ext: h


text:
/*
Pierwszy program stosuje metode Brenta
gsl_root_fsolver_brent
dla rownania
x^2 - 5 = 0

rozwiazaniem jest x = \sqrt 5 = 2.236068...
*/
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>
#include <string.h>

#include "demo_fn.h"

int
main (int argc, char *argv[])
{
  int status;
  int iter = 0, max_iter = 100;
  const gsl_root_fsolver_type *T;
  gsl_root_fsolver *s;
  double r = 0, r_expected = sqrt (5.0);
  double x_lo = 0.0, x_hi = 5.0;
  gsl_function F;
  struct quadratic_params params = {1.0, -2.0, 1.0};

  F.function = &quadratic;
  F.params = &params;

  if (argc < 2) {
      errno = EINVAL;
      perror("");
      return errno;
  }
  char *method = argv[1];
  if (0 == strcmp("bisection", method)) {
      T = gsl_root_fsolver_bisection;
  } else if (0 == strcmp("secant", method)) {
      T = gsl_root_fsolver_brent;
  } else if (0 == strcmp("steffenson", method)) {
      T = gsl_root_fsolver_falsepos;
  } else {
      errno = EINVAL;
      perror("");
      return errno;
  }

  T = gsl_root_fsolver_bisection;
  s = gsl_root_fsolver_alloc (T);
  gsl_root_fsolver_set (s, &F, x_lo, x_hi);

  printf ("using %s method\n",
          gsl_root_fsolver_name (s));

  printf ("%5s [%9s, %9s] %9s %10s %9s\n",
          "iter", "lower", "upper", "root",
          "err", "err(est)");

  do
    {
      iter++;
      status = gsl_root_fsolver_iterate (s);
      r = gsl_root_fsolver_root (s);
      x_lo = gsl_root_fsolver_x_lower (s);
      x_hi = gsl_root_fsolver_x_upper (s);
      status = gsl_root_test_interval (x_lo, x_hi,
                                       0, 0.001);

      if (status == GSL_SUCCESS)
        printf ("Converged:\n");

      printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n",
              iter, x_lo, x_hi,
              r, r - r_expected,
              x_hi - x_lo);
    }
  while (status == GSL_CONTINUE && iter < max_iter);
  return status;
}

max_stars_repo_name: mistyfiky/agh-mownit

max_issues_repo_name: mistyfiky/agh-mownit

max_forks_repo_name: mistyfiky/agh-mownit


text:
//-*-C++-*-
/***************************************************************************
 *
 *   Copyright (C) 2008 by Paul Demorest
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/

#ifndef __Pulsar_Wavelet_h
#define __Pulsar_Wavelet_h

#include "ReferenceAble.h"

#include <gsl/gsl_wavelet.h>

namespace Pulsar {

  class Profile;

  //! Performs 1-D discrete wavelet transforms (DWT)
  /*! This class is basically a wrapper for the GSL wavelet implementation.
   * Refer to GSL documentation for more info.
   */
  class WaveletTransform : public Reference::Able {

    public:

      //! Constructor
      WaveletTransform();

      //! Destructor
      ~WaveletTransform();

      //! Set wavelet type
      /*! Possible types (as of GSL 1.11) are:
       *

       *
 gsl_wavelet_daubechies
       *
 gsl_wavelet_haar
       *
 gsl_wavelet_bspline
       *

       *  All types can have _centered appended to their name to
       *  get centered versions.
       */
      void set_type(const gsl_wavelet_type *t) { type=t; };

      //! Set wavelet order
      /*! Implemented values (as of GSL 1.11) are:
       *

       *
 daubechies: order=4,6,...,20 (even values only).
       *
 haar: order=2 only.
       *
 bspline: order=103,105,202,204,206,208,301,303,305,307,309.
       *

       */
      void set_order(int o) { order=o; };

      //! Set wavelet type and order using short string notation
      void set_wavelet(std::string s);

      //! Set whether or not to decimate during transform
      void set_decimate(bool _dec=true) { decimate = _dec; }

      //! Perform forward transform
      void transform(unsigned npts, const float *in);

      //! Perform forward transform
      void transform(const std::vector<float>& in);

      //! Perform forward transform
      void transform(const Pulsar::Profile* in);

      //! Perform inverse transform in-place
      void invert();

      //! Possible states of the data
      enum State { Wavelet, Time, Empty };

      //! Get decimation setting
      const bool get_decimate() const { return decimate; }

      //! Get current state
      State get_state() { return state; }

      //! Get wavelet data
      double get_data(int level, int k);

      //! Get data element
      double get_data(int n);

      //! Get pointer to data array
      double* get_data();

      //! Get const pointer to data array
      const double* get_data() const;

      //! Get number of pts currently in use
      const int get_npts() const { return npts; };

      //! Get log2 of number of pts currently in use
      const int get_log2_npts() const { return log2_npts; };

      //! Get total number of wavelet coeffs computed
      const int get_ncoeff() const { return decimate ? npts : log2_npts*npts; }

      //! Get number of coeffs at a given level
      const int get_ncoeff(int level) const;

    protected:

      //! Wavelet type
      const gsl_wavelet_type *type;

      //! Use decimated (standard) or undecimated transform
      bool decimate;

      //! Wavelet order
      int order;

      //! Wavelet coeffs
      gsl_wavelet *wave;

      //! Number of points
      int npts;

      //! Log2 of npts
      int log2_npts;

      //! Data storage
      double *data;

      //! Mean of data
      double mean;

      //! Current state of the data
      State state;

      //! GSL workspace
      gsl_wavelet_workspace *work;

    private:

      //! Free internal memory
      void free_mem();

      //! Allocate internal memory, set up wavelet coeffs
      void init_mem();

  };

}

#endif


text:
/**
 *
 * @file core_dgbtype1cb.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 Azzam Haidar
 * @date 2012-12-15
 * @generated d Tue Jan  7 11:44:49 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

//#define AU(m,n) &(A[(m) + LDA*(n)])
//#define AL(m,n) &(A[(m) + LDA*(n)])
#define AL(m_, n_) (A + NB + LDA * (n_) + ((m_)-(n_)))
#define AU(m_, n_) (A + NB + LDA * (n_) + ((m_)-(n_)+NB))
#define VQ(m)     (VQ + (m))
#define VP(m)     (VP + (m))
#define TAUQ(m)   (TAUQ + (m))
#define TAUP(m)   (TAUP + (m))

/***************************************************************************//**
 *
 * @ingroup CORE_double
 *
 *  CORE_dgbtype1cb is a kernel that will operate on a region (triangle) of data
 *  bounded by st and ed. This kernel eliminate a column by an column-wise
 *  annihiliation, then it apply a left+right update on the hermitian triangle.
 *  Note that the column to be eliminated is located at st-1.
 *
 *  All detail are available on technical report or SC11 paper.
 *  Azzam Haidar, Hatem Ltaief, and Jack Dongarra. 2011.
 *  Parallel reduction to condensed forms for symmetric eigenvalue problems
 *  using aggregated fine-grained and memory-aware kernels. In Proceedings
 *  of 2011 International Conference for High Performance Computing,
 *  Networking, Storage and Analysis (SC '11). ACM, New York, NY, USA, ,
 *  Article 8 , 11 pages.
 *  http://doi.acm.org/10.1145/2063384.2063394
 *
 *******************************************************************************
 *
 * @param[in] N
 *          The order of the matrix A.
 *
 * @param[in] NB
 *          The size of the band.
 *
 * @param[in, out] A
 *          A pointer to the matrix A of size (3*NB+1)-by-N.
 *
 * @param[in] LDA
 *          The leading dimension of the matrix A. LDA >= max(1,3*NB+1)
 *
 * @param[out] V
 *          double array, dimension N if eigenvalue only
 *          requested or (LDV*blkcnt*Vblksiz) if Eigenvectors requested
 *          The Householder reflectors are stored in this array.
 *
 * @param[out] TAU
 *          double array, dimension (N).
 *          The scalar factors of the Householder reflectors are stored
 *          in this array.
 *
 * @param[in] st
 *          A pointer to the start index where this kernel will operate.
 *
 * @param[in] ed
 *          A pointer to the end index where this kernel will operate.
 *
 * @param[in] sweep
 *          The sweep number that is eliminated. it serve to calculate the
 *          pointer to the position where to store the Vs and Ts.
 *
 * @param[in] Vblksiz
 *          constant which correspond to the blocking used when applying the Vs.
 *          it serve to calculate the pointer to the position where to store the
 *          Vs and Ts.
 *
 * @param[in] WANTZ
 *          constant which indicate if Eigenvalue are requested or both
 *          Eigenvalue/Eigenvectors.
 *
 * @param[in] WORK
 *          Workspace of size nb.
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 ******************************************************************************/
/***************************************************************************
 *          TYPE 1-BAND-bidiag Lower/Upper columnwise-Householder
 ***************************************************************************/
void
CORE_dgbtype1cb(PLASMA_enum uplo, int N, int NB,
                double *A, int LDA,
                double *VQ, double *TAUQ,
                double *VP, double *TAUP,
                int st, int ed, int sweep, int Vblksiz, int WANTZ,
                double *WORK)
{
    double ctmp;
    int i, len, LDX, lenj;
    int blkid, vpos, taupos, tpos;
    /* find the pointer to the Vs and Ts as stored by the bulgechasing
     * note that in case no eigenvector required V and T are stored
     * on a vector of size N
     * */
     if( WANTZ == 0 ) {
         vpos   = ((sweep+1)%2)*N + st;
         taupos = ((sweep+1)%2)*N + st;
     } else {
         findVTpos(N, NB, Vblksiz, sweep, st,
                   &vpos, &taupos, &tpos, &blkid);
     }

    LDX = LDA-1;
    len = ed-st+1;

    if( uplo == PlasmaUpper ) {
        /* ========================
         *       UPPER CASE
         * ========================*/
        /* Eliminate the row  at st-1 */
        *VP(vpos) = 1.;
        for(i=1; i<len; i++){
            *VP(vpos+i)     = (*AU(st-1, st+i));
            *AU(st-1, st+i) = 0.;
        }
        ctmp = (*AU(st-1, st));
        LAPACKE_dlarfg_work(len, &ctmp, VP(vpos+1), 1, TAUP(taupos) );
        *AU(st-1, st) = ctmp;
        /* Apply right on A(st:ed,st:ed) */
        ctmp = *TAUP(taupos);
        LAPACKE_dlarfx_work(LAPACK_COL_MAJOR, lapack_const(PlasmaRight),
                            len, len, VP(vpos), ctmp, AU(st, st), LDX, WORK);

        /* Eliminate the created col at st */
        *VQ(vpos) = 1.;
        memcpy( VQ(vpos+1), AU(st+1, st), (len-1)*sizeof(double) );
        memset( AU(st+1, st), 0, (len-1)*sizeof(double) );
        LAPACKE_dlarfg_work(len, AU(st, st), VQ(vpos+1), 1, TAUQ(taupos) );
        lenj = len-1;
        ctmp = (*TAUQ(taupos));
        LAPACKE_dlarfx_work(LAPACK_COL_MAJOR, lapack_const(PlasmaLeft),
                            len, lenj, VQ(vpos), ctmp, AU(st, st+1), LDX, WORK);
    }else{
        /* ========================
         *       LOWER CASE
         * ========================*/
        /* Eliminate the col  at st-1 */
        *VQ(vpos) = 1.;
        memcpy( VQ(vpos+1), AL(st+1, st-1), (len-1)*sizeof(double) );
        memset( AL(st+1, st-1), 0, (len-1)*sizeof(double) );
        LAPACKE_dlarfg_work(len, AL(st, st-1), VQ(vpos+1), 1, TAUQ(taupos) );
        /* Apply left on A(st:ed,st:ed) */
        ctmp = (*TAUQ(taupos));
        LAPACKE_dlarfx_work(LAPACK_COL_MAJOR, lapack_const(PlasmaLeft),
                            len, len, VQ(vpos), ctmp, AL(st, st), LDX, WORK);
        /* Eliminate the created row at st */
        *VP(vpos) = 1.;
        for(i=1; i<len; i++){
            *VP(vpos+i)     = (*AL(st, st+i));
            *AL(st, st+i)   = 0.;
        }
        ctmp = (*AL(st, st));
        LAPACKE_dlarfg_work(len, &ctmp, VP(vpos+1), 1, TAUP(taupos) );
        *AL(st, st) = ctmp;
        lenj = len-1;
        ctmp = (*TAUP(taupos));
        LAPACKE_dlarfx_work(LAPACK_COL_MAJOR, lapack_const(PlasmaRight),
                            lenj, len, VP(vpos), ctmp, AL(st+1, st), LDX, WORK);
    }
    /* end of uplo case */
    return;
}
/***************************************************************************/
#undef AU
#undef AL
#undef VQ
#undef VP
#undef TAUQ
#undef TAUP


text:
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpRNG.h>
#include <gbpMCMC.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_interp.h>

void compute_MCMC_chain_stats(double **x,
                              int      n_x,
                              int      n_avg,
                              double * x_min,
                              double * x_bar_in,
                              double * x_max,
                              double * x_sigma,
                              double **auto_cor,
                              double * slopes,
                              double * dP_sub,
                              double * ln_Pr_in,
                              double * ln_Pr_min,
                              double * ln_Pr_avg,
                              double * ln_Pr_max) {
    int     i_P;
    int     i_avg;
    double *x_bar;
    double *s_temp;
    double  c0, c1, cov00, cov01, cov11, sumsq;
    double  covar;
    double  var;

    // Allocate a temporary array for the mean if it hasn't been given/isn't requested but is needed
    if((x_sigma != NULL || auto_cor != NULL || slopes != NULL) && x_bar_in == NULL)
        x_bar = (double *)SID_malloc(sizeof(double) * n_x);
    else
        x_bar = x_bar_in;

    // Compute the mean (if it's requested/needed)
    if(x_bar != NULL) {
        for(i_P = 0; i_P < n_x; i_P++) {
            x_bar[i_P] = 0.;
            for(i_avg = 0; i_avg < n_avg; i_avg++)
                x_bar[i_P] += x[i_P][i_avg];
            x_bar[i_P] /= (double)n_avg;
        }
    }

    // Compute the std deviation (if it's requested/needed)
    if(x_sigma != NULL) {
        for(i_P = 0; i_P < n_x; i_P++) {
            x_sigma[i_P] = 0.;
            for(i_avg = 0; i_avg < n_avg; i_avg++)
                x_sigma[i_P] += pow(x[i_P][i_avg] - x_bar[i_P], 2.);
            x_sigma[i_P] = sqrt(x_sigma[i_P] / (double)n_avg);
        }
    }

    // Compute the slope across the interval (if it's requested/needed)
    if(slopes != NULL) {
        s_temp = (double *)SID_malloc(sizeof(double) * n_avg);
        for(i_avg = 0; i_avg < n_avg; i_avg++)
            s_temp[i_avg] = (double)i_avg;
        for(i_P = 0; i_P < n_x; i_P++) {
            gsl_fit_linear(s_temp, 1, x[i_P], 1, n_avg, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
            slopes[i_P] = c1;
            if(dP_sub != NULL)
                dP_sub[i_P] = sqrt(sumsq / (double)n_avg);
        }
        SID_free(SID_FARG s_temp);
    }

    // Compute the auto correlation function (if it's requested/needed)
    if(auto_cor != NULL) {
        for(i_P = 0; i_P < n_x; i_P++) {
            for(i_P = 1; i_P < n_avg; i_P++) {
                for(i_avg = 0, var = 0., covar = 0.; i_avg < (n_avg - i_P); i_avg++) {
                    covar += (x[i_P][i_avg] - x_bar[i_P]) * (x[i_P][i_avg + i_P] - x_bar[i_P]);
                    var += (x[i_P][i_avg] - x_bar[i_P]) * (x[i_P][i_avg] - x_bar[i_P]);
                }
                auto_cor[i_P][i_P - 1] = covar / var;
            }
        }
    }

    // Compute ln_Pr statistics
    if(ln_Pr_in != NULL) {
        (*ln_Pr_min) = ln_Pr_in[0];
        (*ln_Pr_avg) = ln_Pr_in[0];
        (*ln_Pr_max) = ln_Pr_in[0];
        for(i_avg = 1; i_avg < n_avg; i_avg++) {
            (*ln_Pr_min) = GBP_MIN((*ln_Pr_min), ln_Pr_in[i_avg]);
            (*ln_Pr_avg) += ln_Pr_in[i_avg];
            (*ln_Pr_max) = GBP_MAX((*ln_Pr_max), ln_Pr_in[i_avg]);
        }
        (*ln_Pr_avg) /= (double)n_avg;
    }

    if(x_bar != x_bar_in)
        SID_free(SID_FARG x_bar);
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
/* Copyright (C) 2021 Atsushi Togo */
/* All rights reserved. */

/* This file is part of phonopy. */

/* Redistribution and use in source and binary forms, with or without */
/* modification, are permitted provided that the following conditions */
/* are met: */

/* * Redistributions of source code must retain the above copyright */
/*   notice, this list of conditions and the following disclaimer. */

/* * Redistributions in binary form must reproduce the above copyright */
/*   notice, this list of conditions and the following disclaimer in */
/*   the documentation and/or other materials provided with the */
/*   distribution. */

/* * Neither the name of the phonopy project nor the names of its */
/*   contributors may be used to endorse or promote products derived */
/*   from this software without specific prior written permission. */

/* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS */
/* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT */
/* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS */
/* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE */
/* COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, */
/* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, */
/* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */
/* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER */
/* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT */
/* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN */
/* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */
/* POSSIBILITY OF SUCH DAMAGE. */

#ifndef __phono3py_H__
#define __phono3py_H__

#if defined(MKL_LAPACKE) || defined(SCIPY_MKL_H)
#include <mkl.h>
#else
#include <lapacke.h>
#endif
#include "phonoc_array.h"

long ph3py_get_interaction(
    Darray *fc3_normal_squared, const char *g_zero, const Darray *frequencies,
    const lapack_complex_double *eigenvectors, const long (*triplets)[3],
    const long num_triplets, const long (*bz_grid_addresses)[3],
    const long D_diag[3], const long Q[3][3], const double *fc3,
    const long is_compact_fc3, const double (*svecs)[3],
    const long multi_dims[2], const long (*multi)[2], const double *masses,
    const long *p2s_map, const long *s2p_map, const long *band_indices,
    const long symmetrize_fc3_q, const double cutoff_frequency);
long ph3py_get_pp_collision(
    double *imag_self_energy,
    const long relative_grid_address[24][4][3], /* thm */
    const double *frequencies, const lapack_complex_double *eigenvectors,
    const long (*triplets)[3], const long num_triplets,
    const long *triplet_weights, const long (*bz_grid_addresses)[3], /* thm */
    const long *bz_map,                                              /* thm */
    const long bz_grid_type, const long D_diag[3], const long Q[3][3],
    const double *fc3, const long is_compact_fc3, const double (*svecs)[3],
    const long multi_dims[2], const long (*multi)[2], const double *masses,
    const long *p2s_map, const long *s2p_map, const Larray *band_indices,
    const Darray *temperatures, const long is_NU, const long symmetrize_fc3_q,
    const double cutoff_frequency);
long ph3py_get_pp_collision_with_sigma(
    double *imag_self_energy, const double sigma, const double sigma_cutoff,
    const double *frequencies, const lapack_complex_double *eigenvectors,
    const long (*triplets)[3], const long num_triplets,
    const long *triplet_weights, const long (*bz_grid_addresses)[3],
    const long D_diag[3], const long Q[3][3], const double *fc3,
    const long is_compact_fc3, const double (*svecs)[3],
    const long multi_dims[2], const long (*multi)[2], const double *masses,
    const long *p2s_map, const long *s2p_map, const Larray *band_indices,
    const Darray *temperatures, const long is_NU, const long symmetrize_fc3_q,
    const double cutoff_frequency);
void ph3py_get_imag_self_energy_at_bands_with_g(
    double *imag_self_energy, const Darray *fc3_normal_squared,
    const double *frequencies, const long (*triplets)[3],
    const long *triplet_weights, const double *g, const char *g_zero,
    const double temperature, const double cutoff_frequency,
    const long num_frequency_points, const long frequency_point_index);
void ph3py_get_detailed_imag_self_energy_at_bands_with_g(
    double *detailed_imag_self_energy, double *imag_self_energy_N,
    double *imag_self_energy_U, const Darray *fc3_normal_squared,
    const double *frequencies, const long (*triplets)[3],
    const long *triplet_weights, const long (*bz_grid_addresses)[3],
    const double *g, const char *g_zero, const double temperature,
    const double cutoff_frequency);
void ph3py_get_real_self_energy_at_bands(
    double *real_self_energy, const Darray *fc3_normal_squared,
    const long *band_indices, const double *frequencies,
    const long (*triplets)[3], const long *triplet_weights,
    const double epsilon, const double temperature,
    const double unit_conversion_factor, const double cutoff_frequency);
void ph3py_get_real_self_energy_at_frequency_point(
    double *real_self_energy, const double frequency_point,
    const Darray *fc3_normal_squared, const long *band_indices,
    const double *frequencies, const long (*triplets)[3],
    const long *triplet_weights, const double epsilon, const double temperature,
    const double unit_conversion_factor, const double cutoff_frequency);
void ph3py_get_collision_matrix(
    double *collision_matrix, const Darray *fc3_normal_squared,
    const double *frequencies, const long (*triplets)[3],
    const long *triplets_map, const long *map_q,
    const long *rotated_grid_points, const double *rotations_cartesian,
    const double *g, const long num_ir_gp, const long num_gp,
    const long num_rot, const double temperature,
    const double unit_conversion_factor, const double cutoff_frequency);
void ph3py_get_reducible_collision_matrix(
    double *collision_matrix, const Darray *fc3_normal_squared,
    const double *frequencies, const long (*triplets)[3],
    const long *triplets_map, const long *map_q, const double *g,
    const long num_gp, const double temperature,
    const double unit_conversion_factor, const double cutoff_frequency);
void ph3py_get_isotope_scattering_strength(
    double *gamma, const long grid_point, const double *mass_variances,
    const double *frequencies, const lapack_complex_double *eigenvectors,
    const long num_grid_points, const long *band_indices, const long num_band,
    const long num_band0, const double sigma, const double cutoff_frequency);
void ph3py_get_thm_isotope_scattering_strength(
    double *gamma, const long grid_point, const long *ir_grid_points,
    const long *weights, const double *mass_variances,
    const double *frequencies, const lapack_complex_double *eigenvectors,
    const long num_ir_grid_points, const long *band_indices,
    const long num_band, const long num_band0,
    const double *integration_weights, const double cutoff_frequency);
void ph3py_distribute_fc3(double *fc3, const long target, const long source,
                          const long *atom_mapping, const long num_atom,
                          const double *rot_cart);
void ph3py_rotate_delta_fc2(double (*fc3)[3][3][3],
                            const double (*delta_fc2s)[3][3],
                            const double *inv_U,
                            const double (*site_sym_cart)[3][3],
                            const long *rot_map_syms, const long num_atom,
                            const long num_site_sym, const long num_disp);
void ph3py_get_permutation_symmetry_fc3(double *fc3, const long num_atom);
void ph3py_get_permutation_symmetry_compact_fc3(
    double *fc3, const long p2s[], const long s2pp[], const long nsym_list[],
    const long perms[], const long n_satom, const long n_patom);
void ph3py_transpose_compact_fc3(double *fc3, const long p2s[],
                                 const long s2pp[], const long nsym_list[],
                                 const long perms[], const long n_satom,
                                 const long n_patom, const long t_type);
long ph3py_get_triplets_reciprocal_mesh_at_q(
    long *map_triplets, long *map_q, const long grid_point, const long mesh[3],
    const long is_time_reversal, const long num_rot,
    const long (*rec_rotations)[3][3], const long swappable);
long ph3py_get_BZ_triplets_at_q(long (*triplets)[3], const long grid_point,
                                const long (*bz_grid_addresses)[3],
                                const long *bz_map, const long *map_triplets,
                                const long num_map_triplets,
                                const long D_diag[3], const long Q[3][3],
                                const long bz_grid_type);
long ph3py_get_integration_weight(
    double *iw, char *iw_zero, const double *frequency_points,
    const long num_band0, const long relative_grid_address[24][4][3],
    const long mesh[3], const long (*triplets)[3], const long num_triplets,
    const long (*bz_grid_addresses)[3], const long *bz_map,
    const long bz_grid_type, const double *frequencies1, const long num_band1,
    const double *frequencies2, const long num_band2, const long tp_type,
    const long openmp_per_triplets, const long openmp_per_bands);
void ph3py_get_integration_weight_with_sigma(
    double *iw, char *iw_zero, const double sigma, const double sigma_cutoff,
    const double *frequency_points, const long num_band0,
    const long (*triplets)[3], const long num_triplets,
    const double *frequencies, const long num_band, const long tp_type);
long ph3py_get_grid_index_from_address(const long address[3],
                                       const long mesh[3]);
void ph3py_get_gr_grid_addresses(long gr_grid_addresses[][3],
                                 const long D_diag[3]);
long ph3py_get_reciprocal_rotations(long rec_rotations[48][3][3],
                                    const long (*rotations)[3][3],
                                    const long num_rot,
                                    const long is_time_reversal);
long ph3py_transform_rotations(long (*transformed_rots)[3][3],
                               const long (*rotations)[3][3],
                               const long num_rot, const long D_diag[3],
                               const long Q[3][3]);
long ph3py_get_snf3x3(long D_diag[3], long P[3][3], long Q[3][3],
                      const long A[3][3]);
long ph3py_transform_rotations(long (*transformed_rots)[3][3],
                               const long (*rotations)[3][3],
                               const long num_rot, const long D_diag[3],
                               const long Q[3][3]);
long ph3py_get_ir_grid_map(long *ir_grid_map, const long D_diag[3],
                           const long PS[3], const long (*grg_rotations)[3][3],
                           const long num_rot);
long ph3py_get_bz_grid_addresses(long (*bz_grid_addresses)[3], long *bz_map,
                                 long *bzg2grg, const long D_diag[3],
                                 const long Q[3][3], const long PS[3],
                                 const double rec_lattice[3][3],
                                 const long type);
long ph3py_rotate_bz_grid_index(const long bz_grid_index,
                                const long rotation[3][3],
                                const long (*bz_grid_addresses)[3],
                                const long *bz_map, const long D_diag[3],
                                const long PS[3], const long bz_grid_type);
void ph3py_symmetrize_collision_matrix(double *collision_matrix,
                                       const long num_column,
                                       const long num_temp,
                                       const long num_sigma);
void ph3py_expand_collision_matrix(double *collision_matrix,
                                   const long *rot_grid_points,
                                   const long *ir_grid_points,
                                   const long num_ir_gp,
                                   const long num_grid_points,
                                   const long num_rot, const long num_sigma,
                                   const long num_temp, const long num_band);
long ph3py_get_neighboring_gird_points(
    long *relative_grid_points, const long *grid_points,
    const long (*relative_grid_address)[3], const long mesh[3],
    const long (*bz_grid_addresses)[3], const long *bz_map,
    const long bz_grid_type, const long num_grid_points,
    const long num_relative_grid_address);
long ph3py_get_thm_integration_weights_at_grid_points(
    double *iw, const double *frequency_points, const long num_band0,
    const long num_band, const long num_gp,
    const long (*relative_grid_address)[4][3], const long D_diag[3],
    const long *grid_points, const long (*bz_grid_addresses)[3],
    const long *bz_map, const long bz_grid_type, const double *frequencies,
    const long *gp2irgp_map, const char function);

#endif


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

/* Author:  G. Jungman */

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

gsl_odeiv_control *
gsl_odeiv_control_alloc(const gsl_odeiv_control_type * T)
{
  gsl_odeiv_control * c =
    (gsl_odeiv_control *) malloc(sizeof(gsl_odeiv_control));

  if(c == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for control struct",
                      GSL_ENOMEM);
    };

  c->type = T;
  c->state = c->type->alloc();

  if (c->state == 0)
    {
      free (c);		/* exception in constructor, avoid memory leak */

      GSL_ERROR_NULL ("failed to allocate space for control state",
                      GSL_ENOMEM);
    };

  return c;
}

int
gsl_odeiv_control_init(gsl_odeiv_control * c,
                       double eps_abs, double eps_rel,
                       double a_y, double a_dydt)
{
  return c->type->init (c->state, eps_abs, eps_rel, a_y, a_dydt);
}

void
gsl_odeiv_control_free(gsl_odeiv_control * c)
{
  c->type->free(c->state);
  free(c);
}

const char *
gsl_odeiv_control_name(const gsl_odeiv_control * c)
{
  return c->type->name;
}

int
gsl_odeiv_control_hadjust (gsl_odeiv_control * c, gsl_odeiv_step * s, const double y0[], const double yerr[], const double dydt[], double * h)
{
  return c->type->hadjust(c->state, s->dimension, s->type->order(s->state),
                          y0, yerr, dydt, h);
}


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_symv (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int order = 101;
   int uplo = 121;
   float alpha = 1.0f;
   float beta = -1.0f;
   int N = 1;
   int lda = 1;
   float A[] = { -0.428f };
   float X[] = { -0.34f };
   int incX = -1;
   float Y[] = { -0.888f };
   int incY = -1;
   float y_expected[] = { 1.03352f };
   cblas_ssymv(order, uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "ssymv(case 1054)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   float alpha = 1.0f;
   float beta = -1.0f;
   int N = 1;
   int lda = 1;
   float A[] = { -0.428f };
   float X[] = { -0.34f };
   int incX = -1;
   float Y[] = { -0.888f };
   int incY = -1;
   float y_expected[] = { 1.03352f };
   cblas_ssymv(order, uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "ssymv(case 1055)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   float alpha = 1.0f;
   float beta = -1.0f;
   int N = 1;
   int lda = 1;
   float A[] = { -0.428f };
   float X[] = { -0.34f };
   int incX = -1;
   float Y[] = { -0.888f };
   int incY = -1;
   float y_expected[] = { 1.03352f };
   cblas_ssymv(order, uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "ssymv(case 1056)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   float alpha = 1.0f;
   float beta = -1.0f;
   int N = 1;
   int lda = 1;
   float A[] = { -0.428f };
   float X[] = { -0.34f };
   int incX = -1;
   float Y[] = { -0.888f };
   int incY = -1;
   float y_expected[] = { 1.03352f };
   cblas_ssymv(order, uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "ssymv(case 1057)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   float alpha = 1.0f;
   float beta = -1.0f;
   int N = 1;
   int lda = 1;
   float A[] = { -0.428f };
   float X[] = { -0.34f };
   int incX = -1;
   float Y[] = { -0.888f };
   int incY = -1;
   float y_expected[] = { 1.03352f };
   cblas_ssymv(order, uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "ssymv(case 1058)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   float alpha = 1.0f;
   float beta = -1.0f;
   int N = 1;
   int lda = 1;
   float A[] = { -0.428f };
   float X[] = { -0.34f };
   int incX = -1;
   float Y[] = { -0.888f };
   int incY = -1;
   float y_expected[] = { 1.03352f };
   cblas_ssymv(order, uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "ssymv(case 1059)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   float alpha = 1.0f;
   float beta = -1.0f;
   int N = 1;
   int lda = 1;
   float A[] = { -0.428f };
   float X[] = { -0.34f };
   int incX = -1;
   float Y[] = { -0.888f };
   int incY = -1;
   float y_expected[] = { 1.03352f };
   cblas_ssymv(order, uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "ssymv(case 1060)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   float alpha = 1.0f;
   float beta = -1.0f;
   int N = 1;
   int lda = 1;
   float A[] = { -0.428f };
   float X[] = { -0.34f };
   int incX = -1;
   float Y[] = { -0.888f };
   int incY = -1;
   float y_expected[] = { 1.03352f };
   cblas_ssymv(order, uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "ssymv(case 1061)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   double alpha = 0;
   double beta = -0.3;
   int N = 1;
   int lda = 1;
   double A[] = { 0.544 };
   double X[] = { -0.601 };
   int incX = -1;
   double Y[] = { -0.852 };
   int incY = -1;
   double y_expected[] = { 0.2556 };
   cblas_dsymv(order, uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dsymv(case 1062)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   double alpha = 0;
   double beta = -0.3;
   int N = 1;
   int lda = 1;
   double A[] = { 0.544 };
   double X[] = { -0.601 };
   int incX = -1;
   double Y[] = { -0.852 };
   int incY = -1;
   double y_expected[] = { 0.2556 };
   cblas_dsymv(order, uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dsymv(case 1063)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   double alpha = 0;
   double beta = -0.3;
   int N = 1;
   int lda = 1;
   double A[] = { 0.544 };
   double X[] = { -0.601 };
   int incX = -1;
   double Y[] = { -0.852 };
   int incY = -1;
   double y_expected[] = { 0.2556 };
   cblas_dsymv(order, uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dsymv(case 1064)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   double alpha = 0;
   double beta = -0.3;
   int N = 1;
   int lda = 1;
   double A[] = { 0.544 };
   double X[] = { -0.601 };
   int incX = -1;
   double Y[] = { -0.852 };
   int incY = -1;
   double y_expected[] = { 0.2556 };
   cblas_dsymv(order, uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dsymv(case 1065)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   double alpha = 0;
   double beta = -0.3;
   int N = 1;
   int lda = 1;
   double A[] = { 0.544 };
   double X[] = { -0.601 };
   int incX = -1;
   double Y[] = { -0.852 };
   int incY = -1;
   double y_expected[] = { 0.2556 };
   cblas_dsymv(order, uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dsymv(case 1066)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   double alpha = 0;
   double beta = -0.3;
   int N = 1;
   int lda = 1;
   double A[] = { 0.544 };
   double X[] = { -0.601 };
   int incX = -1;
   double Y[] = { -0.852 };
   int incY = -1;
   double y_expected[] = { 0.2556 };
   cblas_dsymv(order, uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dsymv(case 1067)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   double alpha = 0;
   double beta = -0.3;
   int N = 1;
   int lda = 1;
   double A[] = { 0.544 };
   double X[] = { -0.601 };
   int incX = -1;
   double Y[] = { -0.852 };
   int incY = -1;
   double y_expected[] = { 0.2556 };
   cblas_dsymv(order, uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dsymv(case 1068)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   double alpha = 0;
   double beta = -0.3;
   int N = 1;
   int lda = 1;
   double A[] = { 0.544 };
   double X[] = { -0.601 };
   int incX = -1;
   double Y[] = { -0.852 };
   int incY = -1;
   double y_expected[] = { 0.2556 };
   cblas_dsymv(order, uplo, N, alpha, A, lda, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dsymv(case 1069)");
     }
   };
  };

}


text:
#include <math.h>
#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <Python.h>
int sum(int* npyArray3D, int npyLength1D, int npyLength2D, int npyLength3D)
{
  int i, j, k;
  int sum = 0;
  double dummy = sqrt(9.0);
  int MAT_DIM = 6;

  gsl_matrix *m = gsl_matrix_alloc(MAT_DIM, MAT_DIM);
  for (i=0; i < MAT_DIM*MAT_DIM; i++) {
    gsl_matrix_set (m, i%MAT_DIM, i/MAT_DIM, (double) npyArray3D[i]);
  }

  for (i=0; i < MAT_DIM*MAT_DIM; i++) {
    printf("%f, ", gsl_matrix_get (m, i%MAT_DIM, i/MAT_DIM));
  }
  printf("\n");

  for (i=0;i<npyLength1D;i++)
    for (j=0;j<npyLength2D;j++)
      for (k=0;k<npyLength3D;k++)
        sum += npyArray3D[i*npyLength3D*npyLength2D + k*npyLength2D + j];

  printf("%f\n",dummy);
  return sum;
}

double get_det(PyObject *A)
{
  int MAT_DIM = 6;
  int i, signum;
  double det;
  int nInts = PyList_Size(A);

  gsl_matrix *m = gsl_matrix_alloc(MAT_DIM, MAT_DIM);
  gsl_permutation *p;

  p = gsl_permutation_alloc(m->size1);

  for (i=0; i<nInts; i++)
  {
    PyObject *oo = PyList_GetItem(A, i);
    gsl_matrix_set (m, i%MAT_DIM, i/MAT_DIM, PyFloat_AS_DOUBLE(oo));
  }

  gsl_linalg_LU_decomp(m, p, &signum);
  det = gsl_linalg_LU_det(m, signum);

  return det;
}

max_stars_repo_name: tcrundall/chronostar

max_issues_repo_name: tcrundall/chronostar

max_forks_repo_name: tcrundall/chronostar


text:
#ifndef GSL_INTERFACE_H
#define GSL_INTERFACE_H

#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>

#ifdef __cplusplus
extern "C" {
#endif

/*****************************
 * Structure definitions
 *****************************
 */

/**
 * Encapsulation of the integrator;
*/
typedef struct cgsl_integrator{

  const gsl_odeiv2_step_type * step_type; /* time step allocation */
  gsl_odeiv2_step      * step;		  /* time step control */
  gsl_odeiv2_control   * control;	  /* definition of the system */
  gsl_odeiv2_system      system;	  /* definition of the driver */
  gsl_odeiv2_evolve    * evolution;	  /* emove forward in time*/
  gsl_odeiv2_driver    * driver;	  /* high level interface */

} cgsl_integrator;

//see <gsl/gsl_odeiv2.h> for an explanation of these
typedef int (* ode_function_ptr ) (double t, const double y[], double dydt[], void * params);
typedef int (* ode_jacobian_ptr ) (double t, const double y[], double * dfdy, double dfdt[], void * params);

/** Pre/post step callbacks
 *      t: Start of the current time step (communicationPoint)
 *     dt: Length of the current time step (communicationStepSize)
 *      n: Size of system
 *      y: For pre, the values of the variables before stepping. For post, after.
 * params: Opaque pointer
 */
typedef int (* pre_post_step_ptr ) (double t, double dt, int n, const double y[], void * params);

enum cgsl_callback_actions{
  CGSL_RESTART = 1                     // multistep integrators need a reset when
                                // used with the epce filter.
};

/**
 *
 *  This is intended to contain everything needed to integrate only one
 *  little bit at a time *and* to be able to backtrack when needed, and to
 *  have multiple instances of the model.
 *
 *  Here we use this for inheritance by aggregation so that each model is
 *  declared as
 *  typedef struct my_model{
 *     cgs_model   m;
 *         .
 *         .
 *     extra stuff
 *         .
 *         .
 *  }  my_model ;
*/

typedef struct cgsl_model{
  int n_variables;
  double *x;  	        /** state variables */
  double *x_backup;     /** for get/set FMU state */
  void * parameters;

  /** Definition of the dynamical system: this assumes an *explicit* ODE */
  ode_function_ptr function;
  /** Jacobian */
  ode_jacobian_ptr jacobian;

  /** Pre/post step functions */
  pre_post_step_ptr pre_step;
  pre_post_step_ptr post_step;

  /** Get/set FMU state
   * Used to copy/retreive internal state to/from temporary storage inside the model
   * params: Opaque pointer
   */
  void (* get_state) (struct cgsl_model *model);
  void (* set_state) (struct cgsl_model *model);

  /** Destructor */
  void (* free) (struct cgsl_model * model);

  /** Needed by ModelExchange to store content from parameters */
  void* (*get_model_parameters)(const struct cgsl_model *model);
} cgsl_model;

/**
 * Useful function for allocating the most common type of model
 */
cgsl_model* cgsl_model_default_alloc(
        int n_variables,            /** Number of variables */
        const double *x0,           /** Initial values. If NULL, initialize model->x to all zeroes instead */
        void *parameters,           /** User pointer */
        ode_function_ptr function,  /** ODE function */
        ode_jacobian_ptr jacobian,  /** Jacobian */
        pre_post_step_ptr pre_step, /** Pre-step function */
        pre_post_step_ptr post_step,/** Post-step function */
        size_t sz                   /** If sz > sizeof(cgsl_model) then allocate sz bytes instead.
                                        Useful for the my_model case described earlier in this file */
);

/** Default destructor. Frees model->x and the model itself */
void cgsl_model_default_free(cgsl_model *model);

/**
 * Finally we can put everything in a bag.  The spefic model only has to
 * fill in these fields.
 */
typedef struct cgsl_simulation {

  cgsl_model * model;         /** this contains the state variables */
  cgsl_integrator i;
  double t;		      /** current time */
  double t1;		      /** stop time */
  double h;		      /** first stepsize and current value */

  int   fixed_step;	      /** whether or not we move at fixed step */
  FILE * file;
  int store_data;	      /** whether or not we save the data in an array */
  double * data;              /** store variables as integration proceeds */
  int buffer_size;            /** size of data storage */
  int n;		      /** number of time steps taken */
  int save;		      /** persistence to file */
  int print;		      /** verbose on stderr */

  int iterations;             /** Number of iterations done by last call to cgsl_step_to() */
} cgsl_simulation;

/*****************************
 * Enum  definitions
 *****************************
 */

/**
 * List of available time integration methods in the gsl_odeiv2 library
 */
enum cgsl_integrator_ids
{
  rk2,		/* 0 */
  rk4,		/* 1 */
  rkf45,	/* 2 */
  rkck,		/* 3 */
  rk8pd,	/* 4 */
  rk1imp,	/* 5 */
  rk2imp,	/* 6 */
  rk4imp,	/* 7 */
  bsimp,	/* 8 */
  msadams,	/* 9 */
  msbdf		/*10 */
};

/* this will flag integrators which need restart */
  extern int restart_integrator;

/**************
 * Function declarations.
 **************
 */

/**
 * Essentially the constructor for the cgsl_simulation object.
 * The dynamical model is defined by a function, its Jacobian, an
 * opaque parameter struct, and a set of initial values.
 * All variables are assumed to be continuous.
 *
 *  \TODO: fix semantics for saving.  Use a filename
 *   instead of descriptor?  Open and close file automatically?
 */
cgsl_simulation cgsl_init_simulation_tolerance(
  cgsl_model * model, /** the model we work on */
        enum cgsl_integrator_ids integrator, /** Integrator ID   */
        double h,                    /** Initial time-step: must be non-zero, even  with variable step*/
        int fixed_step,		     /** Boolean */
        int save,		     /** Boolean */
        int print,		     /** Boolean  */
  FILE *f,		             /** File descriptor if 'save' is enabled  */
  double reltol, double abstol
  );

cgsl_simulation cgsl_init_simulation(
  cgsl_model * model, /** the model we work on */
        enum cgsl_integrator_ids integrator, /** Integrator ID   */
        double h,                    /** Initial time-step: must be non-zero, even  with variable step*/
        int fixed_step,		     /** Boolean */
        int save,		     /** Boolean */
        int print,		     /** Boolean  */
        FILE *f		             /** File descriptor if 'save' is enabled  */
  );

/**
 *  Memory deallocation.
 */
void  cgsl_free_simulation( cgsl_simulation sim );

/**  Step from current time to next communication point. */
int cgsl_step_to(void * _s,  double comm_point, double comm_step ) ;

/** Accessor */
const gsl_odeiv2_step_type * cgsl_get_integrator( int  i ) ;

/**  Commit to file. */
void cgsl_save_data( struct cgsl_simulation * sim );

/** Mutators for fixed step*/
/** \TODO: make this a toggle */
void cgsl_simulation_set_fixed_step( cgsl_simulation * s, double h );
void cgsl_simulation_set_variable_step( cgsl_simulation * s );

/** Get/set FMU state */
void cgsl_simulation_get( cgsl_simulation *s );
void cgsl_simulation_set( cgsl_simulation *s );

#ifdef __cplusplus
}
#endif

#endif

max_stars_repo_name: Tjoppen/fmigo

max_issues_repo_name: Tjoppen/fmigo

max_forks_repo_name: Tjoppen/fmigo


text:
/****************************************
 * MIT License
 *
 * Copyright (c) 2020 Miguel Ramos Pernas
 ****************************************/

#ifdef USE_CPU
#include "Python.h"
#include "math.h"
#include <gsl/gsl_integration.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_vegas.h>

#if NDIM == 1
/// Function proxy to integrate using 1-dimensional methods
double function_proxy(double x, void *vparams) {

  double *params = (double *)vparams;

  return FUNCTION(x, FWD_PARAMS(params)); // Definitions in "evaluators.c"
}
#endif

/// Function proxy to integrate using Monte Carlo methods
double monte_function_proxy(double *data, size_t, void *vparams) {

  double *params = (double *)vparams;

  return FUNCTION(DATA(data),
                  FWD_PARAMS(params)); // Definitions in "evaluators.c"
}

extern "C" {

#if NDIM == 1
/// Exposed function to integrate using the QNG method
PyObject *integrate_qng(double lb, double ub, PyObject *config,
                        double *params) {

  gsl_function func = {&function_proxy, params};

  double atol, rtol;
  PyArg_ParseTuple(config, "dd", &atol, &rtol);

  double res, err;
  size_t neval;
  gsl_integration_qng(&func, lb, ub, atol, rtol, &res, &err, &neval);

  return Py_BuildValue("(ddi)", res, err, neval);
}

/// Exposed function to integrate using the QAG method
PyObject *integrate_qag(double lb, double ub, PyObject *config,
                        double *params) {

  gsl_function func = {&function_proxy, params};

  double atol, rtol;
  int limit, key, workspace_size;
  PyArg_ParseTuple(config, "ddiii", &atol, &rtol, &limit, &key,
                   &workspace_size);

  gsl_integration_workspace *w =
      gsl_integration_workspace_alloc(workspace_size);

  double res, aerr;
  gsl_integration_qag(&func, lb, ub, atol, rtol, limit, key, w, &res, &aerr);

  gsl_integration_workspace_free(w);

  return Py_BuildValue("(dd)", res, aerr);
}

/// Exposed function to integrate using the CQUAD method
PyObject *integrate_cquad(double lb, double ub, PyObject *config,
                          double *params) {

  gsl_function func = {&function_proxy, params};

  double atol, rtol;
  int workspace_size;
  PyArg_ParseTuple(config, "ddi", &atol, &rtol, &workspace_size);

  gsl_integration_cquad_workspace *w =
      gsl_integration_cquad_workspace_alloc(workspace_size);

  double res, err;
  size_t neval;
  gsl_integration_cquad(&func, lb, ub, atol, rtol, w, &res, &err, &neval);

  gsl_integration_cquad_workspace_free(w);

  return Py_BuildValue("(ddi)", res, err, neval);
}
#endif

/// Exposed function to integrate using plain MonteCarlo
PyObject *integrate_plain(double *lb, double *ub, PyObject *config,
                          double *params) {

  gsl_monte_function func = {&monte_function_proxy, NDIM, params};

  double res, err;

  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);

  // Define the state
  gsl_monte_plain_state *s = gsl_monte_plain_alloc(NDIM);

  int calls;
  PyArg_ParseTuple(config, "i", &calls);

  // Calculate the integral
  gsl_monte_plain_integrate(&func, lb, ub, NDIM, calls, r, s, &res, &err);

  gsl_monte_plain_free(s);

  gsl_rng_free(r);

  return Py_BuildValue("(dd)", res, err);
}

/// Exposed function to integrate using the MISER method
PyObject *integrate_miser(double *lb, double *ub, PyObject *config,
                          double *params) {

  gsl_monte_function func = {&monte_function_proxy, NDIM, params};

  double res, err;

  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);

  // Define the state
  gsl_monte_miser_state *s = gsl_monte_miser_alloc(NDIM);

  int calls;
  double estimate_frac;
  int min_calls;
  int min_calls_per_bisection;
  double alpha;
  double dither;

  PyArg_ParseTuple(config, "idiidd", &calls, &estimate_frac, &min_calls,
                   &min_calls_per_bisection, &alpha, &dither);

  gsl_monte_miser_params p;
  gsl_monte_miser_params_get(s, &p);
  p.estimate_frac = estimate_frac;
  p.min_calls = min_calls;
  p.min_calls_per_bisection = min_calls_per_bisection;
  p.alpha = alpha;
  p.dither = dither;
  gsl_monte_miser_params_set(s, &p);

  // Calculate the integral
  gsl_monte_miser_integrate(&func, lb, ub, NDIM, calls, r, s, &res, &err);

  gsl_monte_miser_free(s);

  gsl_rng_free(r);

  return Py_BuildValue("(dd)", res, err);
}

/// Exposed function to integrate using the VEGAS method
PyObject *integrate_vegas(double *lb, double *ub, PyObject *config,
                          double *params) {

  gsl_monte_function func = {&monte_function_proxy, NDIM, params};

  double res, err;

  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);

  // Define the state
  gsl_monte_vegas_state *s = gsl_monte_vegas_alloc(NDIM);

  int calls;
  double alpha;
  int iterations;
  int mode;

  PyArg_ParseTuple(config, "idii", &calls, &alpha, &iterations, &mode);

  gsl_monte_vegas_params p;
  gsl_monte_vegas_params_get(s, &p);
  p.alpha = alpha;
  p.iterations = iterations;
  p.mode = mode;
  gsl_monte_vegas_params_set(s, &p);

  // Calculate the integral
  gsl_monte_vegas_integrate(&func, lb, ub, NDIM, calls, r, s, &res, &err);

  do {
    gsl_monte_vegas_integrate(&func, lb, ub, NDIM, calls, r, s, &res, &err);

  } while (fabs(gsl_monte_vegas_chisq(s) - 1.0) > 0.5);

  gsl_monte_vegas_free(s);

  gsl_rng_free(r);

  return Py_BuildValue("(dd)", res, err);
}
}

#endif // USE_CPU

max_stars_repo_name: mramospe/minkit

max_issues_repo_name: mramospe/minkit

max_forks_repo_name: mramospe/minkit


text:
// This file is part of playd.
// playd is licensed under the MIT licence: see LICENSE.txt.

/**
 * @file
 * Declaration of the Basic_audio, Null_audio and Audio classes.
 * @see audio/audio.cpp
 */

#ifndef PLAYD_AUDIO_H
#define PLAYD_AUDIO_H

#include <chrono>
#include <cstdint>
#include <memory>
#include <string>
#include <utility>
#include <vector>

#undef max
#include <gsl/gsl>

#include "../response.h"
#include "sink.h"
#include "source.h"

namespace Playd::Audio
{
/**
 * An audio item.
 *
 * Audio abstractly represents an audio item that can be played, stopped,
 * and queried for its position and path (or equivalent).
 *
 * Audio is a virtual interface implemented concretely by PipeAudio, and
 * also by mock implementations for testing purposes.
 *
 * @see PipeAudio
 */
class Audio
{
public:
	/// Enumeration of possible states for Audio.
	using State = Sink::State;

	/// Virtual, empty destructor for Audio.
	virtual ~Audio() = default;

	//
	// Control interface
	//

	/**
	 * Performs an update cycle on this Audio.
	 *
	 * Depending on the Audio implementation, this may do actions such as
	 * performing a decoding round, checking for end-of-file, transferring
	 * frames, and so on.
	 *
	 * @return The state of the Audio after updating.
	 * @see State
	 */
	virtual State Update() = 0;

	/**
	 * Sets whether this Audio should be playing or not.
	 * @param playing True for playing; false for stopped.
	 * @exception NoAudioError if the current state is NONE.
	 */
	virtual void SetPlaying(bool playing) = 0;

	/**
	 * Attempts to seek to the given position.
	 * @param position The position to seek to, in microseconds.
	 * @exception NoAudioError if the current state is NONE.
	 * @see Position
	 */
	virtual void SetPosition(std::chrono::microseconds position) = 0;

	//
	// Property access
	//

	/**
	 * This Audio's current file.
	 * @return The filename of this current file.
	 * @exception NoAudioError if the current state is NONE.
	 */
	virtual std::string_view File() const = 0;

	/**
	 * The state of this Audio.
	 * @return this Audio's current state.
	 */
	virtual State CurrentState() const = 0;

	/**
	 * This Audio's current position.
	 *
	 * As this may be executing whilst the playing callback is running,
	 * do not expect it to be highly accurate.
	 *
	 * @return The current position, in microseconds.
	 * @exception NoAudioError if the current state is NONE.
	 * @see Seek
	 */
	virtual std::chrono::microseconds Position() const = 0;

	/**
	 * This Audio's length.
	 *
	 * @return The length, in microseconds.
	 * @exception NoAudioError if the current state is NONE.
	 * @see Seek
	 */
	virtual std::chrono::microseconds Length() const = 0;
};

/**
 * A dummy Audio implementation representing a lack of file.
 *
 * NoAudio throws exceptions if any attempt is made to change, start, or stop
 * the audio, and returns Audio::State::NONE during any attempt to Update.
 * If asked to emit the audio file, NoAudio does nothing.
 *
 * @see Audio
 */
class NullAudio : public Audio
{
public:
	Audio::State Update() override;

	Audio::State CurrentState() const override;

	// The following all raise an exception:

	void SetPlaying(bool playing) override;

	void SetPosition(std::chrono::microseconds position) override;

	std::chrono::microseconds Position() const override;

	std::chrono::microseconds Length() const override;

	std::string_view File() const override;
};

/**
 * A concrete implementation of Audio as a 'pipe'.
 *
 * Basic_audio is comprised of a 'source', which decodes frames from a
 * file, and a 'sink', which plays out the decoded frames.  Updating
 * consists of shifting frames from the source to the sink.
 *
 * @see Audio
 * @see Sink
 * @see Source
 */
class BasicAudio : public Audio
{
public:
	/**
	 * Constructs audio from a source and a sink.
	 * @param src The source of decoded audio frames.
	 * @param sink The target of decoded audio frames.
	 * @see AudioSystem::Load
	 */
	BasicAudio(std::unique_ptr<Source> src, std::unique_ptr<Sink> sink);

	Audio::State Update() override;

	std::string_view File() const override;

	void SetPlaying(bool playing) override;

	Audio::State CurrentState() const override;

	void SetPosition(std::chrono::microseconds position) override;

	std::chrono::microseconds Position() const override;

	std::chrono::microseconds Length() const override;

private:
	/// The source of audio data.
	std::unique_ptr<Source> src;

	/// The sink to which audio data is sent.
	std::unique_ptr<Sink> sink;

	/// The current decoded frame.
	Source::DecodeVector frame;

	/// A span representing the unclaimed part of the decoded frame.
	gsl::span<const std::byte> frame_span;

	/// Clears the current frame and its iterator.
	void ClearFrame();

	/**
	 * Decodes a new frame, if the current frame is empty.
	 * @return True if more frames are available to decode; false
	 *   otherwise.
	 */
	bool DecodeIfFrameEmpty();

	/**
	 * Returns whether the current frame has been finished.
	 * If this is true, then either the frame is empty, or all of the
	 * samples in the frame have been fed to the ringbuffer.
	 * @return True if the frame is finished; false otherwise.
	 */
	bool FrameFinished() const;

	/// Transfers as much of the current frame as possible to the sink.
	void TransferFrame();
};

} // namespace Playd::Audio

#endif // PLAYD_AUDIO_H

max_stars_repo_name: UniversityRadioYork/ury-playd

max_issues_repo_name: UniversityRadioYork/ury-playd

max_forks_repo_name: UniversityRadioYork/ury-playd


text:
#ifndef KGS_GSLQR_H
#define KGS_GSLQR_H

#include "math/QR.h"
#include <gsl/gsl_matrix.h>

class QRGSL: public QR {
 public:
  QRGSL(gsl_matrix* M): QR(M){}

 protected:

  void updateFromMatrix() override;
};

#endif //KGS_GSLQR_H