text:
/*  @file mvn.c
 *  @title multivariate normal random variables
 *  @author Carl Boettiger, <cboettig@gmail.com>
 *
 *  Based on the R function rmvnorm, from the mvtnorm package
 *  by Friedrich Leisch and Fabian Scheipl, implemented
 *  using the GSL libraries
 */

#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>

int mvn(gsl_rng * rng, const gsl_vector * mean, gsl_matrix * covar, gsl_vector * ANS)
{
	int i;
	size_t n = mean->size;

	/* Calculate eigenvalues and eigenvectors of covar matrix */
	gsl_vector *eval = gsl_vector_alloc (n);
	gsl_matrix *evec = gsl_matrix_alloc (n, n);
	gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (n);
	gsl_eigen_symmv (covar, eval, evec, w);
	gsl_eigen_symmv_free (w);
//	gsl_eigen_symmv_sort (eval, evec, GSL_EIGEN_SORT_ABS_DESC);

	/* Setup for: evec * matrix(diag(eval)) * transpose(evec)  */
	gsl_matrix *eval_mx = gsl_matrix_calloc (n, n);
	gsl_matrix * x_M = gsl_matrix_alloc (n,n);
	gsl_matrix * x_M_x = gsl_matrix_alloc (n,n);

	gsl_vector_view diagonal = gsl_matrix_diagonal(eval_mx);
	gsl_vector_memcpy(&diagonal.vector, eval);
	for(i=0;i<n;i++)
	{
		gsl_vector_set( &diagonal.vector,
						i,
						sqrt( gsl_vector_get(&diagonal.vector, i) )
					  );
	}

	/* evec * matrix(diag(eval)) * transpose(evec)  */
//	gsl_blas_dsymm (CblasLeft, CblasUpper,
//					1.0, evec, eval_mx, 0.0, x_M);

	gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
					1.0, evec, eval_mx, 0.0, x_M);
	gsl_blas_dgemm (CblasNoTrans, CblasTrans,
					1.0, x_M, evec, 0.0, x_M_x);

	gsl_matrix_free(x_M);
	gsl_matrix_free(eval_mx);
	gsl_matrix_free(evec);
	gsl_vector_free(eval);

	gsl_vector * rnorms = gsl_vector_alloc(n);
	for(i=0;i<n;i++)
	{
		gsl_vector_set
			( rnorms, i,
			  gsl_ran_gaussian_ziggurat(rng, 1)
			);
	}

	gsl_blas_dgemv( CblasTrans, 1.0, x_M_x, rnorms, 0, ANS);
	gsl_vector_add(ANS, mean);
	gsl_matrix_free(x_M_x);
  gsl_vector_free(rnorms);

	return 0;
	/* answer provided through pass by reference */
}

/*
int main()
{
	const int n_tips = 2;
	gsl_rng * rng = gsl_rng_alloc(gsl_rng_default);
	gsl_vector * simdata = gsl_vector_alloc(n_tips);
	gsl_matrix * I = gsl_matrix_alloc(n_tips, n_tips);
	gsl_matrix_set_identity(I);
	gsl_vector * A = gsl_vector_calloc(n_tips);
	gsl_vector_set_all(A, 5.);
	mvn(rng, A, I, simdata);
	gsl_vector_fprintf(stdout, A, "%g");
	return 0;
}
*/

/*
#define N 3
int main(void)
{

	double mean[] = {1,2,5};
	double sigma[] = {1, 0, 0,
					  0, 1, .9,
					  0, .9, 1};
	double out[N] = {0,0,0};

	gsl_vector_view mean_view = gsl_vector_view_array(mean, N);
	gsl_matrix_view sigma_view = gsl_matrix_view_array(sigma, N,N);
	gsl_vector_view out_view = gsl_vector_view_array(out, N);

	gsl_rng * rng = gsl_rng_alloc(gsl_rng_default);

	int i;
	for(i=0;i<10000;i++) {
		mvn(rng, &mean_view.vector, &sigma_view.matrix, &out_view.vector);
		printf("%g %g %g\n", out[0], out[1], out[2]);
	}
	return 0;
}
*/

max_stars_repo_name: cboettig/wrightscape

max_issues_repo_name: cboettig/wrightscape

max_forks_repo_name: cboettig/wrightscape


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

/* Author:  G. Jungman, Modified: B. Gough. */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>

/* The eigen_sort below is not very good, but it is simple and
 * self-contained. We can always implement an improved sort later.  */

int
gsl_eigen_symmv_sort (gsl_vector * eval, gsl_matrix * evec,
                      gsl_eigen_sort_t sort_type)
{
  if (evec->size1 != evec->size2)
    {
      GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
    }
  else if (eval->size != evec->size1)
    {
      GSL_ERROR ("eigenvalues must match eigenvector matrix", GSL_EBADLEN);
    }
  else
    {
      const size_t N = eval->size;
      size_t i;

      for (i = 0; i < N - 1; i++)
        {
          size_t j;
          size_t k = i;

          double ek = gsl_vector_get (eval, i);

          /* search for something to swap */
          for (j = i + 1; j < N; j++)
            {
              int test;
              const double ej = gsl_vector_get (eval, j);

              switch (sort_type)
                {
                case GSL_EIGEN_SORT_VAL_ASC:
                  test = (ej < ek);
                  break;
                case GSL_EIGEN_SORT_VAL_DESC:
                  test = (ej > ek);
                  break;
                case GSL_EIGEN_SORT_ABS_ASC:
                  test = (fabs (ej) < fabs (ek));
                  break;
                case GSL_EIGEN_SORT_ABS_DESC:
                  test = (fabs (ej) > fabs (ek));
                  break;
                default:
                  GSL_ERROR ("unrecognized sort type", GSL_EINVAL);
                }

              if (test)
                {
                  k = j;
                  ek = ej;
                }
            }

          if (k != i)
            {
              /* swap eigenvalues */
              gsl_vector_swap_elements (eval, i, k);

              /* swap eigenvectors */
              gsl_matrix_swap_columns (evec, i, k);
            }
        }

      return GSL_SUCCESS;
    }
}

int
gsl_eigen_hermv_sort (gsl_vector * eval, gsl_matrix_complex * evec,
                      gsl_eigen_sort_t sort_type)
{
  if (evec->size1 != evec->size2)
    {
      GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
    }
  else if (eval->size != evec->size1)
    {
      GSL_ERROR ("eigenvalues must match eigenvector matrix", GSL_EBADLEN);
    }
  else
    {
      const size_t N = eval->size;
      size_t i;

      for (i = 0; i < N - 1; i++)
        {
          size_t j;
          size_t k = i;

          double ek = gsl_vector_get (eval, i);

          /* search for something to swap */
          for (j = i + 1; j < N; j++)
            {
              int test;
              const double ej = gsl_vector_get (eval, j);

              switch (sort_type)
                {
                case GSL_EIGEN_SORT_VAL_ASC:
                  test = (ej < ek);
                  break;
                case GSL_EIGEN_SORT_VAL_DESC:
                  test = (ej > ek);
                  break;
                case GSL_EIGEN_SORT_ABS_ASC:
                  test = (fabs (ej) < fabs (ek));
                  break;
                case GSL_EIGEN_SORT_ABS_DESC:
                  test = (fabs (ej) > fabs (ek));
                  break;
                default:
                  GSL_ERROR ("unrecognized sort type", GSL_EINVAL);
                }

              if (test)
                {
                  k = j;
                  ek = ej;
                }
            }

          if (k != i)
            {
              /* swap eigenvalues */
              gsl_vector_swap_elements (eval, i, k);

              /* swap eigenvectors */
              gsl_matrix_complex_swap_columns (evec, i, k);
            }
        }

      return GSL_SUCCESS;
    }
}

int
gsl_eigen_nonsymmv_sort (gsl_vector_complex * eval,
                         gsl_matrix_complex * evec,
                         gsl_eigen_sort_t sort_type)
{
  if (evec->size1 != evec->size2)
    {
      GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
    }
  else if (eval->size != evec->size1)
    {
      GSL_ERROR ("eigenvalues must match eigenvector matrix", GSL_EBADLEN);
    }
  else
    {
      const size_t N = eval->size;
      size_t i;

      for (i = 0; i < N - 1; i++)
        {
          size_t j;
          size_t k = i;

          gsl_complex ek = gsl_vector_complex_get (eval, i);

          /* search for something to swap */
          for (j = i + 1; j < N; j++)
            {
              int test;
              const gsl_complex ej = gsl_vector_complex_get (eval, j);

              switch (sort_type)
                {
                case GSL_EIGEN_SORT_ABS_ASC:
                  test = (gsl_complex_abs (ej) < gsl_complex_abs (ek));
                  break;
                case GSL_EIGEN_SORT_ABS_DESC:
                  test = (gsl_complex_abs (ej) > gsl_complex_abs (ek));
                  break;
                case GSL_EIGEN_SORT_VAL_ASC:
                case GSL_EIGEN_SORT_VAL_DESC:
                default:
                  GSL_ERROR ("invalid sort type", GSL_EINVAL);
                }

              if (test)
                {
                  k = j;
                  ek = ej;
                }
            }

          if (k != i)
            {
              /* swap eigenvalues */
              gsl_vector_complex_swap_elements (eval, i, k);

              /* swap eigenvectors */
              gsl_matrix_complex_swap_columns (evec, i, k);
            }
        }

      return GSL_SUCCESS;
    }
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/**
 *
 * @precisions normal z -> c d s
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <cblas.h>
#include <lapacke.h>
#include <plasma.h>
#include <core_blas.h>
#include "auxiliary.h"

/*-------------------------------------------------------------------
 * Check the orthogonality of Q
 */

int z_check_orthogonality(int M, int N, int LDQ, PLASMA_Complex64_t *Q)
{
    double alpha, beta;
    double normQ;
    int info_ortho;
    int i;
    int minMN = min(M, N);
    double eps;
    double *work = (double *)malloc(minMN*sizeof(double));

    eps = LAPACKE_dlamch_work('e');
    alpha = 1.0;
    beta  = -1.0;

    /* Build the idendity matrix USE DLASET?*/
    PLASMA_Complex64_t *Id = (PLASMA_Complex64_t *) malloc(minMN*minMN*sizeof(PLASMA_Complex64_t));
    memset((void*)Id, 0, minMN*minMN*sizeof(PLASMA_Complex64_t));
    for (i = 0; i < minMN; i++)
        Id[i*minMN+i] = (PLASMA_Complex64_t)1.0;

    /* Perform Id - Q'Q */
    if (M >= N)
        cblas_zherk(CblasColMajor, CblasUpper, CblasConjTrans, N, M, alpha, Q, LDQ, beta, Id, N);
    else
        cblas_zherk(CblasColMajor, CblasUpper, CblasNoTrans, M, N, alpha, Q, LDQ, beta, Id, M);

    normQ = LAPACKE_zlansy_work(LAPACK_COL_MAJOR, 'i', 'u', minMN, Id, minMN, work);

    printf("============\n");
    printf("Checking the orthogonality of Q \n");
    printf("||Id-Q'*Q||_oo / (N*eps) = %e \n",normQ/(minMN*eps));

    if ( isnan(normQ / (minMN * eps)) || (normQ / (minMN * eps) > 10.0) ) {
        printf("-- Orthogonality is suspicious ! \n");
        info_ortho=1;
    }
    else {
        printf("-- Orthogonality is CORRECT ! \n");
        info_ortho=0;
    }

    free(work); free(Id);

    return info_ortho;
}

/*------------------------------------------------------------
 *  Check the factorization QR
 */

int z_check_QRfactorization(int M, int N, PLASMA_Complex64_t *A1, PLASMA_Complex64_t *A2, int LDA, PLASMA_Complex64_t *Q)
{
    double Anorm, Rnorm;
    PLASMA_Complex64_t alpha, beta;
    int info_factorization;
    int i,j;
    double eps;

    eps = LAPACKE_dlamch_work('e');

    PLASMA_Complex64_t *Ql       = (PLASMA_Complex64_t *)malloc(M*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *Residual = (PLASMA_Complex64_t *)malloc(M*N*sizeof(PLASMA_Complex64_t));
    double *work              = (double *)malloc(max(M,N)*sizeof(double));

    alpha=1.0;
    beta=0.0;

    if (M >= N) {
        /* Extract the R */
        PLASMA_Complex64_t *R = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));
        memset((void*)R, 0, N*N*sizeof(PLASMA_Complex64_t));
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'u', M, N, A2, LDA, R, N);

        /* Perform Ql=Q*R */
        memset((void*)Ql, 0, M*N*sizeof(PLASMA_Complex64_t));
        cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, N, CBLAS_SADDR(alpha), Q, LDA, R, N, CBLAS_SADDR(beta), Ql, M);
        free(R);
    }
    else {
        /* Extract the L */
        PLASMA_Complex64_t *L = (PLASMA_Complex64_t *)malloc(M*M*sizeof(PLASMA_Complex64_t));
        memset((void*)L, 0, M*M*sizeof(PLASMA_Complex64_t));
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'l', M, N, A2, LDA, L, M);

    /* Perform Ql=LQ */
        memset((void*)Ql, 0, M*N*sizeof(PLASMA_Complex64_t));
        cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, M, CBLAS_SADDR(alpha), L, M, Q, LDA, CBLAS_SADDR(beta), Ql, M);
        free(L);
    }

    /* Compute the Residual */
    for (i = 0; i < M; i++)
        for (j = 0 ; j < N; j++)
            Residual[j*M+i] = A1[j*LDA+i]-Ql[j*M+i];

    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'i', M, N, Residual, M, work);
    Anorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'i', M, N, A2, LDA, work);

    if (M >= N) {
        printf("============\n");
        printf("Checking the QR Factorization \n");
        printf("-- ||A-QR||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
    }
    else {
        printf("============\n");
        printf("Checking the LQ Factorization \n");
        printf("-- ||A-LQ||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
    }

    if (isnan(Rnorm / (Anorm * N *eps)) || (Rnorm / (Anorm * N * eps) > 10.0) ) {
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else {
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    free(work); free(Ql); free(Residual);

    return info_factorization;
}

/*------------------------------------------------------------------------
 *  Check the factorization of the matrix A2
 */

int z_check_LLTfactorization(int N, PLASMA_Complex64_t *A1, PLASMA_Complex64_t *A2, int LDA, int uplo)
{
    double Anorm, Rnorm;
    PLASMA_Complex64_t alpha;
    int info_factorization;
    int i,j;
    double eps;

    eps = LAPACKE_dlamch_work('e');

    PLASMA_Complex64_t *Residual = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *L1       = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *L2       = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));
    double *work              = (double *)malloc(N*sizeof(double));

    memset((void*)L1, 0, N*N*sizeof(PLASMA_Complex64_t));
    memset((void*)L2, 0, N*N*sizeof(PLASMA_Complex64_t));

    alpha= 1.0;

    LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,' ', N, N, A1, LDA, Residual, N);

    /* Dealing with L'L or U'U  */
    if (uplo == PlasmaUpper){
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L1, N);
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L2, N);
        cblas_ztrmm(CblasColMajor, CblasLeft, CblasUpper, CblasConjTrans, CblasNonUnit, N, N, CBLAS_SADDR(alpha), L1, N, L2, N);
    }
    else{
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L1, N);
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L2, N);
        cblas_ztrmm(CblasColMajor, CblasRight, CblasLower, CblasConjTrans, CblasNonUnit, N, N, CBLAS_SADDR(alpha), L1, N, L2, N);
    }

    /* Compute the Residual || A -L'L|| */
    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
           Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i];

    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'i', N, N, Residual, N, work);
    Anorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'i', N, N, A1, LDA, work);

    printf("============\n");
    printf("Checking the Cholesky Factorization \n");
    printf("-- ||L'L-A||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));

    if ( isnan(Rnorm/(Anorm*N*eps)) || (Rnorm/(Anorm*N*eps) > 10.0) ){
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else{
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    free(Residual); free(L1); free(L2); free(work);

    return info_factorization;
}

/*--------------------------------------------------------------
 * Check the gemm
 */
double z_check_gemm(PLASMA_enum transA, PLASMA_enum transB, int M, int N, int K,
                   PLASMA_Complex64_t alpha, PLASMA_Complex64_t *A, int LDA,
                   PLASMA_Complex64_t *B, int LDB,
                   PLASMA_Complex64_t beta, PLASMA_Complex64_t *Cplasma,
                   PLASMA_Complex64_t *Cref, int LDC,
                   double *Cinitnorm, double *Cplasmanorm, double *Clapacknorm )
{
    PLASMA_Complex64_t beta_const = -1.0;
    double Rnorm;
    double *work = (double *)malloc(max(K,max(M, N))* sizeof(double));

    *Cinitnorm   = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'i', M, N, Cref,    LDC, work);
    *Cplasmanorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'i', M, N, Cplasma, LDC, work);

    cblas_zgemm(CblasColMajor, (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB, M, N, K,
                CBLAS_SADDR(alpha), A, LDA, B, LDB, CBLAS_SADDR(beta), Cref, LDC);

    *Clapacknorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'i', M, N, Cref, LDC, work);

    cblas_zaxpy(LDC * N, CBLAS_SADDR(beta_const), Cplasma, 1, Cref, 1);

    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'i', M, N, Cref, LDC, work);

    free(work);

    return Rnorm;
}

/*--------------------------------------------------------------
 * Check the trsm
 */
double z_check_trsm(PLASMA_enum side, PLASMA_enum uplo, PLASMA_enum trans, PLASMA_enum diag,
                   int M, int NRHS, PLASMA_Complex64_t alpha,
                   PLASMA_Complex64_t *A, int LDA,
                   PLASMA_Complex64_t *Bplasma, PLASMA_Complex64_t *Bref, int LDB,
                   double *Binitnorm, double *Bplasmanorm, double *Blapacknorm )
{
    PLASMA_Complex64_t beta_const = -1.0;
    double Rnorm;
    double *work = (double *)malloc(max(M, NRHS)* sizeof(double));
    /*double eps = LAPACKE_dlamch_work('e');*/

    *Binitnorm   = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'i', M, NRHS, Bref,    LDB, work);
    *Bplasmanorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'm', M, NRHS, Bplasma, LDB, work);

    cblas_ztrsm(CblasColMajor, (CBLAS_SIDE)side, (CBLAS_UPLO)uplo,
                (CBLAS_TRANSPOSE)trans, (CBLAS_DIAG)diag, M, NRHS,
                CBLAS_SADDR(alpha), A, LDA, Bref, LDB);

    *Blapacknorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'm', M, NRHS, Bref, LDB, work);

    cblas_zaxpy(LDB * NRHS, CBLAS_SADDR(beta_const), Bplasma, 1, Bref, 1);

    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'm', M, NRHS, Bref, LDB, work);
    Rnorm = Rnorm / *Blapacknorm;
    /* max(M,NRHS) * eps);*/

    free(work);

    return Rnorm;
}

/*--------------------------------------------------------------
 * Check the solution
 */

double z_check_solution(int M, int N, int NRHS, PLASMA_Complex64_t *A, int LDA,
                      PLASMA_Complex64_t *B,  PLASMA_Complex64_t *X, int LDB,
                      double *anorm, double *bnorm, double *xnorm )
{
/*     int info_solution; */
    double Rnorm = -1.00;
    PLASMA_Complex64_t zone  =  1.0;
    PLASMA_Complex64_t mzone = -1.0;
    double *work = (double *)malloc(max(M, N)* sizeof(double));

    *anorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'i', M, N,    A, LDA, work);
    *xnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'i', M, NRHS, X, LDB, work);
    *bnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'i', N, NRHS, B, LDB, work);

    cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, NRHS, N, CBLAS_SADDR(zone), A, LDA, X, LDB, CBLAS_SADDR(mzone), B, LDB);

    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'i', N, NRHS, B, LDB, work);

    free(work);

    return Rnorm;
}


text:
/*************************************************************************************
 * Copyright (C) 2017, Trung Pham <me@trungbpham.com>                                *
 * All rights reserved.                                                              *
 *                                                                                *
 * Redistribution and use in source and binary forms, with or without                *
 * modification, are permitted provided that the following conditions are met:       *
 *                                                                                *
 *       * Redistributions of source code must retain the above copyright notice,    *
 *       this list of conditions and the following disclaimer.                       *
 *                                                                                *
 *       * Redistributions in binary form must reproduce the above copyright notice, *
 *       this list of conditions and the following disclaimer in the documentation   *
 *       and/or other materials provided with the distribution.                      *
 *                                                                                *
 *       * Neither the name of cfeap nor the names of its contributors may be used   *
 *       to endorse or promote products derived from this software without specific  *
 *       prior written permission.                                                   *
 *                                                                                *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"       *
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE         *
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE    *
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE      *
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL        *
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR        *
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER        *
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,     *
 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE     *
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.              *
 *************************************************************************************/

#ifndef __CBLAS_WRAPPER_H__
#define __CBLAS_WRAPPER_H__

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

int cblw_scalar_vec_mul_x(size_t dim, double a, double *x);

int cblw_scalar_mat_mul_x(size_t col, size_t row, double a, double *X);

int cblw_scalar_vec_mul_y(size_t dim, double a, double *x, double *y);

int cblw_scalar_mat_mul_y(size_t col, size_t row, double a, double *X, double *Y);

int cblw_mat_vec_mul(size_t col, size_t row, double *A, double *x, double *y);

int cblw_mat_mat_mul(size_t rowA, size_t colB, size_t colA_rowB, double *A, double *B, double *C);

#endif //__CBLAS_WRAPPER_H__

ext: h


text:
/* vector/gsl_vector_complex_long_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_VECTOR_COMPLEX_LONG_DOUBLE_H__
#define __GSL_VECTOR_COMPLEX_LONG_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_complex.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_long_double.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_block_complex_long_double.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;
  long double *data;
  gsl_block_complex_long_double *block;
  int owner;
} gsl_vector_complex_long_double;

typedef struct
{
  gsl_vector_complex_long_double vector;
} _gsl_vector_complex_long_double_view;

typedef _gsl_vector_complex_long_double_view gsl_vector_complex_long_double_view;

typedef struct
{
  gsl_vector_complex_long_double vector;
} _gsl_vector_complex_long_double_const_view;

typedef const _gsl_vector_complex_long_double_const_view gsl_vector_complex_long_double_const_view;

/* Allocation */

GSL_FUN gsl_vector_complex_long_double *gsl_vector_complex_long_double_alloc (const size_t n);
GSL_FUN gsl_vector_complex_long_double *gsl_vector_complex_long_double_calloc (const size_t n);

GSL_FUN gsl_vector_complex_long_double *
gsl_vector_complex_long_double_alloc_from_block (gsl_block_complex_long_double * b,
                                           const size_t offset,
                                           const size_t n,
                                           const size_t stride);

GSL_FUN gsl_vector_complex_long_double *
gsl_vector_complex_long_double_alloc_from_vector (gsl_vector_complex_long_double * v,
                                             const size_t offset,
                                             const size_t n,
                                             const size_t stride);

GSL_FUN void gsl_vector_complex_long_double_free (gsl_vector_complex_long_double * v);

/* Views */

GSL_FUN _gsl_vector_complex_long_double_view
gsl_vector_complex_long_double_view_array (long double *base,
                                     size_t n);

GSL_FUN _gsl_vector_complex_long_double_view
gsl_vector_complex_long_double_view_array_with_stride (long double *base,
                                                 size_t stride,
                                                 size_t n);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_vector_complex_long_double_const_view_array (const long double *base,
                                           size_t n);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_vector_complex_long_double_const_view_array_with_stride (const long double *base,
                                                       size_t stride,
                                                       size_t n);

GSL_FUN _gsl_vector_complex_long_double_view
gsl_vector_complex_long_double_subvector (gsl_vector_complex_long_double *base,
                                         size_t i,
                                         size_t n);

GSL_FUN _gsl_vector_complex_long_double_view
gsl_vector_complex_long_double_subvector_with_stride (gsl_vector_complex_long_double *v,
                                                size_t i,
                                                size_t stride,
                                                size_t n);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_vector_complex_long_double_const_subvector (const gsl_vector_complex_long_double *base,
                                               size_t i,
                                               size_t n);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_vector_complex_long_double_const_subvector_with_stride (const gsl_vector_complex_long_double *v,
                                                      size_t i,
                                                      size_t stride,
                                                      size_t n);

GSL_FUN _gsl_vector_long_double_view
gsl_vector_complex_long_double_real (gsl_vector_complex_long_double *v);

GSL_FUN _gsl_vector_long_double_view
gsl_vector_complex_long_double_imag (gsl_vector_complex_long_double *v);

GSL_FUN _gsl_vector_long_double_const_view
gsl_vector_complex_long_double_const_real (const gsl_vector_complex_long_double *v);

GSL_FUN _gsl_vector_long_double_const_view
gsl_vector_complex_long_double_const_imag (const gsl_vector_complex_long_double *v);

/* Operations */

GSL_FUN void gsl_vector_complex_long_double_set_zero (gsl_vector_complex_long_double * v);
GSL_FUN void gsl_vector_complex_long_double_set_all (gsl_vector_complex_long_double * v,
                                       gsl_complex_long_double z);
GSL_FUN int gsl_vector_complex_long_double_set_basis (gsl_vector_complex_long_double * v, size_t i);

GSL_FUN int gsl_vector_complex_long_double_fread (FILE * stream,
                                    gsl_vector_complex_long_double * v);
GSL_FUN int gsl_vector_complex_long_double_fwrite (FILE * stream,
                                     const gsl_vector_complex_long_double * v);
GSL_FUN int gsl_vector_complex_long_double_fscanf (FILE * stream,
                                     gsl_vector_complex_long_double * v);
GSL_FUN int gsl_vector_complex_long_double_fprintf (FILE * stream,
                                      const gsl_vector_complex_long_double * v,
                                      const char *format);

GSL_FUN int gsl_vector_complex_long_double_memcpy (gsl_vector_complex_long_double * dest, const gsl_vector_complex_long_double * src);

GSL_FUN int gsl_vector_complex_long_double_reverse (gsl_vector_complex_long_double * v);

GSL_FUN int gsl_vector_complex_long_double_swap (gsl_vector_complex_long_double * v, gsl_vector_complex_long_double * w);
GSL_FUN int gsl_vector_complex_long_double_swap_elements (gsl_vector_complex_long_double * v, const size_t i, const size_t j);

GSL_FUN int gsl_vector_complex_long_double_equal (const gsl_vector_complex_long_double * u,
                                    const gsl_vector_complex_long_double * v);

GSL_FUN int gsl_vector_complex_long_double_isnull (const gsl_vector_complex_long_double * v);
GSL_FUN int gsl_vector_complex_long_double_ispos (const gsl_vector_complex_long_double * v);
GSL_FUN int gsl_vector_complex_long_double_isneg (const gsl_vector_complex_long_double * v);
GSL_FUN int gsl_vector_complex_long_double_isnonneg (const gsl_vector_complex_long_double * v);

GSL_FUN int gsl_vector_complex_long_double_add (gsl_vector_complex_long_double * a, const gsl_vector_complex_long_double * b);
GSL_FUN int gsl_vector_complex_long_double_sub (gsl_vector_complex_long_double * a, const gsl_vector_complex_long_double * b);
GSL_FUN int gsl_vector_complex_long_double_mul (gsl_vector_complex_long_double * a, const gsl_vector_complex_long_double * b);
GSL_FUN int gsl_vector_complex_long_double_div (gsl_vector_complex_long_double * a, const gsl_vector_complex_long_double * b);
GSL_FUN int gsl_vector_complex_long_double_scale (gsl_vector_complex_long_double * a, const gsl_complex_long_double x);
GSL_FUN int gsl_vector_complex_long_double_add_constant (gsl_vector_complex_long_double * a, const gsl_complex_long_double x);
GSL_FUN int gsl_vector_complex_long_double_axpby (const gsl_complex_long_double alpha, const gsl_vector_complex_long_double * x, const gsl_complex_long_double beta, gsl_vector_complex_long_double * y);

GSL_FUN INLINE_DECL gsl_complex_long_double gsl_vector_complex_long_double_get (const gsl_vector_complex_long_double * v, const size_t i);
GSL_FUN INLINE_DECL void gsl_vector_complex_long_double_set (gsl_vector_complex_long_double * v, const size_t i, gsl_complex_long_double z);
GSL_FUN INLINE_DECL gsl_complex_long_double *gsl_vector_complex_long_double_ptr (gsl_vector_complex_long_double * v, const size_t i);
GSL_FUN INLINE_DECL const gsl_complex_long_double *gsl_vector_complex_long_double_const_ptr (const gsl_vector_complex_long_double * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
gsl_complex_long_double
gsl_vector_complex_long_double_get (const gsl_vector_complex_long_double * v,
                              const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      gsl_complex_long_double zero = {{0, 0}};
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, zero);
    }
#endif
  return *GSL_COMPLEX_LONG_DOUBLE_AT (v, i);
}

INLINE_FUN
void
gsl_vector_complex_long_double_set (gsl_vector_complex_long_double * v,
                              const size_t i, gsl_complex_long_double z)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  *GSL_COMPLEX_LONG_DOUBLE_AT (v, i) = z;
}

INLINE_FUN
gsl_complex_long_double *
gsl_vector_complex_long_double_ptr (gsl_vector_complex_long_double * v,
                              const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return GSL_COMPLEX_LONG_DOUBLE_AT (v, i);
}

INLINE_FUN
const gsl_complex_long_double *
gsl_vector_complex_long_double_const_ptr (const gsl_vector_complex_long_double * v,
                                    const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return GSL_COMPLEX_LONG_DOUBLE_AT (v, i);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_COMPLEX_LONG_DOUBLE_H__ */

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


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

/* Runge-Kutta 8(9), Prince-Dormand
 *
 * High Order Embedded Runge-Kutta Formulae
 * P.J. Prince and J.R. Dormand
 * J. Comp. Appl. Math.,7, pp. 67-75, 1981
 */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>

#include "odeiv_util.h"

/* Prince-Dormand constants */

static const double Abar[] = {
  14005451.0 / 335480064.0,
  0.0,
  0.0,
  0.0,
  0.0,
  -59238493.0 / 1068277825.0,
  181606767.0 / 758867731.0,
  561292985.0 / 797845732.0,
  -1041891430.0 / 1371343529.0,
  760417239.0 / 1151165299.0,
  118820643.0 / 751138087.0,
  -528747749.0 / 2220607170.0,
  1.0 / 4.0
};

static const double A[] = {
  13451932.0 / 455176623.0,
  0.0,
  0.0,
  0.0,
  0.0,
  -808719846.0 / 976000145.0,
  1757004468.0 / 5645159321.0,
  656045339.0 / 265891186.0,
  -3867574721.0 / 1518517206.0,
  465885868.0 / 322736535.0,
  53011238.0 / 667516719.0,
  2.0 / 45.0
};

static const double ah[] = {
  1.0 / 18.0,
  1.0 / 12.0,
  1.0 / 8.0,
  5.0 / 16.0,
  3.0 / 8.0,
  59.0 / 400.0,
  93.0 / 200.0,
  5490023248.0 / 9719169821.0,
  13.0 / 20.0,
  1201146811.0 / 1299019798.0
};

static const double b21 = 1.0 / 18.0;
static const double b3[] = { 1.0 / 48.0, 1.0 / 16.0 };
static const double b4[] = { 1.0 / 32.0, 0.0, 3.0 / 32.0 };
static const double b5[] = { 5.0 / 16.0, 0.0, -75.0 / 64.0, 75.0 / 64.0 };
static const double b6[] = { 3.0 / 80.0, 0.0, 0.0, 3.0 / 16.0, 3.0 / 20.0 };
static const double b7[] = {
  29443841.0 / 614563906.0,
  0.0,
  0.0,
  77736538.0 / 692538347.0,
  -28693883.0 / 1125000000.0,
  23124283.0 / 1800000000.0
};
static const double b8[] = {
  16016141.0 / 946692911.0,
  0.0,
  0.0,
  61564180.0 / 158732637.0,
  22789713.0 / 633445777.0,
  545815736.0 / 2771057229.0,
  -180193667.0 / 1043307555.0
};
static const double b9[] = {
  39632708.0 / 573591083.0,
  0.0,
  0.0,
  -433636366.0 / 683701615.0,
  -421739975.0 / 2616292301.0,
  100302831.0 / 723423059.0,
  790204164.0 / 839813087.0,
  800635310.0 / 3783071287.0
};
static const double b10[] = {
  246121993.0 / 1340847787.0,
  0.0,
  0.0,
  -37695042795.0 / 15268766246.0,
  -309121744.0 / 1061227803.0,
  -12992083.0 / 490766935.0,
  6005943493.0 / 2108947869.0,
  393006217.0 / 1396673457.0,
  123872331.0 / 1001029789.0
};
static const double b11[] = {
  -1028468189.0 / 846180014.0,
  0.0,
  0.0,
  8478235783.0 / 508512852.0,
  1311729495.0 / 1432422823.0,
  -10304129995.0 / 1701304382.0,
  -48777925059.0 / 3047939560.0,
  15336726248.0 / 1032824649.0,
  -45442868181.0 / 3398467696.0,
  3065993473.0 / 597172653.0
};
static const double b12[] = {
  185892177.0 / 718116043.0,
  0.0,
  0.0,
  -3185094517.0 / 667107341.0,
  -477755414.0 / 1098053517.0,
  -703635378.0 / 230739211.0,
  5731566787.0 / 1027545527.0,
  5232866602.0 / 850066563.0,
  -4093664535.0 / 808688257.0,
  3962137247.0 / 1805957418.0,
  65686358.0 / 487910083.0
};
static const double b13[] = {
  403863854.0 / 491063109.0,
  0.0,
  0.0,
  -5068492393.0 / 434740067.0,
  -411421997.0 / 543043805.0,
  652783627.0 / 914296604.0,
  11173962825.0 / 925320556.0,
  -13158990841.0 / 6184727034.0,
  3936647629.0 / 1978049680.0,
  -160528059.0 / 685178525.0,
  248638103.0 / 1413531060.0,
  0.0
};

typedef struct
{
  double *k[13];
  double *ytmp;
  double *y0;
}
rk8pd_state_t;

static void *
rk8pd_alloc (size_t dim)
{
  rk8pd_state_t *state = (rk8pd_state_t *) malloc (sizeof (rk8pd_state_t));
  int i, j;

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for rk8pd_state", GSL_ENOMEM);
    }

  state->ytmp = (double *) malloc (dim * sizeof (double));

  if (state->ytmp == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
    }

  state->y0 = (double *) malloc (dim * sizeof (double));

  if (state->y0 == 0)
    {
      free (state->ytmp);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
    }

  for (i = 0; i < 13; i++)
    {
      state->k[i] = (double *) malloc (dim * sizeof (double));

      if (state->k[i] == 0)
        {
          for (j = 0; j < i; j++)
            {
              free (state->k[j]);
            }
          free (state->y0);
          free (state->ytmp);
          free (state);
          GSL_ERROR_NULL ("failed to allocate space for k's", GSL_ENOMEM);
        }
    }

  return state;
}

static int
rk8pd_apply (void *vstate,
             size_t dim,
             double t,
             double h,
             double y[],
             double yerr[],
             const double dydt_in[],
             double dydt_out[], const gsl_odeiv_system * sys)
{
  rk8pd_state_t *state = (rk8pd_state_t *) vstate;

  size_t i;

  double *const ytmp = state->ytmp;
  double *const y0 = state->y0;
  /* Note that k1 is stored in state->k[0] due to zero-based indexing */
  double *const k1 = state->k[0];
  double *const k2 = state->k[1];
  double *const k3 = state->k[2];
  double *const k4 = state->k[3];
  double *const k5 = state->k[4];
  double *const k6 = state->k[5];
  double *const k7 = state->k[6];
  double *const k8 = state->k[7];
  double *const k9 = state->k[8];
  double *const k10 = state->k[9];
  double *const k11 = state->k[10];
  double *const k12 = state->k[11];
  double *const k13 = state->k[12];

  DBL_MEMCPY (y0, y, dim);

  /* k1 step */
  if (dydt_in != NULL)
    {
      DBL_MEMCPY (k1, dydt_in, dim);
    }
  else
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t, y, k1);

      if (s != GSL_SUCCESS)
	{
	  return s;
	}
    }

  for (i = 0; i < dim; i++)
    ytmp[i] = y[i] + b21 * h * k1[i];

  /* k2 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[0] * h, ytmp, k2);

    if (s != GSL_SUCCESS)
      {
	return s;
      }
  }

  for (i = 0; i < dim; i++)
    ytmp[i] = y[i] + h * (b3[0] * k1[i] + b3[1] * k2[i]);

  /* k3 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[1] * h, ytmp, k3);

    if (s != GSL_SUCCESS)
      {
	return s;
      }
  }

  for (i = 0; i < dim; i++)
    ytmp[i] = y[i] + h * (b4[0] * k1[i] + b4[2] * k3[i]);

  /* k4 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[2] * h, ytmp, k4);

    if (s != GSL_SUCCESS)
      {
	return s;
      }
  }

  for (i = 0; i < dim; i++)
    ytmp[i] = y[i] + h * (b5[0] * k1[i] + b5[2] * k3[i] + b5[3] * k4[i]);

  /* k5 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[3] * h, ytmp, k5);

    if (s != GSL_SUCCESS)
      {
	return s;
      }
  }

  for (i = 0; i < dim; i++)
    ytmp[i] = y[i] + h * (b6[0] * k1[i] + b6[3] * k4[i] + b6[4] * k5[i]);

  /* k6 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[4] * h, ytmp, k6);

    if (s != GSL_SUCCESS)
      {
	return s;
      }
  }

  for (i = 0; i < dim; i++)
    ytmp[i] =
      y[i] + h * (b7[0] * k1[i] + b7[3] * k4[i] + b7[4] * k5[i] +
                  b7[5] * k6[i]);

  /* k7 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[5] * h, ytmp, k7);

    if (s != GSL_SUCCESS)
      {
	return s;
      }
  }

  for (i = 0; i < dim; i++)
    ytmp[i] =
      y[i] + h * (b8[0] * k1[i] + b8[3] * k4[i] + b8[4] * k5[i] +
                  b8[5] * k6[i] + b8[6] * k7[i]);

  /* k8 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[6] * h, ytmp, k8);

    if (s != GSL_SUCCESS)
      {
	return s;
      }
  }

  for (i = 0; i < dim; i++)
    ytmp[i] =
      y[i] + h * (b9[0] * k1[i] + b9[3] * k4[i] + b9[4] * k5[i] +
                  b9[5] * k6[i] + b9[6] * k7[i] + b9[7] * k8[i]);

  /* k9 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[7] * h, ytmp, k9);

    if (s != GSL_SUCCESS)
      {
	return s;
      }
  }

  for (i = 0; i < dim; i++)
    ytmp[i] =
      y[i] + h * (b10[0] * k1[i] + b10[3] * k4[i] + b10[4] * k5[i] +
                  b10[5] * k6[i] + b10[6] * k7[i] + b10[7] * k8[i] +
                  b10[8] * k9[i]);

  /* k10 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[8] * h, ytmp, k10);

    if (s != GSL_SUCCESS)
      {
	return s;
      }
  }

  for (i = 0; i < dim; i++)
    ytmp[i] =
      y[i] + h * (b11[0] * k1[i] + b11[3] * k4[i] + b11[4] * k5[i] +
                  b11[5] * k6[i] + b11[6] * k7[i] + b11[7] * k8[i] +
                  b11[8] * k9[i] + b11[9] * k10[i]);

  /* k11 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + ah[9] * h, ytmp, k11);

    if (s != GSL_SUCCESS)
      {
	return s;
      }
  }

  for (i = 0; i < dim; i++)
    ytmp[i] =
      y[i] + h * (b12[0] * k1[i] + b12[3] * k4[i] + b12[4] * k5[i] +
                  b12[5] * k6[i] + b12[6] * k7[i] + b12[7] * k8[i] +
                  b12[8] * k9[i] + b12[9] * k10[i] + b12[10] * k11[i]);

  /* k12 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + h, ytmp, k12);

    if (s != GSL_SUCCESS)
      {
	return s;
      }
  }

  for (i = 0; i < dim; i++)
    ytmp[i] =
      y[i] + h * (b13[0] * k1[i] + b13[3] * k4[i] + b13[4] * k5[i] +
                  b13[5] * k6[i] + b13[6] * k7[i] + b13[7] * k8[i] +
                  b13[8] * k9[i] + b13[9] * k10[i] + b13[10] * k11[i] +
                  b13[11] * k12[i]);

  /* k13 step */
  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + h, ytmp, k13);

    if (s != GSL_SUCCESS)
      {
	return s;
      }
  }

  /* final sum  */
  for (i = 0; i < dim; i++)
    {
      const double ksum8 =
        Abar[0] * k1[i] + Abar[5] * k6[i] + Abar[6] * k7[i] +
        Abar[7] * k8[i] + Abar[8] * k9[i] + Abar[9] * k10[i] +
        Abar[10] * k11[i] + Abar[11] * k12[i] + Abar[12] * k13[i];
      y[i] += h * ksum8;
    }

  /* Evaluate dydt_out[]. */

  if (dydt_out != NULL)
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t + h, y, dydt_out);

      if (s != GSL_SUCCESS)
	{
	  /* Restore initial values */
	  DBL_MEMCPY (y, y0, dim);
	  return s;
	}
    }

  /* error estimate */
  for (i = 0; i < dim; i++)
    {
      const double ksum8 =
        Abar[0] * k1[i] + Abar[5] * k6[i] + Abar[6] * k7[i] +
        Abar[7] * k8[i] + Abar[8] * k9[i] + Abar[9] * k10[i] +
        Abar[10] * k11[i] + Abar[11] * k12[i] + Abar[12] * k13[i];
      const double ksum7 =
        A[0] * k1[i] + A[5] * k6[i] + A[6] * k7[i] + A[7] * k8[i] +
        A[8] * k9[i] + A[9] * k10[i] + A[10] * k11[i] + A[11] * k12[i];
      yerr[i] = h * (ksum7 - ksum8);
    }

  return GSL_SUCCESS;
}

static int
rk8pd_reset (void *vstate, size_t dim)
{
  rk8pd_state_t *state = (rk8pd_state_t *) vstate;

  int i;

  for (i = 0; i < 13; i++)
    {
      DBL_ZERO_MEMSET (state->k[i], dim);
    }

  DBL_ZERO_MEMSET (state->y0, dim);
  DBL_ZERO_MEMSET (state->ytmp, dim);

  return GSL_SUCCESS;
}

static unsigned int
rk8pd_order (void *vstate)
{
  rk8pd_state_t *state = (rk8pd_state_t *) vstate;
  state = 0; /* prevent warnings about unused parameters */
  return 8;
}

static void
rk8pd_free (void *vstate)
{
  rk8pd_state_t *state = (rk8pd_state_t *) vstate;
  int i;

  for (i = 0; i < 13; i++)
    {
      free (state->k[i]);
    }
  free (state->y0);
  free (state->ytmp);
  free (state);
}

static const gsl_odeiv_step_type rk8pd_type = { "rk8pd",        /* name */
  1,                            /* can use dydt_in */
  1,                            /* gives exact dydt_out */
  &rk8pd_alloc,
  &rk8pd_apply,
  &rk8pd_reset,
  &rk8pd_order,
  &rk8pd_free
};

const gsl_odeiv_step_type *gsl_odeiv_step_rk8pd = &rk8pd_type;


text:
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <signal.h>
#include <string.h>
#include <math.h>
#include "cblas.h"
#include <lapacke.h>
#include <unistd.h>
#include <sys/types.h>

#include "spidriver_host.h"
#include "adcdriver_host.h"
#include "matrix_utils.h"

#define PI 3.1415926535

// Length of data buffer
#define NUMPTS 128

// Number of signal vectors
#define PSIG 2

// Sampling frequency.  Must match the sampling frequency commanded
// to the A/D.
#define FSAMP 15625

// Parameters used in searching for peak
#define MAXRECURSIONS 5
#define NGRID 25

//===========================================================
// Utility functions for MUSIC

//-----------------------------------------------------
void find_bracket(int N, float *u, int *ileft, int *iright) {
  // Given input vector y, this finds the max element,
  // then returns the indices of the vectors to its
  // left and right.

  int i;

  i = maxeltf(N, u);
  if (i == 0) {
    *ileft = 0;
    *iright = 1;
  } else if (i == (N-1)) {
    *ileft = N-2;
    *iright = N-1;
  } else {
    *ileft = i-1;
    *iright = i+1;
  }
  return;
}

//-----------------------------------------------------
void extract_noise_vectors(float *A, int m, int n, int c, float *E) {
  // This fcn takes input matrix A of size mxn.  It extracts the vectors
  // of A to the right, starting at col c, and puts the extracted
  // vectors into E.  E has size [m, n-c]

  //printf("Entered extract_noise_vectors, m = %d, n = %d, c = %d\n", m, n, c);
  int i, j, l;
  for (i = 0; i < m; i++) {
    for (j = c; j < n; j++) {
      // printf("i = %d, j = %d, A[i,j] = %f\n", i, j, MATRIX_ELEMENT(A, m, n, i, j));
      l = lindex(m, n-c, i, j-c);
      // printf("lindex(m, n-c, i,j-c) = %d\n", l);
      E[l] = MATRIX_ELEMENT(A, m, n, i, j);
    }
  }
}

//-----------------------------------------------------
float music_sum(float f, float *v, int Mr, int Mc) {
  // This performs the sum over noise space vectors.
  // Since complex numbers are not supported, I split
  // the computation into two parts, real and imag.

  int i, j;

  float *er;
  float *ei;
  float *vi;
  float s, tr, ti;

  // printf("--> Entered music_sum, [row, col] = [%d, %d]\n", Mr, Mc);

  er = (float*) malloc(Mr*sizeof(float));
  ei = (float*) malloc(Mr*sizeof(float));
  vi = (float*) malloc(Mr*sizeof(float));

  // Create e vector
  for(i=0; i<Mr; i++) {
    er[i] = cos(2*PI*i*f);
    ei[i] = sin(2*PI*i*f);
  }

  // Compute denominator
  s = 0;
  for(i=0; i<Mc; i++) {     // iterate over columns.

    // v is a matrix, so extract noise vector vi as column vector.
    for(j=0; j<Mr; j++) {   // iterate over rows
      vi[j] = v[lindex(Mr, Mc, j, i)];
    }
    tr = cblas_sdot(Mr, er, 1, vi, 1);
    ti = cblas_sdot(Mr, ei, 1, vi, 1);
    s = s + tr*tr + ti*ti;
  }

  free(er);
  free(ei);
  free(vi);

  // Must guard against returning inf or nan.
  return (1.0f/s);
}

//-----------------------------------------------------
// Called when Ctrl+C is pressed - triggers the program to stop.
void stopHandler(int sig) {
  adc_quit();
  exit(0);
}

//==========================================================
// This is the main program.  It runs a loop, takes a buffer
// of data from the A/D, then uses the MUSIC algorithm to
// compute the frequency of the input sine wave.
int main (void)
{
  // Loop variables
  uint32_t i, j;

  // Buffers for tx and rx data from A/D registers.
  uint32_t tx_buf[3];
  uint32_t rx_buf[4];

  // Used in LAPACK computations
  int m = NUMPTS;
  int info;
  float superb[NUMPTS-1];

  // Measured voltages from A/D
  float v[NUMPTS];           // Vector of measurements
  float Rxx[NUMPTS*NUMPTS];  // Covariance matrix.
  float U[NUMPTS * NUMPTS];
  float S[NUMPTS];
  float VT[NUMPTS * NUMPTS];

  float Nu[NUMPTS * (NUMPTS-PSIG)];  // Vector of noises

  // Used in finding peak corresponding to dominant frequency
  float f[NGRID];
  float Pmu[NGRID];
  int ileft, iright;
  float fleft, fright, fpeak;

  // Stuff used with "hit return when ready..."
  char dummy[8];

  printf("------------   Starting main.....   -------------\n");

  // Run until Ctrl+C pressed:
  signal(SIGINT, stopHandler);

  // Sanity check user.
  if(getuid()!=0){
     printf("You must run this program as root. Exiting.\n");
     exit(EXIT_FAILURE);
  }

  // Initialize A/D converter
  adc_config();
  adc_set_samplerate(SAMP_RATE_15625);
  adc_set_chan0();

  // Now loop forever, read buffer, and compute frequency.
  // printf("--------------------------------------------------\n");
  while(1) {

    // Do A/D read to fill buffer with NUMPTS measurements.
    adc_read_multiple(NUMPTS, v);
    //printf("Values read = \n");
    //for (i=0; i<NUMPTS; i++) {
    //  printf("i = %d, v = %e\n", i, v[i]);
    //}

    // Zero out Rxx prior to filling it using sger
    zeros(m, m, Rxx);

    // Create covariance matrix by doing outer product of v with
    // itself.
    cblas_sger(CblasRowMajor,    /* Row-major storage */
               m,                /* Row count for Rxx  */
               m,                /* Col count for Rxx  */
               1.0f,             /* scale factor to apply to v*v' */
               v,
               1,                /* stride between elements of v. */
               v,
               1,                /* stride between elements of v. */
               Rxx,
               m);               /* leading dimension of matrix Rxx. */
    //printf("\nMatrix Rxx (%d x %d) =\n", m, m);
    //print_matrix(Rxx, m, m);

    // Now compute SVD of Rxx.
    info = LAPACKE_sgesvd(LAPACK_ROW_MAJOR, 'A', 'A',
                  m, m, Rxx,
                  m, S, U, m,
                  VT, m, superb);
    if (info != 0)  {
      fprintf(stderr, "Error: dgesvd returned with a non-zero status (info = %d)\n", info);
      return(-1);
    }

    //printf("\nMatrix U (%d x %d) is:\n", m, m);
    //print_matrix(U, m, m);

    //printf("\nVector S (%d x %d) is:\n", m, 1);
    //print_matrix(S, m, 1);

    //printf("\nMatrix VT (%d x %d) is:\n", m, m);
    //print_matrix(VT, m, m);

    // Extract noise vectors here.  The noise vectors are held in Nu
    extract_noise_vectors(U, m, m, PSIG, Nu);
    //printf("\nMatrix Nu (%d x %d) is:\n", m, NUMPTS-PSIG);
    //print_matrix(Nu, m, NUMPTS-PSIG);

    // Now find max freq.  Set up initial grid endpoints.  Freqs are
    // in units of Hz.
    fleft = 0.0f;
    fright = FSAMP/2.0f;

    for (j=0; j<MAXRECURSIONS; j++) {
      // printf("fleft = %f, fright = %f\n", fleft, fright);

      // Set up search grid
      linspace(fleft, fright, NGRID, f);
      //printf("\nVector f =\n");
      //print_matrix(f, NGRID, 1);

      // Compute vector of amplitudes Pmu on grid.  music_sum wants normalized
      // frequencies
      for (i = 0; i < NGRID; i++) {
        Pmu[i] = music_sum(f[i]/FSAMP, Nu, NUMPTS, NUMPTS-PSIG);
      }
      //printf("\nVector Pmu =\n");
      //print_matrix(Pmu, NGRID, 1);

      find_bracket(NGRID, Pmu, &ileft, &iright);
      fleft = f[ileft];
      fright = f[iright];
    }
    fpeak = (fleft+fright)/2.0f;   // Assume peak is average of fleft and fright
    printf("Peak frequency found at f = %f Hz\n", fpeak);

    // usleep(500000);   // delay 1/2 sec.
  }

}

max_stars_repo_name: RobertPHeller/MusicPlaypen

max_issues_repo_name: RobertPHeller/MusicPlaypen

max_forks_repo_name: RobertPHeller/MusicPlaypen


text:
#define _GNU_SOURCE
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <glib.h>
#if HAVE_MPI
#include <mpi.h>
#endif
#include "genetic/genetic.h"

#define N_SIMULATIONS 10000
#define SEED 707l
int ntasks = 1;
unsigned int nthreads = 1;
GMutex mutex[1];
GeneticVariable v[2];

double
evaluate (Entity * entity)
{
  double x, y, e1, e2;
  x = genetic_get_variable (entity, v);
  y = genetic_get_variable (entity, v + 1);
  e1 = x + y - 3.;
  e2 = x - y - 1.;
  e1 = e1 * e1 + e2 * e2;
  return e1;
}

int
main (int argn, char **argc)
{
  FILE *file;
  double *best_variables;
  char *best_genome;
  double xgenerations, mutation_ratio, reproduction_ratio, adaptation_ratio,
    evolution_ratio, best_objective;
  int rank;
  unsigned int ngenerations;
#if HAVE_MPI
  MPI_Init (&argn, &argc);
  MPI_Comm_size (MPI_COMM_WORLD, &ntasks);
  MPI_Comm_rank (MPI_COMM_WORLD, &rank);
#else
  rank = 0;
#endif
  nthreads = 4;
  v[0].maximum = 10.;
  v[0].minimum = -10.;
  v[0].nbits = 30;
  v[1].maximum = 10.;
  v[1].minimum = -10.;
  v[1].nbits = 30;
  file = fopen (argc[1], "r");
  if (fscanf (file, "%*s%lf%*s%lf%*s%lf%*s%lf", &xgenerations,
              &mutation_ratio, &reproduction_ratio, &adaptation_ratio) != 4)
    return 1;
  fclose (file);
  ngenerations = xgenerations;
  evolution_ratio = mutation_ratio + reproduction_ratio + adaptation_ratio;
  genetic_algorithm_default (2, v,
                             N_SIMULATIONS
                             / (1 + (ngenerations - 1) * evolution_ratio),
                             ngenerations, mutation_ratio, reproduction_ratio,
                             adaptation_ratio, SEED, 0., &evaluate,
                             &best_genome, &best_variables, &best_objective);
  if (rank == 0)
    {
      file = fopen (argc[2], "w");
      fprintf (file, "%.14le", best_objective);
      printf ("objective=%.14le\n", best_objective);
      fclose (file);
      g_free (best_genome);
      g_free (best_variables);
    }
#if HAVE_MPI
  MPI_Finalize ();
#endif
  return 0;
}

max_stars_repo_name: jburguete/mpcotool

max_issues_repo_name: jburguete/mpcotool

max_forks_repo_name: jburguete/mpcotool


text:
/* rng/mrg.c
 *
 * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

/* This is a fifth-order multiple recursive generator. The sequence is,

   x_n = (a_1 x_{n-1} + a_5 x_{n-5}) mod m

   with a_1 = 107374182, a_2 = a_3 = a_4 = 0, a_5 = 104480 and m = 2^31-1.

   We initialize the generator with x_n = s_n MOD m for n = 1..5,
   where s_n = (69069 * s_{n-1}) mod 2^32, and s_0 = s is the
   user-supplied seed.

   NOTE: According to the paper the seeds must lie in the range [0,
   2^31 - 2] with at least one non-zero value -- our seeding procedure
   satisfies these constraints.

   We then use 6 iterations of the generator to "warm up" the internal
   state.

   With this initialization procedure the theoretical value of
   z_{10006} is 2064828650 for s = 1. The subscript 10006 means (1)
   seed the generator with s = 1, (2) do the 6 warm-up iterations
   that are part of the seeding process, (3) then do 10000 actual
   iterations.

   The period of this generator is about 2^155.

   From: P. L'Ecuyer, F. Blouin, and R. Coutre, "A search for good
   multiple recursive random number generators", ACM Transactions on
   Modeling and Computer Simulation 3, 87-98 (1993). */

static inline unsigned long int mrg_get (void *vstate);
static double mrg_get_double (void *vstate);
static void mrg_set (void *state, unsigned long int s);

static const long int m = 2147483647;
static const long int a1 = 107374182, q1 = 20, r1 = 7;
static const long int a5 = 104480, q5 = 20554, r5 = 1727;

typedef struct
  {
    long int x1, x2, x3, x4, x5;
  }
mrg_state_t;

static inline unsigned long int
mrg_get (void *vstate)
{
  mrg_state_t *state = (mrg_state_t *) vstate;

  long int p1, h1, p5, h5;

  h5 = state->x5 / q5;
  p5 = a5 * (state->x5 - h5 * q5) - h5 * r5;
  if (p5 > 0)
    p5 -= m;

  h1 = state->x1 / q1;
  p1 = a1 * (state->x1 - h1 * q1) - h1 * r1;
  if (p1 < 0)
    p1 += m;

  state->x5 = state->x4;
  state->x4 = state->x3;
  state->x3 = state->x2;
  state->x2 = state->x1;

  state->x1 = p1 + p5;

  if (state->x1 < 0)
    state->x1 += m;

  return state->x1;
}

static double
mrg_get_double (void *vstate)
{
  return mrg_get (vstate) / 2147483647.0 ;
}

static void
mrg_set (void *vstate, unsigned long int s)
{
  /* An entirely adhoc way of seeding! This does **not** come from
     L'Ecuyer et al */

  mrg_state_t *state = (mrg_state_t *) vstate;

  if (s == 0)
    s = 1;	/* default seed is 1 */

#define LCG(n) ((69069 * n) & 0xffffffffUL)
  s = LCG (s);
  state->x1 = s % m;
  s = LCG (s);
  state->x2 = s % m;
  s = LCG (s);
  state->x3 = s % m;
  s = LCG (s);
  state->x4 = s % m;
  s = LCG (s);
  state->x5 = s % m;

  /* "warm it up" with at least 5 calls to go through
     all the x values */

  mrg_get (state);
  mrg_get (state);
  mrg_get (state);
  mrg_get (state);
  mrg_get (state);
  mrg_get (state);

  return;
}

static const gsl_rng_type mrg_type =
{"mrg",				/* name */
 2147483646,			/* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (mrg_state_t),
 &mrg_set,
 &mrg_get,
 &mrg_get_double};

const gsl_rng_type *gsl_rng_mrg = &mrg_type;


text:
#pragma once

#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include <gsl\gsl>

namespace Library
{
	class BlendStates final
	{
	public:
		inline static winrt::com_ptr<ID3D11BlendState> AlphaBlending;
		inline static winrt::com_ptr<ID3D11BlendState> MultiplicativeBlending;

		static void Initialize(gsl::not_null<ID3D11Device*> direct3DDevice);
		static void Shutdown();

		BlendStates() = delete;
		BlendStates(const BlendStates&) = delete;
		BlendStates& operator=(const BlendStates&) = delete;
		BlendStates(BlendStates&&) = delete;
		BlendStates& operator=(BlendStates&&) = delete;
		~BlendStates() = default;
	};
}


text:
/***************************************************************************/
/*                                                                         */
/* vector.c - Vector class for mruby                                       */
/* Copyright (C) 2015 Paolo Bosetti                                        */
/* paolo[dot]bosetti[at]unitn.it                                           */
/* Department of Industrial Engineering, University of Trento              */
/*                                                                         */
/* This library is free software.  You can redistribute it and/or          */
/* modify it under the terms of the GNU GENERAL PUBLIC LICENSE 2.0.        */
/*                                                                         */
/* This library is distributed in the hope that it will be useful,         */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of          */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the           */
/* Artistic License 2.0 for more details.                                  */
/*                                                                         */
/* See the file LICENSE                                                    */
/*                                                                         */
/***************************************************************************/

#include <gsl/gsl_blas.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_rng.h>
#include "vector.h"

#pragma mark -
#pragma mark • Utilities

// Garbage collector handler, for play_data struct
// if play_data contains other dynamic data, free it too!
// Check it with GC.start
void vector_destructor(mrb_state *mrb, void *p_) {
  gsl_vector *v = (gsl_vector *)p_;
  gsl_vector_free(v);
};

// Creating data type and reference for GC, in a const struct
const struct mrb_data_type vector_data_type = {"vector_data",
                                               vector_destructor};

// Utility function for getting the struct out of the wrapping IV @data
void mrb_vector_get_data(mrb_state *mrb, mrb_value self, gsl_vector **data) {
  mrb_value data_value;
  data_value = mrb_iv_get(mrb, self, mrb_intern_lit(mrb, "@data"));

  // Loading data from data_value into p_data:
  Data_Get_Struct(mrb, data_value, &vector_data_type, *data);
  if (!*data)
    mrb_raise(mrb, E_RUNTIME_ERROR, "Could not access @data");
}

#pragma mark -
#pragma mark • Init and accessing

// Data Initializer C function (not exposed!)
static void mrb_vector_init(mrb_state *mrb, mrb_value self, mrb_int n) {
  mrb_value data_value; // this IV holds the data
  gsl_vector *p_data;   // pointer to the C struct

  data_value = mrb_iv_get(mrb, self, mrb_intern_lit(mrb, "@data"));

  // if @data already exists, free its content:
  if (!mrb_nil_p(data_value)) {
    Data_Get_Struct(mrb, data_value, &vector_data_type, p_data);
    free(p_data);
  }
  // Allocate and zero-out the data struct:
  p_data = gsl_vector_calloc(n);
  if (!p_data)
    mrb_raise(mrb, E_RUNTIME_ERROR, "Could not allocate @data");

  // Wrap struct into @data:
  mrb_iv_set(
      mrb, self, mrb_intern_lit(mrb, "@data"), // set @data
      mrb_obj_value(                           // with value hold in struct
          Data_Wrap_Struct(mrb, mrb->object_class, &vector_data_type, p_data)));
}

static mrb_value mrb_vector_initialize(mrb_state *mrb, mrb_value self) {
  mrb_int n;
  mrb_get_args(mrb, "i", &n);

  // Call strcut initializer:
  mrb_vector_init(mrb, self, n);
  mrb_iv_set(mrb, self, mrb_intern_lit(mrb, "@length"), mrb_fixnum_value(n));
  mrb_iv_set(mrb, self, mrb_intern_lit(mrb, "@format"),
             mrb_str_new_cstr(mrb, "%10.3f"));
  return mrb_nil_value();
}

static mrb_value mrb_vector_rnd_fill(mrb_state *mrb, mrb_value self) {
  gsl_vector *p_vec = NULL;
  const gsl_rng_type *T;
  gsl_rng *r;
  mrb_int h;

  mrb_vector_get_data(mrb, self, &p_vec);

  gsl_rng_env_setup();
  T = gsl_rng_default;
  r = gsl_rng_alloc(T);

  for (h = 0; h < p_vec->size; h++) {
    gsl_vector_set(p_vec, h, gsl_rng_uniform(r));
  }
  return self;
}

static mrb_value mrb_vector_dup(mrb_state *mrb, mrb_value self) {
  mrb_value other;
  gsl_vector *p_vec = NULL, *p_vec_other = NULL;
  mrb_value args[1];

  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  args[0] = mrb_fixnum_value(p_vec->size);
  other = mrb_obj_new(mrb, mrb_class_get(mrb, "Vector"), 1, args);
  mrb_vector_get_data(mrb, other, &p_vec_other);
  gsl_vector_memcpy(p_vec_other, p_vec);
  return other;
}

static mrb_value mrb_vector_all(mrb_state *mrb, mrb_value self) {
  mrb_float v;
  gsl_vector *p_vec = NULL;

  mrb_get_args(mrb, "f", &v);
  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  gsl_vector_set_all(p_vec, v);
  return self;
}

static mrb_value mrb_vector_zero(mrb_state *mrb, mrb_value self) {
  gsl_vector *p_vec = NULL;

  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  gsl_vector_set_zero(p_vec);
  return self;
}

static mrb_value mrb_vector_basis(mrb_state *mrb, mrb_value self) {
  mrb_int i;
  gsl_vector *p_vec = NULL;

  mrb_get_args(mrb, "i", &i);
  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  gsl_vector_set_basis(p_vec, i);
  return self;
}

#pragma mark -
#pragma mark • Tests

static mrb_value mrb_vector_equal(mrb_state *mrb, mrb_value self) {
  mrb_value other;
  gsl_vector *p_vec, *p_vec_other;
  mrb_get_args(mrb, "o", &other);

  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  mrb_vector_get_data(mrb, other, &p_vec_other);
  if (1 == gsl_vector_equal(p_vec, p_vec_other))
    return mrb_true_value();
  else
    return mrb_false_value();
}

#pragma mark -
#pragma mark • Accessors

static mrb_value mrb_vector_get_i(mrb_state *mrb, mrb_value self) {
  mrb_int i = 0;
  gsl_vector *p_vec = NULL;

  mrb_get_args(mrb, "i", &i);
  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  if (i >= p_vec->size) {
    mrb_raise(mrb, E_VECTOR_ERROR, "Vector index out of range!");
  }
  return mrb_float_value(mrb, gsl_vector_get(p_vec, (size_t)i));
}

static mrb_value mrb_vector_set_i(mrb_state *mrb, mrb_value self) {
  mrb_int i = 0;
  mrb_float f;
  gsl_vector *p_vec = NULL;

  mrb_get_args(mrb, "if", &i, &f);

  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  if (i >= p_vec->size) {
    mrb_raise(mrb, E_VECTOR_ERROR, "Vector index out of range!");
  }
  gsl_vector_set(p_vec, (size_t)i, (double)f);
  return mrb_float_value(mrb, f);
}

static mrb_value mrb_vector_to_a(mrb_state *mrb, mrb_value self) {
  int i;
  mrb_value ary = mrb_nil_value();
  gsl_vector *p_vec = NULL;
  mrb_float e;
  mrb_vector_get_data(mrb, self, &p_vec);
  ary = mrb_ary_new_capa(mrb, p_vec->size);
  for (i = 0; i < p_vec->size; i++) {
    e = *(p_vec->data + i * p_vec->stride);
    mrb_ary_set(mrb, ary, i, mrb_float_value(mrb, e));
  }
  return ary;
}

#pragma mark -
#pragma mark • Properties

static mrb_value mrb_vector_max(mrb_state *mrb, mrb_value self) {
  gsl_vector *p_vec = NULL;
  mrb_vector_get_data(mrb, self, &p_vec);
  return mrb_float_value(mrb, gsl_vector_max(p_vec));
}

static mrb_value mrb_vector_min(mrb_state *mrb, mrb_value self) {
  gsl_vector *p_vec = NULL;
  mrb_vector_get_data(mrb, self, &p_vec);
  return mrb_float_value(mrb, gsl_vector_min(p_vec));
}

static mrb_value mrb_vector_max_index(mrb_state *mrb, mrb_value self) {
  gsl_vector *p_vec = NULL;
  mrb_vector_get_data(mrb, self, &p_vec);
  return mrb_fixnum_value(gsl_vector_max_index(p_vec));
}

static mrb_value mrb_vector_min_index(mrb_state *mrb, mrb_value self) {
  gsl_vector *p_vec = NULL;
  mrb_vector_get_data(mrb, self, &p_vec);
  return mrb_fixnum_value(gsl_vector_min_index(p_vec));
}

#pragma mark -
#pragma mark • Operations

static mrb_value mrb_vector_add(mrb_state *mrb, mrb_value self) {
  mrb_value other;
  gsl_vector *p_vec, *p_vec_other;
  mrb_get_args(mrb, "o", &other);

  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);

  if (mrb_obj_is_kind_of(mrb, other, mrb_class_get(mrb, "Vector"))) {
    mrb_vector_get_data(mrb, other, &p_vec_other);
    if (p_vec->size != p_vec_other->size) {
      mrb_raise(mrb, E_VECTOR_ERROR, "Vector indexes don't match!");
    }
    gsl_vector_add(p_vec, p_vec_other);
  } else if (mrb_obj_is_kind_of(mrb, other, mrb_class_get(mrb, "Numeric"))) {
    gsl_vector_add_constant(p_vec, mrb_to_flo(mrb, other));
  }
  return self;
}

static mrb_value mrb_vector_sub(mrb_state *mrb, mrb_value self) {
  mrb_value other;
  gsl_vector *p_vec, *p_vec_other;
  mrb_get_args(mrb, "o", &other);

  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  mrb_vector_get_data(mrb, other, &p_vec_other);
  if (p_vec->size != p_vec_other->size) {
    mrb_raise(mrb, E_VECTOR_ERROR, "Vector dimensions don't match!");
  }
  gsl_vector_sub(p_vec, p_vec_other);
  return self;
}

static mrb_value mrb_vector_mul(mrb_state *mrb, mrb_value self) {
  mrb_value other;
  gsl_vector *p_vec, *p_vec_other;
  mrb_get_args(mrb, "o", &other);

  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  if (mrb_obj_is_kind_of(mrb, other, mrb_class_get(mrb, "Vector"))) {
    mrb_vector_get_data(mrb, other, &p_vec_other);
    if (p_vec->size != p_vec_other->size) {
      mrb_raise(mrb, E_VECTOR_ERROR, "Vector indexes don't match!");
    }
    gsl_vector_mul(p_vec, p_vec_other);
  } else if (mrb_obj_is_kind_of(mrb, other, mrb_class_get(mrb, "Numeric"))) {
    gsl_vector_scale(p_vec, mrb_to_flo(mrb, other));
  }
  return self;
}

static mrb_value mrb_vector_div(mrb_state *mrb, mrb_value self) {
  mrb_value other;
  gsl_vector *p_vec, *p_vec_other;
  mrb_get_args(mrb, "o", &other);

  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  mrb_vector_get_data(mrb, other, &p_vec_other);
  if (p_vec->size != p_vec_other->size) {
    mrb_raise(mrb, E_VECTOR_ERROR, "Vector indexes don't match!");
  }
  gsl_vector_div(p_vec, p_vec_other);
  return self;
}

static mrb_value mrb_vector_prod(mrb_state *mrb, mrb_value self) {
  mrb_value other;
  gsl_vector *p_vec, *p_vec_other;
  mrb_float res;
  mrb_get_args(mrb, "o", &other);

  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  mrb_vector_get_data(mrb, other, &p_vec_other);
  if (!mrb_obj_is_kind_of(mrb, other, mrb_class_get(mrb, "Vector"))) {
    mrb_raise(mrb, E_ARGUMENT_ERROR, "Need a Vector!");
  }
  if (p_vec->size != p_vec_other->size) {
    mrb_raise(mrb, E_VECTOR_ERROR, "Vector indexes don't match!");
  }
  if (gsl_blas_ddot(p_vec, p_vec_other, &res)) {
    mrb_raise(mrb, E_VECTOR_ERROR, "Cannot multiply");
  }
  return mrb_float_value(mrb, res);
}

static mrb_value mrb_vector_norm(mrb_state *mrb, mrb_value self) {
  gsl_vector *p_vec;

  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  return mrb_float_value(mrb, gsl_blas_dnrm2(p_vec));
}

static mrb_value mrb_vector_sum(mrb_state *mrb, mrb_value self) {
  gsl_vector *p_vec;

  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  return mrb_float_value(mrb, gsl_blas_dasum(p_vec));
}

static mrb_value mrb_vector_swap(mrb_state *mrb, mrb_value self) {
  gsl_vector *p_vec;
  mrb_int i, j;
  mrb_get_args(mrb, "ii", &i, &j);

  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  if (gsl_vector_swap_elements(p_vec, i, j)) {
    mrb_raise(mrb, E_VECTOR_ERROR, "Cannot swap");
  }
  return self;
}

static mrb_value mrb_vector_reverse(mrb_state *mrb, mrb_value self) {
  gsl_vector *p_vec;

  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  if (gsl_vector_reverse(p_vec)) {
    mrb_raise(mrb, E_VECTOR_ERROR, "Cannot reverse");
  }
  return self;
}

#pragma mark -
#pragma mark • Statistics

static mrb_value mrb_vector_mean(mrb_state *mrb, mrb_value self) {
  gsl_vector *p_vec;
  mrb_float result;
  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  result = gsl_stats_mean(p_vec->data, p_vec->stride, p_vec->size);
  return mrb_float_value(mrb, result);
}

static mrb_value mrb_vector_variance(mrb_state *mrb, mrb_value self) {
  gsl_vector *p_vec;
  mrb_float result;
  mrb_float m;
  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  if (mrb_get_args(mrb, "|f", &m) == 1) {
    result = gsl_stats_variance_m(p_vec->data, p_vec->stride, p_vec->size, m);
  } else {
    result = gsl_stats_variance(p_vec->data, p_vec->stride, p_vec->size);
  }
  return mrb_float_value(mrb, result);
}

static mrb_value mrb_vector_sd(mrb_state *mrb, mrb_value self) {
  gsl_vector *p_vec;
  mrb_float result;
  mrb_float m;
  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  if (mrb_get_args(mrb, "|f", &m) == 1) {
    result = gsl_stats_sd_m(p_vec->data, p_vec->stride, p_vec->size, m);
  } else {
    result = gsl_stats_sd(p_vec->data, p_vec->stride, p_vec->size);
  }
  return mrb_float_value(mrb, result);
}

static mrb_value mrb_vector_absdev(mrb_state *mrb, mrb_value self) {
  gsl_vector *p_vec;
  mrb_float result;
  mrb_float m;
  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  if (mrb_get_args(mrb, "|f", &m) == 1) {
    result = gsl_stats_absdev_m(p_vec->data, p_vec->stride, p_vec->size, m);
  } else {
    result = gsl_stats_absdev(p_vec->data, p_vec->stride, p_vec->size);
  }
  return mrb_float_value(mrb, result);
}

static mrb_value mrb_vector_quantile(mrb_state *mrb, mrb_value self) {
  gsl_vector *p_vec = NULL, *p_sort_vec = NULL;
  mrb_float result;
  mrb_float f;
  mrb_int n;
  // call utility for unwrapping @data into p_data:
  mrb_vector_get_data(mrb, self, &p_vec);
  n = mrb_get_args(mrb, "|f", &f);
  if (f < 0 || f > 1) {
    mrb_raise(mrb, E_VECTOR_ERROR, "Quantile must be in [0,1]");
  }
  p_sort_vec = gsl_vector_calloc(p_vec->size);
  if (gsl_vector_memcpy(p_sort_vec, p_vec)) {
    mrb_raise(mrb, E_VECTOR_ERROR, "Cannot copy vector");
  }
  gsl_sort_vector(p_sort_vec);
  if (n == 1) {
    result = gsl_stats_quantile_from_sorted_data(
        p_sort_vec->data, p_sort_vec->stride, p_sort_vec->size, f);
  } else {
    result = gsl_stats_quantile_from_sorted_data(
        p_sort_vec->data, p_sort_vec->stride, p_sort_vec->size, 0.5);
  }
  return mrb_float_value(mrb, result);
}

#pragma mark -
#pragma mark • Gem setup

void mrb_gsl_vector_init(mrb_state *mrb) {
  struct RClass *gsl;

  mrb_load_string(mrb, "class VectorError < Exception; end");

  gsl = mrb_define_class(mrb, "Vector", mrb->object_class);
  mrb_define_method(mrb, gsl, "all", mrb_vector_all, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, gsl, "zero", mrb_vector_zero, MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "basis", mrb_vector_basis, MRB_ARGS_REQ(1));

  mrb_define_method(mrb, gsl, "initialize", mrb_vector_initialize,
                    MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "rnd_fill", mrb_vector_rnd_fill,
                    MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "dup", mrb_vector_dup, MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "===", mrb_vector_equal, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, gsl, "[]", mrb_vector_get_i, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, gsl, "[]=", mrb_vector_set_i, MRB_ARGS_REQ(2));
  mrb_define_method(mrb, gsl, "to_a", mrb_vector_to_a, MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "max", mrb_vector_max, MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "min", mrb_vector_min, MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "max_index", mrb_vector_max_index,
                    MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "min_index", mrb_vector_min_index,
                    MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "add!", mrb_vector_add, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, gsl, "sub!", mrb_vector_sub, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, gsl, "mul!", mrb_vector_mul, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, gsl, "div!", mrb_vector_div, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, gsl, "^", mrb_vector_prod, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, gsl, "norm", mrb_vector_norm, MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "sum", mrb_vector_sum, MRB_ARGS_NONE());
  mrb_define_method(mrb, gsl, "swap!", mrb_vector_swap, MRB_ARGS_REQ(2));
  mrb_define_method(mrb, gsl, "reverse!", mrb_vector_reverse, MRB_ARGS_NONE());

  mrb_define_method(mrb, gsl, "mean", mrb_vector_mean, MRB_ARGS_OPT(1));
  mrb_define_method(mrb, gsl, "variance", mrb_vector_variance, MRB_ARGS_OPT(1));
  mrb_define_method(mrb, gsl, "sd", mrb_vector_sd, MRB_ARGS_OPT(1));
  mrb_define_method(mrb, gsl, "absdev", mrb_vector_absdev, MRB_ARGS_OPT(1));
  mrb_define_method(mrb, gsl, "quantile", mrb_vector_quantile, MRB_ARGS_OPT(1));
}

max_stars_repo_name: UniTN-Mechatronics/mruby-gsl

max_issues_repo_name: UniTN-Mechatronics/mruby-gsl

max_forks_repo_name: UniTN-Mechatronics/mruby-gsl


text:
/* gsl_sf_hermite.h
 *
 * Copyright (C) 2011-2014 Konrad Griessinger
 *
 * 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.
 */

/*----------------------------------------------------------------------*
 * (konradg(at)gmx.net)                                                 *
 *----------------------------------------------------------------------*/

#ifndef __GSL_SF_HERMITE_H__
#define __GSL_SF_HERMITE_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

int gsl_sf_hermite_prob_e(const int n, const double x, gsl_sf_result * result);
double gsl_sf_hermite_prob(const int n, const double x);
int gsl_sf_hermite_prob_deriv_e(const int m, const int n, const double x, gsl_sf_result * result);
double gsl_sf_hermite_prob_deriv(const int m, const int n, const double x);
int gsl_sf_hermite_e(const int n, const double x, gsl_sf_result * result);
double gsl_sf_hermite(const int n, const double x);
int gsl_sf_hermite_deriv_e(const int m, const int n, const double x, gsl_sf_result * result);
double gsl_sf_hermite_deriv(const int m, const int n, const double x);
int gsl_sf_hermite_func_e(const int n, const double x, gsl_sf_result * result);
double gsl_sf_hermite_func(const int n, const double x);
int gsl_sf_hermite_func_fast_e(const int n, const double x, gsl_sf_result * result);
double gsl_sf_hermite_func_fast(const int n, const double x);
int gsl_sf_hermite_prob_array(const int nmax, const double x, double * result_array);
int gsl_sf_hermite_prob_array_deriv(const int m, const int nmax, const double x, double * result_array);
int gsl_sf_hermite_prob_deriv_array(const int mmax, const int n, const double x, double * result_array);
int gsl_sf_hermite_prob_series_e(const int n, const double x, const double * a, gsl_sf_result * result);
double gsl_sf_hermite_prob_series(const int n, const double x, const double * a);
int gsl_sf_hermite_array(const int nmax, const double x, double * result_array);
int gsl_sf_hermite_array_deriv(const int m, const int nmax, const double x, double * result_array);
int gsl_sf_hermite_deriv_array(const int mmax, const int n, const double x, double * result_array);
int gsl_sf_hermite_series_e(const int n, const double x, const double * a, gsl_sf_result * result);
double gsl_sf_hermite_series(const int n, const double x, const double * a);
int gsl_sf_hermite_func_array(const int nmax, const double x, double * result_array);
int gsl_sf_hermite_func_series_e(const int n, const double x, const double * a, gsl_sf_result * result);
double gsl_sf_hermite_func_series(const int n, const double x, const double * a);
int gsl_sf_hermite_func_der_e(const int m, const int n, const double x, gsl_sf_result * result);
double gsl_sf_hermite_func_der(const int m, const int n, const double x);
int gsl_sf_hermite_prob_zero_e(const int n, const int s, gsl_sf_result * result);
double gsl_sf_hermite_prob_zero(const int n, const int s);
int gsl_sf_hermite_zero_e(const int n, const int s, gsl_sf_result * result);
double gsl_sf_hermite_zero(const int n, const int s);
int gsl_sf_hermite_func_zero_e(const int n, const int s, gsl_sf_result * result);
double gsl_sf_hermite_func_zero(const int n, const int s);

#ifndef GSL_DISABLE_DEPRECATED

int gsl_sf_hermite_phys_e(const int n, const double x, gsl_sf_result * result);
double gsl_sf_hermite_phys(const int n, const double x);
int gsl_sf_hermite_phys_der_e(const int m, const int n, const double x, gsl_sf_result * result);
double gsl_sf_hermite_phys_der(const int m, const int n, const double x);
int gsl_sf_hermite_phys_array(const int nmax, const double x, double * result_array);
int gsl_sf_hermite_phys_series_e(const int n, const double x, const double * a, gsl_sf_result * result);
double gsl_sf_hermite_phys_series(const int n, const double x, const double * a);
int gsl_sf_hermite_phys_array_der(const int m, const int nmax, const double x, double * result_array);
int gsl_sf_hermite_phys_der_array(const int mmax, const int n, const double x, double * result_array);
int gsl_sf_hermite_phys_zero_e(const int n, const int s, gsl_sf_result * result);
double gsl_sf_hermite_phys_zero(const int n, const int s);

int gsl_sf_hermite_prob_array_der(const int m, const int nmax, const double x, double * result_array);
int gsl_sf_hermite_prob_der_array(const int mmax, const int n, const double x, double * result_array);
int gsl_sf_hermite_prob_der_e(const int m, const int n, const double x, gsl_sf_result * result);
double gsl_sf_hermite_prob_der(const int m, const int n, const double x);

#endif /* !GSL_DISABLE_DEPRECATED */

__END_DECLS

#endif /* __GSL_SF_HERMITE_H__ */

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: ielomariala/Hex-Game


text:
// matrix/kaldi-blas.h

// Copyright 2009-2011  Ondrej Glembek;  Microsoft Corporation

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

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

// THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
// KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
// WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
// MERCHANTABLITY OR NON-INFRINGEMENT.
// See the Apache 2 License for the specific language governing permissions and
// limitations under the License.
#ifndef KALDI_MATRIX_KALDI_BLAS_H_
#define KALDI_MATRIX_KALDI_BLAS_H_

// This file handles the #includes for BLAS, LAPACK and so on.
// It manipulates the declarations into a common format that kaldi can handle.
// However, the kaldi code will check whether HAVE_ATLAS is defined as that
// code is called a bit differently from CLAPACK that comes from other sources.

// There are three alternatives:
//   (i) you have ATLAS, which includes the ATLAS implementation of CBLAS
//   plus a subset of CLAPACK (but with clapack_ in the function declarations).
//   In this case, define HAVE_ATLAS and make sure the relevant directories are
//   in the include path.

//   (ii) you have CBLAS (some implementation thereof) plus CLAPACK.
//   In this case, define HAVE_CLAPACK.
//   [Since CLAPACK depends on BLAS, the presence of BLAS is implicit].

//  (iii) you have the MKL library, which includes CLAPACK and CBLAS.

// Note that if we are using ATLAS, no Svd implementation is supplied,
// so we define HAVE_Svd to be zero and this directs our implementation to
// supply its own "by hand" implementation which is based on TNT code.

#if (defined(HAVE_CLAPACK) && (defined(HAVE_ATLAS) || defined(HAVE_MKL))) \
    || (defined(HAVE_ATLAS) && defined(HAVE_MKL))
#error "Do not define more than one of HAVE_CLAPACK, HAVE_ATLAS and HAVE_MKL"
#endif

#ifdef HAVE_ATLAS
  extern "C" {
    #include <cblas.h>
    #include <clapack.h>
  }
#elif defined(HAVE_CLAPACK)
  #ifdef __APPLE__
    #include <Accelerate/Accelerate.h>
    typedef __CLPK_integer          integer;
    typedef __CLPK_logical          logical;
    typedef __CLPK_real             real;
    typedef __CLPK_doublereal       doublereal;
    typedef __CLPK_complex          complex;
    typedef __CLPK_doublecomplex    doublecomplex;
    typedef __CLPK_ftnlen           ftnlen;
  #else
    extern "C" {
      // May be in /usr/[local]/include if installed; else this uses the one
      // from the external/CLAPACK-* directory.
      #include <cblas.h>
      #include <f2c.h>
      #include <clapack.h>

      // get rid of macros from f2c.h -- these are dangerous.
      #undef abs
      #undef dabs
      #undef min
      #undef max
      #undef dmin
      #undef dmax
      #undef bit_test
      #undef bit_clear
      #undef bit_set
    }
  #endif
#elif defined(HAVE_MKL)
  extern "C" {
    #include <mkl.h>
  }
#else
  #error "You need to define (using the preprocessor) either HAVE_CLAPACK or HAVE_ATLAS or HAVE_MKL (but not more than one)"
#endif

#ifdef HAVE_CLAPACK
typedef integer KaldiBlasInt;
#endif
#ifdef HAVE_MKL
typedef MKL_INT KaldiBlasInt;
#endif

#ifdef HAVE_ATLAS
// in this case there is no need for KaldiBlasInt-- this typedef is only needed
// for Svd code which is not included in ATLAS (we re-implement it).
#endif

#endif  // KALDI_MATRIX_KALDI_BLAS_H_

max_stars_repo_name: hihihippp/Kaldi

max_issues_repo_name: UdyanSachdev/kaldi

max_forks_repo_name: UdyanSachdev/kaldi


text:
/* rng/ran3.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007, 2010 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 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_rng.h>

/* This is an implementation of the algorithm used in Knuths's
   subtractive generator, with the Numerical Recipe's ran3 paramters.
   It is a subtractive lagged fibonnaci generator. */

static inline unsigned long int ran3_get (void *vstate);
static double ran3_get_double (void *vstate);
static void ran3_set (void *state, unsigned long int s);

#define M_BIG 1000000000
#define M_SEED 161803398

typedef struct
  {
    unsigned int x;
    unsigned int y;
    unsigned long int buffer[56];
  }
ran3_state_t;

static inline unsigned long int
ran3_get (void *vstate)
{
  ran3_state_t *state = (ran3_state_t *) vstate;
  long int j;

  state->x++;

  if (state->x == 56)
    state->x = 1;

  state->y++;

  if (state->y == 56)
    state->y = 1;

  j = state->buffer[state->x] - state->buffer[state->y];

  if (j < 0)
    j += M_BIG;

  state->buffer[state->x] = j;

  return j;
}

static double
ran3_get_double (void *vstate)
{
  return ran3_get (vstate) / (double) M_BIG ;
}

static void
ran3_set (void *vstate, unsigned long int s)
{
  ran3_state_t *state = (ran3_state_t *) vstate;
  int i, i1;
  long int j, k;

  if (s == 0)
    s = 1;      /* default seed is 1 */

  j = (M_SEED - s) % M_BIG;

  /* Avoid potential problem with negative values */
  if (j < 0) j += M_BIG;

  /* the zeroth element is never used, but we initialize it for
     consistency between states */

  state->buffer[0] = 0;

  state->buffer[55] = j;

  k = 1;
  for (i = 1; i < 55; i++)
    {
      int n = (21 * i) % 55;
      state->buffer[n] = k;
      k = j - k;
      if (k < 0)
        k += M_BIG;
      j = state->buffer[n];

    }

  for (i1 = 0; i1 < 4; i1++)
    {
      for (i = 1; i < 56; i++)
        {
          long int t = state->buffer[i] - state->buffer[1 + (i + 30) % 55];
          if (t < 0)
            t += M_BIG;
          state->buffer[i] = t;
        }
    }

  state->x = 0;
  state->y = 31;

  return;
}

static const gsl_rng_type ran3_type =
{"ran3",                        /* name */
 M_BIG,                         /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (ran3_state_t),
 &ran3_set,
 &ran3_get,
 &ran3_get_double};

const gsl_rng_type *gsl_rng_ran3 = &ran3_type;

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#ifndef ENVIRONMENT_INCLUDE
#define ENVIRONMENT_INCLUDE

#include <map>
#include <string>
#include <vector>

#include <gsl/string_span>

namespace execHelper {
namespace config {
using EnvArg = std::string;
using EnvArgs = std::vector<EnvArg>;

using EnvironmentCollection = std::map<std::string, EnvArg>;
using EnvironmentValue = std::pair<std::string, EnvArg>;

static const gsl::czstring<> ENVIRONMENT_KEY = "environment";

} // namespace config
} // namespace execHelper

#endif /* ENVIRONMENT_INCLUDE */

max_stars_repo_name: exec-helper/source

max_issues_repo_name: exec-helper/source

max_forks_repo_name: exec-helper/source