text:
/* linalg/test_cholesky.c
 *
 * Copyright (C) 2016 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_choleskyc_decomp_eps(const gsl_matrix_complex * m, const double eps, const char * desc);
static int test_choleskyc_decomp(gsl_rng * r);
static int test_choleskyc_invert(gsl_rng * r);

static int
test_choleskyc_decomp_eps(const gsl_matrix_complex * m, const double eps, const char * desc)
{
  int s = 0;
  const size_t M = m->size1;
  const size_t N = m->size2;
  size_t i, j;

  gsl_matrix_complex * v  = gsl_matrix_complex_alloc(M,N);
  gsl_matrix_complex * a  = gsl_matrix_complex_alloc(M,N);
  gsl_matrix_complex * l  = gsl_matrix_complex_alloc(M,N);
  gsl_matrix_complex * lh  = gsl_matrix_complex_alloc(N,N);

  gsl_matrix_complex_memcpy(v, m);
  gsl_matrix_complex_set_zero(l);
  gsl_matrix_complex_set_zero(lh);

  s += gsl_linalg_complex_cholesky_decomp(v);

  /* Compute L L^H */

  for (i = 0; i < N ; i++)
    {
      for (j = 0; j <= i; j++)
        {
          gsl_complex vij = gsl_matrix_complex_get(v, i, j);
          gsl_matrix_complex_set (l, i, j, vij);
          gsl_matrix_complex_set (lh, j, i, gsl_complex_conjugate(vij));
        }
    }

  /* compute a = l lh */
  gsl_blas_zgemm (CblasNoTrans,
                  CblasNoTrans,
                  GSL_COMPLEX_ONE,
                  l,
                  lh,
                  GSL_COMPLEX_ZERO,
                  a);

  for (i = 0; i < M; i++)
    {
      for (j = 0; j < N; j++)
        {
          gsl_complex aij = gsl_matrix_complex_get(a, i, j);
          gsl_complex mij = gsl_matrix_complex_get(m, i, j);

          gsl_test_rel(GSL_REAL(aij), GSL_REAL(mij), eps,
                       "%s: real (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, N, i, j, GSL_REAL(aij), GSL_REAL(mij));

          gsl_test_rel(GSL_IMAG(aij), GSL_IMAG(mij), eps,
                       "%s: imag (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, N, i, j, GSL_IMAG(aij), GSL_IMAG(mij));
        }
    }

  gsl_matrix_complex_free(v);
  gsl_matrix_complex_free(a);
  gsl_matrix_complex_free(l);
  gsl_matrix_complex_free(lh);

  return s;
}

static int
test_choleskyc_decomp(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 50;
  size_t N;

  for (N = 1; N <= N_max; ++N)
    {
      gsl_matrix_complex * m = gsl_matrix_complex_alloc(N, N);

      create_posdef_complex_matrix(m, r);
      test_choleskyc_decomp_eps(m, 1.0e6 * N * GSL_DBL_EPSILON, "cholesky_complex_decomp random");

      gsl_matrix_complex_free(m);
    }

  return s;
}

static int
test_choleskyc_invert_eps(const gsl_matrix_complex * m, const double eps, const char * desc)
{
  int s = 0;
  const size_t N = m->size1;
  size_t i, j;
  gsl_matrix_complex * v  = gsl_matrix_complex_alloc(N, N);
  gsl_matrix_complex * c  = gsl_matrix_complex_alloc(N, N);

  gsl_matrix_complex_memcpy(v, m);

  s += gsl_linalg_complex_cholesky_decomp(v);
  s += gsl_linalg_complex_cholesky_invert(v);

  gsl_blas_zhemm(CblasLeft, CblasUpper, GSL_COMPLEX_ONE, m, v, GSL_COMPLEX_ZERO, c);

  /* c should be the identity matrix */
  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          gsl_complex cij = gsl_matrix_complex_get(c, i, j);
          double expected = (i == j) ? 1.0 : 0.0;

          /* check real part */
          gsl_test_rel(GSL_REAL(cij), expected, eps,
                       "%s: real (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, N, i, j, GSL_REAL(cij), expected);

          /* check imaginary part */
          gsl_test_rel(GSL_IMAG(cij), 0.0, eps,
                       "%s: imag (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, N, i, j, GSL_IMAG(cij), 0.0);
        }
    }

  gsl_matrix_complex_free(v);
  gsl_matrix_complex_free(c);

  return s;
}

static int
test_choleskyc_invert(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 50;
  size_t N;

  for (N = 1; N <= N_max; ++N)
    {
      gsl_matrix_complex * m = gsl_matrix_complex_alloc(N, N);

      create_posdef_complex_matrix(m, r);
      test_choleskyc_invert_eps(m, 32.0 * N * GSL_DBL_EPSILON, "cholesky_complex_invert random");

      gsl_matrix_complex_free(m);
    }

  return s;
}


text:
// SPDX-License-Identifier: MIT
// The MIT License (MIT)
//
// Copyright (c) 2014-2018, Institute for Software & Systems Engineering
// Copyright (c) 2018-2019, Johannes Leupolz
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
// THE SOFTWARE.

#ifndef PEMC_EXECUTABLE_MODEL_ABSTRACT_MODEL_H_
#define PEMC_EXECUTABLE_MODEL_ABSTRACT_MODEL_H_

#include <vector>
#include <gsl/span>
#include <gsl/gsl_byte>
#include <cstdint>
#include <atomic>
#include <stack>
#include <limits>
#include <functional>

#include "pemc/basic/tsc_index.h"
#include "pemc/basic/label.h"
#include "pemc/basic/model_capacity.h"
#include "pemc/basic/probability.h"
#include "pemc/basic/raw_memory.h"
#include "pemc/formula/formula.h"
#include "pemc/generic_traverser/i_transitions_calculator.h"
#include "pemc/generic_traverser/i_pre_state_storage_modifier.h"
#include "pemc/generic_traverser/i_post_state_storage_modifier.h"
#include "pemc/executable_model/i_choice_resolver.h"

namespace pemc {

  class AbstractModel {
  protected:
      IChoiceResolver* choiceResolver = nullptr;
      std::vector<std::function<bool()>> formulaEvaluators;
  public:
      AbstractModel() {}
      virtual ~AbstractModel() {}

      //AbstractModel(const AbstractModel&) = delete;

      void setChoiceResolver(IChoiceResolver* _choiceResolver);

      virtual void serialize(gsl::span<gsl::byte> position) {}

      virtual void deserialize(gsl::span<gsl::byte> position) {}

      virtual void setFormulasForLabel(const std::vector<std::shared_ptr<Formula>>& _formulas) {}

      Label calculateLabel();

      size_t choose(const gsl::span<Probability>& choices);

      size_t choose(size_t numberOfChoices);

      virtual void resetToInitialState() {}

      virtual void step() {}

      virtual int32_t getStateVectorSize() {return 0;}
  };

}

#endif  // PEMC_EXECUTABLE_MODEL_ABSTRACT_MODEL_H_

ext: h

max_stars_repo_name: joleuger/pemc

max_issues_repo_name: joleuger/pemc

max_forks_repo_name: joleuger/pemc


text:
#ifndef BIO_MATRIX_DEPENDENCIES_H_
#define BIO_MATRIX_DEPENDENCIES_H_

#include "bio/defs.h"
#include "bio/matrix.h"
#include "bio/contingency_homogeneity.h"

#include <boost/array.hpp>

#include <gsl/gsl_math.h>

#include <vector>
#include <numeric>
#include <set>

BIO_NS_START

/** Indexes a pair of bases for comparison. */
struct BasePair
{
	/** The index of the base we wish to analyse. */
	unsigned observed_index;

	/** The index of the base we condition on. */
	unsigned conditioned_index;

	BasePair(unsigned observed_index, unsigned conditioned_index);

	bool operator<(const BasePair & rhs) const;
};

std::ostream &
operator<<(std::ostream & os, const BasePair & base_pair);

/** A contingency table to compare a pair of bases. */
typedef ContingencyTable<4, 4> base_pair_contingency_table_t;

/** The contingency tables for all base pairs in a pssm. */
typedef std::map<BasePair, base_pair_contingency_table_t> pssm_contingency_tables_t;

/** The results of analysing the dependencies between a pair of bases. */
struct BasePairDependencyResults
{
	double homogeneity_bayes_factor;
};

/** The results of the dependencies between all base pairs in a pssm. */
typedef std::map< BasePair, BasePairDependencyResults > pssm_dependency_results_t;

/** An iterator pointing to a result. */
typedef pssm_dependency_results_t::const_iterator pssm_dependency_result_it;

/** A less than operator that lets us sort by scores. */
struct HomogeneityLessThen
{
	bool operator()(pssm_dependency_result_it lhs, pssm_dependency_result_it rhs) const;
};

/** An ordered set of results. */
typedef std::set< pssm_dependency_result_it, HomogeneityLessThen > pssm_dependency_result_set;

/** Contains all the sequences that were used to compose a Pssm. Indexed first by position then by sequence number. */
typedef std::vector< std::vector<char> > pssm_source_matrix_t;
typedef boost::shared_ptr< pssm_source_matrix_t > pssm_source_matrix_ptr_t;
typedef std::map< TableLink, pssm_source_matrix_ptr_t > pssm_source_map_t;

std::ostream &
operator<<(std::ostream & os, pssm_source_map_t::value_type value);

void
build_all_pssm_sources(pssm_source_map_t & pssm_source_map);

/** Take a Transfac matrix and generate a matrix from the sequences that defined it. */
void
build_pssm_source(
	const Matrix * matrix,
	pssm_source_matrix_t & pssm_source);

/** Adds a sequence to a pssm source. */
void
add_sequence_to_pssm_source(
	const seq_t & sequence,
	pssm_source_matrix_t & pssm_source);

/** Count the conditional frequencies in a matrix and generate contingency tables. */
void
count_frequencies(
	const pssm_source_matrix_t & pssm_source,
	pssm_contingency_tables_t & contingency_tables);

/** Calculate the dependencies in all the contingency tables. */
void
calculate_dependencies(
	pssm_contingency_tables_t & contingency_tables,
	pssm_dependency_results_t & results);

/** Build an ordered set of results. */
void
build_result_set(
	pssm_dependency_results_t & results,
	pssm_dependency_result_set & result_set);

struct MatrixDependencies
{
	pssm_source_matrix_t pssm_source;
	pssm_contingency_tables_t contingency_tables;
	pssm_dependency_results_t results;
	pssm_dependency_result_set result_set;

	MatrixDependencies(const pssm_source_matrix_t & pssm_source);
	MatrixDependencies(const Matrix * matrix);
	MatrixDependencies();
};

typedef std::map<const Matrix *, MatrixDependencies> matrix_dependencies_map_t;

BIO_NS_END

#endif //BIO_MATRIX_DEPENDENCIES_H_

max_stars_repo_name: JohnReid/biopsy

max_issues_repo_name: JohnReid/biopsy

max_forks_repo_name: JohnReid/biopsy


text:
/* Copyright 2016 by

  Laboratoire d'Informatique de Paris 6 - Équipe PEQUAN
  Sorbonne Universités
  UPMC Univ Paris 06
  UMR 7606, LIP6
  4, place Jussieu
  F-75252 Paris Cedex 05
  France

  Laboratoire d'Informatique de Paris 6, equipe PEQUAN,
  UPMC Universite Paris 06 - CNRS - UMR 7606 - LIP6, Paris, France

  Contributors:
  Anastasia Volkova anastasia.volkova@lip6.fr

  This software is a mathematical library whose purpose is to provide
  functions to compute the Worst-Case Peak Gain measure of Linear
  Time-Invariant Digital Filters.

  This software is governed by the CeCILL-C license under French law and
  abiding by the rules of distribution of free software.  You can  use,
  modify and/ or redistribute the software under the terms of the CeCILL-C
  license as circulated by CEA, CNRS and INRIA at the following URL
  "http://www.cecill.info".

  As a counterpart to the access to the source code and  rights to copy,
  modify and redistribute granted by the license, users are provided only
  with a limited warranty  and the software's author,  the holder of the
  economic rights,  and the successive licensors  have only  limited
  liability.

  In this respect, the user's attention is drawn to the risks associated
  with loading,  using,  modifying and/or developing or reproducing the
  software by the user in light of its specific status of free software,
  that may mean  that it is complicated to manipulate,  and  that  also
  therefore means  that it is reserved for developers  and  experienced
  professionals having in-depth computer knowledge. Users are therefore
  encouraged to load and test the software's suitability as regards their
  requirements in conditions enabling the security of their systems and/or
  data to be ensured and,  more generally, to use and operate it in the
  same conditions as regards security.

  The fact that you are presently reading this means that you have had
  knowledge of the CeCILL-C license and that you accept its terms.

  This program is distributed WITHOUT ANY WARRANTY; without even the
  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

 */

#ifndef SRC_LIB_CLAPACK_FUNCTIONS_CONFIG_H_
#define SRC_LIB_CLAPACK_FUNCTIONS_CONFIG_H_

#ifdef __cplusplus
extern "C" {
#endif

#include <stdlib.h>
#include "config.h"
#include "aux_funcs.h"

#include <f2c.h>

#ifdef CLAPACK_HEADER
	#include <clapack.h>
#endif

#ifdef LAPACK_HEADER
	#include <lapack.h>
#endif

#ifdef LAPACKE_HEADER
	#include <lapacke.h>
#endif

#ifdef BLAS_HEADER
	#include <blas.h>
#endif

#ifdef CBLAS_HEADER
	#include <cblas.h>
#endif

#ifdef BLASWRAP_HEADER
	#include <blaswrap.h>
#endif

#include <complex.h>
typedef struct complexdouble
{
	double r;
	double i;
}complexdouble;

#define doublereal double
#define lapack_int int
#define integer long int
#define lapacke_complex double _Complex

void my_zgetrf(int *m, int *n, complexdouble *A,
 	int *lda, int *ipiv, int *info);

void my_zgetri(int *n, complexdouble *A, int *lda,
		int *ipiv, complexdouble *work, int *lwork, int *info);
void my_dgeevx(int* n, double *A, int* lda, double *wr, double *wi, double *vl,
	       int* ldvl, double *vr, int* ldvr, int *ilo, int *ihi,
	       double *scale, double *abnrm, double *rconde, double *rcondv,
	       double *work, int *lwork, long int *iwork, int *info);

#ifdef __cplusplus
}
#endif

#endif /* SRC_LIB_CLAPACK_FUNCTIONS_CONFIG_H_ */

max_stars_repo_name: remi-garcia/WCPG

max_issues_repo_name: remi-garcia/WCPG

max_forks_repo_name: remi-garcia/WCPG


text:
/* linalg/qr_ur.c
 *
 * Copyright (C) 2019, 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 <string.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

/*
 * this module contains routines for the QR factorization of a matrix
 * using the recursive Level 3 BLAS algorithm of Elmroth and Gustavson with
 * additional modifications courtesy of Julien Langou.
 */

static double qrtr_householder_transform (double *v0, gsl_vector * v);

/*
gsl_linalg_QR_UR_decomp()
  Compute the QR decomposition of the "triangle on top of rectangle" matrix

  [ S ] = Q [ R ]
  [ A ]     [ 0 ]

where S is N-by-N upper triangular and A is M-by-N dense.

Inputs: S   - on input, upper triangular N-by-N matrix
              on output, R factor in upper triangle
        A   - on input, dense M-by-N matrix
              on output, Householder matrix V
        T   - (output) block reflector matrix, N-by-N

Notes:
1) Based on the Elmroth/Gustavson algorithm, taking into account the
sparse structure of the S matrix

2) The Householder matrix V has the special form:

      N
V = [ I  ] N
    [ V~ ] M

The matrix V~ is stored in A on output; the identity is not stored

3) The orthogonal matrix is

Q = I - V T V^T
*/

int
gsl_linalg_QR_UR_decomp (gsl_matrix * S, gsl_matrix * A, gsl_matrix * T)
{
  const size_t M = A->size1;
  const size_t N = S->size1;

  if (N != S->size2)
    {
      GSL_ERROR ("S matrix must be square", GSL_ENOTSQR);
    }
  else if (N != A->size2)
    {
      GSL_ERROR ("S and A have different number of columns", GSL_EBADLEN);
    }
  else if (T->size1 != N || T->size2 != N)
    {
      GSL_ERROR ("T matrix has wrong dimensions", GSL_EBADLEN);
    }
  else if (N == 1)
    {
      /* base case, compute Householder transform for single column matrix */
      double * T00 = gsl_matrix_ptr(T, 0, 0);
      double * S00 = gsl_matrix_ptr(S, 0, 0);
      gsl_vector_view v = gsl_matrix_column(A, 0);
      *T00 = qrtr_householder_transform(S00, &v.vector);
      return GSL_SUCCESS;
    }
  else
    {
      /*
       * partition matrices:
       *
       *       N1  N2              N1  N2
       * N1 [ S11 S12 ] and  N1 [ T11 T12 ]
       * N2 [  0  S22 ]      N2 [  0  T22 ]
       * M  [  A1  A2 ]
       */
      int status;
      const size_t N1 = N / 2;
      const size_t N2 = N - N1;

      gsl_matrix_view S11 = gsl_matrix_submatrix(S, 0, 0, N1, N1);
      gsl_matrix_view S12 = gsl_matrix_submatrix(S, 0, N1, N1, N2);
      gsl_matrix_view S22 = gsl_matrix_submatrix(S, N1, N1, N2, N2);

      gsl_matrix_view A1 = gsl_matrix_submatrix(A, 0, 0, M, N1);
      gsl_matrix_view A2 = gsl_matrix_submatrix(A, 0, N1, M, N2);

      gsl_matrix_view T11 = gsl_matrix_submatrix(T, 0, 0, N1, N1);
      gsl_matrix_view T12 = gsl_matrix_submatrix(T, 0, N1, N1, N2);
      gsl_matrix_view T22 = gsl_matrix_submatrix(T, N1, N1, N2, N2);

      /*
       * Eq. 2: recursively factor
       *
       *       N1           N1
       * N1 [ S11 ] = Q1 [ R11 ] N1
       * N2 [  0  ]      [  0  ] N2
       * M  [  A1 ]      [  0  ] M
       */
      status = gsl_linalg_QR_UR_decomp(&S11.matrix, &A1.matrix, &T11.matrix);
      if (status)
        return status;

      /*
       * Eq. 3:
       *
       *      N2              N2            N2
       * N1 [ R12  ] = Q1^T [ S12 ] = [   S12 - W  ] N1
       * N2 [ S22~ ]        [ S22 ]   [     S22    ] N2
       * M  [  A2~ ]        [  A2 ]   [ A2 - V1~ W ] M
       *
       * where W = T11^T ( S12 + V1~^T A2 ), using T12 as temporary storage, and
       *
       *        N1
       * V1 = [  I  ] N1
       *      [  0  ] N2
       *      [ V1~ ] M
       */
      gsl_matrix_memcpy(&T12.matrix, &S12.matrix);                                                    /* W := S12 */
      gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, &A1.matrix, &A2.matrix, 1.0, &T12.matrix);        /* W := S12 + V1~^T A2 */
      gsl_blas_dtrmm(CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, 1.0, &T11.matrix, &T12.matrix); /* W := T11^T W */
      gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1.0, &A1.matrix, &T12.matrix, 1.0, &A2.matrix);     /* A2 := A2 - V1~ W */
      gsl_matrix_sub(&S12.matrix, &T12.matrix);                                                       /* R12 := S12 - W */

      /*
       * Eq. 4: recursively factor
       *
       * [ S22~ ] = Q2~ [ R22 ]
       * [  A2~ ]       [  0  ]
       */
      status = gsl_linalg_QR_UR_decomp(&S22.matrix, &A2.matrix, &T22.matrix);
      if (status)
        return status;

      /*
       * Eq. 13: update T12 := -T11 * V1^T * V2 * T22
       *
       * where:
       *
       *        N1                N2
       * V1 = [  I  ] N1   V2 = [  0  ] N1
       *      [  0  ] N2        [  I  ] N2
       *      [ V1~ ] M         [ V2~ ] M
       *
       * Note: V1^T V2 = V1~^T V2~
       */

      gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, &A1.matrix, &A2.matrix, 0.0, &T12.matrix);           /* T12 := V1~^T * V2~ */
      gsl_blas_dtrmm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, -1.0, &T11.matrix, &T12.matrix); /* T12 := -T11 * T12 */
      gsl_blas_dtrmm(CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0, &T22.matrix, &T12.matrix); /* T12 := T12 * T22 */

      return GSL_SUCCESS;
    }
}

/*
qrtr_householder_transform()
  This routine is an optimized version of
gsl_linalg_householder_transform(), designed for the QR
decomposition of M-by-N matrices of the form:

B = [ S ]
    [ A ]

where S is N-by-N upper triangular, and A is M-by-N dense.
This routine computes a householder transformation (tau,v) of a
x so that P x = [ I - tau*v*v' ] x annihilates x(1:n-1). x will
be a subcolumn of the matrix B, and so its structure will be:

x = [ x0 ] <- 1 nonzero value for the diagonal element of S
    [ 0  ] <- N - j - 1 zeros, where j is column of matrix in [0,N-1]
    [ x  ] <- M nonzero values for the dense part A

Inputs: v0 - pointer to diagonal element of S
             on input, v0 = x0;
        v  - on input, x vector
             on output, householder vector v
*/

static double
qrtr_householder_transform (double *v0, gsl_vector * v)
{
  /* replace v[0:M-1] with a householder vector (v[0:M-1]) and
     coefficient tau that annihilate v[1:M-1] */

  double alpha, beta, tau ;

  /* compute xnorm = || [ 0 ; v ] ||, ignoring zero part of vector */
  double xnorm = gsl_blas_dnrm2(v);

  if (xnorm == 0)
    {
      return 0.0; /* tau = 0 */
    }

  alpha = *v0;
  beta = - GSL_SIGN(alpha) * hypot(alpha, xnorm) ;
  tau = (beta - alpha) / beta ;

  {
    double s = (alpha - beta);

    if (fabs(s) > GSL_DBL_MIN)
      {
        gsl_blas_dscal (1.0 / s, v);
        *v0 = beta;
      }
    else
      {
        gsl_blas_dscal (GSL_DBL_EPSILON / s, v);
        gsl_blas_dscal (1.0 / GSL_DBL_EPSILON, v);
        *v0 = beta;
      }
  }

  return tau;
}

ext: c

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#ifndef __GSL_PERMUTE_VECTOR_H__
#define __GSL_PERMUTE_VECTOR_H__

#include <gsl/gsl_permute_vector_complex_long_double.h>
#include <gsl/gsl_permute_vector_complex_double.h>
#include <gsl/gsl_permute_vector_complex_float.h>

#include <gsl/gsl_permute_vector_long_double.h>
#include <gsl/gsl_permute_vector_double.h>
#include <gsl/gsl_permute_vector_float.h>

#include <gsl/gsl_permute_vector_ulong.h>
#include <gsl/gsl_permute_vector_long.h>

#include <gsl/gsl_permute_vector_uint.h>
#include <gsl/gsl_permute_vector_int.h>

#include <gsl/gsl_permute_vector_ushort.h>
#include <gsl/gsl_permute_vector_short.h>

#include <gsl/gsl_permute_vector_uchar.h>
#include <gsl/gsl_permute_vector_char.h>

#endif /* __GSL_PERMUTE_VECTOR_H__ */

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: ielomariala/Hex-Game


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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_errno.h>

void
gsl_ieee_env_setup (void)
{
  const char * p = getenv("GSL_IEEE_MODE") ;

  int precision = 0, rounding = 0, exception_mask = 0 ;

  int comma = 0 ;

  if (p == 0)  /* GSL_IEEE_MODE environment variable is not set */
    return ;

  if (*p == '\0') /* GSL_IEEE_MODE environment variable is empty */
    return ;

  gsl_ieee_read_mode_string (p, &precision, &rounding, &exception_mask) ;

  gsl_ieee_set_mode (precision, rounding, exception_mask) ;

  fprintf(stderr, "GSL_IEEE_MODE=\"") ;

  /* Print string with a preceeding comma if the list has already begun */

#define PRINTC(x) do {if(comma) fprintf(stderr,","); fprintf(stderr,x); comma++ ;} while(0)

  switch (precision)
    {
    case GSL_IEEE_SINGLE_PRECISION:
      PRINTC("single-precision") ;
      break ;
    case GSL_IEEE_DOUBLE_PRECISION:
      PRINTC("double-precision") ;
      break ;
    case GSL_IEEE_EXTENDED_PRECISION:
      PRINTC("extended-precision") ;
      break ;
    }

  switch (rounding)
    {
    case GSL_IEEE_ROUND_TO_NEAREST:
      PRINTC("round-to-nearest") ;
      break ;
    case GSL_IEEE_ROUND_DOWN:
      PRINTC("round-down") ;
      break ;
    case GSL_IEEE_ROUND_UP:
      PRINTC("round-up") ;
      break ;
    case GSL_IEEE_ROUND_TO_ZERO:
      PRINTC("round-to-zero") ;
      break ;
    }

  if ((exception_mask & GSL_IEEE_MASK_ALL) == GSL_IEEE_MASK_ALL)
    {
      PRINTC("mask-all") ;
    }
  else if ((exception_mask & GSL_IEEE_MASK_ALL) == 0)
    {
      PRINTC("trap-common") ;
    }
  else
    {
      if (exception_mask & GSL_IEEE_MASK_INVALID)
	PRINTC("mask-invalid") ;

      if (exception_mask & GSL_IEEE_MASK_DENORMALIZED)
	PRINTC("mask-denormalized") ;

      if (exception_mask & GSL_IEEE_MASK_DIVISION_BY_ZERO)
	PRINTC("mask-division-by-zero") ;

      if (exception_mask & GSL_IEEE_MASK_OVERFLOW)
	PRINTC("mask-overflow") ;

      if (exception_mask & GSL_IEEE_MASK_UNDERFLOW)
	PRINTC("mask-underflow") ;
    }

  if (exception_mask & GSL_IEEE_TRAP_INEXACT)
    PRINTC("trap-inexact") ;

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


text:
/* specfunc/bessel_Yn.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"
#include "bessel_amp_phase.h"
#include "bessel_olver.h"

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

/* assumes n >= 1 */
static int bessel_Yn_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);
  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 = -exp(ln_pre1) * sum1 / M_PI;

  pre2 = -exp(n*ln_x_2) / M_PI;
  if(fabs(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./k;
      psi_npkp1 += 1./(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 = pre2 * sum2;
  }
  else {
    term2 = 0.0;
  }

  result->val  = term1 + term2;
  result->err  = 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_Yn_e(int n, const double x, gsl_sf_result * result)
{
  int sign = 1;

  if(n < 0) {
    /* reduce to case n >= 0 */
    n = -n;
    if(GSL_IS_ODD(n)) sign = -1;
  }

  /* CHECK_POINTER(result) */

  if(n == 0) {
    int status = gsl_sf_bessel_Y0_e(x, result);
    result->val *= sign;
    return status;
  }
  else if(n == 1) {
    int status = gsl_sf_bessel_Y1_e(x, result);
    result->val *= sign;
    return status;
  }
  else {
    if(x <= 0.0) {
      DOMAIN_ERROR(result);
    }
    if(x < 5.0) {
      int status = bessel_Yn_small_x(n, x, result);
      result->val *= sign;
      return status;
    }
    else if(GSL_ROOT3_DBL_EPSILON * x > (n*n + 1.0)) {
      int status = gsl_sf_bessel_Ynu_asympx_e((double)n, x, result);
      result->val *= sign;
      return status;
    }
    else if(n > 50) {
      int status = gsl_sf_bessel_Ynu_asymp_Olver_e((double)n, x, result);
      result->val *= sign;
      return status;
    }
    else {
      double two_over_x = 2.0/x;
      gsl_sf_result r_by;
      gsl_sf_result r_bym;
      int stat_1 = gsl_sf_bessel_Y1_e(x, &r_by);
      int stat_0 = gsl_sf_bessel_Y0_e(x, &r_bym);
      double bym = r_bym.val;
      double by  = r_by.val;
      double byp;
      int j;

      for(j=1; j<n; j++) {
        byp = j*two_over_x*by - bym;
        bym = by;
        by  = byp;
      }
      result->val  = sign * by;
      result->err  = fabs(result->val) * (fabs(r_by.err/r_by.val) + fabs(r_bym.err/r_bym.val));
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

      return GSL_ERROR_SELECT_2(stat_1, stat_0);
    }
  }
}

int
gsl_sf_bessel_Yn_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 ("error", GSL_EDOM);
  }
  else {
    gsl_sf_result r_Ynm1;
    gsl_sf_result r_Yn;
    int stat_nm1 = gsl_sf_bessel_Yn_e(nmin,   x, &r_Ynm1);
    int stat_n   = gsl_sf_bessel_Yn_e(nmin+1, x, &r_Yn);
    double Ynp1;
    double Yn   = r_Yn.val;
    double Ynm1 = r_Ynm1.val;
    int n;

    int stat = GSL_ERROR_SELECT_2(stat_nm1, stat_n);

    if(stat == GSL_SUCCESS) {
      for(n=nmin+1; n<=nmax+1; n++) {
        result_array[n-nmin-1] = Ynm1;
        Ynp1 = -Ynm1 + 2.0*n/x * Yn;
        Ynm1 = Yn;
        Yn   = Ynp1;
      }
    }
    else {
      for(n=nmin; n<=nmax; n++) {
        result_array[n-nmin] = 0.0;
      }
    }

    return stat;
  }
}

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

#include "eval.h"

double gsl_sf_bessel_Yn(const int n, const double x)
{
  EVAL_RESULT(gsl_sf_bessel_Yn_e(n, x, &result));
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

#ifndef __GSL_VECTOR_UCHAR_H__
#define __GSL_VECTOR_UCHAR_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_uchar.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 char *data;
  gsl_block_uchar *block;
  int owner;
}
gsl_vector_uchar;

typedef struct
{
  gsl_vector_uchar vector;
} _gsl_vector_uchar_view;

typedef _gsl_vector_uchar_view gsl_vector_uchar_view;

typedef struct
{
  gsl_vector_uchar vector;
} _gsl_vector_uchar_const_view;

typedef const _gsl_vector_uchar_const_view gsl_vector_uchar_const_view;

/* Allocation */

GSL_EXPORT gsl_vector_uchar *gsl_vector_uchar_alloc (const size_t n);
GSL_EXPORT gsl_vector_uchar *gsl_vector_uchar_calloc (const size_t n);

GSL_EXPORT gsl_vector_uchar *gsl_vector_uchar_alloc_from_block (gsl_block_uchar * b,
                                                                const size_t offset,
                                                                const size_t n,
                                                                const size_t stride);

GSL_EXPORT gsl_vector_uchar *gsl_vector_uchar_alloc_from_vector (gsl_vector_uchar * v,
                                                                 const size_t offset,
                                                                 const size_t n,
                                                                 const size_t stride);

GSL_EXPORT void gsl_vector_uchar_free (gsl_vector_uchar * v);

/* Views */

GSL_EXPORT
_gsl_vector_uchar_view
gsl_vector_uchar_view_array (unsigned char *v, size_t n);

GSL_EXPORT
_gsl_vector_uchar_view
gsl_vector_uchar_view_array_with_stride (unsigned char *base,
                                         size_t stride,
                                         size_t n);

GSL_EXPORT
_gsl_vector_uchar_const_view
gsl_vector_uchar_const_view_array (const unsigned char *v, size_t n);

GSL_EXPORT
_gsl_vector_uchar_const_view
gsl_vector_uchar_const_view_array_with_stride (const unsigned char *base,
                                               size_t stride,
                                               size_t n);

GSL_EXPORT
_gsl_vector_uchar_view
gsl_vector_uchar_subvector (gsl_vector_uchar *v,
                            size_t i,
                            size_t n);

GSL_EXPORT
_gsl_vector_uchar_view
gsl_vector_uchar_subvector_with_stride (gsl_vector_uchar *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_EXPORT
_gsl_vector_uchar_const_view
gsl_vector_uchar_const_subvector (const gsl_vector_uchar *v,
                                  size_t i,
                                  size_t n);

GSL_EXPORT
_gsl_vector_uchar_const_view
gsl_vector_uchar_const_subvector_with_stride (const gsl_vector_uchar *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_EXPORT unsigned char gsl_vector_uchar_get (const gsl_vector_uchar * v, const size_t i);
GSL_EXPORT void gsl_vector_uchar_set (gsl_vector_uchar * v, const size_t i, unsigned char x);

GSL_EXPORT unsigned char *gsl_vector_uchar_ptr (gsl_vector_uchar * v, const size_t i);
GSL_EXPORT const unsigned char *gsl_vector_uchar_const_ptr (const gsl_vector_uchar * v, const size_t i);

GSL_EXPORT void gsl_vector_uchar_set_zero (gsl_vector_uchar * v);
GSL_EXPORT void gsl_vector_uchar_set_all (gsl_vector_uchar * v, unsigned char x);
GSL_EXPORT int gsl_vector_uchar_set_basis (gsl_vector_uchar * v, size_t i);

GSL_EXPORT int gsl_vector_uchar_fread (FILE * stream, gsl_vector_uchar * v);
GSL_EXPORT int gsl_vector_uchar_fwrite (FILE * stream, const gsl_vector_uchar * v);
GSL_EXPORT int gsl_vector_uchar_fscanf (FILE * stream, gsl_vector_uchar * v);
GSL_EXPORT int gsl_vector_uchar_fprintf (FILE * stream, const gsl_vector_uchar * v,
                                         const char *format);

GSL_EXPORT int gsl_vector_uchar_memcpy (gsl_vector_uchar * dest, const gsl_vector_uchar * src);

GSL_EXPORT int gsl_vector_uchar_reverse (gsl_vector_uchar * v);

GSL_EXPORT int gsl_vector_uchar_swap (gsl_vector_uchar * v, gsl_vector_uchar * w);
GSL_EXPORT int gsl_vector_uchar_swap_elements (gsl_vector_uchar * v, const size_t i, const size_t j);

GSL_EXPORT unsigned char gsl_vector_uchar_max (const gsl_vector_uchar * v);
GSL_EXPORT unsigned char gsl_vector_uchar_min (const gsl_vector_uchar * v);
GSL_EXPORT void gsl_vector_uchar_minmax (const gsl_vector_uchar * v, unsigned char * min_out, unsigned char * max_out);

GSL_EXPORT size_t gsl_vector_uchar_max_index (const gsl_vector_uchar * v);
GSL_EXPORT size_t gsl_vector_uchar_min_index (const gsl_vector_uchar * v);
GSL_EXPORT void gsl_vector_uchar_minmax_index (const gsl_vector_uchar * v, size_t * imin, size_t * imax);

GSL_EXPORT int gsl_vector_uchar_add (gsl_vector_uchar * a, const gsl_vector_uchar * b);
GSL_EXPORT int gsl_vector_uchar_sub (gsl_vector_uchar * a, const gsl_vector_uchar * b);
GSL_EXPORT int gsl_vector_uchar_mul (gsl_vector_uchar * a, const gsl_vector_uchar * b);
GSL_EXPORT int gsl_vector_uchar_div (gsl_vector_uchar * a, const gsl_vector_uchar * b);
GSL_EXPORT int gsl_vector_uchar_scale (gsl_vector_uchar * a, const double x);
GSL_EXPORT int gsl_vector_uchar_add_constant (gsl_vector_uchar * a, const double x);

GSL_EXPORT int gsl_vector_uchar_isnull (const gsl_vector_uchar * v);

#ifdef HAVE_INLINE

extern inline
unsigned char
gsl_vector_uchar_get (const gsl_vector_uchar * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

extern inline
void
gsl_vector_uchar_set (gsl_vector_uchar * v, const size_t i, unsigned char x)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

extern inline
unsigned char *
gsl_vector_uchar_ptr (gsl_vector_uchar * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (unsigned char *) (v->data + i * v->stride);
}

extern inline
const unsigned char *
gsl_vector_uchar_const_ptr (const gsl_vector_uchar * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const unsigned char *) (v->data + i * v->stride);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_UCHAR_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
#ifndef SEARCHPATTERNBUILDER_HHHHH
#define SEARCHPATTERNBUILDER_HHHHH

#include "detail/SearchPatternBase.h"
#include "SearchPattern.h"

#include <vector>
#include <string>
#include <memory>

#include <gsl/string_span>

namespace MPSig {

    class SearchPatternBuilder {
        std::vector<std::unique_ptr<detail::SearchPatternBase>> m_patternSteps;
    public:
        SearchPatternBuilder();

        SearchPatternBuilder& HasHexPattern(gsl::cstring_span<-1> pattern);
        SearchPatternBuilder& ReferencesBSTR(gsl::cwstring_span<-1> bstrText);

        SearchPattern Compile();
    };
}

#endif


text:
#pragma once
#include <string>
#include <ntdll.h>
#include <gsl/span>

namespace pe
{
  class module : public HINSTANCE__
  {
  public:
    module() = delete;
    uintptr_t handle() const;
    template <class T>
    inline auto rva_to(uint32_t rva)
    {
      return reinterpret_cast<T *>(reinterpret_cast<uintptr_t>(this) + rva);
    }
    template <class T>
    inline auto rva_to(uint32_t rva) const
    {
      return reinterpret_cast<const T *>(reinterpret_cast<uintptr_t>(this) + rva);
    }
    std::wstring base_name() const;
    std::wstring full_name() const;
    IMAGE_DOS_HEADER *dos_header();
    const IMAGE_DOS_HEADER *dos_header() const;
    IMAGE_NT_HEADERS *nt_header();
    const IMAGE_NT_HEADERS *nt_header() const;
    size_t size() const;
    gsl::span<class segment> segments();
    gsl::span<const class segment> segments() const;
    class segment *segment(const char* namename);
    const class segment *segment(const char* name) const;
    class export_directory *export_directory();
    const class export_directory *export_directory() const;
    void *find_function(const char *name) const;
    void *find_function(uint32_t num) const;
  };
  class module *get_module(const wchar_t *name = nullptr);
  class module *get_module_from_address(void *pc);
  const class module *get_module_from_address(const void *pc);
  class module *instance_module();
}

#include "module.inl"

max_stars_repo_name: bnsmodpolice/loginhelper

max_issues_repo_name: bnsmodpolice/loginhelper

max_forks_repo_name: bnsmodpolice/loginhelper


text:
// Copyright (C) 2015 Vincent Lejeune
// For conditions of distribution and use, see copyright notice in License.txt
#pragma once
#include <vulkan\vulkan.h>
#include <array>
#include <gsl/gsl>

#define CHECK_VKRESULT(cmd) { VkResult res = (cmd); if (res != VK_SUCCESS) throw; }

namespace structures
{
	constexpr VkAttachmentDescription attachment_description(
		VkFormat format,
		VkAttachmentLoadOp load_op,
		VkAttachmentStoreOp store_op,
		VkImageLayout initial_layout,
		VkImageLayout final_layout,
		VkSampleCountFlagBits sample_count = VK_SAMPLE_COUNT_1_BIT,
		VkAttachmentLoadOp stencil_load_op = VK_ATTACHMENT_LOAD_OP_DONT_CARE,
		VkAttachmentStoreOp stencil_store_op = VK_ATTACHMENT_STORE_OP_DONT_CARE,
		VkAttachmentDescriptionFlags flag = 0)
	{
		return{ flag, format, sample_count, load_op, store_op, stencil_load_op, stencil_store_op, initial_layout, final_layout };
	}
}

max_stars_repo_name: vlj/YAGF

max_issues_repo_name: vlj/YAGF

max_forks_repo_name: vlj/YAGF


text:
#pragma once
#include <gsl/gsl_multifit_nlinear.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_cdf.h>

#include <vector>
#include <functional>
#include <algorithm>
#include <tuple>

// this "fdf_Helper" helper class is used for create the polynomial to be fit in the runtime.
// Since GSL would like a static function, have to do some trick to make it happy
// see https://stackoverflow.com/questions/13074756/how-to-avoid-static-member-function-when-using-gsl-with-c?noredirect=1&lq=1
// https://stackoverflow.com/questions/6610046/stdfunction-and-stdbind-what-are-they-and-when-should-they-be-used
// https://stackoverflow.com/questions/7582546/using-generic-stdfunction-objects-with-member-functions-in-one-class

class PolynomialFit
{
public:
    struct data
    {
        double* x;
        double* y;
        size_t n;
        std::function<int(const gsl_vector* p, void* datafit, gsl_vector* f)> f = NULL;
        std::function<int(const gsl_vector* p, void* datafit, gsl_matrix* J)> df = NULL;
        std::function<int(const gsl_vector* p, const gsl_vector* v, void* datafit, gsl_vector* fvv)> fvv = NULL;
    };
    class fdf_Helper
    {
    public:
        static int invoke_f(const gsl_vector* p, void* datafit, gsl_vector* f) {
            PolynomialFit::data* dataf = static_cast<PolynomialFit::data*>(datafit);
            return dataf->f(p, datafit, f);
        }
        static int invoke_df(const gsl_vector* p, void* datafit, gsl_matrix* J) {
            PolynomialFit::data* dataf = static_cast<PolynomialFit::data*>(datafit);
            return dataf->df(p, datafit, J);
        }
        static int invoke_fvv(const gsl_vector* p, const gsl_vector* v, void* datafit, gsl_vector* fvv) {
            PolynomialFit::data* dataf = static_cast<PolynomialFit::data*>(datafit);
            return dataf->fvv(p, v, datafit, fvv);
        }
    };

private:

    gsl_vector* p0;     /* initial fitting parameters: a0,...,an, (b0,b1)*/
    gsl_vector* p;      /* fitting parameters */
    gsl_vector* f;      /* residual yi-y(xi) */
    gsl_matrix* covar;  /*covariance matrix, not yet multiplied by sigma, so is not variance-covariance matrix*/

    gsl_multifit_nlinear_fdf         fdf;
    gsl_multifit_nlinear_parameters  fdf_params;
    gsl_multifit_nlinear_workspace*  work;

    data                             fit_data;

    size_t                           max_iter;
    double                           ptol; /* tolerance on fitting parameter p */
    double                           gtol; /* tolerance on gradient */
    double                           ftol;

    int                              info; /* fiting stop reason, for debug*/
    double                           rss;  /* mean of sum of residual square, dof corrected*/
    gsl_vector*                      confid95;

public:

    /* model function: a0 + a1*x + a2*x^2 +...+ b0*sqrt(x-b1) */

    PolynomialFit() {};
    PolynomialFit(size_t n, double* datax, double* datay, unsigned order, bool useSqrt,
        std::vector<double> init_para,
        size_t max_iter = 200, double ptol = 1.0e-8, double gtol = 1.0e-8, double ftol = 1.0e-8);
    double polynomial(std::vector<double> para, double t);
    ~PolynomialFit();
    void solve_system();
    void set_data(double* datay);
    void set_data(size_t n, double* datax, double* datay, bool free = true);
    void set_initialP(std::vector<double> para);

    std::vector<double> calcFittedGaussian();
    std::vector<double> fittedPara() const;
    std::vector<double> confidence95Interval() const;

private:

    bool useSqrt;
    unsigned order;
    unsigned numOfPara;

    void init_dataFunc();

    //static int funcf_Wrapper/*next line is the argument of this "funcf_Wrapper"*/
    //(const gsl_vector* p, void* datafit, gsl_vector* f);

    //static int (*funcf_Wrapper/*next line is the argument of this "funcf_Wrapper"*/
    //(std::function<int(const gsl_vector* p, void* datafit, gsl_vector* f)> func_f))
    //    (const gsl_vector* p, void* datafit, gsl_vector* f);

    int func_f(const gsl_vector* p, void* datafit, gsl_vector* f);
    int func_df(const gsl_vector* p, void* datafit, gsl_matrix* J);
    int func_fvv(const gsl_vector* p, const gsl_vector* v, void* datafit, gsl_vector* fvv);

};


text:
/**
 *
 * @file core_cgetf2_nopiv.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Omar Zenati
 * @author Mathieu Faverge
 * @date 2013-02-01
 * @generated c Tue Jan  7 11:44:50 2014
 *
 **/
#include "common.h"
#include <math.h>
#include <lapacke.h>

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex32_t
 *
 *  CORE_cgetf2_nopiv computes an LU factorization of a general diagonal
 *  dominant M-by-N matrix A witout no pivoting and no blocking. It is the
 *  internal function called by CORE_cgetrf_nopiv().
 *
 *  The factorization has the form
 *     A = L * U
 *  where L is lower triangular with unit
 *  diagonal elements (lower trapezoidal if m > n), and U is upper
 *  triangular (upper trapezoidal if m < n).
 *
 *  This is the right-looking Level 3 BLAS version of the algorithm.
 *
 *******************************************************************************
 *
 *  @param[in] M
 *          The number of rows of the matrix A.  M >= 0.
 *
 *  @param[in] N
 *          The number of columns of the matrix A.  N >= 0.
 *
 *  @param[in,out] A
 *          On entry, the M-by-N matrix to be factored.
 *          On exit, the factors L and U from the factorization
 *          A = P*L*U; the unit diagonal elements of L are not stored.
 *
 *  @param[in] LDA
 *          The leading dimension of the array A.  LDA >= max(1,M).
 *
 *******************************************************************************
 *
 * @return
 *         \retval PLASMA_SUCCESS successful exit
 *         \retval <0 if INFO = -k, the k-th argument had an illegal value
 *         \retval >0 if INFO = k, U(k,k) is exactly zero. The factorization
 *              has been completed, but the factor U is exactly
 *              singular, and division by zero will occur if it is used
 *              to solve a system of equations.
 *
 ******************************************************************************/
int
CORE_cgetf2_nopiv(int M, int N,
                  PLASMA_Complex32_t *A, int LDA)
{
    PLASMA_Complex32_t mzone = (PLASMA_Complex32_t)-1.0;
    PLASMA_Complex32_t alpha;
    float sfmin;
    int i, j, k;
    int info;

    /* Check input arguments */
    info = 0;
    if (M < 0) {
        coreblas_error(1, "Illegal value of M");
        return -1;
    }
    if (N < 0) {
        coreblas_error(2, "Illegal value of N");
        return -2;
    }
    if ((LDA < max(1,M)) && (M > 0)) {
        coreblas_error(5, "Illegal value of LDA");
        return -5;
    }

    /* Quick return */
    if ( (M == 0) || (N == 0) )
        return PLASMA_SUCCESS;

    sfmin = LAPACKE_slamch_work('S');
    k = min(M, N);
    for(i=0 ; i < k; i++) {
        alpha = A[i*LDA+i];
        if ( alpha != (PLASMA_Complex32_t)0.0 ) {
            /* Compute elements J+1:M of J-th column. */
            if (i < M) {
                if ( cabsf(alpha) > sfmin ) {
                    alpha = 1.0 / alpha;
                    cblas_cscal( M-i-1, CBLAS_SADDR(alpha), &(A[i*LDA+i+1]), 1);
                } else {
                    for(j=i+1; j<M; j++)
                        A[LDA*i+j] = A[LDA*i+j] / alpha;
                }
            }
        } else if ( info == 0 ) {
            info = i;
            goto end;
        }

        if ( i < k ) {
            /* Update trailing submatrix */
            cblas_cgeru(CblasColMajor,
                        M-i-1, N-i-1, CBLAS_SADDR(mzone),
                        &A[LDA* i   +i+1], 1,
                        &A[LDA*(i+1)+i  ], LDA,
                        &A[LDA*(i+1)+i+1], LDA);
        }
    }

 end:
    return info;
}


text:
/* matrix/gsl_matrix_double.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_MATRIX_DOUBLE_H__
#define __GSL_MATRIX_DOUBLE_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_vector_double.h>
#include <gsl/gsl_blas_types.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 size1;
  size_t size2;
  size_t tda;
  double * data;
  gsl_block * block;
  int owner;
} gsl_matrix;

typedef struct
{
  gsl_matrix matrix;
} _gsl_matrix_view;

typedef _gsl_matrix_view gsl_matrix_view;

typedef struct
{
  gsl_matrix matrix;
} _gsl_matrix_const_view;

typedef const _gsl_matrix_const_view gsl_matrix_const_view;

/* Allocation */

GSL_FUN gsl_matrix *
gsl_matrix_alloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix *
gsl_matrix_calloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix *
gsl_matrix_alloc_from_block (gsl_block * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_FUN gsl_matrix *
gsl_matrix_alloc_from_matrix (gsl_matrix * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN gsl_vector *
gsl_vector_alloc_row_from_matrix (gsl_matrix * m,
                                        const size_t i);

GSL_FUN gsl_vector *
gsl_vector_alloc_col_from_matrix (gsl_matrix * m,
                                        const size_t j);

GSL_FUN void gsl_matrix_free (gsl_matrix * m);

/* Views */

GSL_FUN _gsl_matrix_view
gsl_matrix_submatrix (gsl_matrix * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_view
gsl_matrix_row (gsl_matrix * m, const size_t i);

GSL_FUN _gsl_vector_view
gsl_matrix_column (gsl_matrix * m, const size_t j);

GSL_FUN _gsl_vector_view
gsl_matrix_diagonal (gsl_matrix * m);

GSL_FUN _gsl_vector_view
gsl_matrix_subdiagonal (gsl_matrix * m, const size_t k);

GSL_FUN _gsl_vector_view
gsl_matrix_superdiagonal (gsl_matrix * m, const size_t k);

GSL_FUN _gsl_vector_view
gsl_matrix_subrow (gsl_matrix * m, const size_t i,
                         const size_t offset, const size_t n);

GSL_FUN _gsl_vector_view
gsl_matrix_subcolumn (gsl_matrix * m, const size_t j,
                            const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_view
gsl_matrix_view_array (double * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUN _gsl_matrix_view
gsl_matrix_view_array_with_tda (double * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUN _gsl_matrix_view
gsl_matrix_view_vector (gsl_vector * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUN _gsl_matrix_view
gsl_matrix_view_vector_with_tda (gsl_vector * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUN _gsl_matrix_const_view
gsl_matrix_const_submatrix (const gsl_matrix * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_const_view
gsl_matrix_const_row (const gsl_matrix * m,
                            const size_t i);

GSL_FUN _gsl_vector_const_view
gsl_matrix_const_column (const gsl_matrix * m,
                               const size_t j);

GSL_FUN _gsl_vector_const_view
gsl_matrix_const_diagonal (const gsl_matrix * m);

GSL_FUN _gsl_vector_const_view
gsl_matrix_const_subdiagonal (const gsl_matrix * m,
                                    const size_t k);

GSL_FUN _gsl_vector_const_view
gsl_matrix_const_superdiagonal (const gsl_matrix * m,
                                      const size_t k);

GSL_FUN _gsl_vector_const_view
gsl_matrix_const_subrow (const gsl_matrix * m, const size_t i,
                               const size_t offset, const size_t n);

GSL_FUN _gsl_vector_const_view
gsl_matrix_const_subcolumn (const gsl_matrix * m, const size_t j,
                                  const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_const_view
gsl_matrix_const_view_array (const double * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUN _gsl_matrix_const_view
gsl_matrix_const_view_array_with_tda (const double * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUN _gsl_matrix_const_view
gsl_matrix_const_view_vector (const gsl_vector * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN _gsl_matrix_const_view
gsl_matrix_const_view_vector_with_tda (const gsl_vector * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUN void gsl_matrix_set_zero (gsl_matrix * m);
GSL_FUN void gsl_matrix_set_identity (gsl_matrix * m);
GSL_FUN void gsl_matrix_set_all (gsl_matrix * m, double x);

GSL_FUN int gsl_matrix_fread (FILE * stream, gsl_matrix * m) ;
GSL_FUN int gsl_matrix_fwrite (FILE * stream, const gsl_matrix * m) ;
GSL_FUN int gsl_matrix_fscanf (FILE * stream, gsl_matrix * m);
GSL_FUN int gsl_matrix_fprintf (FILE * stream, const gsl_matrix * m, const char * format);

GSL_FUN int gsl_matrix_memcpy(gsl_matrix * dest, const gsl_matrix * src);
GSL_FUN int gsl_matrix_swap(gsl_matrix * m1, gsl_matrix * m2);
GSL_FUN int gsl_matrix_tricpy(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix * dest, const gsl_matrix * src);

GSL_FUN int gsl_matrix_swap_rows(gsl_matrix * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_swap_columns(gsl_matrix * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_swap_rowcol(gsl_matrix * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_transpose (gsl_matrix * m);
GSL_FUN int gsl_matrix_transpose_memcpy (gsl_matrix * dest, const gsl_matrix * src);
GSL_FUN int gsl_matrix_transpose_tricpy (CBLAS_UPLO_t Uplo_src, CBLAS_DIAG_t Diag, gsl_matrix * dest, const gsl_matrix * src);

GSL_FUN double gsl_matrix_max (const gsl_matrix * m);
GSL_FUN double gsl_matrix_min (const gsl_matrix * m);
GSL_FUN void gsl_matrix_minmax (const gsl_matrix * m, double * min_out, double * max_out);

GSL_FUN void gsl_matrix_max_index (const gsl_matrix * m, size_t * imax, size_t *jmax);
GSL_FUN void gsl_matrix_min_index (const gsl_matrix * m, size_t * imin, size_t *jmin);
GSL_FUN void gsl_matrix_minmax_index (const gsl_matrix * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_FUN int gsl_matrix_equal (const gsl_matrix * a, const gsl_matrix * b);

GSL_FUN int gsl_matrix_isnull (const gsl_matrix * m);
GSL_FUN int gsl_matrix_ispos (const gsl_matrix * m);
GSL_FUN int gsl_matrix_isneg (const gsl_matrix * m);
GSL_FUN int gsl_matrix_isnonneg (const gsl_matrix * m);

GSL_FUN int gsl_matrix_add (gsl_matrix * a, const gsl_matrix * b);
GSL_FUN int gsl_matrix_sub (gsl_matrix * a, const gsl_matrix * b);
GSL_FUN int gsl_matrix_mul_elements (gsl_matrix * a, const gsl_matrix * b);
GSL_FUN int gsl_matrix_div_elements (gsl_matrix * a, const gsl_matrix * b);
GSL_FUN int gsl_matrix_scale (gsl_matrix * a, const double x);
GSL_FUN int gsl_matrix_scale_rows (gsl_matrix * a, const gsl_vector * x);
GSL_FUN int gsl_matrix_scale_columns (gsl_matrix * a, const gsl_vector * x);
GSL_FUN int gsl_matrix_add_constant (gsl_matrix * a, const double x);
GSL_FUN int gsl_matrix_add_diagonal (gsl_matrix * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUN int gsl_matrix_get_row(gsl_vector * v, const gsl_matrix * m, const size_t i);
GSL_FUN int gsl_matrix_get_col(gsl_vector * v, const gsl_matrix * m, const size_t j);
GSL_FUN int gsl_matrix_set_row(gsl_matrix * m, const size_t i, const gsl_vector * v);
GSL_FUN int gsl_matrix_set_col(gsl_matrix * m, const size_t j, const gsl_vector * v);
/***********************************************************************/

/* inline functions if you are using GCC */

GSL_FUN INLINE_DECL double   gsl_matrix_get(const gsl_matrix * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL void    gsl_matrix_set(gsl_matrix * m, const size_t i, const size_t j, const double x);
GSL_FUN INLINE_DECL double * gsl_matrix_ptr(gsl_matrix * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL const double * gsl_matrix_const_ptr(const gsl_matrix * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE
INLINE_FUN
double
gsl_matrix_get(const gsl_matrix * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
        }
    }
#endif
  return m->data[i * m->tda + j] ;
}

INLINE_FUN
void
gsl_matrix_set(gsl_matrix * m, const size_t i, const size_t j, const double x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  m->data[i * m->tda + j] = x ;
}

INLINE_FUN
double *
gsl_matrix_ptr(gsl_matrix * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (double *) (m->data + (i * m->tda + j)) ;
}

INLINE_FUN
const double *
gsl_matrix_const_ptr(const gsl_matrix * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const double *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_DOUBLE_H__ */


text:
/* specfunc/gsl_sf_legendre.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2004 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_LEGENDRE_H__
#define __GSL_SF_LEGENDRE_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

/* P_l(x)   l >= 0; |x| <= 1
 *
 * exceptions: GSL_EDOM
 */
int     gsl_sf_legendre_Pl_e(const int l, const double x, gsl_sf_result * result);
double  gsl_sf_legendre_Pl(const int l, const double x);

/* P_l(x) for l=0,...,lmax; |x| <= 1
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_legendre_Pl_array(
  const int lmax, const double x,
  double * result_array
  );

/* P_l(x) and P_l'(x) for l=0,...,lmax; |x| <= 1
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_legendre_Pl_deriv_array(
  const int lmax, const double x,
  double * result_array,
  double * result_deriv_array
  );

/* P_l(x), l=1,2,3
 *
 * exceptions: none
 */
int gsl_sf_legendre_P1_e(double x, gsl_sf_result * result);
int gsl_sf_legendre_P2_e(double x, gsl_sf_result * result);
int gsl_sf_legendre_P3_e(double x, gsl_sf_result * result);
double gsl_sf_legendre_P1(const double x);
double gsl_sf_legendre_P2(const double x);
double gsl_sf_legendre_P3(const double x);

/* Q_0(x), x > -1, x != 1
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_legendre_Q0_e(const double x, gsl_sf_result * result);
double gsl_sf_legendre_Q0(const double x);

/* Q_1(x), x > -1, x != 1
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_legendre_Q1_e(const double x, gsl_sf_result * result);
double gsl_sf_legendre_Q1(const double x);

/* Q_l(x), x > -1, x != 1, l >= 0
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_legendre_Ql_e(const int l, const double x, gsl_sf_result * result);
double gsl_sf_legendre_Ql(const int l, const double x);

/* P_l^m(x)  m >= 0; l >= m; |x| <= 1.0
 *
 * Note that this function grows combinatorially with l.
 * Therefore we can easily generate an overflow for l larger
 * than about 150.
 *
 * There is no trouble for small m, but when m and l are both large,
 * then there will be trouble. Rather than allow overflows, these
 * functions refuse to calculate when they can sense that l and m are
 * too big.
 *
 * If you really want to calculate a spherical harmonic, then DO NOT
 * use this. Instead use legendre_sphPlm() below, which  uses a similar
 * recursion, but with the normalized functions.
 *
 * exceptions: GSL_EDOM, GSL_EOVRFLW
 */
int     gsl_sf_legendre_Plm_e(const int l, const int m, const double x, gsl_sf_result * result);
double  gsl_sf_legendre_Plm(const int l, const int m, const double x);

/* P_l^m(x)  m >= 0; l >= m; |x| <= 1.0
 * l=|m|,...,lmax
 *
 * exceptions: GSL_EDOM, GSL_EOVRFLW
 */
int gsl_sf_legendre_Plm_array(
  const int lmax, const int m, const double x,
  double * result_array
  );

/* P_l^m(x)  and d(P_l^m(x))/dx;  m >= 0; lmax >= m; |x| <= 1.0
 * l=|m|,...,lmax
 *
 * exceptions: GSL_EDOM, GSL_EOVRFLW
 */
int gsl_sf_legendre_Plm_deriv_array(
  const int lmax, const int m, const double x,
  double * result_array,
  double * result_deriv_array
  );

/* P_l^m(x), normalized properly for use in spherical harmonics
 * m >= 0; l >= m; |x| <= 1.0
 *
 * There is no overflow problem, as there is for the
 * standard normalization of P_l^m(x).
 *
 * Specifically, it returns:
 *
 *        sqrt((2l+1)/(4pi)) sqrt((l-m)!/(l+m)!) P_l^m(x)
 *
 * exceptions: GSL_EDOM
 */
int     gsl_sf_legendre_sphPlm_e(const int l, int m, const double x, gsl_sf_result * result);
double  gsl_sf_legendre_sphPlm(const int l, const int m, const double x);

/* sphPlm(l,m,x) values
 * m >= 0; l >= m; |x| <= 1.0
 * l=|m|,...,lmax
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_legendre_sphPlm_array(
  const int lmax, int m, const double x,
  double * result_array
  );

/* sphPlm(l,m,x) and d(sphPlm(l,m,x))/dx values
 * m >= 0; l >= m; |x| <= 1.0
 * l=|m|,...,lmax
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_legendre_sphPlm_deriv_array(
  const int lmax, const int m, const double x,
  double * result_array,
  double * result_deriv_array
  );

/* size of result_array[] needed for the array versions of Plm
 * (lmax - m + 1)
 */
int gsl_sf_legendre_array_size(const int lmax, const int m);

/* Irregular Spherical Conical Function
 * P^{1/2}_{-1/2 + I lambda}(x)
 *
 * x > -1.0
 * exceptions: GSL_EDOM
 */
int gsl_sf_conicalP_half_e(const double lambda, const double x, gsl_sf_result * result);
double gsl_sf_conicalP_half(const double lambda, const double x);

/* Regular Spherical Conical Function
 * P^{-1/2}_{-1/2 + I lambda}(x)
 *
 * x > -1.0
 * exceptions: GSL_EDOM
 */
int gsl_sf_conicalP_mhalf_e(const double lambda, const double x, gsl_sf_result * result);
double gsl_sf_conicalP_mhalf(const double lambda, const double x);

/* Conical Function
 * P^{0}_{-1/2 + I lambda}(x)
 *
 * x > -1.0
 * exceptions: GSL_EDOM
 */
int gsl_sf_conicalP_0_e(const double lambda, const double x, gsl_sf_result * result);
double gsl_sf_conicalP_0(const double lambda, const double x);

/* Conical Function
 * P^{1}_{-1/2 + I lambda}(x)
 *
 * x > -1.0
 * exceptions: GSL_EDOM
 */
int gsl_sf_conicalP_1_e(const double lambda, const double x, gsl_sf_result * result);
double gsl_sf_conicalP_1(const double lambda, const double x);

/* Regular Spherical Conical Function
 * P^{-1/2-l}_{-1/2 + I lambda}(x)
 *
 * x > -1.0, l >= -1
 * exceptions: GSL_EDOM
 */
int gsl_sf_conicalP_sph_reg_e(const int l, const double lambda, const double x, gsl_sf_result * result);
double gsl_sf_conicalP_sph_reg(const int l, const double lambda, const double x);

/* Regular Cylindrical Conical Function
 * P^{-m}_{-1/2 + I lambda}(x)
 *
 * x > -1.0, m >= -1
 * exceptions: GSL_EDOM
 */
int gsl_sf_conicalP_cyl_reg_e(const int m, const double lambda, const double x, gsl_sf_result * result);
double gsl_sf_conicalP_cyl_reg(const int m, const double lambda, const double x);

/* The following spherical functions are specializations
 * of Legendre functions which give the regular eigenfunctions
 * of the Laplacian on a 3-dimensional hyperbolic space.
 * Of particular interest is the flat limit, which is
 * Flat-Lim := {lambda->Inf, eta->0, lambda*eta fixed}.
 */

/* Zeroth radial eigenfunction of the Laplacian on the
 * 3-dimensional hyperbolic space.
 *
 * legendre_H3d_0(lambda,eta) := sin(lambda*eta)/(lambda*sinh(eta))
 *
 * Normalization:
 * Flat-Lim legendre_H3d_0(lambda,eta) = j_0(lambda*eta)
 *
 * eta >= 0.0
 * exceptions: GSL_EDOM
 */
int gsl_sf_legendre_H3d_0_e(const double lambda, const double eta, gsl_sf_result * result);
double gsl_sf_legendre_H3d_0(const double lambda, const double eta);

/* First radial eigenfunction of the Laplacian on the
 * 3-dimensional hyperbolic space.
 *
 * legendre_H3d_1(lambda,eta) :=
 *    1/sqrt(lambda^2 + 1) sin(lam eta)/(lam sinh(eta))
 *    (coth(eta) - lambda cot(lambda*eta))
 *
 * Normalization:
 * Flat-Lim legendre_H3d_1(lambda,eta) = j_1(lambda*eta)
 *
 * eta >= 0.0
 * exceptions: GSL_EDOM
 */
int gsl_sf_legendre_H3d_1_e(const double lambda, const double eta, gsl_sf_result * result);
double gsl_sf_legendre_H3d_1(const double lambda, const double eta);

/* l'th radial eigenfunction of the Laplacian on the
 * 3-dimensional hyperbolic space.
 *
 * Normalization:
 * Flat-Lim legendre_H3d_l(l,lambda,eta) = j_l(lambda*eta)
 *
 * eta >= 0.0, l >= 0
 * exceptions: GSL_EDOM
 */
int gsl_sf_legendre_H3d_e(const int l, const double lambda, const double eta, gsl_sf_result * result);
double gsl_sf_legendre_H3d(const int l, const double lambda, const double eta);

/* Array of H3d(ell),  0 <= ell <= lmax
 */
int gsl_sf_legendre_H3d_array(const int lmax, const double lambda, const double eta, double * result_array);

/* associated legendre P_{lm} routines */

typedef enum
{
  GSL_SF_LEGENDRE_SCHMIDT,
  GSL_SF_LEGENDRE_SPHARM,
  GSL_SF_LEGENDRE_FULL,
  GSL_SF_LEGENDRE_NONE
} gsl_sf_legendre_t;

int gsl_sf_legendre_array(const gsl_sf_legendre_t norm,
                          const size_t lmax, const double x,
                          double result_array[]);
int gsl_sf_legendre_array_e(const gsl_sf_legendre_t norm,
                            const size_t lmax, const double x,
                            const double csphase,
                            double result_array[]);
int gsl_sf_legendre_deriv_array(const gsl_sf_legendre_t norm,
                                const size_t lmax, const double x,
                                double result_array[],
                                double result_deriv_array[]);
int gsl_sf_legendre_deriv_array_e(const gsl_sf_legendre_t norm,
                                  const size_t lmax, const double x,
                                  const double csphase,
                                  double result_array[],
                                  double result_deriv_array[]);
int gsl_sf_legendre_deriv_alt_array(const gsl_sf_legendre_t norm,
                                    const size_t lmax, const double x,
                                    double result_array[],
                                    double result_deriv_array[]);
int gsl_sf_legendre_deriv_alt_array_e(const gsl_sf_legendre_t norm,
                                      const size_t lmax, const double x,
                                      const double csphase,
                                      double result_array[],
                                      double result_deriv_array[]);
int gsl_sf_legendre_deriv2_array(const gsl_sf_legendre_t norm,
                                 const size_t lmax, const double x,
                                 double result_array[],
                                 double result_deriv_array[],
                                 double result_deriv2_array[]);
int gsl_sf_legendre_deriv2_array_e(const gsl_sf_legendre_t norm,
                                   const size_t lmax, const double x,
                                   const double csphase,
                                   double result_array[],
                                   double result_deriv_array[],
                                   double result_deriv2_array[]);
int gsl_sf_legendre_deriv2_alt_array(const gsl_sf_legendre_t norm,
                                     const size_t lmax, const double x,
                                     double result_array[],
                                     double result_deriv_array[],
                                     double result_deriv2_array[]);
int gsl_sf_legendre_deriv2_alt_array_e(const gsl_sf_legendre_t norm,
                                       const size_t lmax, const double x,
                                       const double csphase,
                                       double result_array[],
                                       double result_deriv_array[],
                                       double result_deriv2_array[]);
size_t gsl_sf_legendre_array_n(const size_t lmax);
size_t gsl_sf_legendre_array_index(const size_t l, const size_t m);
size_t gsl_sf_legendre_nlm(const size_t lmax);

__END_DECLS

#endif /* __GSL_SF_LEGENDRE_H__ */

max_stars_repo_name: snipekill/FPGen

max_issues_repo_name: snipekill/FPGen

max_forks_repo_name: snipekill/FPGen


text:
// Author: Samuel D. Relton
#include <stdio.h>
#include <cblas.h>
#include "sgemm_batched.h"

float A1[] = {
3, 1, 3,
1, 5, 9,
2, 6, 5
};

float B1[] = {
3, 1, 3,
1, 5, 9,
2, 6, 5
};

float C1[] = {
0, 0, 0,
0, 0, 0,
0, 0, 0
};

float A2[] = {
3, 1, 3,
1, 5, 9,
2, 6, 5
};

float B2[] = {
3, 1, 3,
1, 5, 9,
2, 6, 5
};

float C2[] = {
0, 0, 0,
0, 0, 0,
0, 0, 0
};

float ALPHA[] = {1.0, 1.0};
float BETA[] = {0.0, 0.0};

int LDA[] = {3, 3};
int LDB[] = {3, 3};
int LDC[] = {3, 3};

bblas_trans_t  TRANSA[] = {CblasNoTrans, CblasNoTrans};
bblas_trans_t  TRANSB[] = {CblasNoTrans, CblasNoTrans};

int M[] = {3, 3};
int N[] = {3, 3};
int K[] = {3, 3};

int BATCHCOUNT = 2;
bblas_batch_type_t BATCH_TYPE = BBLASFixed;

int INFO[] = {0, 0};

float *arrayA[] = {A1, A2};
float *arrayB[] = {B1, B2};
float *arrayC[] = {C1, C2};

int main()
{
	int i, j, k;

	sgemm_batched(TRANSA, TRANSB, M, N, K, ALPHA, arrayA,
		      LDA, arrayB, LDB, BETA, arrayC, LDC, BATCHCOUNT,
		      BATCH_TYPE, INFO);

	for (k = 0; k < BATCHCOUNT; k++) {
		putchar('\n');
		for (i=0; i<3; ++i) {
			for (j=0; j<3; ++j) {
				printf("%5.1f", arrayC[k][i*3+j]);

			}
			putchar('\n');
		}
	}

    return 0;
}

max_stars_repo_name: NLAFET/BBLAS-ref

max_issues_repo_name: NLAFET/BBLAS-ref

max_forks_repo_name: NLAFET/BBLAS-ref


text:
/*
** smooth variance computations
**
** M. Kuhlmann, MPI-KYB, Sept 2015
*/

#include <viaio/Vlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_combination.h>

VImage variance_smoothing(VImage *src[], float width)
{
  size_t nslices, nrows, ncols;

  /* image dims */
  nslices = VImageNBands(src[0]);
  nrows  = VImageNRows(src[0]);
  ncols  = VImageNColumns(src[0]);

  VImage var = VCreateImage(nslices,nrows,ncols,VFloatRepn);
  VFillImage(var,VAllBands,0);

  /*double *data,int n,double *a,double *v)*/

  int j;
  double ave,var,nx,s,u;

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

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

  /* smooth variance */

  return var;
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


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

#ifndef __GSL_RANDIST_H__
#define __GSL_RANDIST_H__
#include <gsl/gsl_rng.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

unsigned int gsl_ran_bernoulli (const gsl_rng * r, double p);
double gsl_ran_bernoulli_pdf (const unsigned int k, double p);

double gsl_ran_beta (const gsl_rng * r, const double a, const double b);
double gsl_ran_beta_pdf (const double x, const double a, const double b);

unsigned int gsl_ran_binomial (const gsl_rng * r, double p, unsigned int n);
unsigned int gsl_ran_binomial_knuth (const gsl_rng * r, double p, unsigned int n);
unsigned int gsl_ran_binomial_tpe (const gsl_rng * r, double p, unsigned int n);
double gsl_ran_binomial_pdf (const unsigned int k, const double p, const unsigned int n);

double gsl_ran_exponential (const gsl_rng * r, const double mu);
double gsl_ran_exponential_pdf (const double x, const double mu);

double gsl_ran_exppow (const gsl_rng * r, const double a, const double b);
double gsl_ran_exppow_pdf (const double x, const double a, const double b);

double gsl_ran_cauchy (const gsl_rng * r, const double a);
double gsl_ran_cauchy_pdf (const double x, const double a);

double gsl_ran_chisq (const gsl_rng * r, const double nu);
double gsl_ran_chisq_pdf (const double x, const double nu);

void gsl_ran_dirichlet (const gsl_rng * r, const size_t K, const double alpha[], double theta[]);
double gsl_ran_dirichlet_pdf (const size_t K, const double alpha[], const double theta[]);
double gsl_ran_dirichlet_lnpdf (const size_t K, const double alpha[], const double theta[]);

double gsl_ran_erlang (const gsl_rng * r, const double a, const double n);
double gsl_ran_erlang_pdf (const double x, const double a, const double n);

double gsl_ran_fdist (const gsl_rng * r, const double nu1, const double nu2);
double gsl_ran_fdist_pdf (const double x, const double nu1, const double nu2);

double gsl_ran_flat (const gsl_rng * r, const double a, const double b);
double gsl_ran_flat_pdf (double x, const double a, const double b);

double gsl_ran_gamma (const gsl_rng * r, const double a, const double b);
double gsl_ran_gamma_int (const gsl_rng * r, const unsigned int a);
double gsl_ran_gamma_pdf (const double x, const double a, const double b);
double gsl_ran_gamma_mt (const gsl_rng * r, const double a, const double b);
double gsl_ran_gamma_knuth (const gsl_rng * r, const double a, const double b);

double gsl_ran_gaussian (const gsl_rng * r, const double sigma);
double gsl_ran_gaussian_ratio_method (const gsl_rng * r, const double sigma);
double gsl_ran_gaussian_ziggurat (const gsl_rng * r, const double sigma);
double gsl_ran_gaussian_pdf (const double x, const double sigma);

double gsl_ran_ugaussian (const gsl_rng * r);
double gsl_ran_ugaussian_ratio_method (const gsl_rng * r);
double gsl_ran_ugaussian_pdf (const double x);

double gsl_ran_gaussian_tail (const gsl_rng * r, const double a, const double sigma);
double gsl_ran_gaussian_tail_pdf (const double x, const double a, const double sigma);

double gsl_ran_ugaussian_tail (const gsl_rng * r, const double a);
double gsl_ran_ugaussian_tail_pdf (const double x, const double a);

void gsl_ran_bivariate_gaussian (const gsl_rng * r, double sigma_x, double sigma_y, double rho, double *x, double *y);
double gsl_ran_bivariate_gaussian_pdf (const double x, const double y, const double sigma_x, const double sigma_y, const double rho);

double gsl_ran_landau (const gsl_rng * r);
double gsl_ran_landau_pdf (const double x);

unsigned int gsl_ran_geometric (const gsl_rng * r, const double p);
double gsl_ran_geometric_pdf (const unsigned int k, const double p);

unsigned int gsl_ran_hypergeometric (const gsl_rng * r, unsigned int n1, unsigned int n2, unsigned int t);
double gsl_ran_hypergeometric_pdf (const unsigned int k, const unsigned int n1, const unsigned int n2, unsigned int t);

double gsl_ran_gumbel1 (const gsl_rng * r, const double a, const double b);
double gsl_ran_gumbel1_pdf (const double x, const double a, const double b);

double gsl_ran_gumbel2 (const gsl_rng * r, const double a, const double b);
double gsl_ran_gumbel2_pdf (const double x, const double a, const double b);

double gsl_ran_logistic (const gsl_rng * r, const double a);
double gsl_ran_logistic_pdf (const double x, const double a);

double gsl_ran_lognormal (const gsl_rng * r, const double zeta, const double sigma);
double gsl_ran_lognormal_pdf (const double x, const double zeta, const double sigma);

unsigned int gsl_ran_logarithmic (const gsl_rng * r, const double p);
double gsl_ran_logarithmic_pdf (const unsigned int k, const double p);

void gsl_ran_multinomial (const gsl_rng * r, const size_t K,
                          const unsigned int N, const double p[],
                          unsigned int n[] );
double gsl_ran_multinomial_pdf (const size_t K,
                                const double p[], const unsigned int n[] );
double gsl_ran_multinomial_lnpdf (const size_t K,
                           const double p[], const unsigned int n[] );

unsigned int gsl_ran_negative_binomial (const gsl_rng * r, double p, double n);
double gsl_ran_negative_binomial_pdf (const unsigned int k, const double p, double n);

unsigned int gsl_ran_pascal (const gsl_rng * r, double p, unsigned int n);
double gsl_ran_pascal_pdf (const unsigned int k, const double p, unsigned int n);

double gsl_ran_pareto (const gsl_rng * r, double a, const double b);
double gsl_ran_pareto_pdf (const double x, const double a, const double b);

unsigned int gsl_ran_poisson (const gsl_rng * r, double mu);
void gsl_ran_poisson_array (const gsl_rng * r, size_t n, unsigned int array[],
                            double mu);
double gsl_ran_poisson_pdf (const unsigned int k, const double mu);

double gsl_ran_rayleigh (const gsl_rng * r, const double sigma);
double gsl_ran_rayleigh_pdf (const double x, const double sigma);

double gsl_ran_rayleigh_tail (const gsl_rng * r, const double a, const double sigma);
double gsl_ran_rayleigh_tail_pdf (const double x, const double a, const double sigma);

double gsl_ran_tdist (const gsl_rng * r, const double nu);
double gsl_ran_tdist_pdf (const double x, const double nu);

double gsl_ran_laplace (const gsl_rng * r, const double a);
double gsl_ran_laplace_pdf (const double x, const double a);

double gsl_ran_levy (const gsl_rng * r, const double c, const double alpha);
double gsl_ran_levy_skew (const gsl_rng * r, const double c, const double alpha, const double beta);

double gsl_ran_weibull (const gsl_rng * r, const double a, const double b);
double gsl_ran_weibull_pdf (const double x, const double a, const double b);

void gsl_ran_dir_2d (const gsl_rng * r, double * x, double * y);
void gsl_ran_dir_2d_trig_method (const gsl_rng * r, double * x, double * y);
void gsl_ran_dir_3d (const gsl_rng * r, double * x, double * y, double * z);
void gsl_ran_dir_nd (const gsl_rng * r, size_t n, double * x);

void gsl_ran_shuffle (const gsl_rng * r, void * base, size_t nmembm, size_t size);
int gsl_ran_choose (const gsl_rng * r, void * dest, size_t k, void * src, size_t n, size_t size) ;
void gsl_ran_sample (const gsl_rng * r, void * dest, size_t k, void * src, size_t n, size_t size) ;

typedef struct {                /* struct for Walker algorithm */
    size_t K;
    size_t *A;
    double *F;
} gsl_ran_discrete_t;

gsl_ran_discrete_t * gsl_ran_discrete_preproc (size_t K, const double *P);
void gsl_ran_discrete_free(gsl_ran_discrete_t *g);
size_t gsl_ran_discrete (const gsl_rng *r, const gsl_ran_discrete_t *g);
double gsl_ran_discrete_pdf (size_t k, const gsl_ran_discrete_t *g);

__END_DECLS

#endif /* __GSL_RANDIST_H__ */

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


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

#include "../../runtime/pcp.h"
#include "driver.h"

static double random_double(void)
{
    return (double) rand() / RAND_MAX;
}

static void usage(char *prog)
{
    printf("Usage: %s n m b p q\n", prog);
    printf("\n");
    printf("n: The size of the triangular matrix\n");
    printf("m: The number of right-hand sides\n");
    printf("b: The tile size\n");
    printf("p: The number of threads (or 0 for one per core)\n");
    printf("q: The size of the reserved set (or 0 to use the regular mode)\n");
    exit(EXIT_FAILURE);
}

int main(int argc, char *argv[])
{
    // Verify the number of command line options.
    if (argc != 6 && argc != 7) {
        usage(argv[0]);
    }

    // Parse command line.
    int n = atoi(argv[1]);
    int m = atoi(argv[2]);
    int blksz = atoi(argv[3]);
    int num_threads = atoi(argv[4]);
    int reserved_set_size = atoi(argv[5]);
    int verify = 1;
    if (argc == 7) {
        verify = 0;
    }

    // Verify options.
    if (n < 1 ||
        m < 1 ||
        blksz < 1 ||
        num_threads < 0 ||
        reserved_set_size < 0 ||
        (reserved_set_size >= num_threads && num_threads > 0)) {
        usage(argv[0]);
    }

    // Start the runtime system.
    pcp_start(num_threads);

    // Set the execution mode.
    if (reserved_set_size == 0) {
        pcp_set_mode(PCP_REGULAR);
    } else {
        pcp_set_mode(PCP_FIXED);
        pcp_set_reserved_set_size(reserved_set_size);
    }

    // Allocate matrices.
    int ldL, ldX, ldB;
    double *L, *X, *B;
    ldL = ldX = ldB = n;
    L = (double*) malloc(sizeof(double) * ldL * n);
    X = (double*) malloc(sizeof(double) * ldX * m);
    B = (double*) malloc(sizeof(double) * ldB * m);

#define L(i,j) L[(i) + (j) * ldL]
#define X(i,j) X[(i) + (j) * ldX]
#define B(i,j) B[(i) + (j) * ldB]

    // Generate L.
    printf("Generating lower triangular matrix L...\n");
    for (int j = 0; j < n; ++j) {
        for (int i = 0; i < j; ++i) {
            L(i,j) = 0;
        }
        for (int i = j; i < n; ++i) {
            L(i,j) = random_double();
        }
        L(j,j) = n;
    }

    // Generate B.
    printf("Generating right-hand sides B...\n");
    for (int j = 0; j < m; ++j) {
        for (int i = 0; i < n; ++i) {
            B(i,j) = random_double();
        }
    }

    // Copy B to X.
    printf("Copying B to X...\n");
    for (int j = 0; j < m; ++j) {
        for (int i = 0; i < n; ++i) {
            X(i,j) = B(i,j);
        }
    }

    // Call the driver.
    printf("Calling the driver routine...\n");
    parallel_forward_substitution(n, m, blksz, L, ldL, X, ldX);

    if (verify) {
        // Verify the solution.
        printf("Verifying the solution... ");
        fflush(stdout);
        cblas_dtrmm(CblasColMajor, CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit,
                    n, m,
                    1.0, L, ldL,
                    X, ldX);
        double err = 0;
        for (int j = 0; j < m; ++j) {
            for (int i = 0; i < n; ++i) {
                double e = X(i,j) - B(i,j);
                err += e * e;
            }
        }
        err = sqrt(err);
        if (err / n < 1e-14) {
            printf("PASSED\n");
        } else {
            printf("FAILED\n");
        }
    } else {
        printf("Verification disabled\n");
    }

    // Save trace.
    printf("Saving the trace...\n");
    pcp_view_trace_tikz();

    // Print statistics.
    printf("Printing statistics...\n");
    pcp_view_statistics_stdout();

    // Stop the runtime system.
    pcp_stop();

#undef L
#undef X
#undef B

    // Clean up.
    free(L);
    free(X);
    free(B);

    return EXIT_SUCCESS;
}

max_stars_repo_name: NLAFET/pcp-runtime

max_issues_repo_name: NLAFET/pcp-runtime

max_forks_repo_name: NLAFET/pcp-runtime


text:
#ifndef UTIL_GSL_CONVERTER_H
#define UTIL_GSL_CONVERTER_H

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

namespace gsl {

inline
gsl_vector* convert_vec(const arma::vec & vector) {

  gsl_vector * gsl_vec = gsl_vector_alloc(vector.n_elem);

#pragma omp parallel for
  for(arma::uword i=0; i<vector.n_elem; i++) {
    gsl_vector_set(gsl_vec, i, vector(i));
  }

  return gsl_vec;
}

inline
arma::vec convert_vec(const gsl_vector * gsl_vec) {

  arma::vec arma_vec(gsl_vec->size);

#pragma omp parallel for
  for(arma::uword i=0; i<arma_vec.n_elem; i++) {
    arma_vec(i) = gsl_vector_get(gsl_vec, i);
  }

  return arma_vec;
}

inline
gsl_matrix* convert_mat(const arma::mat & mat) {

  gsl_matrix * gsl_mat = gsl_matrix_alloc(mat.n_rows, mat.n_cols);

#pragma omp parallel for
  for(arma::uword i=0; i<mat.n_rows; i++) {
    for(arma::uword j=0; j<mat.n_cols; j++) {
     gsl_matrix_set(gsl_mat, i, j, mat(i,j));
    }
  }

  return gsl_mat;
}

inline
arma::mat convert_mat(const gsl_matrix * gsl_mat) {

  arma::mat arma_mat(gsl_mat->size1, gsl_mat->size2);

#pragma omp parallel for
  for(arma::uword i=0; i<arma_mat.n_rows; i++) {
    for(arma::uword j=0; j<arma_mat.n_cols; j++) {
      arma_mat(i,j) = gsl_matrix_get(gsl_mat, i, j);
    }
  }

  return arma_mat;
}

}

#endif //UTIL_GSL_CONVERTER_H

max_stars_repo_name: Walter-Feng/Quartz

max_issues_repo_name: Walter-Feng/Quartz

max_forks_repo_name: Walter-Feng/Quartz


text:
/**
 *  Copyright (C) 2021 FISCO BCOS.
 *  SPDX-License-Identifier: Apache-2.0
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 * @brief interface for Ledger
 * @file LedgerInterface.h
 * @author: kyonRay
 * @date: 2021-04-07
 */

#pragma once

#include "../../interfaces/crypto/CommonType.h"
#include "../../interfaces/protocol/Block.h"
#include "../../interfaces/protocol/BlockHeader.h"
#include "../../interfaces/protocol/Transaction.h"
#include "../../interfaces/protocol/TransactionReceipt.h"
#include "../../interfaces/storage/StorageInterface.h"
#include "LedgerConfig.h"
#include "LedgerTypeDef.h"
#include "bcos-utilities/Error.h"
#include <gsl/span>
#include <map>

namespace bcos::ledger
{
class LedgerInterface
{
public:
    using Ptr = std::shared_ptr<LedgerInterface>;
    LedgerInterface() = default;
    virtual ~LedgerInterface() {}

    /**
     * @brief async prewrite a block in scheduler module
     * @param block the block to commit
     * @param callback trigger this callback when write is finished
     */
    virtual void asyncPrewriteBlock(bcos::storage::StorageInterface::Ptr storage,
        bcos::protocol::Block::ConstPtr block, std::function<void(Error::Ptr&&)> callback) = 0;

    /**
     * @brief async store txs in block when tx pool verify
     * @param _txToStore tx bytes data list
     * @param _txHashList tx hash list
     * @param _onTxsStored callback
     */
    virtual void asyncStoreTransactions(std::shared_ptr<std::vector<bytesConstPtr>> _txToStore,
        crypto::HashListPtr _txHashList, std::function<void(Error::Ptr)> _onTxStored) = 0;

    /**
     * @brief async get block by blockNumber
     * @param _blockNumber number of block
     * @param _blockFlag flag bit of what the block be callback contains,
     *                   you can checkout all flags in LedgerTypeDef.h
     * @param _onGetBlock
     *
     * @example
     * asyncGetBlockDataByNumber(10, HEADER|TRANSACTIONS, [](error, block){ doSomething(); });
     */
    virtual void asyncGetBlockDataByNumber(protocol::BlockNumber _blockNumber, int32_t _blockFlag,
        std::function<void(Error::Ptr, protocol::Block::Ptr)> _onGetBlock) = 0;

    /**
     * @brief async get latest block number
     * @param _onGetBlock
     */
    virtual void asyncGetBlockNumber(
        std::function<void(Error::Ptr, protocol::BlockNumber)> _onGetBlock) = 0;

    /**
     * @brief async get block hash by block number
     * @param _blockNumber the number of block to get
     * @param _onGetBlock
     */
    virtual void asyncGetBlockHashByNumber(protocol::BlockNumber _blockNumber,
        std::function<void(Error::Ptr, crypto::HashType)> _onGetBlock) = 0;

    /**
     * @brief async get block number by block hash
     * @param _blockHash the hash of block to get
     * @param _onGetBlock
     */
    virtual void asyncGetBlockNumberByHash(crypto::HashType const& _blockHash,
        std::function<void(Error::Ptr, protocol::BlockNumber)> _onGetBlock) = 0;

    /**
     * @brief async get a batch of transaction by transaction hash list
     * @param _txHashList transaction hash list, hash should be hex
     * @param _withProof if true then it will callback MerkleProofPtr map in _onGetTx
     *                   if false then MerkleProofPtr map will be nullptr
     * @param _onGetTx return <error, [tx data in bytes], map<txHash, merkleProof>
     */
    virtual void asyncGetBatchTxsByHashList(crypto::HashListPtr _txHashList, bool _withProof,
        std::function<void(Error::Ptr, bcos::protocol::TransactionsPtr,
            std::shared_ptr<std::map<std::string, MerkleProofPtr>>)>
            _onGetTx) = 0;

    /**
     * @brief async get a transaction receipt by tx hash
     * @param _txHash hash of transaction
     * @param _withProof if true then it will callback MerkleProofPtr in _onGetTx
     *                   if false then MerkleProofPtr will be nullptr
     * @param _onGetTx
     */
    virtual void asyncGetTransactionReceiptByHash(crypto::HashType const& _txHash, bool _withProof,
        std::function<void(Error::Ptr, protocol::TransactionReceipt::ConstPtr, MerkleProofPtr)>
            _onGetTx) = 0;

    /**
     * @brief async get total transaction count and latest block number
     * @param _callback callback totalTxCount, totalFailedTxCount, and latest block number
     */
    virtual void asyncGetTotalTransactionCount(std::function<void(Error::Ptr, int64_t _totalTxCount,
            int64_t _failedTxCount, protocol::BlockNumber _latestBlockNumber)>
            _callback) = 0;

    /**
     * @brief async get system config by table key
     * @param _key the key of row, you can checkout all key in LedgerTypeDef.h
     * @param _onGetConfig callback when get config, <value, latest block number>
     */
    virtual void asyncGetSystemConfigByKey(std::string const& _key,
        std::function<void(Error::Ptr, std::string, protocol::BlockNumber)> _onGetConfig) = 0;

    /**
     * @brief async get node list by type, can be sealer or observer
     * @param _type the type of node, CONSENSUS_SEALER or CONSENSUS_OBSERVER
     * @param _onGetConfig
     */
    virtual void asyncGetNodeListByType(std::string const& _type,
        std::function<void(Error::Ptr, consensus::ConsensusNodeListPtr)> _onGetConfig) = 0;

    /**
     * @brief async get a batch of nonce lists in blocks
     * @param _startNumber start block number
     * @param _offset batch offset, if batch is 0, then callback nonce list in start block number;
     * if (_startNumber + _offset) > latest block number, then callback nonce lists in
     * [_startNumber, latest number]
     * @param _onGetList
     */
    virtual void asyncGetNonceList(protocol::BlockNumber _startNumber, int64_t _offset,
        std::function<void(
            Error::Ptr, std::shared_ptr<std::map<protocol::BlockNumber, protocol::NonceListPtr>>)>
            _onGetList) = 0;
};
}  // namespace bcos::ledger