text:
/* vector/gsl_vector_ulong.h
 *
 * 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.
 */

#ifndef __GSL_VECTOR_ULONG_H__
#define __GSL_VECTOR_ULONG_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_ulong.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  unsigned long *data;
  gsl_block_ulong *block;
  int owner;
}
gsl_vector_ulong;

typedef struct
{
  gsl_vector_ulong vector;
} _gsl_vector_ulong_view;

typedef _gsl_vector_ulong_view gsl_vector_ulong_view;

typedef struct
{
  gsl_vector_ulong vector;
} _gsl_vector_ulong_const_view;

typedef const _gsl_vector_ulong_const_view gsl_vector_ulong_const_view;

/* Allocation */

GSL_FUN gsl_vector_ulong *gsl_vector_ulong_alloc (const size_t n);
GSL_FUN gsl_vector_ulong *gsl_vector_ulong_calloc (const size_t n);

GSL_FUN gsl_vector_ulong *gsl_vector_ulong_alloc_from_block (gsl_block_ulong * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

GSL_FUN gsl_vector_ulong *gsl_vector_ulong_alloc_from_vector (gsl_vector_ulong * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

GSL_FUN void gsl_vector_ulong_free (gsl_vector_ulong * v);

/* Views */

GSL_FUN _gsl_vector_ulong_view
gsl_vector_ulong_view_array (unsigned long *v, size_t n);

GSL_FUN _gsl_vector_ulong_view
gsl_vector_ulong_view_array_with_stride (unsigned long *base,
                                         size_t stride,
                                         size_t n);

GSL_FUN _gsl_vector_ulong_const_view
gsl_vector_ulong_const_view_array (const unsigned long *v, size_t n);

GSL_FUN _gsl_vector_ulong_const_view
gsl_vector_ulong_const_view_array_with_stride (const unsigned long *base,
                                               size_t stride,
                                               size_t n);

GSL_FUN _gsl_vector_ulong_view
gsl_vector_ulong_subvector (gsl_vector_ulong *v,
                            size_t i,
                            size_t n);

GSL_FUN _gsl_vector_ulong_view
gsl_vector_ulong_subvector_with_stride (gsl_vector_ulong *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_FUN _gsl_vector_ulong_const_view
gsl_vector_ulong_const_subvector (const gsl_vector_ulong *v,
                                  size_t i,
                                  size_t n);

GSL_FUN _gsl_vector_ulong_const_view
gsl_vector_ulong_const_subvector_with_stride (const gsl_vector_ulong *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_FUN void gsl_vector_ulong_set_zero (gsl_vector_ulong * v);
GSL_FUN void gsl_vector_ulong_set_all (gsl_vector_ulong * v, unsigned long x);
GSL_FUN int gsl_vector_ulong_set_basis (gsl_vector_ulong * v, size_t i);

GSL_FUN int gsl_vector_ulong_fread (FILE * stream, gsl_vector_ulong * v);
GSL_FUN int gsl_vector_ulong_fwrite (FILE * stream, const gsl_vector_ulong * v);
GSL_FUN int gsl_vector_ulong_fscanf (FILE * stream, gsl_vector_ulong * v);
GSL_FUN int gsl_vector_ulong_fprintf (FILE * stream, const gsl_vector_ulong * v,
                              const char *format);

GSL_FUN int gsl_vector_ulong_memcpy (gsl_vector_ulong * dest, const gsl_vector_ulong * src);

GSL_FUN int gsl_vector_ulong_reverse (gsl_vector_ulong * v);

GSL_FUN int gsl_vector_ulong_swap (gsl_vector_ulong * v, gsl_vector_ulong * w);
GSL_FUN int gsl_vector_ulong_swap_elements (gsl_vector_ulong * v, const size_t i, const size_t j);

GSL_FUN unsigned long gsl_vector_ulong_max (const gsl_vector_ulong * v);
GSL_FUN unsigned long gsl_vector_ulong_min (const gsl_vector_ulong * v);
GSL_FUN void gsl_vector_ulong_minmax (const gsl_vector_ulong * v, unsigned long * min_out, unsigned long * max_out);

GSL_FUN size_t gsl_vector_ulong_max_index (const gsl_vector_ulong * v);
GSL_FUN size_t gsl_vector_ulong_min_index (const gsl_vector_ulong * v);
GSL_FUN void gsl_vector_ulong_minmax_index (const gsl_vector_ulong * v, size_t * imin, size_t * imax);

GSL_FUN int gsl_vector_ulong_add (gsl_vector_ulong * a, const gsl_vector_ulong * b);
GSL_FUN int gsl_vector_ulong_sub (gsl_vector_ulong * a, const gsl_vector_ulong * b);
GSL_FUN int gsl_vector_ulong_mul (gsl_vector_ulong * a, const gsl_vector_ulong * b);
GSL_FUN int gsl_vector_ulong_div (gsl_vector_ulong * a, const gsl_vector_ulong * b);
GSL_FUN int gsl_vector_ulong_scale (gsl_vector_ulong * a, const unsigned long x);
GSL_FUN int gsl_vector_ulong_add_constant (gsl_vector_ulong * a, const double x);
GSL_FUN int gsl_vector_ulong_axpby (const unsigned long alpha, const gsl_vector_ulong * x, const unsigned long beta, gsl_vector_ulong * y);
GSL_FUN unsigned long gsl_vector_ulong_sum (const gsl_vector_ulong * a);

GSL_FUN int gsl_vector_ulong_equal (const gsl_vector_ulong * u,
                            const gsl_vector_ulong * v);

GSL_FUN int gsl_vector_ulong_isnull (const gsl_vector_ulong * v);
GSL_FUN int gsl_vector_ulong_ispos (const gsl_vector_ulong * v);
GSL_FUN int gsl_vector_ulong_isneg (const gsl_vector_ulong * v);
GSL_FUN int gsl_vector_ulong_isnonneg (const gsl_vector_ulong * v);

GSL_FUN INLINE_DECL unsigned long gsl_vector_ulong_get (const gsl_vector_ulong * v, const size_t i);
GSL_FUN INLINE_DECL void gsl_vector_ulong_set (gsl_vector_ulong * v, const size_t i, unsigned long x);
GSL_FUN INLINE_DECL unsigned long * gsl_vector_ulong_ptr (gsl_vector_ulong * v, const size_t i);
GSL_FUN INLINE_DECL const unsigned long * gsl_vector_ulong_const_ptr (const gsl_vector_ulong * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
unsigned long
gsl_vector_ulong_get (const gsl_vector_ulong * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

INLINE_FUN
void
gsl_vector_ulong_set (gsl_vector_ulong * v, const size_t i, unsigned long x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

INLINE_FUN
unsigned long *
gsl_vector_ulong_ptr (gsl_vector_ulong * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (unsigned long *) (v->data + i * v->stride);
}

INLINE_FUN
const unsigned long *
gsl_vector_ulong_const_ptr (const gsl_vector_ulong * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const unsigned long *) (v->data + i * v->stride);
}
#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_ULONG_H__ */

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/* specfunc/gsl_sf_bessel.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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 */

#ifndef __GSL_SF_BESSEL_H__
#define __GSL_SF_BESSEL_H__

#include <stdlib.h>
#include <gsl/gsl_mode.h>
#include <gsl/gsl_precision.h>
#include <gsl/gsl_sf_result.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

/* Regular Bessel Function J_0(x)
 *
 * exceptions: none
 */
int gsl_sf_bessel_J0_e(const double x,  gsl_sf_result * result);
double gsl_sf_bessel_J0(const double x);

/* Regular Bessel Function J_1(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
int gsl_sf_bessel_J1_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_J1(const double x);

/* Regular Bessel Function J_n(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
int gsl_sf_bessel_Jn_e(int n, double x, gsl_sf_result * result);
double gsl_sf_bessel_Jn(const int n, const double x);

/* Regular Bessel Function J_n(x),  nmin <= n <= nmax
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
int gsl_sf_bessel_Jn_array(int nmin, int nmax, double x, double * result_array);

/* Irregular Bessel function Y_0(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
int gsl_sf_bessel_Y0_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_Y0(const double x);

/* Irregular Bessel function Y_1(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
int gsl_sf_bessel_Y1_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_Y1(const double x);

/* Irregular Bessel function Y_n(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
int gsl_sf_bessel_Yn_e(int n,const double x, gsl_sf_result * result);
double gsl_sf_bessel_Yn(const int n,const double x);

/* Irregular Bessel function Y_n(x), nmin <= n <= nmax
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
int gsl_sf_bessel_Yn_array(const int nmin, const int nmax, const double x, double * result_array);

/* Regular modified Bessel function I_0(x)
 *
 * exceptions: GSL_EOVRFLW
 */
int gsl_sf_bessel_I0_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_I0(const double x);

/* Regular modified Bessel function I_1(x)
 *
 * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW
 */
int gsl_sf_bessel_I1_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_I1(const double x);

/* Regular modified Bessel function I_n(x)
 *
 * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW
 */
int gsl_sf_bessel_In_e(const int n, const double x, gsl_sf_result * result);
double gsl_sf_bessel_In(const int n, const double x);

/* Regular modified Bessel function  I_n(x) for n=nmin,...,nmax
 *
 * nmin >=0, nmax >= nmin
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
int gsl_sf_bessel_In_array(const int nmin, const int nmax, const double x, double * result_array);

/* Scaled regular modified Bessel function
 *  exp(-|x|) I_0(x)
 *
 * exceptions: none
 */
int gsl_sf_bessel_I0_scaled_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_I0_scaled(const double x);

/* Scaled regular modified Bessel function
 *  exp(-|x|) I_1(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
int gsl_sf_bessel_I1_scaled_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_I1_scaled(const double x);

/* Scaled regular modified Bessel function
 *  exp(-|x|) I_n(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
int gsl_sf_bessel_In_scaled_e(int n, const double x, gsl_sf_result * result);
double gsl_sf_bessel_In_scaled(const int n, const double x);

/* Scaled regular modified Bessel function
 *  exp(-|x|) I_n(x)  for n=nmin,...,nmax
 *
 * nmin >=0, nmax >= nmin
 * exceptions: GSL_EUNDRFLW
 */
int gsl_sf_bessel_In_scaled_array(const int nmin, const int nmax, const double x, double * result_array);

/* Irregular modified Bessel function K_0(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
int gsl_sf_bessel_K0_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_K0(const double x);

/* Irregular modified Bessel function K_1(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
int gsl_sf_bessel_K1_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_K1(const double x);

/* Irregular modified Bessel function K_n(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
int gsl_sf_bessel_Kn_e(const int n, const double x, gsl_sf_result * result);
double gsl_sf_bessel_Kn(const int n, const double x);

/* Irregular modified Bessel function  K_n(x)  for n=nmin,...,nmax
 *
 * x > 0.0, nmin >=0, nmax >= nmin
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
int gsl_sf_bessel_Kn_array(const int nmin, const int nmax, const double x, double * result_array);

/* Scaled irregular modified Bessel function
 *  exp(x) K_0(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM
 */
int gsl_sf_bessel_K0_scaled_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_K0_scaled(const double x);

/* Scaled irregular modified Bessel function
 *  exp(x) K_1(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
int gsl_sf_bessel_K1_scaled_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_K1_scaled(const double x);

/* Scaled irregular modified Bessel function
 *  exp(x) K_n(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
int gsl_sf_bessel_Kn_scaled_e(int n, const double x, gsl_sf_result * result);
double gsl_sf_bessel_Kn_scaled(const int n, const double x);

/* Scaled irregular modified Bessel function  exp(x) K_n(x)  for n=nmin,...,nmax
 *
 * x > 0.0, nmin >=0, nmax >= nmin
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
int gsl_sf_bessel_Kn_scaled_array(const int nmin, const int nmax, const double x, double * result_array);

/* Regular spherical Bessel function j_0(x) = sin(x)/x
 *
 * exceptions: none
 */
int gsl_sf_bessel_j0_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_j0(const double x);

/* Regular spherical Bessel function j_1(x) = (sin(x)/x - cos(x))/x
 *
 * exceptions: GSL_EUNDRFLW
 */
int gsl_sf_bessel_j1_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_j1(const double x);

/* Regular spherical Bessel function j_2(x) = ((3/x^2 - 1)sin(x) - 3cos(x)/x)/x
 *
 * exceptions: GSL_EUNDRFLW
 */
int gsl_sf_bessel_j2_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_j2(const double x);

/* Regular spherical Bessel function j_l(x)
 *
 * l >= 0, x >= 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
int gsl_sf_bessel_jl_e(const int l, const double x, gsl_sf_result * result);
double gsl_sf_bessel_jl(const int l, const double x);

/* Regular spherical Bessel function j_l(x) for l=0,1,...,lmax
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
int gsl_sf_bessel_jl_array(const int lmax, const double x, double * result_array);

/* Regular spherical Bessel function j_l(x) for l=0,1,...,lmax
 * Uses Steed's method.
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
int gsl_sf_bessel_jl_steed_array(const int lmax, const double x, double * jl_x_array);

/* Irregular spherical Bessel function y_0(x)
 *
 * exceptions: none
 */
int gsl_sf_bessel_y0_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_y0(const double x);

/* Irregular spherical Bessel function y_1(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
int gsl_sf_bessel_y1_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_y1(const double x);

/* Irregular spherical Bessel function y_2(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
int gsl_sf_bessel_y2_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_y2(const double x);

/* Irregular spherical Bessel function y_l(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
int gsl_sf_bessel_yl_e(int l, const double x, gsl_sf_result * result);
double gsl_sf_bessel_yl(const int l, const double x);

/* Irregular spherical Bessel function y_l(x) for l=0,1,...,lmax
 *
 * exceptions: GSL_EUNDRFLW
 */
int gsl_sf_bessel_yl_array(const int lmax, const double x, double * result_array);

/* Regular scaled modified spherical Bessel function
 *
 * Exp[-|x|] i_0(x)
 *
 * exceptions: none
 */
int gsl_sf_bessel_i0_scaled_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_i0_scaled(const double x);

/* Regular scaled modified spherical Bessel function
 *
 * Exp[-|x|] i_1(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
int gsl_sf_bessel_i1_scaled_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_i1_scaled(const double x);

/* Regular scaled modified spherical Bessel function
 *
 * Exp[-|x|] i_2(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
int gsl_sf_bessel_i2_scaled_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_i2_scaled(const double x);

/* Regular scaled modified spherical Bessel functions
 *
 * Exp[-|x|] i_l(x)
 *
 * i_l(x) = Sqrt[Pi/(2x)] BesselI[l+1/2,x]
 *
 * l >= 0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
int gsl_sf_bessel_il_scaled_e(const int l, double x, gsl_sf_result * result);
double gsl_sf_bessel_il_scaled(const int l, const double x);

/* Regular scaled modified spherical Bessel functions
 *
 * Exp[-|x|] i_l(x)
 * for l=0,1,...,lmax
 *
 * exceptions: GSL_EUNDRFLW
 */
int gsl_sf_bessel_il_scaled_array(const int lmax, const double x, double * result_array);

/* Irregular scaled modified spherical Bessel function
 * Exp[x] k_0(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
int gsl_sf_bessel_k0_scaled_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_k0_scaled(const double x);

/* Irregular modified spherical Bessel function
 * Exp[x] k_1(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW, GSL_EOVRFLW
 */
int gsl_sf_bessel_k1_scaled_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_k1_scaled(const double x);

/* Irregular modified spherical Bessel function
 * Exp[x] k_2(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW, GSL_EOVRFLW
 */
int gsl_sf_bessel_k2_scaled_e(const double x, gsl_sf_result * result);
double gsl_sf_bessel_k2_scaled(const double x);

/* Irregular modified spherical Bessel function
 * Exp[x] k_l[x]
 *
 * k_l(x) = Sqrt[Pi/(2x)] BesselK[l+1/2,x]
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
int gsl_sf_bessel_kl_scaled_e(int l, const double x, gsl_sf_result * result);
double gsl_sf_bessel_kl_scaled(const int l, const double x);

/* Irregular scaled modified spherical Bessel function
 * Exp[x] k_l(x)
 *
 * for l=0,1,...,lmax
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
int gsl_sf_bessel_kl_scaled_array(const int lmax, const double x, double * result_array);

/* Regular cylindrical Bessel function J_nu(x)
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
int gsl_sf_bessel_Jnu_e(const double nu, const double x, gsl_sf_result * result);
double gsl_sf_bessel_Jnu(const double nu, const double x);

/* Irregular cylindrical Bessel function Y_nu(x)
 *
 * exceptions:
 */
int gsl_sf_bessel_Ynu_e(double nu, double x, gsl_sf_result * result);
double gsl_sf_bessel_Ynu(const double nu, const double x);

/* Regular cylindrical Bessel function J_nu(x)
 * evaluated at a series of x values. The array
 * contains the x values. They are assumed to be
 * strictly ordered and positive. The array is
 * over-written with the values of J_nu(x_i).
 *
 * exceptions: GSL_EDOM, GSL_EINVAL
 */
int gsl_sf_bessel_sequence_Jnu_e(double nu, gsl_mode_t mode, size_t size, double * v);

/* Scaled modified cylindrical Bessel functions
 *
 * Exp[-|x|] BesselI[nu, x]
 * x >= 0, nu >= 0
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_bessel_Inu_scaled_e(double nu, double x, gsl_sf_result * result);
double gsl_sf_bessel_Inu_scaled(double nu, double x);

/* Modified cylindrical Bessel functions
 *
 * BesselI[nu, x]
 * x >= 0, nu >= 0
 *
 * exceptions: GSL_EDOM, GSL_EOVRFLW
 */
int gsl_sf_bessel_Inu_e(double nu, double x, gsl_sf_result * result);
double gsl_sf_bessel_Inu(double nu, double x);

/* Scaled modified cylindrical Bessel functions
 *
 * Exp[+|x|] BesselK[nu, x]
 * x > 0, nu >= 0
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_bessel_Knu_scaled_e(const double nu, const double x, gsl_sf_result * result);
double gsl_sf_bessel_Knu_scaled(const double nu, const double x);

/* Modified cylindrical Bessel functions
 *
 * BesselK[nu, x]
 * x > 0, nu >= 0
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
int gsl_sf_bessel_Knu_e(const double nu, const double x, gsl_sf_result * result);
double gsl_sf_bessel_Knu(const double nu, const double x);

/* Logarithm of modified cylindrical Bessel functions.
 *
 * Log[BesselK[nu, x]]
 * x > 0, nu >= 0
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_bessel_lnKnu_e(const double nu, const double x, gsl_sf_result * result);
double gsl_sf_bessel_lnKnu(const double nu, const double x);

/* s'th positive zero of the Bessel function J_0(x).
 *
 * exceptions:
 */
int gsl_sf_bessel_zero_J0_e(unsigned int s, gsl_sf_result * result);
double gsl_sf_bessel_zero_J0(unsigned int s);

/* s'th positive zero of the Bessel function J_1(x).
 *
 * exceptions:
 */
int gsl_sf_bessel_zero_J1_e(unsigned int s, gsl_sf_result * result);
double gsl_sf_bessel_zero_J1(unsigned int s);

/* s'th positive zero of the Bessel function J_nu(x).
 *
 * exceptions:
 */
int gsl_sf_bessel_zero_Jnu_e(double nu, unsigned int s, gsl_sf_result * result);
double gsl_sf_bessel_zero_Jnu(double nu, unsigned int s);

__END_DECLS

#endif /* __GSL_SF_BESSEL_H__ */

ext: h

max_forks_repo_name: gersteinlab/LESSeq


text:
/* eigen/jacobi.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
 */
/* Simple linear algebra operations, operating directly
 * on the gsl_vector and gsl_matrix objects. These are
 * meant for "generic" and "small" systems. Anyone
 * interested in large systems will want to use more
 * sophisticated methods, presumably involving native
 * BLAS operations, specialized data representations,
 * or other optimizations.
 */
#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include "gsl_eigen.h"

#define REAL double

inline static void
jac_rotate(gsl_matrix * a,
           unsigned int i, unsigned int j, unsigned int k, unsigned int l,
           double * g, double * h,
           double s, double tau)
{
  *g = gsl_matrix_get(a, i, j);
  *h = gsl_matrix_get(a, k, l);
  gsl_matrix_set(a, i, j, (*g) - s*((*h) + (*g)*tau));
  gsl_matrix_set(a, k, l, (*h) + s*((*g) - (*h)*tau));
}

int
gsl_eigen_jacobi(gsl_matrix * a,
                         gsl_vector * eval,
                         gsl_matrix * evec,
                         unsigned int max_rot,
                         unsigned int * nrot)
{
  if(a->size1 != a->size2) {
     GSL_ERROR ("eigenproblem requires square matrix", GSL_ENOTSQR);
  }
  else if(a->size1 != evec->size1 || a->size1 != evec->size2) {
     GSL_ERROR ("eigenvector matrix must match input matrix", GSL_EBADLEN);
  }
  else if(a->size1 != eval->size) {
    GSL_ERROR ("eigenvalue vector must match input matrix", GSL_EBADLEN);
  }
  else {
    const unsigned int n = a->size1;
    unsigned int i, j, iq, ip;
    double t, s;

    REAL * b = (REAL *) malloc(n * sizeof(REAL));
    REAL * z = (REAL *) malloc(n * sizeof(REAL));
    if(b == 0 || z == 0) {
      if(b != 0) free(b);
      if(z != 0) free(z);
      GSL_ERROR ("could not allocate memory for workspace", GSL_ENOMEM);
    }

    /* Set eigenvectors to coordinate basis. */
    for(ip=0; ip<n; ip++) {
      for(iq=0; iq<n; iq++) {
        gsl_matrix_set(evec, ip, iq, 0.0);
      }
      gsl_matrix_set(evec, ip, ip, 1.0);
    }

    /* Initialize eigenvalues and workspace. */
    for(ip=0; ip<n; ip++) {
      REAL a_ipip = gsl_matrix_get(a, ip, ip);
      z[ip] = 0.0;
      b[ip] = a_ipip;
      gsl_vector_set(eval, ip, a_ipip);
    }

    *nrot = 0;

    for(i=1; i<=max_rot; i++) {
      REAL thresh;
      REAL tau;
      REAL g, h, c;
      REAL sm = 0.0;
      for(ip=0; ip<n-1; ip++) {
        for(iq=ip+1; iq<n; iq++) {
          sm += fabs(gsl_matrix_get(a, ip, iq));
        }
      }
      if(sm == 0.0) {
        free(z);
        free(b);
        return GSL_SUCCESS;
      }

      if(i < 4)
        thresh = 0.2*sm/(n*n);
      else
        thresh = 0.0;

      for(ip=0; ip<n-1; ip++) {
        for(iq=ip+1; iq<n; iq++) {
          const REAL d_ip = gsl_vector_get(eval, ip);
          const REAL d_iq = gsl_vector_get(eval, iq);
	  const REAL a_ipiq = gsl_matrix_get(a, ip, iq);
          g = 100.0 * fabs(a_ipiq);
          if(   i > 4
             && fabs(d_ip)+g == fabs(d_ip)
             && fabs(d_iq)+g == fabs(d_iq)
	     ) {
            gsl_matrix_set(a, ip, iq, 0.0);
          }
          else if(fabs(a_ipiq) > thresh) {
            h = d_iq - d_ip;
            if(fabs(h) + g == fabs(h)) {
              t = a_ipiq/h;
            }
            else {
              REAL theta = 0.5*h/a_ipiq;
              t = 1.0/(fabs(theta) + sqrt(1.0 + theta*theta));
              if(theta < 0.0) t = -t;
            }

            c   = 1.0/sqrt(1.0+t*t);
            s   = t*c;
            tau = s/(1.0+c);
            h   = t * a_ipiq;
            z[ip] -= h;
            z[iq] += h;
	    gsl_vector_set(eval, ip, d_ip - h);
	    gsl_vector_set(eval, iq, d_iq + h);
	    gsl_matrix_set(a, ip, iq, 0.0);

            for(j=0; j<ip; j++){
              jac_rotate(a, j, ip, j, iq, &g, &h, s, tau);
            }
            for(j=ip+1; j<iq; j++){
              jac_rotate(a, ip, j, j, iq, &g, &h, s, tau);
            }
            for(j=iq+1; j<n; j++){
              jac_rotate(a, ip, j, iq, j, &g, &h, s, tau);
            }
            for (j=0; j<n; j++){
              jac_rotate(evec, j, ip, j, iq, &g, &h, s, tau);
            }
            ++(*nrot);
          }
        }
      }
      for (ip=0; ip<n; ip++) {
        b[ip] += z[ip];
        z[ip]  = 0.0;
	gsl_vector_set(eval, ip, b[ip]);
      }

      /* continue iteration */
    }

    return GSL_EMAXITER;
  }
}

int
gsl_eigen_invert_jacobi(const gsl_matrix * a,
                             gsl_matrix * ainv,
                             unsigned int max_rot)
{
  if(a->size1 != a->size2 || ainv->size1 != ainv->size2) {
    return GSL_ENOTSQR;
  }
  else if(a->size1 != ainv->size2) {
    return GSL_EBADLEN;
  }
  else {
    const unsigned int n = a->size1;
    unsigned int nrot;
    unsigned int i,j,k,l;

    /* This is annoying because I do not want
     * the error handling in these functions.
     * But there are no "impl"-like versions
     * of these allocators... sigh.
     */
    gsl_vector * eval = gsl_vector_alloc(n);
    gsl_matrix * evec = gsl_matrix_alloc(n, n);
    gsl_matrix * inv_diag = gsl_matrix_alloc(n, n);

    if(eval == 0 || evec == 0 || inv_diag == 0) {
      if(eval != 0) gsl_vector_free(eval);
      if(evec != 0) gsl_matrix_free(evec);
      if(inv_diag != 0) gsl_matrix_free(inv_diag);
      return GSL_ENOMEM;
    }

    memcpy(ainv->data, a->data, n*n*sizeof(REAL));

    gsl_eigen_jacobi(ainv, eval, evec, max_rot, &nrot);

    for(i=0; i<n; i++) {
      if(fabs(gsl_vector_get(eval, i)) < 100.0 * GSL_DBL_EPSILON) {
        /* apparent singularity */
        gsl_vector_free(eval);
        gsl_matrix_free(evec);
        gsl_matrix_free(inv_diag);
        return GSL_ESING;
      }
    }

    /* Invert the diagonalized matrix. */
    for(i=0; i<n; i++) {
      for(j=0; j<n; j++) {
        gsl_matrix_set(inv_diag, i, j, 0.0);
      }
      gsl_matrix_set(inv_diag, i, i, 1.0/gsl_vector_get(eval, i));
    }

    for(i=0; i<n; i++) {
      for(j=0; j<n; j++) {
        gsl_matrix_set(ainv, i, j, 0.0);
        for(k=0; k<n; k++) {
          for(l=0; l<n; l++) {
	    REAL ainv_ij = gsl_matrix_get(ainv, i, j);
	    REAL evec_il = gsl_matrix_get(evec, i, l);
	    REAL evec_jk = gsl_matrix_get(evec, j, k);
	    REAL inv_diag_lk = gsl_matrix_get(inv_diag, l, k);
	    REAL delta = evec_il * inv_diag_lk * evec_jk;
	    gsl_matrix_set(ainv, i, j, ainv_ij + delta);
	  }
	}
      }
    }

    gsl_vector_free(eval);
    gsl_matrix_free(evec);
    gsl_matrix_free(inv_diag);
    return GSL_SUCCESS;
  }
}

ext: c


text:
//Gets deltas (first order differences) of X.
//Only the deltas are output in Y.

//I implement this just like FIR for speed and shorter code,
//except non-causal and mid-sample of B is 0,
//so I don't explicitly make B (e.g., B[n] just equals sc*n).

//It seems that Kaldi and others use the same method for edge samples,
//which is numpy.edge, so that is used here.
//Comment out those lines to set out-of-range samps to 0.

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

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

int get_deltas_s (float *Y, const float *X, const int iscolmajor, const int R, const int C, const int dim, const int N);
int get_deltas_d (double *Y, const double *X, const int iscolmajor, const int R, const int C, const int dim, const int N);

int get_deltas_s (float *Y, const float *X, const int iscolmajor, const int R, const int C, const int dim, const int N)
{
    const float z = 0.0f;
    int r, c, n;
    float sc = 1.0f;

    //Checks
    if (R<1) { fprintf(stderr,"error in get_deltas_s: R (nrows Y) must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in get_deltas_s: C (ncols Y) must be positive\n"); return 1; }
    if (N<1) { fprintf(stderr,"error in get_deltas_s: N (delta winlength) must be positive\n"); return 1; }

    //Get sc (normalizer)
    for (n=2; n<=N; n++) { sc += n*n; }
    sc = 0.5f/sc;

    //Initialize Y
    cblas_scopy(R*C,&z,0,&Y[0],1);

    if (dim==0)
    {
        if (iscolmajor)
        {
            for (c=0; c<C; c++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[c*R],0,&Y[c*R],1);        //beg edge samps
                    cblas_saxpy(R-n,-sc*n,&X[c*R],1,&Y[c*R+n],1);    //past samps
                    cblas_saxpy(R-n,sc*n,&X[c*R+n],1,&Y[c*R],1);     //future samps
                    cblas_saxpy(n,sc*n,&X[c*R+R-1],0,&Y[c*R+R-n],1); //end edge samps
                }
            }
        }
        else
        {
            for (n=1; n<=N; n++)
            {
                cblas_saxpy(C*(R-n),-sc*n,&X[0],1,&Y[n*C],1); //past samps
                cblas_saxpy(C*(R-n),sc*n,&X[n*C],1,&Y[0],1);  //future samps
                for (c=0; c<C; c++)
                {
                    cblas_saxpy(n,-sc*n,&X[c],0,&Y[c],C);                //beg edge samps
                    cblas_saxpy(n,sc*n,&X[c+C*(R-1)],0,&Y[c+C*(R-n)],C); //end edge samps
                }
            }
        }
    }
    else if (dim==1)
    {
        if (iscolmajor)
        {
            for (n=1; n<=N; n++)
            {
                cblas_saxpy(R*(C-n),-sc*n,&X[0],1,&Y[n*R],1); //past samps
                cblas_saxpy(R*(C-n),sc*n,&X[n*R],1,&Y[0],1);  //future samps
                for (r=0; r<R; r++)
                {
                    cblas_saxpy(n,-sc*n,&X[r],0,&Y[r],R);                //beg edge samps
                    cblas_saxpy(n,sc*n,&X[r+R*(C-1)],0,&Y[r+R*(C-n)],R); //end edge samps
                }
            }
        }
        else
        {
            for (r=0; r<R; r++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[r*C],0,&Y[r*C],1);        //beg edge samps
                    cblas_saxpy(C-n,-sc*n,&X[r*C],1,&Y[r*C+n],1);    //past samps
                    cblas_saxpy(C-n,sc*n,&X[r*C+n],1,&Y[r*C],1);     //future samps
                    cblas_saxpy(n,sc*n,&X[r*C+C-1],0,&Y[r*C+C-n],1); //end edge samps
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in get_deltas_s: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int get_deltas_d (double *Y, const double *X, const int iscolmajor, const int R, const int C, const int dim, const int N)
{
    const double z = 0.0;
    int r, c, n;
    double sc = 1.0;

    //Checks
    if (R<1) { fprintf(stderr,"error in get_deltas_d: R (nrows Y) must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in get_deltas_d: C (ncols Y) must be positive\n"); return 1; }
    if (N<1) { fprintf(stderr,"error in get_deltas_d: N (delta winlength) must be positive\n"); return 1; }

    //Get sc (normalizer)
    for (n=2; n<=N; n++) { sc += n*n; }
    sc = 0.5/sc;

    //Initialize Y
    cblas_dcopy(R*C,&z,0,&Y[0],1);

    if (dim==0)
    {
        if (iscolmajor)
        {
            for (c=0; c<C; c++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[c*R],0,&Y[c*R],1);        //beg edge samps
                    cblas_daxpy(R-n,-sc*n,&X[c*R],1,&Y[c*R+n],1);    //past samps
                    cblas_daxpy(R-n,sc*n,&X[c*R+n],1,&Y[c*R],1);     //future samps
                    cblas_daxpy(n,sc*n,&X[c*R+R-1],0,&Y[c*R+R-n],1); //end edge samps
                }
            }
        }
        else
        {
            for (n=1; n<=N; n++)
            {
                cblas_daxpy(C*(R-n),-sc*n,&X[0],1,&Y[n*C],1); //past samps
                cblas_daxpy(C*(R-n),sc*n,&X[n*C],1,&Y[0],1);  //future samps
                for (c=0; c<C; c++)
                {
                    cblas_daxpy(n,-sc*n,&X[c],0,&Y[c],C);                //beg edge samps
                    cblas_daxpy(n,sc*n,&X[c+C*(R-1)],0,&Y[c+C*(R-n)],C); //end edge samps
                }
            }
        }
    }
    else if (dim==1)
    {
        if (iscolmajor)
        {
            for (n=1; n<=N; n++)
            {
                cblas_daxpy(R*(C-n),-sc*n,&X[0],1,&Y[n*R],1); //past samps
                cblas_daxpy(R*(C-n),sc*n,&X[n*R],1,&Y[0],1);  //future samps
                for (r=0; r<R; r++)
                {
                    cblas_daxpy(n,-sc*n,&X[r],0,&Y[r],R);                //beg edge samps
                    cblas_daxpy(n,sc*n,&X[r+R*(C-1)],0,&Y[r+R*(C-n)],R); //end edge samps
                }
            }
        }
        else
        {
            for (r=0; r<R; r++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[r*C],0,&Y[r*C],1);        //beg edge samps
                    cblas_daxpy(C-n,-sc*n,&X[r*C],1,&Y[r*C+n],1);    //past samps
                    cblas_daxpy(C-n,sc*n,&X[r*C+n],1,&Y[r*C],1);     //future samps
                    cblas_daxpy(n,sc*n,&X[r*C+C-1],0,&Y[r*C+C-n],1); //end edge samps
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in get_deltas_d: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif


text:
#ifndef __GSL_WRAPPER_H__
#define __GSL_WRAPPER_H__

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

#include <iostream>
#include <cassert>

namespace dvrlib{

struct gsl_exception {
  const char* reason;
  const char* file;
  int line;
  int gsl_errno;
};

void gsl_enable_exceptions();

class vector_view;

class vector {
  gsl_vector* v;
  vector();
public:
  vector(int n);
  vector(int n, double x);
  vector(int n, const double* x);
  vector(const vector& src);
  ~vector();

  gsl_vector* gsl_internal();
  const gsl_vector* gsl_internal() const;

  int size() const;
  void set(int i, double val);
  double get(int i) const;
  double operator[](int i);

  vector& operator=(const vector& src);
  vector& operator+=(const vector& src);
  vector operator+(const vector& src) const;
  vector& operator-=(const vector& src);
  vector operator-(const vector& src) const;
  vector operator-() const;

  vector& operator*=(double d);
  vector operator*(double d) const;
  double norm1() const;
  double norm2() const;

  vector_view subvector(int k, int n);
  const vector_view subvector(int k, int n) const;
  friend class vector_view;
  friend class matrix;
};

vector operator*(double d, const vector& src);
std::ostream& operator<<(std::ostream& out, const vector& vec);

class vector_view: public vector {
  gsl_vector_view vv;
  vector_view(gsl_vector_view vv);

public:
  vector_view(const vector_view& src);
  gsl_vector_view* gsl_internal();
  vector_view& operator=(const vector& src);

  friend class vector;
  friend class matrix;
};

class matrix_view;

class matrix {
  gsl_matrix* m;
  matrix();
public:
  matrix(int n1, int n2, bool id=false, const double* diag=0);
  matrix(int n1, int n2, const double* x);
  template<int n>
    matrix(int n1, int n2, const double (*x)[n]);

  matrix(const matrix& src);
  matrix(const gsl_matrix* src);
  ~matrix();

  gsl_matrix* gsl_internal();
  const gsl_matrix* gsl_internal() const;

  int size1() const;
  int size2() const;
  void set(int i, int j, double val);
  double get(int i, int j) const;
  vector_view operator[](int i);

  matrix& operator=(const matrix& src);

  matrix operator+(const matrix& src) const;
  matrix operator+=(const matrix& src) const;
  matrix operator-(const matrix& src) const;
  matrix operator-=(const matrix& src) const;
  matrix operator-() const;

  vector operator*(const vector& src) const;
  matrix operator*(const matrix& src) const;
  matrix operator*(double d) const;
  matrix operator*=(double d) const;

  matrix transpose() const;
  matrix inverse() const;
  vector linsolve(const vector& b) const;
  void svd(matrix& U, matrix& V, vector& s) const;
  vector svd() const;

  matrix_view submatrix(int k1, int k2, int n1, int n2);
  const matrix_view submatrix(int k1, int k2, int n1, int n2) const;
  friend class matrix_view;
};

template<int n>
matrix::matrix(int n1, int n2, const double (*x)[n]) {
  assert(n==n2);
  m = gsl_matrix_alloc(n1, n2);
  gsl_matrix_const_view  src = gsl_matrix_const_view_array(x[0], n1, n2);
  gsl_matrix_memcpy(m, &src.matrix);
}

matrix operator*(double d, const matrix& src);

std::ostream& operator<<(std::ostream& out, const matrix& vec);

class matrix_view: public matrix {
  gsl_matrix_view mv;
  matrix_view(gsl_matrix_view mv);

public:
  matrix_view(const matrix_view& src);
  gsl_matrix_view* gsl_internal();
  matrix_view& operator=(const matrix& src);

  friend class matrix;
};

std::ostream& operator<<(std::ostream& out, const gsl_vector& v);
std::ostream& operator<<(std::ostream& out, const gsl_matrix& m);
std::ostream& operator<<(std::ostream& out, const gsl_matrix_view& mv);

} // namespace dvrlib

#endif // __GSL_WRAPPER_H__

max_stars_repo_name: ezander/dvrlib

max_issues_repo_name: ezander/dvrlib

max_forks_repo_name: ezander/dvrlib


text:
#ifndef S3D_CV_IMAGE_OPERATION_IMAGE_OPERATION_H
#define S3D_CV_IMAGE_OPERATION_IMAGE_OPERATION_H

#include <gsl/gsl>

namespace cv {
  class Mat;
}

namespace s3d {
  struct StanResults;
}

namespace s3d {
namespace image_operation {

class ImageOperation {
public:
  bool applyOnImagesIfEnabled(cv::Mat *leftImage, cv::Mat *rightImage, StanResults* results);

  bool isEnabled();
  void enable();
  void disable();

private:
  virtual bool applyOnImage(cv::Mat *leftImage, cv::Mat *rightImage, StanResults* results) = 0;

  bool isEnabled_{true};
};

} // namespace s3d
} // namespace image_operation

#endif // S3D_CV_IMAGE_OPERATION_IMAGE_OPERATION_H


text:
//
//  symm.h
//  Linear Algebra Template Library
//
//  Created by Rodney James on 1/2/12.
//  Copyright (c) 2012 University of Colorado Denver. All rights reserved.
//

#ifndef _symm_h
#define _symm_h

/// @file symm.h Performs general symmetric matrix-matrix multiplication.

#include <cctype>
#include "latl.h"

namespace LATL
{
   /// @brief  Performs general real symmetric matrix-matrix multiplication.
   ///
   /// For real matrices B and C, symmetric matrix A, and real scalars alpha and beta,
   ///
   ///         C := alpha*A*B + beta*C  or   C := alpha*B*A + beta*C
   /// is computed.
   /// @return 0 if success.
   /// @return -i if the ith argument is invalid.
   /// @tparam real_t Floating point type.
   /// @param side Specifies whether the matrix A appears on the left or right side as follows:
   ///
   ///        if side = 'L' or 'l' then C := alpha*A*B + beta*C,
   ///        if side = 'R' or 'r' then C := alpha*B*A + beta*C.
   /// @param uplo Specifies whether the upper or lower triangular part of the symmetric matrix A
   /// is to be referenced:
   ///
   ///        if uplo = 'U' or 'u' then A is upper triangular,
   ///        if uplo = 'L' or 'l' then A is lower triangular.
   /// @param m Specifies the number of rows of the matrices B and C.  m>=0
   /// @param n Specifies the number of columns of the matrices B and C.  n>=0
   /// @param alpha Real scalar.
   /// @param A Pointer to real symmetric matrix A.  If side = 'L' or 'l', then A is m-by-m;
   /// if side = 'R' or 'r', then A is n-by-n.
   /// @param ldA Column length of the matrix A.  If side = 'L' or 'l', ldA>=m.  If side = 'R' or 'r', ldA>=n.
   /// @param B Pointer to real m-by-n matrix B.
   /// @param ldB Column length of the matrix B.  ldB>=m.
   /// @param beta Real scalar.
   /// @param C Pointer to real m-by-n matrix C.
   /// @param ldC Column length of the matrix C.  ldC>=m.
   /// @ingroup BLAS

   template <typename real_t>
   int SYMM(char side, char uplo, int_t m, int_t n, real_t alpha, real_t *A, int_t ldA, real_t *B, int_t ldB, real_t beta, real_t *C, int_t ldC)
   {
      using std::toupper;

      const real_t zero(0.0);
      const real_t one(1.0);
      int_t i,j,k;
      real_t *a,*b,*c,*at,*bt;
      real_t s,t;

      side=toupper(side);
      uplo=toupper(uplo);

      if((side!='L')&&(side!='R'))
         return -1;
      else if((uplo!='U')&&(uplo!='L'))
         return -2;
      else if(m<0)
         return -3;
      else if(n<0)
         return -4;
      else if(ldA<((side=='L')?m:n))
         return -7;
      else if(ldB<m)
         return -9;
      else if(ldC<m)
         return -12;
      else if((m==0)||(n==0)||((alpha==zero)&&(beta==one)))
         return 0;

      if(alpha==zero)
      {
         c=C;
         for(j=0;j<n;j++)
         {
            for(i=0;i<m;i++)
               c[i]*=beta;
            c+=ldC;
         }
      }
      else if(side=='L')
      {
         if(uplo=='U')
         {
            b=B;
            c=C;
            for(j=0;j<n;j++)
            {
               a=A;
               for(i=0;i<m;i++)
               {
                  t=alpha*b[i];
                  s=zero;
                  for(k=0;k<i;k++)
                  {
                     c[k]+=t*a[k];
                     s+=b[k]*a[k];
                  }
                  c[i]=beta*c[i]+t*a[i]+alpha*s;
                  a+=ldA;
               }
               b+=ldB;
               c+=ldC;
            }
         }
         else
         {
            b=B;
            c=C;
            for(j=0;j<n;j++)
            {
               a=A+m*ldA;
               for(i=m-1;i>=0;i--)
               {
                  a-=ldA;
                  t=alpha*b[i];
                  s=zero;
                  for(k=i+1;k<m;k++)
                  {
                     c[k]+=t*a[k];
                     s+=b[k]*a[k];
                  }
                  c[i]=beta*c[i]+t*a[i]+alpha*s;
               }
               b+=ldB;
               c+=ldC;
            }
         }
      }
      else
      {
         if(uplo=='U')
         {
            a=A;
            c=C;
            b=B;
            for(j=0;j<n;j++)
            {
               t=alpha*a[j];
               for(i=0;i<m;i++)
                  c[i]=c[i]*beta+t*b[i];
               bt=B;
               for(k=0;k<j;k++)
               {
                  t=alpha*a[k];
                  for(i=0;i<m;i++)
                     c[i]+=bt[i]*t;
                  bt+=ldB;
               }
               at=A+(j+1)*ldA;
               bt=B+(j+1)*ldB;
               for(k=j+1;k<n;k++)
               {
                  t=alpha*at[j];
                  for(i=0;i<m;i++)
                     c[i]+=t*bt[i];
                  at+=ldA;
                  bt+=ldB;
               }
               a+=ldA;
               b+=ldB;
               c+=ldC;
            }
         }
         else
         {
            a=A;
            c=C;
            b=B;
            for(j=0;j<n;j++)
            {
               t=alpha*a[j];
               for(i=0;i<m;i++)
                  c[i]=c[i]*beta+t*b[i];
               bt=B;
               at=A;
               for(k=0;k<j;k++)
               {
                  t=alpha*at[j];
                  for(i=0;i<m;i++)
                     c[i]+=bt[i]*t;
                  at+=ldA;
                  bt+=ldB;
               }
               bt=B+(j+1)*ldB;
               for(k=j+1;k<n;k++)
               {
                  t=alpha*at[k];
                  for(i=0;i<m;i++)
                     c[i]+=t*bt[i];
                  bt+=ldB;
               }
               a+=ldA;
               b+=ldB;
               c+=ldC;
            }
         }
      }
      return 0;
   }

   /// @brief  Performs general complex symmetric matrix-matrix multiplication.
   ///
   /// For complex matrices B and C, symmetric matrix A, and complex scalars alpha and beta,
   ///
   ///         C := alpha*A*B + beta*C  or   C := alpha*B*A + beta*C
   /// is computed.
   /// @return 0 if success.
   /// @return -i if the ith argument is invalid.
   /// @tparam real_t Floating point type.
   /// @param side Specifies whether the matrix A appears on the left or right side as follows:
   ///
   ///        if side = 'L' or 'l' then C := alpha*A*B + beta*C,
   ///        if side = 'R' or 'r' then C := alpha*B*A + beta*C.
   /// @param uplo Specifies whether the upper or lower triangular part of the symmetric matrix A
   /// is to be referenced:
   ///
   ///        if uplo = 'U' or 'u' then A is upper triangular,
   ///        if uplo = 'L' or 'l' then A is lower triangular.
   /// @param m Specifies the number of rows of the matrices B and C.  m>=0
   /// @param n Specifies the number of columns of the matrices B and C.  n>=0
   /// @param alpha Complex scalar.
   /// @param A Pointer to complex symmetric matrix A.  If side = 'L' or 'l', then A is m-by-m;
   /// if side = 'R' or 'r', then A is n-by-n.
   /// @param ldA Column length of the matrix A.  If side = 'L' or 'l', ldA>=m.  If side = 'R' or 'r', ldA>=n.
   /// @param B Pointer to complex m-by-n matrix B.
   /// @param ldB Column length of the matrix B.  ldB>=m.
   /// @param beta Complex scalar.
   /// @param C Pointer to complex m-by-n matrix C.
   /// @param ldC Column length of the matrix C.  ldC>=m.
   /// @ingroup BLAS

   template <typename real_t>
   int SYMM(char side, char uplo, int_t m, int_t n, complex<real_t> alpha, complex<real_t> *A, int_t ldA, complex<real_t> *B, int_t ldB, complex<real_t> beta, complex<real_t> *C, int_t ldC)
   {
      return SYMM< complex<real_t> >(side,uplo,m,n,alpha,A,ldA,B,ldB,beta,C,ldC);
   }

#ifdef __latl_cblas
#include <cblas.h>
   template <> int SYMM<float>(char side, char uplo, int_t m, int_t n, float alpha, float *A, int_t ldA, float *B, int_t ldB, float beta, float *C, int_t ldC)
   {
      using std::toupper;
      side=toupper(side);
      uplo=toupper(uplo);
      if((side!='L')&&(side!='R'))
         return -1;
      else if((uplo!='U')&&(uplo!='L'))
         return -2;
      else if(m<0)
         return -3;
      else if(n<0)
         return -4;
      else if(ldA<((side=='L')?m:n))
         return -7;
      else if(ldB<m)
         return -9;
      else if(ldC<m)
         return -12;

      const CBLAS_SIDE Side=(side=='L')?CblasLeft:CblasRight;
      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;

      cblas_ssymm(CblasColMajor,Side,Uplo,m,n,alpha,A,ldA,B,ldB,beta,C,ldC);

      return 0;

   }

   template <> int SYMM<double>(char side, char uplo, int_t m, int_t n, double alpha, double *A, int_t ldA, double *B, int_t ldB, double beta, double *C, int_t ldC)
   {
      using std::toupper;
      side=toupper(side);
      uplo=toupper(uplo);
      if((side!='L')&&(side!='R'))
         return -1;
      else if((uplo!='U')&&(uplo!='L'))
         return -2;
      else if(m<0)
         return -3;
      else if(n<0)
         return -4;
      else if(ldA<((side=='L')?m:n))
         return -7;
      else if(ldB<m)
         return -9;
      else if(ldC<m)
         return -12;

      const CBLAS_SIDE Side=(side=='L')?CblasLeft:CblasRight;
      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;

      cblas_dsymm(CblasColMajor,Side,Uplo,m,n,alpha,A,ldA,B,ldB,beta,C,ldC);

      return 0;

   }
   template <> int SYMM<float>(char side, char uplo, int_t m, int_t n, complex<float> alpha, complex<float> *A, int_t ldA, complex<float> *B, int_t ldB, complex<float> beta, complex<float> *C, int_t ldC)
   {
      using std::toupper;
      side=toupper(side);
      uplo=toupper(uplo);
      if((side!='L')&&(side!='R'))
         return -1;
      else if((uplo!='U')&&(uplo!='L'))
         return -2;
      else if(m<0)
         return -3;
      else if(n<0)
         return -4;
      else if(ldA<((side=='L')?m:n))
         return -7;
      else if(ldB<m)
         return -9;
      else if(ldC<m)
         return -12;

      const CBLAS_SIDE Side=(side=='L')?CblasLeft:CblasRight;
      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;

      cblas_csymm(CblasColMajor,Side,Uplo,m,n,&alpha,A,ldA,B,ldB,&beta,C,ldC);

      return 0;

   }
   template <> int SYMM<double>(char side, char uplo, int_t m, int_t n, complex<double> alpha, complex<double> *A, int_t ldA, complex<double> *B, int_t ldB, complex<double> beta, complex<double> *C, int_t ldC)
   {
      using std::toupper;
      side=toupper(side);
      uplo=toupper(uplo);
      if((side!='L')&&(side!='R'))
         return -1;
      else if((uplo!='U')&&(uplo!='L'))
         return -2;
      else if(m<0)
         return -3;
      else if(n<0)
         return -4;
      else if(ldA<((side=='L')?m:n))
         return -7;
      else if(ldB<m)
         return -9;
      else if(ldC<m)
         return -12;

      const CBLAS_SIDE Side=(side=='L')?CblasLeft:CblasRight;
      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;

      cblas_zsymm(CblasColMajor,Side,Uplo,m,n,&alpha,A,ldA,B,ldB,&beta,C,ldC);

      return 0;

   }
#endif
}
#endif

max_stars_repo_name: langou/latl

max_issues_repo_name: langou/latl

max_forks_repo_name: langou/latl


text:
#include <petsc.h>
#include "um.h"

PetscErrorCode UMInitialize(UM *mesh) {
    mesh->N = 0;
    mesh->K = 0;
    mesh->P = 0;
    mesh->loc = NULL;
    mesh->e = NULL;
    mesh->bf = NULL;
    mesh->ns = NULL;
    return 0;
}

PetscErrorCode UMDestroy(UM *mesh) {
    PetscErrorCode ierr;
    ierr = VecDestroy(&(mesh->loc)); CHKERRQ(ierr);
    ierr = ISDestroy(&(mesh->e)); CHKERRQ(ierr);
    ierr = ISDestroy(&(mesh->bf)); CHKERRQ(ierr);
    ierr = ISDestroy(&(mesh->ns)); CHKERRQ(ierr);
    return 0;
}

PetscErrorCode UMViewASCII(UM *mesh, PetscViewer viewer) {
    PetscErrorCode ierr;
    PetscInt        n, k;
    const Node      *aloc;
    const PetscInt  *ae, *abf, *ans;

    ierr = PetscViewerASCIIPushSynchronized(viewer); CHKERRQ(ierr);
    if (mesh->loc && (mesh->N > 0)) {
        ierr = PetscViewerASCIISynchronizedPrintf(viewer,"%d nodes at (x,y) coordinates:\n",mesh->N); CHKERRQ(ierr);
        ierr = VecGetArrayRead(mesh->loc,(const PetscReal **)&aloc); CHKERRQ(ierr);
        for (n = 0; n < mesh->N; n++) {
            ierr = PetscViewerASCIISynchronizedPrintf(viewer,"    %3d : (%g,%g)\n",
                               n,aloc[n].x,aloc[n].y); CHKERRQ(ierr);
        }
        ierr = VecRestoreArrayRead(mesh->loc,(const PetscReal **)&aloc); CHKERRQ(ierr);
    } else {
        ierr = PetscViewerASCIISynchronizedPrintf(viewer,"node coordinates empty or unallocated\n"); CHKERRQ(ierr);
    }
    if (mesh->e && (mesh->K > 0)) {
        ierr = PetscViewerASCIISynchronizedPrintf(viewer,"%d elements:\n",mesh->K); CHKERRQ(ierr);
        ierr = ISGetIndices(mesh->e,&ae); CHKERRQ(ierr);
        for (k = 0; k < mesh->K; k++) {
            ierr = PetscPrintf(PETSC_COMM_WORLD,"    %3d : %3d %3d %3d\n",
                               k,ae[3*k+0],ae[3*k+1],ae[3*k+2]); CHKERRQ(ierr);
        }
        ierr = ISRestoreIndices(mesh->e,&ae); CHKERRQ(ierr);
    } else {
        ierr = PetscViewerASCIISynchronizedPrintf(viewer,"element index triples empty or unallocated\n"); CHKERRQ(ierr);
    }
    if (mesh->bf && (mesh->N > 0)) {
        ierr = PetscViewerASCIISynchronizedPrintf(viewer,"%d boundary flags at nodes (0 = interior, 1 = boundary, 2 = Dirichlet):\n",mesh->N); CHKERRQ(ierr);
        ierr = ISGetIndices(mesh->bf,&abf); CHKERRQ(ierr);
        for (n = 0; n < mesh->N; n++) {
            ierr = PetscViewerASCIISynchronizedPrintf(viewer,"    %3d : %1d\n",
                               n,abf[n]); CHKERRQ(ierr);
        }
        ierr = ISRestoreIndices(mesh->bf,&abf); CHKERRQ(ierr);
    } else {
        ierr = PetscViewerASCIISynchronizedPrintf(viewer,"boundary flags empty or unallocated\n"); CHKERRQ(ierr);
    }
    if (mesh->ns && (mesh->P > 0)) {
        ierr = PetscViewerASCIISynchronizedPrintf(viewer,"%d Neumann boundary segments:\n",mesh->P); CHKERRQ(ierr);
        ierr = ISGetIndices(mesh->ns,&ans); CHKERRQ(ierr);
        for (n = 0; n < mesh->P; n++) {
            ierr = PetscViewerASCIISynchronizedPrintf(viewer,"    %3d : %3d %3d\n",
                               n,ans[2*n+0],ans[2*n+1]); CHKERRQ(ierr);
        }
        ierr = ISRestoreIndices(mesh->ns,&ans); CHKERRQ(ierr);
    } else {
        ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Neumann boundary segments empty or unallocated\n"); CHKERRQ(ierr);
    }
    ierr = PetscViewerASCIIPopSynchronized(viewer); CHKERRQ(ierr);
    return 0;
}

PetscErrorCode UMViewSolutionBinary(UM *mesh, char *filename, Vec u) {
    PetscErrorCode ierr;
    PetscInt       Nu;
    PetscViewer viewer;
    ierr = VecGetSize(u,&Nu); CHKERRQ(ierr);
    if (Nu != mesh->N) {
        SETERRQ2(PETSC_COMM_SELF,1,
           "incompatible sizes of u (=%d) and number of nodes (=%d)\n",Nu,mesh->N);
    }
    ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer); CHKERRQ(ierr);
    ierr = VecView(u,viewer); CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
    return 0;
}

PetscErrorCode UMReadNodes(UM *mesh, char *filename) {
    PetscErrorCode ierr;
    PetscInt       twoN;
    PetscViewer viewer;
    if (mesh->N > 0) {
        SETERRQ(PETSC_COMM_SELF,1,"nodes already created?\n");
    }
    ierr = VecCreate(PETSC_COMM_WORLD,&mesh->loc); CHKERRQ(ierr);
    ierr = VecSetFromOptions(mesh->loc); CHKERRQ(ierr);
    ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer); CHKERRQ(ierr);
    ierr = VecLoad(mesh->loc,viewer); CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
    ierr = VecGetSize(mesh->loc,&twoN); CHKERRQ(ierr);
    if (twoN % 2 != 0) {
        SETERRQ1(PETSC_COMM_SELF,2,"node locations loaded from %s are not N pairs\n",filename);
    }
    mesh->N = twoN / 2;
    return 0;
}

PetscErrorCode UMCheckElements(UM *mesh) {
    PetscErrorCode ierr;
    const PetscInt  *ae;
    PetscInt        k, m;
    if ((mesh->K == 0) || (mesh->e == NULL)) {
        SETERRQ(PETSC_COMM_SELF,1,
                "number of elements unknown; call UMReadElements() first\n");
    }
    if (mesh->N == 0) {
        SETERRQ(PETSC_COMM_SELF,2,
                "node size unknown so element check impossible; call UMReadNodes() first\n");
    }
    ierr = ISGetIndices(mesh->e,&ae); CHKERRQ(ierr);
    for (k = 0; k < mesh->K; k++) {
        for (m = 0; m < 3; m++) {
            if ((ae[3*k+m] < 0) || (ae[3*k+m] >= mesh->N)) {
                SETERRQ3(PETSC_COMM_SELF,3,
                   "index e[%d]=%d invalid: not between 0 and N-1=%d\n",
                   3*k+m,ae[3*k+m],mesh->N-1);
            }
        }
        // FIXME: could add check for distinct indices
    }
    ierr = ISRestoreIndices(mesh->e,&ae); CHKERRQ(ierr);
    return 0;
}

PetscErrorCode UMCheckBoundaryData(UM *mesh) {
    PetscErrorCode ierr;
    const PetscInt  *ans, *abf;
    PetscInt        n, m;
    if (mesh->N == 0) {
        SETERRQ(PETSC_COMM_SELF,2,
                "node size unknown so boundary flag check impossible; call UMReadNodes() first\n");
    }
    if (mesh->bf == NULL) {
        SETERRQ(PETSC_COMM_SELF,1,
                "boundary flags at nodes not allocated; call UMReadNodes() first\n");
    }
    if ((mesh->P > 0) && (mesh->ns == NULL)) {
        SETERRQ(PETSC_COMM_SELF,3,
                "inconsistent data for Neumann boundary segments\n");
    }
    ierr = ISGetIndices(mesh->bf,&abf); CHKERRQ(ierr);
    for (n = 0; n < mesh->N; n++) {
        switch (abf[n]) {
            case 0 :
            case 1 :
            case 2 :
                break;
            default :
                SETERRQ2(PETSC_COMM_SELF,5,
                   "boundary flag bf[%d]=%d invalid: not in {0,1,2}\n",
                   n,abf[n]);
        }
    }
    ierr = ISRestoreIndices(mesh->bf,&abf); CHKERRQ(ierr);
    if (mesh->P > 0) {
        ierr = ISGetIndices(mesh->ns,&ans); CHKERRQ(ierr);
        for (n = 0; n < mesh->P; n++) {
            for (m = 0; m < 2; m++) {
                if ((ans[2*n+m] < 0) || (ans[2*n+m] >= mesh->N)) {
                    SETERRQ3(PETSC_COMM_SELF,6,
                       "index ns[%d]=%d invalid: not between 0 and N-1=%d\n",
                       2*n+m,ans[3*n+m],mesh->N-1);
                }
            }
        }
        ierr = ISRestoreIndices(mesh->ns,&ans); CHKERRQ(ierr);
    }
    return 0;
}

PetscErrorCode UMReadISs(UM *mesh, char *filename) {
    PetscErrorCode ierr;
    PetscViewer  viewer;
    PetscInt     n_bf;
    if ((!mesh->loc) || (mesh->N == 0)) {
        SETERRQ(PETSC_COMM_SELF,2,
                "node coordinates not created ... do that first ... stopping\n");
    }
    if ((mesh->K > 0) || (mesh->P > 0) || (mesh->e != NULL) || (mesh->bf != NULL) || (mesh->ns != NULL)) {
        SETERRQ(PETSC_COMM_SELF,1,
                "elements, boundary flags, Neumann boundary segments already created? ... stopping\n");
    }
    ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer); CHKERRQ(ierr);
    // create and load e
    ierr = ISCreate(PETSC_COMM_WORLD,&(mesh->e)); CHKERRQ(ierr);
    ierr = ISLoad(mesh->e,viewer); CHKERRQ(ierr);
    ierr = ISGetSize(mesh->e,&(mesh->K)); CHKERRQ(ierr);
    if (mesh->K % 3 != 0) {
        SETERRQ1(PETSC_COMM_SELF,3,
                 "IS e loaded from %s is wrong size for list of element triples\n",filename);
    }
    mesh->K /= 3;
    // create and load bf
    ierr = ISCreate(PETSC_COMM_WORLD,&(mesh->bf)); CHKERRQ(ierr);
    ierr = ISLoad(mesh->bf,viewer); CHKERRQ(ierr);
    ierr = ISGetSize(mesh->bf,&n_bf); CHKERRQ(ierr);
    if (n_bf != mesh->N) {
        SETERRQ1(PETSC_COMM_SELF,4,
                 "IS bf loaded from %s is wrong size for list of boundary flags\n",filename);
    }
    // FIXME  seems there is no way to tell if file is empty at this point
    // create and load ns last ... may *start with a negative value* in which case set P = 0
    const PetscInt *ans;
    ierr = ISCreate(PETSC_COMM_WORLD,&(mesh->ns)); CHKERRQ(ierr);
    ierr = ISLoad(mesh->ns,viewer); CHKERRQ(ierr);
    ierr = ISGetIndices(mesh->ns,&ans); CHKERRQ(ierr);
    if (ans[0] < 0) {
        ISDestroy(&(mesh->ns));
        mesh->ns = NULL;
        mesh->P = 0;
    } else {
        ierr = ISGetSize(mesh->ns,&(mesh->P)); CHKERRQ(ierr);
        if (mesh->P % 2 != 0) {
            SETERRQ1(PETSC_COMM_SELF,4,
                     "IS s loaded from %s is wrong size for list of Neumann boundary segment pairs\n",filename);
        }
        mesh->P /= 2;
    }
    ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);

    // check that mesh is complete now
    ierr = UMCheckElements(mesh); CHKERRQ(ierr);
    ierr = UMCheckBoundaryData(mesh); CHKERRQ(ierr);
    return 0;
}

PetscErrorCode UMStats(UM *mesh, PetscReal *maxh, PetscReal *meanh,
                       PetscReal *maxa, PetscReal *meana) {
    PetscErrorCode ierr;
    const PetscInt *ae;
    const Node     *aloc;
    PetscInt       k;
    PetscReal      x[3], y[3], ax, ay, bx, by, cx, cy, h, a,
                   Maxh = 0.0, Maxa = 0.0, Sumh = 0.0, Suma = 0.0;
    if ((mesh->K == 0) || (mesh->e == NULL)) {
        SETERRQ(PETSC_COMM_SELF,1,
                "number of elements unknown; call UMReadElements() first\n");
    }
    if (mesh->N == 0) {
        SETERRQ(PETSC_COMM_SELF,2,
                "node size unknown so element check impossible; call UMReadNodes() first\n");
    }
    ierr = UMGetNodeCoordArrayRead(mesh,&aloc); CHKERRQ(ierr);
    ierr = ISGetIndices(mesh->e,&ae); CHKERRQ(ierr);
    for (k = 0; k < mesh->K; k++) {
        x[0] = aloc[ae[3*k]].x;
        y[0] = aloc[ae[3*k]].y;
        x[1] = aloc[ae[3*k+1]].x;
        y[1] = aloc[ae[3*k+1]].y;
        x[2] = aloc[ae[3*k+2]].x;
        y[2] = aloc[ae[3*k+2]].y;
        ax = x[1] - x[0];
        ay = y[1] - y[0];
        bx = x[2] - x[0];
        by = y[2] - y[0];
        cx = x[1] - x[2];
        cy = y[1] - y[2];
        h = PetscMax(ax*ax+ay*ay, PetscMax(bx*bx+by*by, cx*cx+cy*cy));
        h = sqrt(h);
        a = 0.5 * PetscAbs(ax*by-ay*bx);
        Maxh = PetscMax(Maxh,h);
        Sumh += h;
        Maxa = PetscMax(Maxa,a);
        Suma += a;
    }
    ierr = ISRestoreIndices(mesh->e,&ae); CHKERRQ(ierr);
    ierr = UMRestoreNodeCoordArrayRead(mesh,&aloc); CHKERRQ(ierr);
    if (maxh)  *maxh = Maxh;
    if (maxa)  *maxa = Maxa;
    if (meanh)  *meanh = Sumh / mesh->K;
    if (meana)  *meana = Suma / mesh->K;
    return 0;
}

PetscErrorCode UMGetNodeCoordArrayRead(UM *mesh, const Node **xy) {
    PetscErrorCode ierr;
    if ((!mesh->loc) || (mesh->N == 0)) {
        SETERRQ(PETSC_COMM_SELF,1,"node coordinates not created ... stopping\n");
    }
    ierr = VecGetArrayRead(mesh->loc,(const PetscReal **)xy); CHKERRQ(ierr);
    return 0;
}

PetscErrorCode UMRestoreNodeCoordArrayRead(UM *mesh, const Node **xy) {
    PetscErrorCode ierr;
    if ((!mesh->loc) || (mesh->N == 0)) {
        SETERRQ(PETSC_COMM_SELF,1,"node coordinates not created ... stopping\n");
    }
    ierr = VecRestoreArrayRead(mesh->loc,(const PetscReal **)xy); CHKERRQ(ierr);
    return 0;
}


text:
#include <stdio.h>
#include <malloc.h>
#include <stdlib.h>
#include <stdbool.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <functions.h>
#define FOR(i,n) for(i=0;i<n;i++)
typedef double DD, LD;
typedef long LL;

void quicksort_mds(double* x, long first, long last)
{
	long pivot, j, i;
	double temp;
	if (first<last)
	{
		pivot = first;
		i = first;
		j = last;
		while (i<j)
		{
			while (x[i]<=x[pivot] && i<last)
				i++;
			while (x[j]>x[pivot])
				j--;
			if (i<j)
			{
				temp = x[i];
				x[i] = x[j];
				x[j] = temp;
			}
		}
		temp = x[pivot];
		x[pivot] = x[j];
		x[j] = temp;
		quicksort(x, first, j - 1);
		quicksort(x, j + 1, last);

	}
}

void find_eigen(double **mat, int n, double*** res_vec, double** res_val)
{
	int i, j;
	double *data = (double*)calloc(((long)(n*n)), sizeof(double));
	double *eig_values = (double*)calloc(n, sizeof(double));
	double** eig_vec = (double**)calloc(n, sizeof(double*));
	for (i = 0; i<n; i++)eig_vec[i] = (double*)calloc(n, sizeof(double));
	for (i = 0; i<n; i++)
	{
		for (j = 0; j<n; j++)
		{
			data[i*n + j] = mat[i][j];
		}
	}
	gsl_matrix_view m= gsl_matrix_view_array(data, n, n);
	gsl_vector *eval = gsl_vector_alloc(n);
	gsl_matrix *evec = gsl_matrix_alloc(n, n);
	gsl_eigen_symmv_workspace * w =gsl_eigen_symmv_alloc(n);
	gsl_eigen_symmv(&m.matrix, eval, evec, w);
	gsl_eigen_symmv_free(w);
	gsl_eigen_symmv_sort(eval, evec,GSL_EIGEN_SORT_VAL_ASC);
	for (i = 0; i<n; i++)
	{
		eig_values[i] = gsl_vector_get(eval, i);
		for (j = 0; j<n; j++)
		{
			eig_vec[i][j] = gsl_matrix_get(evec, i, j);
		}
	}
	*res_val = eig_values;
	*res_vec = eig_vec;
	gsl_vector_free(eval);
	gsl_matrix_free(evec);
	free(data);
}

DD** Multiply_Matrices_mds(DD **Mat1, DD **Mat2, LL M, LL Q, LL P)
{
	LL i, c, d, k;
	DD Sum = 0;
	DD** Temp = (DD**)calloc(M, sizeof(DD*));
	FOR(i, M)Temp[i] = (DD*)calloc(Q, sizeof(DD));
	FOR(c, M)
	{
		FOR(d, Q)
		{
			FOR(k, P)
			{
				Sum += (Mat1[c][k] * Mat2[k][d]);
			}
			Temp[c][d] = Sum;
			Sum = 0;
		}
	}
	return Temp;
}

double ** H_mat(long n)
{
	long i, j;
	DD** k = (DD**)calloc(n, sizeof(DD*));
	FOR(i, n)
	{
		k[i] = (DD*)calloc(n, sizeof(DD));
	}
	FOR(i, n){
		FOR(j, n){ k[i][j] = (i == j) ? (1.00 - 1.00 / n) : (-1.00 / n);
		}
    }
	return k;
}
typedef struct{
	MAT Y;
	DD* eigen;
	LL number;
} Output_mds;

Output_mds MDS(DD** D,long rows,long columns, DD delta)
{
	//      ************** function to calculate Classical Multidimensional Scaling *************
	long  i, j, n;
	DD    **H, **new_D, **K, **new_k;
	new_D = (DD**)calloc(rows, sizeof(DD*));
	n = max(rows, columns);
	FOR(i, rows)
	{
		new_D[i] = (DD*)calloc(columns, sizeof(DD));
	}
	if (delta >= 0)
	{
		H = H_mat(n);
		new_k = (DD**)calloc(n, sizeof(DD*));
		FOR(i, rows)
		{
			new_k[i] = (DD*)calloc(n, sizeof(DD));
		}
		FOR(i, rows)
		{
			FOR(j, columns)
				new_D[i][j] = -0.5*D[i][j];
		}
		K = Multiply_Matrices_mds(Multiply_Matrices_mds(H, new_D, n, columns, n), H, columns, n, n);
        FOR(i,rows)free(new_D[i]);free(new_D);
        FOR(i,n)free(H[i]);free(H);
		DD* idx = (DD*)calloc(n, sizeof(DD));
		DD  ** temp, *val;
		FOR(i, n)
		{
			idx[i] = n - i;
			FOR(j, rows)
				new_k[i][j] = (K[i][j] + K[j][i]) / 2;
		}
		find_eigen(new_k, n, &temp, &val);
		find_eigen(new_k, n, &temp, &val);
		j = n - 1;
		for (i = 0; i < n / 2; i++)
		{
			DD temporary = val[i];
			val[i] = val[j];
			val[j] = temporary;
			j--;
		}
		LL* keep;
		DD max_val = val[0];
		LL magn_keep = 0;
		FOR(i, n){
			if (val[i] > max_val)max_val = val[i];
		}
		max_val = max_val*exp(-10);
		FOR(i, n){
			if (val[i] > max_val)magn_keep++;
		}
		keep = (LL*)calloc(magn_keep, sizeof(LL));
		magn_keep = 0;
		FOR(i, n){
			if (val[i] > max_val){
				keep[magn_keep++] = i;
			}
		}
		DD** V = (DD**)calloc(n, sizeof(DD*));
		DD** mult_V = (DD**)calloc(magn_keep, sizeof(DD*));
		FOR(i, n)V[i] = (DD*)calloc(magn_keep, sizeof(DD));
		FOR(i, magn_keep)
		{
			mult_V[i] = (DD*)calloc(magn_keep, sizeof(DD));
			FOR(j, n)
			{
				if (j < magn_keep)
				{
					if (i != j)mult_V[i][j] = 0;
					else mult_V[i][j] = sqrt(val[keep[i]]);
				}
				V[j][i] = temp[j][(int)idx[keep[i]]-1];
			}
		}
		V = Multiply_Matrices_mds(V, mult_V, n, magn_keep, magn_keep);
		MAT return_MAT;
		return_MAT.matrix = V;
		return_MAT.rows = n;
		return_MAT.columns = magn_keep;
		Output_mds t;
		t.Y = return_MAT; t.eigen = val; t.number=n;
		return t;
	}
	Output_mds t;
    return t;
}


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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_permutation.h>

static int
test_QR_band_decomp_eps(const size_t p, const size_t q, const gsl_matrix * A, const double eps, const char * desc)
{
  int s = 0;
  const size_t M = A->size1;
  const size_t N = A->size2;
  size_t i, j;

  gsl_matrix * AB = gsl_matrix_alloc(N, 2*p + q + 1);
  gsl_vector * tau = gsl_vector_alloc(N);
  gsl_matrix * Q = gsl_matrix_alloc(M, M);
  gsl_matrix * R = gsl_matrix_alloc(M, N);
  gsl_matrix * B = gsl_matrix_alloc(M, N);

  /* convert A to packed banded format */
  gen2band_matrix(p, q, A, AB);

#if 0
  print_octave(A, "A");
  printqrb_octave(M, p, q, AB, "AB");
#endif

  s += gsl_linalg_QR_band_decomp_L2(M, p, q, AB, tau);
  s += gsl_linalg_QR_band_unpack_L2(p, q, AB, tau, Q, R);

#if 0
  printqrb_octave(M, p, q, AB, "QR");
  printv_octave(tau, "tau");
  print_octave(Q, "Q");
  print_octave(R, "R");
#endif

  /* compute B = Q R */
  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, Q, R, 0.0, B);

  for (i = 0; i < M; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          double aij = gsl_matrix_get(A, i, j);
          double bij = gsl_matrix_get(B, i, j);

          gsl_test_rel(bij, aij, eps, "%s (M=%3lu,N=%3lu)(p=%3lu,q=%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, M, N, p, q, i, j, aij, bij);
        }
    }

  gsl_matrix_free(AB);
  gsl_vector_free(tau);
  gsl_matrix_free(Q);
  gsl_matrix_free(R);
  gsl_matrix_free(B);

  return s;
}

static int
test_QR_band_decomp(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 20;
  size_t M, N, p, q;

  /*M = 7; N = 6; p = 2; q = 1;*/
  for (M = 1; M <= N_max; ++M)
    {
      for (N = 1; N <= N_max; ++N)
        {
          gsl_matrix * A = gsl_matrix_alloc(M, N);

          for (p = 0; p < GSL_MIN(M, 10); ++p)
            {
              for (q = 0; q < GSL_MIN(N, 10); ++q)
                {
                  create_band_matrix(p, q, A, r);
                  s += test_QR_band_decomp_eps(p, q, A, 1.0e5 * GSL_MAX(M,N) * GSL_DBL_EPSILON, "QR_band_decomp random");
                }
            }

          gsl_matrix_free(A);
        }
    }

  return s;
}

#if 0
static int
test_QR_band_decomp(void)
{
  int f;
  int s = 0;

  f = test_QR_decomp_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp m(3,5)");
  s += f;

  f = test_QR_decomp_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp m(5,3)");
  s += f;

  f = test_QR_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp hilbert(2)");
  s += f;

  f = test_QR_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp hilbert(3)");
  s += f;

  f = test_QR_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp hilbert(4)");
  s += f;

  f = test_QR_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp hilbert(12)");
  s += f;

  f = test_QR_decomp_dim(vander2, 8.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp vander(2)");
  s += f;

  f = test_QR_decomp_dim(vander3, 64.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp vander(3)");
  s += f;

  f = test_QR_decomp_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp vander(4)");
  s += f;

  f = test_QR_decomp_dim(vander12, 0.0005); /* FIXME: bad accuracy */
  gsl_test(f, "  QR_decomp vander(12)");
  s += f;

  return s;
}
#endif


text:
#pragma once

#include <nextalign/nextalign.h>

#include <gsl/string_span>
#include <string>
#include <vector>

template<typename Letter>
using SequenceSpan = gsl::basic_string_span<Letter, gsl::dynamic_extent>;

using NucleotideSequenceSpan = SequenceSpan<Nucleotide>;

using AminoacidSequenceSpan = SequenceSpan<Aminoacid>;

Nucleotide toNucleotide(char nuc);

Nucleotide stringToNuc(const std::string& nuc);

char nucToChar(Nucleotide nuc);

std::string nucToString(Nucleotide nuc);

Aminoacid charToAa(char aa);

Aminoacid stringToAa(const std::string& aa);

char aaToChar(Aminoacid aa);

std::string aaToString(Aminoacid aa);

template<typename Letter>
struct LetterTag {};

template<typename Letter>
inline Letter stringToLetter(const std::string& str, LetterTag<Letter>);

template<>
inline Nucleotide stringToLetter<Nucleotide>(const std::string& str, LetterTag<Nucleotide>) {
  return stringToNuc(str);
}

template<>
inline Aminoacid stringToLetter<Aminoacid>(const std::string& str, LetterTag<Aminoacid>) {
  return stringToAa(str);
}

template<typename Letter>
inline std::string letterToString(Letter letter);

template<>
inline std::string letterToString(Nucleotide letter) {
  return nucToString(letter);
}

template<>
inline std::string letterToString(Aminoacid letter) {
  return aaToString(letter);
}

std::vector<Insertion> toInsertionsExternal(const std::vector<InsertionInternal<Nucleotide>>& insertions);

std::vector<Peptide> toPeptidesExternal(const std::vector<PeptideInternal>& peptides);

std::vector<RefPeptide> toRefPeptidesExternal(const std::vector<RefPeptideInternal>& peptides);

inline std::ostream& operator<<(std::ostream& os, const Nucleotide& nuc) {
  os << "'" << nucToString(nuc) << "'";
  return os;
}

inline std::ostream& operator<<(std::ostream& os, const NucleotideSequence& seq) {
  os << "\"";
  for (const auto& nuc : seq) {
    os << nucToString(nuc);
  }
  os << "\"";
  return os;
}

inline std::ostream& operator<<(std::ostream& os, const Aminoacid& aa) {
  os << "'" << aaToString(aa) << "'";
  return os;
}

inline std::ostream& operator<<(std::ostream& os, const AminoacidSequence& seq) {
  os << "\"";
  for (const auto& aa : seq) {
    os << aaToString(aa);
  }
  os << "\"";
  return os;
}

inline std::ostream& operator<<(std::ostream& os, const Range& f) {
  os << "{ " << f.begin << ", " << f.end << " }";
  return os;
}

inline std::ostream& operator<<(std::ostream& os, const FrameShiftContext& f) {
  os << "{ "                        //
     << "codon: " << f.codon << ", "//
     << "}"                         //
    ;
  return os;
}

inline std::ostream& operator<<(std::ostream& os, const FrameShiftResult& f) {
  os << "{ "                                      //
     << "geneName: \"" << f.geneName << "\", "    //
     << "nucRel: " << f.nucRel << ", "            //
     << "nucAbs: " << f.nucAbs << ", "            //
     << "codon: " << f.codon << ", "              //
     << "gapsLeading: " << f.gapsLeading << ", "  //
     << "gapsTrailing: " << f.gapsTrailing << ", "//
     << "}"                                       //
    ;
  return os;
}

max_stars_repo_name: neherlab/nextclade

max_issues_repo_name: neherlab/nextclade

max_forks_repo_name: neherlab/nextclade


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

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_psi.h>
#include <gsl/gsl_sf_bessel.h>

#include "error.h"

#include "bessel.h"

/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/

/* [Abramowitz+Stegun, 9.6.11]
 * assumes n >= 1
 */
static
int
bessel_Kn_scaled_small_x(const int n, const double x, gsl_sf_result * result)
{
  int k;
  double y = 0.25 * x * x;
  double ln_x_2 = log(0.5*x);
  double ex = exp(x);
  gsl_sf_result ln_nm1_fact;
  double k_term;
  double term1, sum1, ln_pre1;
  double term2, sum2, pre2;

  gsl_sf_lnfact_e((unsigned int)(n-1), &ln_nm1_fact);

  ln_pre1 = -n*ln_x_2 + ln_nm1_fact.val;
  if(ln_pre1 > GSL_LOG_DBL_MAX - 3.0) GSL_ERROR ("error", GSL_EOVRFLW);

  sum1 = 1.0;
  k_term = 1.0;
  for(k=1; k<=n-1; k++) {
    k_term *= -y/(k * (n-k));
    sum1 += k_term;
  }
  term1 = 0.5 * exp(ln_pre1) * sum1;

  pre2 = 0.5 * exp(n*ln_x_2);
  if(pre2 > 0.0) {
    const int KMAX = 20;
    gsl_sf_result psi_n;
    gsl_sf_result npk_fact;
    double yk = 1.0;
    double k_fact  = 1.0;
    double psi_kp1 = -M_EULER;
    double psi_npkp1;
    gsl_sf_psi_int_e(n, &psi_n);
    gsl_sf_fact_e((unsigned int)n, &npk_fact);
    psi_npkp1 = psi_n.val + 1.0/n;
    sum2 = (psi_kp1 + psi_npkp1 - 2.0*ln_x_2)/npk_fact.val;
    for(k=1; k<KMAX; k++) {
      psi_kp1   += 1.0/k;
      psi_npkp1 += 1.0/(n+k);
      k_fact    *= k;
      npk_fact.val *= n+k;
      yk *= y;
      k_term = yk*(psi_kp1 + psi_npkp1 - 2.0*ln_x_2)/(k_fact*npk_fact.val);
      sum2 += k_term;
    }
    term2 = ( GSL_IS_ODD(n) ? -1.0 : 1.0 ) * pre2 * sum2;
  }
  else {
    term2 = 0.0;
  }

  result->val  = ex * (term1 + term2);
  result->err  = ex * GSL_DBL_EPSILON * (fabs(ln_pre1)*fabs(term1) + fabs(term2));
  result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

  return GSL_SUCCESS;
}

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

int gsl_sf_bessel_Kn_scaled_e(int n, const double x, gsl_sf_result * result)
{
  n = abs(n); /* K(-n, z) = K(n, z) */

  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(n == 0) {
    return gsl_sf_bessel_K0_scaled_e(x, result);
  }
  else if(n == 1) {
    return gsl_sf_bessel_K1_scaled_e(x, result);
  }
  else if(x <= 5.0) {
    return bessel_Kn_scaled_small_x(n, x, result);
  }
  else if(GSL_ROOT3_DBL_EPSILON * x > 0.25 * (n*n + 1)) {
    return gsl_sf_bessel_Knu_scaled_asympx_e((double)n, x, result);
  }
  else if(GSL_MIN(0.29/(n*n), 0.5/(n*n + x*x)) < GSL_ROOT3_DBL_EPSILON) {
    return gsl_sf_bessel_Knu_scaled_asymp_unif_e((double)n, x, result);
  }
  else {
    /* Upward recurrence. [Gradshteyn + Ryzhik, 8.471.1] */
    double two_over_x = 2.0/x;
    gsl_sf_result r_b_jm1;
    gsl_sf_result r_b_j;
    int stat_0 = gsl_sf_bessel_K0_scaled_e(x, &r_b_jm1);
    int stat_1 = gsl_sf_bessel_K1_scaled_e(x, &r_b_j);
    double b_jm1 = r_b_jm1.val;
    double b_j   = r_b_j.val;
    double b_jp1;
    int j;

    for(j=1; j<n; j++) {
      b_jp1 = b_jm1 + j * two_over_x * b_j;
      b_jm1 = b_j;
      b_j   = b_jp1;
    }

    result->val  = b_j;
    result->err  = n * (fabs(b_j) * (fabs(r_b_jm1.err/r_b_jm1.val) + fabs(r_b_j.err/r_b_j.val)));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

    return GSL_ERROR_SELECT_2(stat_0, stat_1);
  }
}

int gsl_sf_bessel_Kn_e(const int n, const double x, gsl_sf_result * result)
{
  const int status = gsl_sf_bessel_Kn_scaled_e(n, x, result);
  const double ex = exp(-x);
  result->val *= ex;
  result->err *= ex;
  result->err += x * GSL_DBL_EPSILON * fabs(result->val);
  return status;
}

int gsl_sf_bessel_Kn_scaled_array(const int nmin, const int nmax, const double x, double * result_array)
{
  /* CHECK_POINTER(result_array) */

  if(nmin < 0 || nmax < nmin || x <= 0.0) {
    int j;
    for(j=0; j<=nmax-nmin; j++) result_array[j] = 0.0;
    GSL_ERROR ("domain error", GSL_EDOM);
  }
  else if(nmax == 0) {
    gsl_sf_result b;
    int stat = gsl_sf_bessel_K0_scaled_e(x, &b);
    result_array[0] = b.val;
    return stat;
  }
  else {
    double two_over_x = 2.0/x;
    gsl_sf_result r_Knm1;
    gsl_sf_result r_Kn;
    int stat_0 = gsl_sf_bessel_Kn_scaled_e(nmin,   x, &r_Knm1);
    int stat_1 = gsl_sf_bessel_Kn_scaled_e(nmin+1, x, &r_Kn);
    int stat = GSL_ERROR_SELECT_2(stat_0, stat_1);
    double Knp1;
    double Kn   = r_Kn.val;
    double Knm1 = r_Knm1.val;
    int n;

    for(n=nmin+1; n<=nmax+1; n++) {
      if(Knm1 < GSL_DBL_MAX) {
        result_array[n-1-nmin] = Knm1;
        Knp1 = Knm1 + n * two_over_x * Kn;
        Knm1 = Kn;
        Kn   = Knp1;
      }
      else {
        /* Overflow. Set the rest of the elements to
         * zero and bug out.
         * FIXME: Note: this relies on the convention
         * that the test x < DBL_MIN fails for x not
         * a number. This may be only an IEEE convention,
         * so the portability is unclear.
         */
        int j;
        for(j=n; j<=nmax+1; j++) result_array[j-1-nmin] = 0.0;
        GSL_ERROR ("overflow", GSL_EOVRFLW);
      }
    }

    return stat;
  }
}

int
gsl_sf_bessel_Kn_array(const int nmin, const int nmax, const double x, double * result_array)
{
  int status = gsl_sf_bessel_Kn_scaled_array(nmin, nmax, x, result_array);
  double ex = exp(-x);
  int i;
  for(i=0; i<=nmax-nmin; i++) result_array[i] *= ex;
  return status;
}

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

#include "eval.h"

double gsl_sf_bessel_Kn_scaled(const int n, const double x)
{
  EVAL_RESULT(gsl_sf_bessel_Kn_scaled_e(n, x, &result));
}

double gsl_sf_bessel_Kn(const int n, const double x)
{
  EVAL_RESULT(gsl_sf_bessel_Kn_e(n, x, &result));
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include "gemm.h"
#if defined(USE_BLAS)
#include <cblas.h>
#else
// see https://github.com/pjreddie/darknet/blob/master/src/gemm.c

static void gemm_nn(int M, int N, int K, float ALPHA,
                    float *A, int lda,
                    float *B, int ldb,
                    float *C, int ldc)
{
    int i;
    #pragma omp parallel for
    for (i = 0; i < M; ++i)
    {
        int j, k;
        for (k = 0; k < K; ++k)
        {
            float A_PART = ALPHA * A[i * lda + k];
            for (j = 0; j < N; ++j)
            {
                C[i * ldc + j] += A_PART * B[k * ldb + j];
            }
        }
    }
}

static void gemm_nt(int M, int N, int K, float ALPHA,
                    float *A, int lda,
                    float *B, int ldb,
                    float *C, int ldc)
{
    int i;
    #pragma omp parallel for
    for (i = 0; i < M; ++i)
    {
        int j, k;
        for (j = 0; j < N; ++j)
        {
            float sum = 0;
            for (k = 0; k < K; ++k)
            {
                sum += ALPHA * A[i * lda + k] * B[j * ldb + k];
            }
            C[i * ldc + j] += sum;
        }
    }
}

static void gemm_tn(int M, int N, int K, float ALPHA,
                    float *A, int lda,
                    float *B, int ldb,
                    float *C, int ldc)
{
    int i;
    #pragma omp parallel for
    for (i = 0; i < M; ++i)
    {
        int j, k;
        for (k = 0; k < K; ++k)
        {
            float A_PART = ALPHA * A[k * lda + i];
            for (j = 0; j < N; ++j)
            {
                C[i * ldc + j] += A_PART * B[k * ldb + j];
            }
        }
    }
}

static void gemm_tt(int M, int N, int K, float ALPHA,
                    float *A, int lda,
                    float *B, int ldb,
                    float *C, int ldc)
{
    int i;
    #pragma omp parallel for
    for (i = 0; i < M; ++i)
    {
        int j, k;
        for (j = 0; j < N; ++j)
        {
            float sum = 0;
            for (k = 0; k < K; ++k)
            {
                sum += ALPHA * A[i + k * lda] * B[k + j * ldb];
            }
            C[i * ldc + j] += sum;
        }
    }
}
#endif

void gemm(int TA, int TB, int M, int N, int K, float ALPHA,
          data_val_t **A, int offa, int lda,
          data_val_t **B, int offb, int ldb,
          float BETA,
          data_val_t **C, int offc, int ldc)
{
#if defined(USE_BLAS)
    cblas_sgemm(
        CblasRowMajor,
        TA ? CblasTrans : CblasNoTrans,
        TB ? CblasTrans : CblasNoTrans,
        M, N, K, ALPHA, *A + offa, lda, *B + offb, ldb, BETA, *C + offc, ldc);
#else
    int i, j;
    for (i = 0; i < M; ++i)
    {
        for (j = 0; j < N; ++j)
        {
            (*C + offc)[i * ldc + j] *= BETA;
        }
    }
    if (!TA && !TB)
        gemm_nn(M, N, K, ALPHA, *A + offa, lda, *B + offb, ldb, *C + offc, ldc);
    else if (TA && !TB)
        gemm_tn(M, N, K, ALPHA, *A + offa, lda, *B + offb, ldb, *C + offc, ldc);
    else if (!TA && TB)
        gemm_nt(M, N, K, ALPHA, *A + offa, lda, *B + offb, ldb, *C + offc, ldc);
    else
        gemm_tt(M, N, K, ALPHA, *A + offa, lda, *B + offb, ldb, *C + offc, ldc);
#endif
}

max_stars_repo_name: yang-le/cnet

max_issues_repo_name: yang-le/cnet

max_forks_repo_name: yang-le/cnet


text:
/*
gsl-sprng.h

Code declaring a new GSL random number type "gsl_rng_sprng20"
which is a thin wrapper over the SPRNG 2.0 parallel random
number generator.

To use, just add the line:
#include "gsl-sprng.h"
immediately after the line:
#include <gsl/gsl_rng.h>
near the start of your code. The new type should now be available.
Make sure you alloc the rng on each processor. If you wish to
set a seed, you should set it to be the same on each processor.

Darren Wilkinson
d.j.wilkinson@ncl.ac.uk
http://www.staff.ncl.ac.uk/d.j.wilkinson/

Last updated: 28/8/2002

*/

#define SIMPLE_SPRNG
#define USE_MPI
#include "sprng.h"

static void sprng_set(void * vstate,unsigned long int s)
{
  init_sprng(DEFAULT_RNG_TYPE,s,SPRNG_DEFAULT);
}

static unsigned long sprng_get(void * vstate)
{
  return( (long) isprng() );
}

static double sprng_get_double(void * vstate)
{
  return( (double) sprng());
}

static const gsl_rng_type sprng_type =
  {"sprng20",        /* name */
   0x7fffffffUL,     /* RAND_MAX */
   0,                /* RAND_MIN */
   0,                /* size of state - not sure about this */
   &sprng_set,          /* initialisation */
   &sprng_get,          /* get integer RN */
   &sprng_get_double};  /* get double RN */

const gsl_rng_type *gsl_rng_sprng20 = &sprng_type;

/* eof */


text:
#include <gsl/gsl_sf_hyperg.h>

#define kl 2
#define ll 3
#define alpha 1

double potential(double x) {
  return 0.5 * pow(alpha, 2) *
         (kl * (kl - 1) / pow(sin(alpha * x), 2) +
          ll * (ll - 1) / pow(cos(alpha * x), 2));
}
double ground_state_energy() { return 0.5 * pow(alpha, 2) * pow(kl + ll, 2); }
double exact(double x) {
  return pow(cos(alpha * x), ll) * pow(sin(alpha * x), kl) /
         sqrt(3.0 * M_PI / 512.0);
}

#undef V0
#undef kl
#undef ll
#undef alpha


text:
/*
 * Library:  ransampl (random number sampling)
 *
 * File:     sampling1.c
 *
 * Contents: Draw samples from a given discrete probability distribution;
 *           specifically: draw representative inhabitants of the nine states
 *
 * Note:     Any modification of this example should be copied to
 *           the manual page source ransampl.pod and to the wiki.
 *
 * Author:   Joachim Wuttke 2013
 *
 * Licence:  see ../COPYING (FreeBSD)
 *
 * Homepage: apps.jcns.fz-juelich.de/ransampl
 */

#include <stdio.h>
#include <gsl/gsl_rng.h>
#include "ransampl.h"

#define RANSAMPL_SHALL_REVEAL_ITS_INTERNALS 1

int main()
{
    const int M=1000000;
    int i, m;

    // Discrete probability distribution example:
    const int n = 9;
    // states of Austria
    const char* names[] = {
        "Wien", "Niederoesterreich", "Oberoesterreich", "Tirol",
        "Kaernten", "Salzburg", "Vorarlberg", "Burgenland", "Steiermark" };
    // inhabitants in millions as of 2011 according to www.statistik.at
    double p[] = { 1.721573, 1.614661, 1.415020, .711161,
                   .558056, .532713, .370833, .285377, .1211506 };

    // Initialize random number generator:
    gsl_rng_env_setup();
    gsl_rng* rng = gsl_rng_alloc( gsl_rng_default );

    // Allocate workspace, and precompute tables:
    printf( "Precomputing tables ...\n" );
    ransampl_ws* ws = ransampl_alloc( n );
    ransampl_set( ws, p );

#ifdef RANSAMPL_SHALL_REVEAL_ITS_INTERNALS
    // Inspect tables:
    printf( "  %-3s  %-3s  %-9s\n", "i", "alias", "prob" );
    for ( int i=0; i<n; ++i )
        printf( "  %3i  %3i  %9.7f\n", i, ws->alias[i], ws->prob[i] );
#endif

    // Draw M random samples; accumulate statistic in histogram 'cumul':
    printf( "Drawing %i samples ...\n", M );
    double cumul[n];
    for ( i=0; i<n; ++i )
        cumul[i] = 0;
    for ( m=0; m<M; ++m ) {
        i = ransampl_draw( ws, gsl_rng_uniform(rng), gsl_rng_uniform(rng) );
        cumul[i] += 1;
    }

    // Print given probability and obtained frequency:
    printf( "Result (input->output):\n");
    double sum = 0;
    for ( int i=0; i<n; ++i )
        sum += p[i];
    printf( "  %-18s  %-9s  %-9s  %-9s\n", "state", "N (Mio.)", "rel", "sim" );
    for ( int i=0; i<n; ++i )
        printf( "  %-18s  %9.7f  %9.7f  %9.7f\n",
                names[i], p[i], p[i]/sum, ((double)cumul[i])/M );

    // Free workspace and terminate:
    ransampl_free( ws );
    return 0;
}

max_stars_repo_name: rmaganza/parsec-improved

max_issues_repo_name: rmaganza/parsec-improved

max_forks_repo_name: rmaganza/parsec-improved


text:
#ifndef MATRIX_H_
#define MATRIX_H_

#include "matrix_funcs.h"
#ifdef NUMPY_INTERFACE
#include <Python.h>
#include <arrayobject.h>
#endif
#include <limits>
#include <assert.h>
#include <stdio.h>
#include <string.h>

#ifdef USE_MKL
#include <mkl.h>
#include <mkl_cblas.h>
#include <mkl_vsl.h>
#include <mkl_vml.h>

#define IS_MKL true

#ifdef DOUBLE_PRECISION
#define MKL_UNIFORM vdRngUniform
#define MKL_NORMAL vdRngGaussian
#define MKL_UNIFORM_RND_METHOD VSL_METHOD_DUNIFORM_STD_ACCURATE
#define MKL_GAUSSIAN_RND_METHOD VSL_METHOD_DGAUSSIAN_BOXMULLER
#define MKL_EXP vdExp
#define MKL_RECIP vdInv
#define MKL_SQUARE vdSqr
#define MKL_TANH vdTanh
#define MKL_LOG vdLn
#define MKL_VECMUL vdMul
#define MKL_VECDIV vdDiv
#else
#define MKL_UNIFORM vsRngUniform
#define MKL_NORMAL vsRngGaussian
#define MKL_UNIFORM_RND_METHOD VSL_METHOD_SUNIFORM_STD_ACCURATE
#define MKL_GAUSSIAN_RND_METHOD VSL_METHOD_SGAUSSIAN_BOXMULLER
#define MKL_EXP vsExp
#define MKL_RECIP vsInv
#define MKL_SQUARE vsSqr
#define MKL_TANH vsTanh
#define MKL_LOG vsLn
#define MKL_VECMUL vsMul
#define MKL_VECDIV vsDiv
#endif /* DOUBLE_PRECISION */

#else
#include <cblas.h>
#define IS_MKL false
#endif /* USE_MKL */

#ifdef DOUBLE_PRECISION
#define CBLAS_GEMM cblas_dgemm
#define CBLAS_SCAL cblas_dscal
#define CBLAS_AXPY cblas_daxpy
#else
#define CBLAS_GEMM cblas_sgemm
#define CBLAS_SCAL cblas_sscal
#define CBLAS_AXPY cblas_saxpy
#endif /* DOUBLE_PRECISION */

#define MTYPE_MAX numeric_limits<MTYPE>::max()

class Matrix {
private:
    MTYPE* _data;
    bool _ownsData;
    int _numRows, _numCols;
    int _numElements;
    int _numDataBytes;
    CBLAS_TRANSPOSE _trans;

    void _init(MTYPE* data, int numRows, int numCols, bool transpose, bool ownsData);
    void _tileTo2(Matrix& target) const;
    void _copyAllTo(Matrix& target) const;
    MTYPE _sum_column(int col) const;
    MTYPE _sum_row(int row) const;
    MTYPE _aggregate(MTYPE(*agg_func)(MTYPE, MTYPE), MTYPE initialValue) const;
    void _aggregate(int axis, Matrix& target, MTYPE(*agg_func)(MTYPE, MTYPE), MTYPE initialValue) const;
    MTYPE _aggregateRow(int row, MTYPE(*agg_func)(MTYPE, MTYPE), MTYPE initialValue) const;
    MTYPE _aggregateCol(int row, MTYPE(*agg_func)(MTYPE, MTYPE), MTYPE initialValue) const;
    void _updateDims(int numRows, int numCols);
    void _applyLoop(MTYPE(*func)(MTYPE));
    void _applyLoop(MTYPE (*func)(MTYPE), Matrix& target);
    void _applyLoop2(const Matrix& a, MTYPE(*func)(MTYPE, MTYPE), Matrix& target) const;
    void _applyLoop2(const Matrix& a, MTYPE (*func)(MTYPE,MTYPE, MTYPE), MTYPE scalar, Matrix& target) const;
    void _applyLoopScalar(const MTYPE scalar, MTYPE(*func)(MTYPE, MTYPE), Matrix& target) const;
    void _checkBounds(int startRow, int endRow, int startCol, int endCol) const;
    void _divideByVector(const Matrix& vec, Matrix& target);
    inline int _getNumColsBackEnd() const {
        return _trans == CblasNoTrans ? _numCols : _numRows;
    }
public:
    enum FUNCTION {
        TANH, RECIPROCAL, SQUARE, ABS, EXP, LOG, ZERO, ONE, LOGISTIC1, LOGISTIC2
    };
    Matrix();
    Matrix(int numRows, int numCols);
#ifdef NUMPY_INTERFACE
    Matrix(const PyArrayObject *src);
#endif
    Matrix(const Matrix &like);
    Matrix(MTYPE* data, int numRows, int numCols);
    Matrix(MTYPE* data, int numRows, int numCols, bool transpose);
    ~Matrix();

    inline MTYPE& getCell(int i, int j) const {
//        assert(i >= 0 && i < _numRows);
//        assert(j >= 0 && j < _numCols);
        if (_trans == CblasTrans) {
            return _data[j * _numRows + i];
        }
        return _data[i * _numCols + j];
    }

    MTYPE& operator()(int i, int j) const {
        return getCell(i, j);
    }

    inline MTYPE* getData() const {
        return _data;
    }

    inline bool isView() const {
        return !_ownsData;
    }

    inline int getNumRows() const {
        return _numRows;
    }

    inline int getNumCols() const {
        return _numCols;
    }

    inline int getNumDataBytes() const {
        return _numDataBytes;
    }

    inline int getNumElements() const {
        return _numElements;
    }

    inline CBLAS_TRANSPOSE getBLASTrans() const {
        return _trans;
    }

    inline bool isSameDims(const Matrix& a) const {
        return a.getNumRows() == getNumRows() && a.getNumCols() == getNumCols();
    }

    inline bool isTrans() const {
        return _trans == CblasTrans;
    }

    /*
     * Only use if you know what you're doing!
     * Does not update any dimensions. Just flips the _trans flag.
     *
     * Use transpose() if you want to get the transpose of this matrix.
     */
    inline void setTrans(bool trans) {
        _trans = trans ? CblasTrans : CblasNoTrans;
    }

    void apply(FUNCTION f);
    void apply(Matrix::FUNCTION f, Matrix& target);
    void subtractFromScalar(MTYPE scalar);
    void subtractFromScalar(MTYPE scalar, Matrix &target) const;
    void biggerThanScalar(MTYPE scalar);
    void smallerThanScalar(MTYPE scalar);
    void equalsScalar(MTYPE scalar);
    void biggerThanScalar(MTYPE scalar, Matrix& target) const;
    void smallerThanScalar(MTYPE scalar, Matrix& target) const;
    void equalsScalar(MTYPE scalar, Matrix& target) const;
    void biggerThan(Matrix& a);
    void biggerThan(Matrix& a, Matrix& target) const;
    void smallerThan(Matrix& a);
    void smallerThan(Matrix& a, Matrix& target) const;
    void minWith(Matrix &a);
    void minWith(Matrix &a, Matrix &target) const;
    void maxWith(Matrix &a);
    void maxWith(Matrix &a, Matrix &target) const;
    void equals(Matrix& a);
    void equals(Matrix& a, Matrix& target) const;
    void notEquals(Matrix& a) ;
    void notEquals(Matrix& a, Matrix& target) const;
    void add(const Matrix &m);
    void add(const Matrix &m, MTYPE scale);
    void add(const Matrix &m, Matrix& target);
    void add(const Matrix &m, MTYPE scale, Matrix& target);
    void subtract(const Matrix &m);
    void subtract(const Matrix &m, Matrix& target);
    void subtract(const Matrix &m, MTYPE scale);
    void subtract(const Matrix &m, MTYPE scale, Matrix& target);
    void addVector(const Matrix& vec, MTYPE scale);
    void addVector(const Matrix& vec, MTYPE scale, Matrix& target);
    void addVector(const Matrix& vec);
    void addVector(const Matrix& vec, Matrix& target);
    void addScalar(MTYPE scalar);
    void addScalar(MTYPE scalar, Matrix& target) const;
    void maxWithScalar(MTYPE scalar);
    void maxWithScalar(MTYPE scalar, Matrix &target) const;
    void minWithScalar(MTYPE scalar);
    void minWithScalar(MTYPE scalar, Matrix &target) const;
    void eltWiseMultByVector(const Matrix& vec);
    void eltWiseMultByVector(const Matrix& vec, Matrix& target);
    void eltWiseDivideByVector(const Matrix& vec);
    void eltWiseDivideByVector(const Matrix& vec, Matrix& target);
    void resize(int newNumRows, int newNumCols);
    void resize(const Matrix& like);
    Matrix& slice(int startRow, int endRow, int startCol, int endCol) const;
    void slice(int startRow, int endRow, int startCol, int endCol, Matrix &target) const;
    Matrix& sliceRows(int startRow, int endRow) const;
    void sliceRows(int startRow, int endRow, Matrix& target) const;
    Matrix& sliceCols(int startCol, int endCol) const;
    void sliceCols(int startCol, int endCol, Matrix& target) const;
    void rightMult(const Matrix &b, MTYPE scale);
    void rightMult(const Matrix &b, Matrix &target) const;
    void rightMult(const Matrix &b);
    void rightMult(const Matrix &b, MTYPE scaleAB, Matrix &target) const;
    void addProduct(const Matrix &a, const Matrix &b, MTYPE scaleAB, MTYPE scaleThis);
    void addProduct(const Matrix& a, const Matrix& b);
    void eltWiseMult(const Matrix& a);
    void eltWiseMult(const Matrix& a, Matrix& target) const;
    void eltWiseDivide(const Matrix& a);
    void eltWiseDivide(const Matrix& a, Matrix &target) const;
    Matrix& transpose() const;
    Matrix& transpose(bool hard) const;
    Matrix& tile(int timesY, int timesX) const;
    void tile(int timesY, int timesX, Matrix& target) const;
    void copy(Matrix &dest, int srcStartRow, int srcEndRow, int srcStartCol, int srcEndCol, int destStartRow, int destStartCol) const;
    Matrix& copy() const;
    void copy(Matrix& target) const;
    Matrix& sum(int axis) const;
    void sum(int axis, Matrix &target) const;
    MTYPE sum() const;
    MTYPE max() const;
    Matrix& max(int axis) const;
    void max(int axis, Matrix& target) const;
    MTYPE min() const;
    Matrix& min(int axis) const;
    void min(int axis, Matrix& target) const;
    void scale(MTYPE scale);
#ifdef USE_MKL
    void randomizeNormal(VSLStreamStatePtr stream, MTYPE mean, MTYPE stdev);
    void randomizeUniform(VSLStreamStatePtr stream);
    void randomizeNormal(VSLStreamStatePtr stream);
#else
    void randomizeNormal(MTYPE mean, MTYPE stdev);
    void randomizeUniform();
    void randomizeNormal();
#endif
    void print() const;
    void print(int startRow,int rows, int startCol,int cols) const;
    void print(int rows, int cols) const;
};

#endif /* MATRIX_H_ */

max_stars_repo_name: barapa/HF-RNN

max_issues_repo_name: barapa/HF-RNN

max_forks_repo_name: barapa/HF-RNN