text:
/**
 *
 * @file testing_dposv.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Bilel Hadri, Hatem Ltaief
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:45:18 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>
#include "testing_dmain.h"

enum blas_order_type {
            blas_rowmajor = 101,
            blas_colmajor = 102 };

enum blas_cmach_type {
            blas_base      = 151,
            blas_t         = 152,
            blas_rnd       = 153,
            blas_ieee      = 154,
            blas_emin      = 155,
            blas_emax      = 156,
            blas_eps       = 157,
            blas_prec      = 158,
            blas_underflow = 159,
            blas_overflow  = 160,
            blas_sfmin     = 161};

enum blas_norm_type {
            blas_one_norm       = 171,
            blas_real_one_norm  = 172,
            blas_two_norm       = 173,
            blas_frobenius_norm = 174,
            blas_inf_norm       = 175,
            blas_real_inf_norm  = 176,
            blas_max_norm       = 177,
            blas_real_max_norm  = 178 };

static void
BLAS_error(char *rname, int err, int val, int x) {
  fprintf( stderr, "%s %d %d %d\n", rname, err, val, x );
  abort();
}

static
void
BLAS_dge_norm(enum blas_order_type order, enum blas_norm_type norm,
  int m, int n, const double *a, int lda, double *res) {
  int i, j; float anorm, v;
  char rname[] = "BLAS_dge_norm";

  if (order != blas_colmajor) BLAS_error( rname, -1, order, 0 );

  if (norm == blas_frobenius_norm) {
    anorm = 0.0f;
    for (j = n; j; --j) {
      for (i = m; i; --i) {
        v = a[0];
        anorm += v * v;
        a++;
      }
      a += lda - m;
    }
    anorm = sqrt( anorm );
  } else if (norm == blas_inf_norm) {
    anorm = 0.0f;
    for (i = 0; i < m; ++i) {
      v = 0.0f;
      for (j = 0; j < n; ++j) {
        v += fabs( a[i + j * lda] );
      }
      if (v > anorm)
        anorm = v;
    }
  } else {
    BLAS_error( rname, -2, norm, 0 );
    return;
  }

  if (res) *res = anorm;
}

static
double
BLAS_dpow_di(double x, int n) {
  double rv = 1.0;

  if (n < 0) {
    n = -n;
    x = 1.0 / x;
  }

  for (; n; n >>= 1, x *= x) {
    if (n & 1)
      rv *= x;
  }

  return rv;
}

static
double
BLAS_dfpinfo(enum blas_cmach_type cmach) {
  double eps = 1.0, r = 1.0, o = 1.0, b = 2.0;
  int t = 53, l = 1024, m = -1021;
  char rname[] = "BLAS_dfpinfo";

  if ((sizeof eps) == sizeof(float)) {
    t = 24;
    l = 128;
    m = -125;
  } else {
    t = 53;
    l = 1024;
    m = -1021;
  }

  /* for (i = 0; i < t; ++i) eps *= half; */
  eps = BLAS_dpow_di( b, -t );
  /* for (i = 0; i >= m; --i) r *= half; */
  r = BLAS_dpow_di( b, m-1 );

  o -= eps;
  /* for (i = 0; i < l; ++i) o *= b; */
  o = (o * BLAS_dpow_di( b, l-1 )) * b;

  switch (cmach) {
    case blas_eps: return eps;
    case blas_sfmin: return r;
    default:
      BLAS_error( rname, -1, cmach, 0 );
      break;
  }
  return 0.0;
}

static int check_factorization(int, double*, double*, int, int , double);
static int check_solution(int, int, double*, int, double*, double*, int, double);
static int check_estimator(PLASMA_enum, int, double *, int,
                           double *, double, double,
                           double);

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

    /* Check for number of arguments*/
    if (argc != 4){
        USAGE("POSV", "N LDA NRHS LDB",
              "   - N    : the size of the matrix\n"
              "   - LDA  : leading dimension of the matrix A\n"
              "   - NRHS : number of RHS\n"
              "   - LDB  : leading dimension of the RHS B\n");
        return -1;
    }

    int N     = atoi(argv[0]);
    int LDA   = atoi(argv[1]);
    int NRHS  = atoi(argv[2]);
    int LDB   = atoi(argv[3]);
    double eps;
    int info_solution, info_factorization;
    int u, trans1, trans2;

    double *A1   = (double *)malloc(LDA*N*sizeof(double));
    double *A2   = (double *)malloc(LDA*N*sizeof(double));
    double *B1   = (double *)malloc(LDB*NRHS*sizeof(double));
    double *B2   = (double *)malloc(LDB*NRHS*sizeof(double));

    /* Check if unable to allocate memory */
    if ((!A1)||(!A2)||(!B1)||(!B2)){
        printf("Out of Memory \n ");
        return -2;
    }

    eps = BLAS_dfpinfo( blas_eps );

    for(u=0; u<2; u++) {

        trans1 = uplo[u] == PlasmaUpper ? PlasmaTrans : PlasmaNoTrans;
        trans2 = uplo[u] == PlasmaUpper ? PlasmaNoTrans : PlasmaTrans;

        /*-------------------------------------------------------------
         *  TESTING DPOSV
         */

        /* Initialize A1 and A2 for Symmetric Positif Matrix */
        PLASMA_dplgsy( (double)N, N, A1, LDA, 51 );
        PLASMA_dlacpy( PlasmaUpperLower, N, N, A1, LDA, A2, LDA );

        /* Initialize B1 and B2 */
        PLASMA_dplrnt( N, NRHS, B1, LDB, 371 );
        PLASMA_dlacpy( PlasmaUpperLower, N, NRHS, B1, LDB, B2, LDB );

        printf("\n");
        printf("------ TESTS FOR PLASMA DPOSV ROUTINE -------  \n");
        printf("            Size of the Matrix %d by %d\n", N, N);
        printf("\n");
        printf(" The matrix A is randomly generated for each test.\n");
        printf("============\n");
        printf(" The relative machine precision (eps) is to be %e \n", eps);
        printf(" Computational tests pass if scaled residuals are less than 60.\n");

        /* PLASMA DPOSV */
        PLASMA_dposv(uplo[u], N, NRHS, A2, LDA, B2, LDB);

        /* Check the factorization and the solution */
        info_factorization = check_factorization( N, A1, A2, LDA, uplo[u], eps);
        info_solution = check_solution(N, NRHS, A1, LDA, B1, B2, LDB, eps);

        if ( (info_solution == 0) && (info_factorization == 0) ) {
            printf("***************************************************\n");
            printf(" ---- TESTING DPOSV(%s) ...................... PASSED !\n", uplostr[u]);
            printf("***************************************************\n");
        }
        else {
            printf("***************************************************\n");
            printf(" - TESTING DPOSV(%s) ... FAILED !\n", uplostr[u]);
            printf("***************************************************\n");
        }

        /*-------------------------------------------------------------
         *  TESTING DPOTRF + DPOTRS
         */

        /* Initialize A1 and A2 for Symmetric Positif Matrix */
        PLASMA_dplgsy( (double)N, N, A1, LDA, 51 );
        PLASMA_dlacpy( PlasmaUpperLower, N, N, A1, LDA, A2, LDA );

        /* Initialize B1 and B2 */
        PLASMA_dplrnt( N, NRHS, B1, LDB, 371 );
        PLASMA_dlacpy( PlasmaUpperLower, N, NRHS, B1, LDB, B2, LDB );

        /* Plasma routines */
        PLASMA_dpotrf(uplo[u], N, A2, LDA);
        PLASMA_dpotrs(uplo[u], N, NRHS, A2, LDA, B2, LDB);

        printf("\n");
        printf("------ TESTS FOR PLASMA DPOTRF + DPOTRS ROUTINE -------  \n");
        printf("            Size of the Matrix %d by %d\n", N, N);
        printf("\n");
        printf(" The matrix A is randomly generated for each test.\n");
        printf("============\n");
        printf(" The relative machine precision (eps) is to be %e \n", eps);
        printf(" Computational tests pass if scaled residuals are less than 60.\n");

        /* Check the factorization and the solution */
        info_factorization = check_factorization( N, A1, A2, LDA, uplo[u], eps);
        info_solution = check_solution(N, NRHS, A1, LDA, B1, B2, LDB, eps);

        if ((info_solution == 0)&(info_factorization == 0)){
            printf("***************************************************\n");
            printf(" ---- TESTING DPOTRF + DPOTRS (%s)............ PASSED !\n", uplostr[u]);
            printf("***************************************************\n");
        }
        else{
            printf("****************************************************\n");
            printf(" - TESTING DPOTRF + DPOTRS (%s)... FAILED !\n", uplostr[u]);
            printf("****************************************************\n");
        }

        /*-------------------------------------------------------------
         *  TESTING DPOTRF + ZPTRSM + DTRSM
         */

        /* Initialize A1 and A2 for Symmetric Positif Matrix */
        PLASMA_dplgsy( (double)N, N, A1, LDA, 51 );
        PLASMA_dlacpy( PlasmaUpperLower, N, N, A1, LDA, A2, LDA );

        /* Initialize B1 and B2 */
        PLASMA_dplrnt( N, NRHS, B1, LDB, 371 );
        PLASMA_dlacpy( PlasmaUpperLower, N, NRHS, B1, LDB, B2, LDB );

        /* PLASMA routines */
        PLASMA_dpotrf(uplo[u], N, A2, LDA);
        PLASMA_dtrsm(PlasmaLeft, uplo[u], trans1, PlasmaNonUnit,
                     N, NRHS, 1.0, A2, LDA, B2, LDB);
        PLASMA_dtrsm(PlasmaLeft, uplo[u], trans2, PlasmaNonUnit,
                     N, NRHS, 1.0, A2, LDA, B2, LDB);

        printf("\n");
        printf("------ TESTS FOR PLASMA DPOTRF + DTRSM + DTRSM  ROUTINE -------  \n");
        printf("            Size of the Matrix %d by %d\n", N, N);
        printf("\n");
        printf(" The matrix A is randomly generated for each test.\n");
        printf("============\n");
        printf(" The relative machine precision (eps) is to be %e \n", eps);
        printf(" Computational tests pass if scaled residuals are less than 60.\n");

        /* Check the factorization and the solution */
        info_factorization = check_factorization( N, A1, A2, LDA, uplo[u], eps);
        info_solution = check_solution(N, NRHS, A1, LDA, B1, B2, LDB, eps);

        if ((info_solution == 0)&(info_factorization == 0)){
            printf("***************************************************\n");
            printf(" ---- TESTING DPOTRF + DTRSM + DTRSM (%s)..... PASSED !\n", uplostr[u]);
            printf("***************************************************\n");
        }
        else{
            printf("***************************************************\n");
            printf(" - TESTING DPOTRF + DTRSM + DTRSM (%s)... FAILED !\n", uplostr[u]);
            printf("***************************************************\n");
        }

        /*-------------------------------------------------------------
         *  TESTING ZPOCON on the last call
         */
        {
            double Anorm = PLASMA_dlansy( PlasmaOneNorm, uplo[u], N, A1, LDA );
            double Acond;

            info_solution = PLASMA_dpocon(uplo[u], N, A2, LDA, Anorm, &Acond);
            if ( info_solution == 0 ) {
                info_solution = check_estimator(uplo[u], N, A1, LDA, A2, Anorm, Acond, eps);
            } else {
                printf(" PLASMA_dpocon returned info = %d\n", info_solution );
            }
            if ((info_solution == 0)){
                printf("***************************************************\n");
                printf(" ---- TESTING DPOTRF + ZPOCON (%s) ........... PASSED !\n", uplostr[u]);
                printf("***************************************************\n");
            }
            else{
                printf("**************************************************\n");
                printf(" - TESTING DPOTRF + ZPOCON (%s) ... FAILED !\n", uplostr[u]);
                printf("**************************************************\n");
            }
        }
    }

    free(A1); free(A2); free(B1); free(B2);

    return 0;
}

/*------------------------------------------------------------------------
 *  Check the factorization of the matrix A2
 */
static int check_factorization(int N, double *A1, double *A2, int LDA, int uplo, double eps)
{
    double Anorm, Rnorm;
    double alpha;
    int info_factorization;
    int i,j;

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

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

    alpha= 1.0;

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

    /* Dealing with L'L or U'U  */
    if (uplo == PlasmaUpper){
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L1, N);
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L2, N);
        cblas_dtrmm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, N, N, (alpha), L1, N, L2, N);
    }
    else{
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L1, N);
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L2, N);
        cblas_dtrmm(CblasColMajor, CblasRight, CblasLower, CblasTrans, CblasNonUnit, N, N, (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];

    BLAS_dge_norm( blas_colmajor, blas_inf_norm, N, N, Residual, N, &Rnorm );
    BLAS_dge_norm( blas_colmajor, blas_inf_norm, N, N, A1, LDA, &Anorm );

    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)) || isinf(Rnorm/(Anorm*N*eps)) || (Rnorm/(Anorm*N*eps) > 60.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 accuracy of the solution of the linear system
 */

static int check_solution(int N, int NRHS, double *A1, int LDA, double *B1, double *B2, int LDB, double eps )
{
    int info_solution;
    double Rnorm, Anorm, Xnorm, Bnorm, result;
    double alpha, beta;
    double *work = (double *)malloc(N*sizeof(double));

    alpha = 1.0;
    beta  = -1.0;

    BLAS_dge_norm( blas_colmajor, blas_inf_norm, N, NRHS, B2, LDB, &Xnorm );
    BLAS_dge_norm( blas_colmajor, blas_inf_norm, N, N,    A1, LDA, &Anorm );
    BLAS_dge_norm( blas_colmajor, blas_inf_norm, N, NRHS, B1, LDB, &Bnorm );

    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, N, NRHS, N, (alpha), A1, LDA, B2, LDB, (beta), B1, LDB);
    BLAS_dge_norm( blas_colmajor, blas_inf_norm, N, NRHS, B1, LDB, &Rnorm );

    if (getenv("PLASMA_TESTING_VERBOSE"))
      printf( "||A||_oo=%f\n||X||_oo=%f\n||B||_oo=%f\n||A X - B||_oo=%e\n", Anorm, Xnorm, Bnorm, Rnorm );

    result = Rnorm / ( (Anorm*Xnorm+Bnorm)*N*eps ) ;
    printf("============\n");
    printf("Checking the Residual of the solution \n");
    printf("-- ||Ax-B||_oo/((||A||_oo||x||_oo+||B||_oo).N.eps) = %e \n", result);

    if (  isnan(Xnorm) || isinf(Xnorm) || isnan(result) || isinf(result) || (result > 60.0) ) {
        printf("-- The solution is suspicious ! \n");
        info_solution = 1;
     }
    else{
        printf("-- The solution is CORRECT ! \n");
        info_solution = 0;
    }

    free(work);

    return info_solution;
}

/*------------------------------------------------------------------------
 *  Check the accuracy of the condition estimator
 */
static int check_estimator(PLASMA_enum uplo, int N, double *A1, int LDA,
                           double *A2, double Anorm, double Acond,
                           double eps)
{
    int info_solution;
    double result, Acond_lapack;
    double invcond, invcond_lapack;

    info_solution = LAPACKE_dpocon(LAPACK_COL_MAJOR, lapack_const(uplo), N, A2, LDA, Anorm, &Acond_lapack);

    if ( info_solution != 0 ) {
        printf(" PLASMA_dgecon returned info = %d\n", info_solution );
        return info_solution;
    }

    invcond_lapack = 1. / ( Acond_lapack );
    invcond        = 1. / ( Acond );

    printf("============\n");
    printf("Checking the condition number \n");
    printf("-- Acond_plasma = %e, Acond_lapack = %e \n"
           "-- Ainvcond_plasma = %e, Ainvcond_lapack = %e \n",
           Acond, Acond_lapack, invcond, invcond_lapack );

    result = fabs( Acond_lapack - Acond ) / eps;
    if ( result > 60. ) {
        printf("-- The solution is suspicious ! \n");
        info_solution = 1;
     }
    else{
        printf("-- The solution is CORRECT ! \n");
        info_solution = 0;
    }
    return info_solution;
}


text:
/////////////////////////////////////////////////////////////////////
// = NMatrix
//
// A linear algebra library for scientific computation in Ruby.
// NMatrix is part of SciRuby.
//
// NMatrix was originally inspired by and derived from NArray, by
// Masahiro Tanaka: http://narray.rubyforge.org
//
// == Copyright Information
//
// SciRuby is Copyright (c) 2010 - 2012, Ruby Science Foundation
// NMatrix is Copyright (c) 2012, Ruby Science Foundation
//
// Please see LICENSE.txt for additional copyright notices.
//
// == Contributing
//
// By contributing source code to SciRuby, you agree to be bound by
// our Contributor Agreement:
//
// * https://github.com/SciRuby/sciruby/wiki/Contributor-Agreement
//
// == cblas.c
//
// Functions in this file allow us to call CBLAS functions using
// arrays of function pointers, by ensuring that each has the same
// signature.

#ifndef CBLAS_C
# define CBLAS_C

#include <cblas.h>
#include "nmatrix.h"

//extern const enum CBLAS_ORDER;
//extern const enum CBLAS_TRANSPOSE;

inline void cblas_r32gemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p) {
  if (Order == CblasColMajor) r32gemm(TransA, TransB, p.M, p.N, p.K, p.alpha.r[0], p.A, p.lda, p.B, p.ldb, p.beta.r[0], p.C, p.ldc);
  else                        r32gemm(TransB, TransA, p.N, p.M, p.K, p.alpha.r[0], p.B, p.ldb, p.A, p.lda, p.beta.r[0], p.C, p.ldc);
}

inline void cblas_r32gemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p) {
  r32gemv(TransA, p.M, p.N, p.alpha.r[0], p.A, p.lda, p.B, p.ldb, p.beta.r[0], p.C, p.ldc);
}

inline void cblas_r64gemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p) {
  if (Order == CblasColMajor) r64gemm(TransA, TransB, p.M, p.N, p.K, p.alpha.ra[0], p.A, p.lda, p.B, p.ldb, p.beta.ra[0], p.C, p.ldc);
  else                        r64gemm(TransB, TransA, p.N, p.M, p.K, p.alpha.ra[0], p.B, p.ldb, p.A, p.lda, p.beta.ra[0], p.C, p.ldc);
}

inline void cblas_r64gemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p) {
  r64gemv(TransA, p.M, p.N, p.alpha.ra[0], p.A, p.lda, p.B, p.ldb, p.beta.ra[0], p.C, p.ldc);
}

inline void cblas_r128gemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p) {
  if (Order == CblasColMajor) r128gemm(TransA, TransB, p.M, p.N, p.K, p.alpha.rat, p.A, p.lda, p.B, p.ldb, p.beta.rat, p.C, p.ldc);
  else                        r128gemm(TransB, TransA, p.N, p.M, p.K, p.alpha.rat, p.B, p.ldb, p.A, p.lda, p.beta.rat, p.C, p.ldc);
}

inline void cblas_r128gemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p) {
  r128gemv(TransA, p.M, p.N, p.alpha.rat, p.A, p.lda, p.B, p.ldb, p.beta.rat, p.C, p.ldc);
}

inline void cblas_bgemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p) {
  bgemv(TransA, p.M, p.N, p.alpha.b[0], p.A, p.lda, p.B, p.ldb, p.beta.b[0], p.C, p.ldc);
}

inline void cblas_bgemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p) {
  if (Order == CblasColMajor) bgemm(TransA, TransB, p.M, p.N, p.K, p.alpha.b[0], p.A, p.lda, p.B, p.ldb, p.beta.b[0], p.C, p.ldc);
  else                        bgemm(TransB, TransA, p.N, p.M, p.K, p.alpha.b[0], p.B, p.ldb, p.A, p.lda, p.beta.b[0], p.C, p.ldc);
}

inline void cblas_i8gemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p) {
  i8gemv(TransA, p.M, p.N, p.alpha.i[0], p.A, p.lda, p.B, p.ldb, p.beta.i[0], p.C, p.ldc);
}

inline void cblas_i8gemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p) {
  if (Order == CblasColMajor) i8gemm(TransA, TransB, p.M, p.N, p.K, p.alpha.i[0], p.A, p.lda, p.B, p.ldb, p.beta.i[0], p.C, p.ldc);
  else                        i8gemm(TransB, TransA, p.N, p.M, p.K, p.alpha.i[0], p.B, p.ldb, p.A, p.lda, p.beta.i[0], p.C, p.ldc);
}

inline void cblas_i16gemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p) {
  i16gemv(TransA, p.M, p.N, p.alpha.i[0], p.A, p.lda, p.B, p.ldb, p.beta.i[0], p.C, p.ldc);
}

inline void cblas_i16gemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p) {
  if (Order == CblasColMajor) i16gemm(TransA, TransB, p.M, p.N, p.K, p.alpha.i[0], p.A, p.lda, p.B, p.ldb, p.beta.i[0], p.C, p.ldc);
  else                        i16gemm(TransB, TransA, p.N, p.M, p.K, p.alpha.i[0], p.B, p.ldb, p.A, p.lda, p.beta.i[0], p.C, p.ldc);
}

inline void cblas_i32gemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p) {
  i32gemv(TransA, p.M, p.N, p.alpha.i[0], p.A, p.lda, p.B, p.ldb, p.beta.i[0], p.C, p.ldc);
}

inline void cblas_i32gemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p) {
  if (Order == CblasColMajor) i32gemm(TransA, TransB, p.M, p.N, p.K, p.alpha.i[0], p.A, p.lda, p.B, p.ldb, p.beta.i[0], p.C, p.ldc);
  else                        i32gemm(TransB, TransA, p.N, p.M, p.K, p.alpha.i[0], p.B, p.ldb, p.A, p.lda, p.beta.i[0], p.C, p.ldc);
}

inline void cblas_i64gemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p) {
  i64gemv(TransA, p.M, p.N, p.alpha.i[0], p.A, p.lda, p.B, p.ldb, p.beta.i[0], p.C, p.ldc);
}

inline void cblas_i64gemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p) {
  if (Order == CblasColMajor) i64gemm(TransA, TransB, p.M, p.N, p.K, p.alpha.i[0], p.A, p.lda, p.B, p.ldb, p.beta.i[0], p.C, p.ldc);
  else                        i64gemm(TransB, TransA, p.N, p.M, p.K, p.alpha.i[0], p.B, p.ldb, p.A, p.lda, p.beta.i[0], p.C, p.ldc);
}

inline void cblas_vgemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p) {
  if (Order == CblasColMajor) vgemm(TransA, TransB, p.M, p.N, p.K, p.alpha.v[0], p.A, p.lda, p.B, p.ldb, p.beta.v[0], p.C, p.ldc);
  else                        vgemm(TransB, TransA, p.N, p.M, p.K, p.alpha.v[0], p.B, p.ldb, p.A, p.lda, p.beta.v[0], p.C, p.ldc);
}

inline void cblas_sgemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p) {
  cblas_sgemv(Order, TransA, p.M, p.N, p.alpha.d[0], p.A, p.lda, p.B, p.ldb, p.beta.d[0], p.C, p.ldc);
}

inline void cblas_sgemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p) {
  cblas_sgemm(Order, TransA, TransB, p.M, p.N, p.K, p.alpha.d[0], p.A, p.lda, p.B, p.ldb, p.beta.d[0], p.C, p.ldc);
}

inline void cblas_dgemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p) {
  cblas_dgemv(Order, TransA, p.M, p.N, p.alpha.d[0], p.A, p.lda, p.B, p.ldb, p.beta.d[0], p.C, p.ldc);
}

inline void cblas_dgemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p) {
  cblas_dgemm(Order, TransA, TransB, p.M, p.N, p.K, p.alpha.d[0], p.A, p.lda, p.B, p.ldb, p.beta.d[0], p.C, p.ldc);
}

inline void cblas_cgemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p) {
  cblas_cgemv(Order, TransA, p.M, p.N, &(p.alpha.c), p.A, p.lda, p.B, p.ldb, &(p.beta.c), p.C, p.ldc);
}

inline void cblas_cgemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p) {
  cblas_cgemm(Order, TransA, TransB, p.M, p.N, p.K, &(p.alpha.c), p.A, p.lda, p.B, p.ldb, &(p.beta.c), p.C, p.ldc);
}

inline void cblas_zgemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p) {
  cblas_zgemv(Order, TransA, p.M, p.N, &(p.alpha.z), p.A, p.lda, p.B, p.ldb, &(p.beta.z), p.C, p.ldc);
}

inline void cblas_zgemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p) {
  cblas_zgemm(Order, TransA, TransB, p.M, p.N, p.K, &(p.alpha.z), p.A, p.lda, p.B, p.ldb, &(p.beta.z), p.C, p.ldc);
}

#endif

ext: c


text:
#ifndef XYWEBSERVER_TEST_UTILS_H_
#define XYWEBSERVER_TEST_UTILS_H_

#include <gtest/gtest.h>

#include <gsl/pointers>

#include "runtime/runtime.h"

#define CO_ASSERT_EQ(val1, val2)               \
  ({                                           \
    auto f = [&]() { ASSERT_EQ(val1, val2); }; \
    f();                                       \
  })

constexpr std::chrono::milliseconds time_deviation(7);

class TestRuntimeCtxGuard {
 public:
  TestRuntimeCtxGuard(
      gsl::owner<std::function<xyco::runtime::Future<void>()> *> co_wrapper,
      bool in_runtime);

  TestRuntimeCtxGuard(const TestRuntimeCtxGuard &guard) = delete;

  TestRuntimeCtxGuard(TestRuntimeCtxGuard &&guard) = delete;

  auto operator=(const TestRuntimeCtxGuard &guard)
      -> TestRuntimeCtxGuard & = delete;

  auto operator=(TestRuntimeCtxGuard &&guard) -> TestRuntimeCtxGuard & = delete;

  ~TestRuntimeCtxGuard();

 private:
  gsl::owner<std::function<xyco::runtime::Future<void>()> *> co_wrapper_;
};

class TestRuntimeCtx {
 public:
  // run until co finished
  template <typename Fn>
  static auto co_run(Fn &&co)
      -> void requires(std::is_invocable_r_v<xyco::runtime::Future<void>, Fn>) {
    // co_outer's lifetime is managed by the caller to avoid being destroyed
    // automatically before resumed.

    std::mutex mutex;
    std::unique_lock<std::mutex> lock_guard(mutex);
    std::condition_variable cv;

    auto co_outer = [&]() -> xyco::runtime::Future<void> {
      try {
        co_await co();
        cv.notify_one();
      } catch (std::exception e) {
        auto f = [&]() { ASSERT_NO_THROW(throw e); };
        f();
        cv.notify_one();
      }
    };
    runtime_->spawn(co_outer());
    cv.wait(lock_guard);
  }

  template <typename Fn>
  static auto co_run_no_wait(Fn &&co) -> TestRuntimeCtxGuard
      requires(std::is_invocable_r_v<xyco::runtime::Future<void>, Fn>) {
    auto *co_outer = gsl::owner<std::function<xyco::runtime::Future<void>()> *>(
        new std::function<xyco::runtime::Future<void>()>(
            [=]() -> xyco::runtime::Future<void> { co_await co(); }));

    return {co_outer, true};
  }

  template <typename Fn>
  static auto co_run_without_runtime(Fn &&co) -> TestRuntimeCtxGuard
      requires(std::is_invocable_r_v<xyco::runtime::Future<void>, Fn>) {
    auto *co_outer = gsl::owner<std::function<xyco::runtime::Future<void>()> *>(
        new std::function<xyco::runtime::Future<void>()>(
            [=]() -> xyco::runtime::Future<void> { co_await co(); }));

    return {co_outer, false};
  }

 private : static std::unique_ptr<xyco::runtime::Runtime> runtime_;
};

#endif  // XYWEBSERVER_TEST_UTILS_H_

ext: h


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

#ifndef __GSL_MULTIFIT_NLIN_H__
#define __GSL_MULTIFIT_NLIN_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.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_multifit_gradient (const gsl_matrix * J, const gsl_vector * f,
                           gsl_vector * g);

int gsl_multifit_covar (const gsl_matrix * J, double epsrel, gsl_matrix * covar);

/* Definition of vector-valued functions with parameters based on gsl_vector */

struct gsl_multifit_function_struct
{
  int (* f) (const gsl_vector * x, void * params, gsl_vector * f);
  size_t n;   /* number of functions */
  size_t p;   /* number of independent variables */
  void * params;
};

typedef struct gsl_multifit_function_struct gsl_multifit_function ;

#define GSL_MULTIFIT_FN_EVAL(F,x,y) (*((F)->f))(x,(F)->params,(y))

typedef struct
  {
    const char *name;
    size_t size;
    int (*alloc) (void *state, size_t n, size_t p);
    int (*set) (void *state, gsl_multifit_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx);
    int (*iterate) (void *state, gsl_multifit_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx);
    void (*free) (void *state);
  }
gsl_multifit_fsolver_type;

typedef struct
  {
    const gsl_multifit_fsolver_type * type;
    gsl_multifit_function * function ;
    gsl_vector * x ;
    gsl_vector * f ;
    gsl_vector * dx ;
    void *state;
  }
gsl_multifit_fsolver;

gsl_multifit_fsolver *
gsl_multifit_fsolver_alloc (const gsl_multifit_fsolver_type * T,
                            size_t n, size_t p);

void gsl_multifit_fsolver_free (gsl_multifit_fsolver * s);

int gsl_multifit_fsolver_set (gsl_multifit_fsolver * s,
                                   gsl_multifit_function * f,
                                   const gsl_vector * x);

int gsl_multifit_fsolver_iterate (gsl_multifit_fsolver * s);

const char * gsl_multifit_fsolver_name (const gsl_multifit_fsolver * s);
gsl_vector * gsl_multifit_fsolver_position (const gsl_multifit_fsolver * s);

/* Definition of vector-valued functions and gradient with parameters
   based on gsl_vector */

struct gsl_multifit_function_fdf_struct
{
  int (* f) (const gsl_vector * x, void * params, gsl_vector * f);
  int (* df) (const gsl_vector * x, void * params, gsl_matrix * df);
  int (* fdf) (const gsl_vector * x, void * params, gsl_vector * f, gsl_matrix *df);
  size_t n;   /* number of functions */
  size_t p;   /* number of independent variables */
  void * params;
};

typedef struct gsl_multifit_function_fdf_struct gsl_multifit_function_fdf ;

#define GSL_MULTIFIT_FN_EVAL_F(F,x,y) ((*((F)->f))(x,(F)->params,(y)))
#define GSL_MULTIFIT_FN_EVAL_DF(F,x,dy) ((*((F)->df))(x,(F)->params,(dy)))
#define GSL_MULTIFIT_FN_EVAL_F_DF(F,x,y,dy) ((*((F)->fdf))(x,(F)->params,(y),(dy)))

typedef struct
  {
    const char *name;
    size_t size;
    int (*alloc) (void *state, size_t n, size_t p);
    int (*set) (void *state, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
    int (*iterate) (void *state, gsl_multifit_function_fdf * fdf, gsl_vector * x, gsl_vector * f, gsl_matrix * J, gsl_vector * dx);
    void (*free) (void *state);
  }
gsl_multifit_fdfsolver_type;

typedef struct
  {
    const gsl_multifit_fdfsolver_type * type;
    gsl_multifit_function_fdf * fdf ;
    gsl_vector * x;
    gsl_vector * f;
    gsl_matrix * J;
    gsl_vector * dx;
    void *state;
  }
gsl_multifit_fdfsolver;

gsl_multifit_fdfsolver *
gsl_multifit_fdfsolver_alloc (const gsl_multifit_fdfsolver_type * T,
                              size_t n, size_t p);

int
gsl_multifit_fdfsolver_set (gsl_multifit_fdfsolver * s,
                                 gsl_multifit_function_fdf * fdf,
                                 const gsl_vector * x);

int
gsl_multifit_fdfsolver_iterate (gsl_multifit_fdfsolver * s);

void
gsl_multifit_fdfsolver_free (gsl_multifit_fdfsolver * s);

const char * gsl_multifit_fdfsolver_name (const gsl_multifit_fdfsolver * s);
gsl_vector * gsl_multifit_fdfsolver_position (const gsl_multifit_fdfsolver * s);

int gsl_multifit_test_delta (const gsl_vector * dx, const gsl_vector * x,
                             double epsabs, double epsrel);

int gsl_multifit_test_gradient (const gsl_vector * g, double epsabs);

/* extern const gsl_multifit_fsolver_type * gsl_multifit_fsolver_gradient; */

GSL_VAR const gsl_multifit_fdfsolver_type * gsl_multifit_fdfsolver_lmder;
GSL_VAR const gsl_multifit_fdfsolver_type * gsl_multifit_fdfsolver_lmsder;

__END_DECLS

#endif /* __GSL_MULTIFIT_NLIN_H__ */

max_stars_repo_name: robpearc/FoldDesign

max_issues_repo_name: robpearc/FoldDesign

max_forks_repo_name: robpearc/FoldDesign


text:
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include <math.h>
#include <libgen.h>

#include <options/options.h>
#include "../csm/csm_all.h"

struct ld_noise_params {
	int seed;
	double discretization;
	double sigma;
	const char* file_input;
	const char* file_output;
	int lambertian;
};

int main(int argc, const char * argv[]) {
	sm_set_program_name(argv[0]);

	options_banner("ld_noise: Adds noise to readings in a scan");

	struct ld_noise_params p;

	struct option* ops = options_allocate(20);
	options_double(ops, "discretization", &p.discretization, 0.0,
		"Size of discretization (disabled if 0)");
	options_double(ops, "sigma", &p.sigma, 0.0,
		"Std deviation of gaussian noise (disabled if 0)");
	options_int(ops, "lambertian", &p.lambertian, 0,
		"Use lambertian model cov = sigma^2 / cos(beta^2) where beta is the incidence. Need have alpha or true_alpha.");
	options_int(ops, "seed", &p.seed, 0,
		"Seed for random number generator (if 0, use GSL_RNG_SEED env. variable).");
	options_string(ops, "in", &p.file_input, "stdin", "Input file ");
	options_string(ops, "out", &p.file_output, "stdout", "Output file ");

	if(!options_parse_args(ops, argc, argv)) {
		fprintf(stderr, "A simple program for adding noise to sensor scans.\n\nUsage:\n");
		options_print_help(ops, stderr);
		return -1;
	}

	FILE * in = open_file_for_reading(p.file_input);
	if(!in) return -3;

	FILE * out = open_file_for_writing(p.file_output);
	if(!out) return -2;

	gsl_rng_env_setup();
	gsl_rng * rng = gsl_rng_alloc (gsl_rng_ranlxs0);
	if(p.seed != 0)
	gsl_rng_set(rng, (unsigned int) p.seed);

	LDP ld; int count = 0;
	while( (ld = ld_from_json_stream(in))) {
		if(!ld_valid_fields(ld))  {
			sm_error("Invalid laser data (#%d in file)\n", count);
			continue;
		}

		int i;
		for(i=0;i<ld->nrays;i++) {
			if(!ld->valid[i]) continue;

			double * reading = ld->readings + i;
			if(p.sigma > 0) {
				double add_sigma = p.sigma;

				if(p.lambertian) {

					int have_alpha = 0;
					double alpha = 0;
					if(!is_nan(ld->true_alpha[i])) {
						alpha = ld->true_alpha[i];
						have_alpha = 1;
					} else if(ld->alpha_valid[i]) {
						alpha = ld->alpha[i];;
						have_alpha = 1;
					} else have_alpha = 0;

					if(have_alpha) {
						/* Recall that alpha points outside the surface */
						double beta = (alpha+M_PI) - ld->theta[i];
					    add_sigma = p.sigma / cos(beta);
					} else {
						sm_error("Because lambertian is active, I need either true_alpha[] or alpha[]");
						ld_write_as_json(ld, stderr);
						return -1;
					}

				}

			   *reading += gsl_ran_gaussian(rng, add_sigma);

				if(is_nan(ld->readings_sigma[i])) {
					ld->readings_sigma[i] = add_sigma;
				} else {
					ld->readings_sigma[i] = sqrt(square(add_sigma) + square(ld->readings_sigma[i]));
				}
			}
			if(p.discretization > 0)
				*reading -= fmod(*reading , p.discretization);
		}

		ld_write_as_json(ld, out);
		ld_free(ld);
	}

	return 0;
}

max_stars_repo_name: ozaslan/basics

max_issues_repo_name: ozaslan/basics

max_forks_repo_name: ozaslan/basics


text:
/* interpolation/bicubic.c
 *
 * Copyright 2012 David Zaslavsky
 *
 * 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 <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_interp2d.h>

#define IDX2D(i, j, w) ((j) * ((w)->xsize) + (i))

typedef struct
{
  double * zx;
  double * zy;
  double * zxy;
  size_t xsize;
  size_t ysize;
} bicubic_state_t;

static void bicubic_free (void * vstate);

static void *
bicubic_alloc(size_t xsize, size_t ysize)
{
  bicubic_state_t *state;

  state = calloc(1, sizeof (bicubic_state_t));

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

  state->zx = (double *) malloc (xsize * ysize * sizeof (double));
  if (state->zx == NULL)
    {
      bicubic_free(state);
      GSL_ERROR_NULL("failed to allocate space for zx", GSL_ENOMEM);
    }

  state->zy = (double *) malloc (xsize * ysize * sizeof (double));
  if (state->zy == NULL)
    {
      bicubic_free(state);
      GSL_ERROR_NULL("failed to allocate space for zy", GSL_ENOMEM);
    }

  state->zxy = (double *) malloc (xsize * ysize * sizeof (double));
  if (state->zxy == NULL)
    {
      bicubic_free(state);
      GSL_ERROR_NULL("failed to allocate space for zxy", GSL_ENOMEM);
    }

  state->xsize = xsize;
  state->ysize = ysize;

  return state;
} /* bicubic_alloc() */

static void
bicubic_free (void * vstate)
{
  bicubic_state_t *state = (bicubic_state_t *) vstate;

  RETURN_IF_NULL(state);

  if (state->zx)
    free (state->zx);

  if (state->zy)
    free (state->zy);

  if (state->zxy)
    free (state->zxy);

  free (state);
} /* bicubic_free() */

static int
bicubic_init(void * vstate, const double xa[], const double ya[],
             const double za[], size_t xsize, size_t ysize)
{
  bicubic_state_t *state = (bicubic_state_t *) vstate;

  gsl_interp_accel *acc = gsl_interp_accel_alloc();
  gsl_spline *spline;
  gsl_vector *x;
  gsl_vector *y;
  size_t i, j;

  x = gsl_vector_alloc(xsize);
  y = gsl_vector_alloc(xsize);
  spline = gsl_spline_alloc(gsl_interp_cspline, xsize);
  for (j = 0; j <= ysize - 1; j++)
    {
      for (i = 0; i <= xsize - 1; i++)
        {
          gsl_vector_set(x, i, xa[i]);
          gsl_vector_set(y, i, za[IDX2D(i, j, state)]);
        }
      gsl_spline_init(spline, x->data, y->data, xsize);
      for (i = 0; i <= xsize - 1; i++)
        {
          state->zx[IDX2D(i, j, state)] = gsl_spline_eval_deriv(spline, xa[i], acc);
        }
    }
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_spline_free(spline);
  gsl_interp_accel_reset(acc);

  x = gsl_vector_alloc(ysize);
  y = gsl_vector_alloc(ysize);
  spline = gsl_spline_alloc(gsl_interp_cspline, ysize);
  for (i = 0; i <= xsize - 1; i++)
    {
      for (j = 0; j <= ysize - 1; j++)
        {
          gsl_vector_set(x, j, ya[j]);
          gsl_vector_set(y, j, za[IDX2D(i, j, state)]);
        }
      gsl_spline_init(spline, x->data, y->data, ysize);
      for (j = 0; j <= ysize - 1; j++)
        {
          state->zy[IDX2D(i, j, state)] = gsl_spline_eval_deriv(spline, ya[j], acc);
        }
    }
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_spline_free(spline);
  gsl_interp_accel_reset(acc);

  x = gsl_vector_alloc(xsize);
  y = gsl_vector_alloc(xsize);
  spline = gsl_spline_alloc(gsl_interp_cspline, xsize);
  for (j = 0; j <= ysize - 1; j++)
    {
      for (i = 0; i <= xsize - 1; i++)
        {
          gsl_vector_set(x, i, xa[i]);
          gsl_vector_set(y, i, state->zy[IDX2D(i, j, state)]);
        }
      gsl_spline_init(spline, x->data, y->data, xsize);
      for (i = 0; i <= xsize - 1; i++)
        {
          state->zxy[IDX2D(i, j, state)] = gsl_spline_eval_deriv(spline, xa[i], acc);
        }
    }
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_spline_free(spline);
  gsl_interp_accel_free(acc);

  return GSL_SUCCESS;
} /* bicubic_init() */

static int
bicubic_eval(const void * vstate, const double xarr[], const double yarr[],
             const double zarr[], size_t xsize, size_t ysize,
             double x, double y, gsl_interp_accel * xa,
             gsl_interp_accel * ya, double * z)
{
  bicubic_state_t *state = (bicubic_state_t *) vstate;

  double xmin, xmax, ymin, ymax;
  double zminmin, zminmax, zmaxmin, zmaxmax;
  double zxminmin, zxminmax, zxmaxmin, zxmaxmax;
  double zyminmin, zyminmax, zymaxmin, zymaxmax;
  double zxyminmin, zxyminmax, zxymaxmin, zxymaxmax;

  double dx, dy;   /* size of the grid cell */
  double dt, du;

  /*
   * t and u are the positions within the grid cell at which we are computing
   * the interpolation, in units of grid cell size
   */
  double t, u;
  double t0, t1, t2, t3, u0, u1, u2, u3;
  double v;
  size_t xi, yi;

  /* first compute the indices into the data arrays where we are interpolating */
  if (xa != NULL)
    xi = gsl_interp_accel_find(xa, xarr, xsize, x);
  else
    xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);

  if (ya != NULL)
    yi = gsl_interp_accel_find(ya, yarr, ysize, y);
  else
    yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);

  /* find the minimum and maximum values on the grid cell in each dimension */
  xmin = xarr[xi];
  xmax = xarr[xi + 1];
  ymin = yarr[yi];
  ymax = yarr[yi + 1];
  zminmin = zarr[IDX2D(xi, yi, state)];
  zminmax = zarr[IDX2D(xi, yi + 1, state)];
  zmaxmin = zarr[IDX2D(xi + 1, yi, state)];
  zmaxmax = zarr[IDX2D(xi + 1, yi + 1, state)];
  /* Get the width and height of the grid cell */
  dx = xmax - xmin;
  dy = ymax - ymin;
  t = (x - xmin)/dx;
  u = (y - ymin)/dy;
  dt = 1./dx; /* partial t / partial x */
  du = 1./dy; /* partial u / partial y */
  zxminmin = state->zx[IDX2D(xi, yi, state)]/dt;
  zxminmax = state->zx[IDX2D(xi, yi + 1, state)]/dt;
  zxmaxmin = state->zx[IDX2D(xi + 1, yi, state)]/dt;
  zxmaxmax = state->zx[IDX2D(xi + 1, yi + 1, state)]/dt;
  zyminmin = state->zy[IDX2D(xi, yi, state)]/du;
  zyminmax = state->zy[IDX2D(xi, yi + 1, state)]/du;
  zymaxmin = state->zy[IDX2D(xi + 1, yi, state)]/du;
  zymaxmax = state->zy[IDX2D(xi + 1, yi + 1, state)]/du;
  zxyminmin = state->zxy[IDX2D(xi, yi, state)]/(dt*du);
  zxyminmax = state->zxy[IDX2D(xi, yi + 1, state)]/(dt*du);
  zxymaxmin = state->zxy[IDX2D(xi + 1, yi, state)]/(dt*du);
  zxymaxmax = state->zxy[IDX2D(xi + 1, yi + 1, state)]/(dt*du);
  t0 = 1;
  t1 = t;
  t2 = t*t;
  t3 = t*t2;
  u0 = 1;
  u1 = u;
  u2 = u*u;
  u3 = u*u2;

  *z = 0;
  v = zminmin;
  *z += v*t0*u0;
  v = zyminmin;
  *z += v*t0*u1;
  v = -3*zminmin + 3*zminmax - 2*zyminmin - zyminmax;
  *z += v*t0*u2;
  v = 2*zminmin - 2*zminmax + zyminmin + zyminmax;
  *z += v*t0*u3;
  v = zxminmin;
  *z += v*t1*u0;
  v = zxyminmin;
  *z += v*t1*u1;
  v = -3*zxminmin + 3*zxminmax - 2*zxyminmin - zxyminmax;
  *z += v*t1*u2;
  v = 2*zxminmin - 2*zxminmax + zxyminmin + zxyminmax;
  *z += v*t1*u3;
  v = -3*zminmin + 3*zmaxmin - 2*zxminmin - zxmaxmin;
  *z += v*t2*u0;
  v = -3*zyminmin + 3*zymaxmin - 2*zxyminmin - zxymaxmin;
  *z += v*t2*u1;
  v = 9*zminmin - 9*zmaxmin + 9*zmaxmax - 9*zminmax + 6*zxminmin + 3*zxmaxmin - 3*zxmaxmax - 6*zxminmax + 6*zyminmin - 6*zymaxmin - 3*zymaxmax + 3*zyminmax + 4*zxyminmin + 2*zxymaxmin + zxymaxmax + 2*zxyminmax;
  *z += v*t2*u2;
  v = -6*zminmin + 6*zmaxmin - 6*zmaxmax + 6*zminmax - 4*zxminmin - 2*zxmaxmin + 2*zxmaxmax + 4*zxminmax - 3*zyminmin + 3*zymaxmin + 3*zymaxmax - 3*zyminmax - 2*zxyminmin - zxymaxmin - zxymaxmax - 2*zxyminmax;
  *z += v*t2*u3;
  v = 2*zminmin - 2*zmaxmin + zxminmin + zxmaxmin;
  *z += v*t3*u0;
  v = 2*zyminmin - 2*zymaxmin + zxyminmin + zxymaxmin;
  *z += v*t3*u1;
  v = -6*zminmin + 6*zmaxmin - 6*zmaxmax + 6*zminmax - 3*zxminmin - 3*zxmaxmin + 3*zxmaxmax + 3*zxminmax - 4*zyminmin + 4*zymaxmin + 2*zymaxmax - 2*zyminmax - 2*zxyminmin - 2*zxymaxmin - zxymaxmax - zxyminmax;
  *z += v*t3*u2;
  v = 4*zminmin - 4*zmaxmin + 4*zmaxmax - 4*zminmax + 2*zxminmin + 2*zxmaxmin - 2*zxmaxmax - 2*zxminmax + 2*zyminmin - 2*zymaxmin - 2*zymaxmax + 2*zyminmax + zxyminmin + zxymaxmin + zxymaxmax + zxyminmax;
  *z += v*t3*u3;

  return GSL_SUCCESS;
} /* bicubic_eval() */

static int
bicubic_deriv_x(const void * vstate, const double xarr[], const double yarr[],
                const double zarr[], size_t xsize, size_t ysize,
                double x, double y,
                gsl_interp_accel * xa, gsl_interp_accel * ya, double * z_p)
{
  bicubic_state_t *state = (bicubic_state_t *) vstate;

  double xmin, xmax, ymin, ymax;
  double zminmin, zminmax, zmaxmin, zmaxmax;
  double zxminmin, zxminmax, zxmaxmin, zxmaxmax;
  double zyminmin, zyminmax, zymaxmin, zymaxmax;
  double zxyminmin, zxyminmax, zxymaxmin, zxymaxmax;
  double dx, dy; /* size of the grid cell */
  double dt, du;

  /*
   * t and u are the positions within the grid cell at which we are computing
   * the interpolation, in units of grid cell size
   */
  double t, u;
  double t0, t1, t2, u0, u1, u2, u3;
  double v;
  size_t xi, yi;

  /* first compute the indices into the data arrays where we are interpolating */
  if (xa != NULL)
    xi = gsl_interp_accel_find(xa, xarr, xsize, x);
  else
    xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);

  if (ya != NULL)
    yi = gsl_interp_accel_find(ya, yarr, ysize, y);
  else
    yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);

  /* find the minimum and maximum values on the grid cell in each dimension */
  xmin = xarr[xi];
  xmax = xarr[xi + 1];
  ymin = yarr[yi];
  ymax = yarr[yi + 1];
  zminmin = zarr[IDX2D(xi, yi, state)];
  zminmax = zarr[IDX2D(xi, yi + 1, state)];
  zmaxmin = zarr[IDX2D(xi + 1, yi, state)];
  zmaxmax = zarr[IDX2D(xi + 1, yi + 1, state)];

  /* get the width and height of the grid cell */
  dx = xmax - xmin;
  dy = ymax - ymin;
  t = (x - xmin)/dx;
  u = (y - ymin)/dy;
  dt = 1./dx; /* partial t / partial x */
  du = 1./dy; /* partial u / partial y */

  zxminmin = state->zx[IDX2D(xi, yi, state)]/dt;
  zxminmax = state->zx[IDX2D(xi, yi + 1, state)]/dt;
  zxmaxmin = state->zx[IDX2D(xi + 1, yi, state)]/dt;
  zxmaxmax = state->zx[IDX2D(xi + 1, yi + 1, state)]/dt;
  zyminmin = state->zy[IDX2D(xi, yi, state)]/du;
  zyminmax = state->zy[IDX2D(xi, yi + 1, state)]/du;
  zymaxmin = state->zy[IDX2D(xi + 1, yi, state)]/du;
  zymaxmax = state->zy[IDX2D(xi + 1, yi + 1, state)]/du;
  zxyminmin = state->zxy[IDX2D(xi, yi, state)]/(dt*du);
  zxyminmax = state->zxy[IDX2D(xi, yi + 1, state)]/(dt*du);
  zxymaxmin = state->zxy[IDX2D(xi + 1, yi, state)]/(dt*du);
  zxymaxmax = state->zxy[IDX2D(xi + 1, yi + 1, state)]/(dt*du);

  t0 = 1;
  t1 = t;
  t2 = t*t;
  u0 = 1;
  u1 = u;
  u2 = u*u;
  u3 = u*u2;

  *z_p = 0;
  v = zxminmin;
  *z_p += v*t0*u0;
  v = zxyminmin;
  *z_p += v*t0*u1;
  v = -3*zxminmin + 3*zxminmax - 2*zxyminmin - zxyminmax;
  *z_p += v*t0*u2;
  v = 2*zxminmin - 2*zxminmax + zxyminmin + zxyminmax;
  *z_p += v*t0*u3;
  v = -3*zminmin + 3*zmaxmin - 2*zxminmin - zxmaxmin;
  *z_p += 2*v*t1*u0;
  v = -3*zyminmin + 3*zymaxmin - 2*zxyminmin - zxymaxmin;
  *z_p += 2*v*t1*u1;
  v = 9*zminmin - 9*zmaxmin + 9*zmaxmax - 9*zminmax + 6*zxminmin + 3*zxmaxmin - 3*zxmaxmax - 6*zxminmax + 6*zyminmin - 6*zymaxmin - 3*zymaxmax + 3*zyminmax + 4*zxyminmin + 2*zxymaxmin + zxymaxmax + 2*zxyminmax;
  *z_p += 2*v*t1*u2;
  v = -6*zminmin + 6*zmaxmin - 6*zmaxmax + 6*zminmax - 4*zxminmin - 2*zxmaxmin + 2*zxmaxmax + 4*zxminmax - 3*zyminmin + 3*zymaxmin + 3*zymaxmax - 3*zyminmax - 2*zxyminmin - zxymaxmin - zxymaxmax - 2*zxyminmax;
  *z_p += 2*v*t1*u3;
  v = 2*zminmin - 2*zmaxmin + zxminmin + zxmaxmin;
  *z_p += 3*v*t2*u0;
  v = 2*zyminmin - 2*zymaxmin + zxyminmin + zxymaxmin;
  *z_p += 3*v*t2*u1;
  v = -6*zminmin + 6*zmaxmin - 6*zmaxmax + 6*zminmax - 3*zxminmin - 3*zxmaxmin + 3*zxmaxmax + 3*zxminmax - 4*zyminmin + 4*zymaxmin + 2*zymaxmax - 2*zyminmax - 2*zxyminmin - 2*zxymaxmin - zxymaxmax - zxyminmax;
  *z_p += 3*v*t2*u2;
  v = 4*zminmin - 4*zmaxmin + 4*zmaxmax - 4*zminmax + 2*zxminmin + 2*zxmaxmin - 2*zxmaxmax - 2*zxminmax + 2*zyminmin - 2*zymaxmin - 2*zymaxmax + 2*zyminmax + zxyminmin + zxymaxmin + zxymaxmax + zxyminmax;
  *z_p += 3*v*t2*u3;
  *z_p *= dt;

  return GSL_SUCCESS;
} /* bicubic_deriv_x() */

static int
bicubic_deriv_y(const void * vstate, const double xarr[], const double yarr[],
                const double zarr[], size_t xsize, size_t ysize,
                double x, double y,
                gsl_interp_accel * xa, gsl_interp_accel * ya, double * z_p)
{
  bicubic_state_t *state = (bicubic_state_t *) vstate;

  double xmin, xmax, ymin, ymax;
  double zminmin, zminmax, zmaxmin, zmaxmax;
  double zxminmin, zxminmax, zxmaxmin, zxmaxmax;
  double zyminmin, zyminmax, zymaxmin, zymaxmax;
  double zxyminmin, zxyminmax, zxymaxmin, zxymaxmax;
  /* dx and dy are the size of the grid cell */
  double dx, dy;
  double dt, du;
  /* t and u are the positions within the grid cell at which we are
   * computing the interpolation, in units of grid cell size */
  double t, u;
  double t0, t1, t2, t3, u0, u1, u2;
  double v;
  size_t xi, yi;

  /* first compute the indices into the data arrays where we are interpolating */
  if (xa != NULL)
    xi = gsl_interp_accel_find(xa, xarr, xsize, x);
  else
    xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);

  if (ya != NULL)
    yi = gsl_interp_accel_find(ya, yarr, ysize, y);
  else
    yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);

  /* find the minimum and maximum values on the grid cell in each dimension */
  xmin = xarr[xi];
  xmax = xarr[xi + 1];
  ymin = yarr[yi];
  ymax = yarr[yi + 1];
  zminmin = zarr[IDX2D(xi, yi, state)];
  zminmax = zarr[IDX2D(xi, yi + 1, state)];
  zmaxmin = zarr[IDX2D(xi + 1, yi, state)];
  zmaxmax = zarr[IDX2D(xi + 1, yi + 1, state)];

  /* get the width and height of the grid cell */
  dx = xmax - xmin;
  dy = ymax - ymin;
  t = (x - xmin)/dx;
  u = (y - ymin)/dy;
  dt = 1./dx; /* partial t / partial x */
  du = 1./dy; /* partial u / partial y */

  zxminmin = state->zx[IDX2D(xi, yi, state)]/dt;
  zxminmax = state->zx[IDX2D(xi, yi + 1, state)]/dt;
  zxmaxmin = state->zx[IDX2D(xi + 1, yi, state)]/dt;
  zxmaxmax = state->zx[IDX2D(xi + 1, yi + 1, state)]/dt;
  zyminmin = state->zy[IDX2D(xi, yi, state)]/du;
  zyminmax = state->zy[IDX2D(xi, yi + 1, state)]/du;
  zymaxmin = state->zy[IDX2D(xi + 1, yi, state)]/du;
  zymaxmax = state->zy[IDX2D(xi + 1, yi + 1, state)]/du;
  zxyminmin = state->zxy[IDX2D(xi, yi, state)]/(dt*du);
  zxyminmax = state->zxy[IDX2D(xi, yi + 1, state)]/(dt*du);
  zxymaxmin = state->zxy[IDX2D(xi + 1, yi, state)]/(dt*du);
  zxymaxmax = state->zxy[IDX2D(xi + 1, yi + 1, state)]/(dt*du);

  t0 = 1;
  t1 = t;
  t2 = t*t;
  t3 = t*t2;
  u0 = 1;
  u1 = u;
  u2 = u*u;

  *z_p = 0;
  v = zyminmin;
  *z_p += v*t0*u0;
  v = -3*zminmin + 3*zminmax - 2*zyminmin - zyminmax;
  *z_p += 2*v*t0*u1;
  v = 2*zminmin-2*zminmax + zyminmin + zyminmax;
  *z_p += 3*v*t0*u2;
  v = zxyminmin;
  *z_p += v*t1*u0;
  v = -3*zxminmin + 3*zxminmax - 2*zxyminmin - zxyminmax;
  *z_p += 2*v*t1*u1;
  v = 2*zxminmin - 2*zxminmax + zxyminmin + zxyminmax;
  *z_p += 3*v*t1*u2;
  v = -3*zyminmin + 3*zymaxmin - 2*zxyminmin - zxymaxmin;
  *z_p += v*t2*u0;
  v = 9*zminmin - 9*zmaxmin + 9*zmaxmax - 9*zminmax + 6*zxminmin + 3*zxmaxmin - 3*zxmaxmax - 6*zxminmax + 6*zyminmin - 6*zymaxmin - 3*zymaxmax + 3*zyminmax + 4*zxyminmin + 2*zxymaxmin + zxymaxmax + 2*zxyminmax;
  *z_p += 2*v*t2*u1;
  v = -6*zminmin + 6*zmaxmin - 6*zmaxmax + 6*zminmax - 4*zxminmin - 2*zxmaxmin + 2*zxmaxmax + 4*zxminmax - 3*zyminmin + 3*zymaxmin + 3*zymaxmax - 3*zyminmax - 2*zxyminmin - zxymaxmin - zxymaxmax - 2*zxyminmax;
  *z_p += 3*v*t2*u2;
  v = 2*zyminmin - 2*zymaxmin + zxyminmin + zxymaxmin;
  *z_p += v*t3*u0;
  v = -6*zminmin + 6*zmaxmin - 6*zmaxmax + 6*zminmax - 3*zxminmin - 3*zxmaxmin + 3*zxmaxmax + 3*zxminmax - 4*zyminmin + 4*zymaxmin + 2*zymaxmax - 2*zyminmax - 2*zxyminmin - 2*zxymaxmin - zxymaxmax - zxyminmax;
  *z_p += 2*v*t3*u1;
  v = 4*zminmin - 4*zmaxmin + 4*zmaxmax - 4*zminmax + 2*zxminmin + 2*zxmaxmin - 2*zxmaxmax - 2*zxminmax + 2*zyminmin - 2*zymaxmin - 2*zymaxmax + 2*zyminmax + zxyminmin + zxymaxmin + zxymaxmax + zxyminmax;
  *z_p += 3*v*t3*u2;
  *z_p *= du;

  return GSL_SUCCESS;
}

static int
bicubic_deriv_xx(const void * vstate, const double xarr[], const double yarr[],
                 const double zarr[], size_t xsize, size_t ysize,
                 double x, double y,
                 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z_pp)
{
  bicubic_state_t *state = (bicubic_state_t *) vstate;

  double xmin, xmax, ymin, ymax;
  double zminmin, zminmax, zmaxmin, zmaxmax;
  double zxminmin, zxminmax, zxmaxmin, zxmaxmax;
  double zyminmin, zyminmax, zymaxmin, zymaxmax;
  double zxyminmin, zxyminmax, zxymaxmin, zxymaxmax;

  double dx, dy; /* size of the grid cell */
  double dt, du;

  /*
   * t and u are the positions within the grid cell at which we are computing
   * the interpolation, in units of grid cell size
   */
  double t, u;
  double t0, t1, u0, u1, u2, u3;
  double v;
  size_t xi, yi;

  /* first compute the indices into the data arrays where we are interpolating */
  if (xa != NULL)
    xi = gsl_interp_accel_find(xa, xarr, xsize, x);
  else
    xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);

  if (ya != NULL)
    yi = gsl_interp_accel_find(ya, yarr, ysize, y);
  else
    yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);

  /* find the minimum and maximum values on the grid cell in each dimension */
  xmin = xarr[xi];
  xmax = xarr[xi + 1];
  ymin = yarr[yi];
  ymax = yarr[yi + 1];
  zminmin = zarr[IDX2D(xi, yi, state)];
  zminmax = zarr[IDX2D(xi, yi + 1, state)];
  zmaxmin = zarr[IDX2D(xi + 1, yi, state)];
  zmaxmax = zarr[IDX2D(xi + 1, yi + 1, state)];

  /* get the width and height of the grid cell */
  dx = xmax - xmin;
  dy = ymax - ymin;
  t = (x - xmin)/dx;
  u = (y - ymin)/dy;
  dt = 1./dx; /* partial t / partial x */
  du = 1./dy; /* partial u / partial y */

  zxminmin = state->zx[IDX2D(xi, yi, state)]/dt;
  zxminmax = state->zx[IDX2D(xi, yi + 1, state)]/dt;
  zxmaxmin = state->zx[IDX2D(xi + 1, yi, state)]/dt;
  zxmaxmax = state->zx[IDX2D(xi + 1, yi + 1, state)]/dt;
  zyminmin = state->zy[IDX2D(xi, yi, state)]/du;
  zyminmax = state->zy[IDX2D(xi, yi + 1, state)]/du;
  zymaxmin = state->zy[IDX2D(xi + 1, yi, state)]/du;
  zymaxmax = state->zy[IDX2D(xi + 1, yi + 1, state)]/du;
  zxyminmin = state->zxy[IDX2D(xi, yi, state)]/(dt*du);
  zxyminmax = state->zxy[IDX2D(xi, yi + 1, state)]/(dt*du);
  zxymaxmin = state->zxy[IDX2D(xi + 1, yi, state)]/(dt*du);
  zxymaxmax = state->zxy[IDX2D(xi + 1, yi + 1, state)]/(dt*du);

  t0 = 1;
  t1 = t;
  u0 = 1;
  u1 = u;
  u2 = u*u;
  u3 = u*u2;

  *z_pp = 0;
  v = -3*zminmin + 3*zmaxmin - 2*zxminmin - zxmaxmin;
  *z_pp += 2*v*t0*u0;
  v = -3*zyminmin + 3*zymaxmin - 2*zxyminmin - zxymaxmin;
  *z_pp += 2*v*t0*u1;
  v = 9*zminmin - 9*zmaxmin + 9*zmaxmax - 9*zminmax + 6*zxminmin + 3*zxmaxmin - 3*zxmaxmax - 6*zxminmax + 6*zyminmin - 6*zymaxmin - 3*zymaxmax + 3*zyminmax + 4*zxyminmin + 2*zxymaxmin + zxymaxmax + 2*zxyminmax;
  *z_pp += 2*v*t0*u2;
  v = -6*zminmin + 6*zmaxmin - 6*zmaxmax + 6*zminmax - 4*zxminmin - 2*zxmaxmin + 2*zxmaxmax + 4*zxminmax - 3*zyminmin + 3*zymaxmin + 3*zymaxmax - 3*zyminmax - 2*zxyminmin - zxymaxmin - zxymaxmax - 2*zxyminmax;
  *z_pp += 2*v*t0*u3;
  v = 2*zminmin - 2*zmaxmin + zxminmin + zxmaxmin;
  *z_pp += 6*v*t1*u0;
  v = 2*zyminmin - 2*zymaxmin + zxyminmin + zxymaxmin;
  *z_pp += 6*v*t1*u1;
  v = -6*zminmin + 6*zmaxmin - 6*zmaxmax + 6*zminmax - 3*zxminmin - 3*zxmaxmin + 3*zxmaxmax + 3*zxminmax - 4*zyminmin + 4*zymaxmin + 2*zymaxmax - 2*zyminmax - 2*zxyminmin - 2*zxymaxmin - zxymaxmax - zxyminmax;
  *z_pp += 6*v*t1*u2;
  v = 4*zminmin - 4*zmaxmin + 4*zmaxmax - 4*zminmax + 2*zxminmin + 2*zxmaxmin - 2*zxmaxmax - 2*zxminmax + 2*zyminmin - 2*zymaxmin - 2*zymaxmax + 2*zyminmax + zxyminmin + zxymaxmin + zxymaxmax + zxyminmax;
  *z_pp += 6*v*t1*u3;
  *z_pp *= dt*dt;

  return GSL_SUCCESS;
}

static int
bicubic_deriv_xy(const void * vstate, const double xarr[], const double yarr[],
                 const double zarr[], size_t xsize, size_t ysize,
                 double x, double y,
                 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z_pp)
{
  bicubic_state_t *state = (bicubic_state_t *) vstate;

  double xmin, xmax, ymin, ymax;
  double zminmin, zminmax, zmaxmin, zmaxmax;
  double zxminmin, zxminmax, zxmaxmin, zxmaxmax;
  double zyminmin, zyminmax, zymaxmin, zymaxmax;
  double zxyminmin, zxyminmax, zxymaxmin, zxymaxmax;

  double dx, dy; /* size of the grid cell */
  double dt, du;

  /*
   * t and u are the positions within the grid cell at which we are computing
   * the interpolation, in units of grid cell size
   */
  double t, u;
  double t0, t1, t2, u0, u1, u2;
  double v;
  size_t xi, yi;

  /* first compute the indices into the data arrays where we are interpolating */
  if (xa != NULL)
    xi = gsl_interp_accel_find(xa, xarr, xsize, x);
  else
    xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);

  if (ya != NULL)
    yi = gsl_interp_accel_find(ya, yarr, ysize, y);
  else
    yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);

  /* find the minimum and maximum values on the grid cell in each dimension */
  xmin = xarr[xi];
  xmax = xarr[xi + 1];
  ymin = yarr[yi];
  ymax = yarr[yi + 1];
  zminmin = zarr[IDX2D(xi, yi, state)];
  zminmax = zarr[IDX2D(xi, yi + 1, state)];
  zmaxmin = zarr[IDX2D(xi + 1, yi, state)];
  zmaxmax = zarr[IDX2D(xi + 1, yi + 1, state)];

  /* get the width and height of the grid cell */
  dx = xmax - xmin;
  dy = ymax - ymin;
  t = (x - xmin)/dx;
  u = (y - ymin)/dy;
  dt = 1./dx; /* partial t / partial x */
  du = 1./dy; /* partial u / partial y */

  zxminmin = state->zx[IDX2D(xi, yi, state)]/dt;
  zxminmax = state->zx[IDX2D(xi, yi + 1, state)]/dt;
  zxmaxmin = state->zx[IDX2D(xi + 1, yi, state)]/dt;
  zxmaxmax = state->zx[IDX2D(xi + 1, yi + 1, state)]/dt;
  zyminmin = state->zy[IDX2D(xi, yi, state)]/du;
  zyminmax = state->zy[IDX2D(xi, yi + 1, state)]/du;
  zymaxmin = state->zy[IDX2D(xi + 1, yi, state)]/du;
  zymaxmax = state->zy[IDX2D(xi + 1, yi + 1, state)]/du;
  zxyminmin = state->zxy[IDX2D(xi, yi, state)]/(dt*du);
  zxyminmax = state->zxy[IDX2D(xi, yi + 1, state)]/(dt*du);
  zxymaxmin = state->zxy[IDX2D(xi + 1, yi, state)]/(dt*du);
  zxymaxmax = state->zxy[IDX2D(xi + 1, yi + 1, state)]/(dt*du);

  t0 = 1;
  t1 = t;
  t2 = t*t;
  u0 = 1;
  u1 = u;
  u2 = u*u;

  *z_pp = 0;
  v = zxyminmin;
  *z_pp += v*t0*u0;
  v = -3*zxminmin + 3*zxminmax - 2*zxyminmin - zxyminmax;
  *z_pp += 2*v*t0*u1;
  v = 2*zxminmin - 2*zxminmax + zxyminmin + zxyminmax;
  *z_pp += 3*v*t0*u2;
  v = -3*zyminmin + 3*zymaxmin - 2*zxyminmin - zxymaxmin;
  *z_pp += 2*v*t1*u0;
  v = 9*zminmin - 9*zmaxmin + 9*zmaxmax - 9*zminmax + 6*zxminmin + 3*zxmaxmin - 3*zxmaxmax - 6*zxminmax + 6*zyminmin - 6*zymaxmin - 3*zymaxmax + 3*zyminmax + 4*zxyminmin + 2*zxymaxmin + zxymaxmax + 2*zxyminmax;
  *z_pp += 4*v*t1*u1;
  v = -6*zminmin + 6*zmaxmin - 6*zmaxmax + 6*zminmax - 4*zxminmin - 2*zxmaxmin + 2*zxmaxmax + 4*zxminmax - 3*zyminmin + 3*zymaxmin + 3*zymaxmax - 3*zyminmax - 2*zxyminmin - zxymaxmin - zxymaxmax - 2*zxyminmax;
  *z_pp += 6*v*t1*u2;
  v = 2*zyminmin - 2*zymaxmin + zxyminmin + zxymaxmin;
  *z_pp += 3*v*t2*u0;
  v = -6*zminmin + 6*zmaxmin - 6*zmaxmax + 6*zminmax - 3*zxminmin - 3*zxmaxmin + 3*zxmaxmax + 3*zxminmax - 4*zyminmin + 4*zymaxmin + 2*zymaxmax - 2*zyminmax - 2*zxyminmin - 2*zxymaxmin - zxymaxmax - zxyminmax;
  *z_pp += 6*v*t2*u1;
  v = 4*zminmin - 4*zmaxmin + 4*zmaxmax - 4*zminmax + 2*zxminmin + 2*zxmaxmin - 2*zxmaxmax - 2*zxminmax + 2*zyminmin - 2*zymaxmin - 2*zymaxmax + 2*zyminmax + zxyminmin + zxymaxmin + zxymaxmax + zxyminmax;
  *z_pp += 9*v*t2*u2;
  *z_pp *= dt*du;

  return GSL_SUCCESS;
}

static int
bicubic_deriv_yy(const void * vstate, const double xarr[], const double yarr[],
                 const double zarr[], size_t xsize, size_t ysize,
                 double x, double y,
                 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z_pp)
{
  bicubic_state_t *state = (bicubic_state_t *) vstate;

  double xmin, xmax, ymin, ymax;
  double zminmin, zminmax, zmaxmin, zmaxmax;
  double zxminmin, zxminmax, zxmaxmin, zxmaxmax;
  double zyminmin, zyminmax, zymaxmin, zymaxmax;
  double zxyminmin, zxyminmax, zxymaxmin, zxymaxmax;

  double dx, dy; /* size of the grid cell */
  double dt, du;

  /*
   * t and u are the positions within the grid cell at which we are computing
   * the interpolation, in units of grid cell size
   */
  double t, u;
  double t0, t1, t2, t3, u0, u1;
  double v;
  size_t xi, yi;

  /* first compute the indices into the data arrays where we are interpolating */
  if (xa != NULL)
    xi = gsl_interp_accel_find(xa, xarr, xsize, x);
  else
    xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);

  if (ya != NULL)
    yi = gsl_interp_accel_find(ya, yarr, ysize, y);
  else
    yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);

  /* find the minimum and maximum values on the grid cell in each dimension */
  xmin = xarr[xi];
  xmax = xarr[xi + 1];
  ymin = yarr[yi];
  ymax = yarr[yi + 1];
  zminmin = zarr[IDX2D(xi, yi, state)];
  zminmax = zarr[IDX2D(xi, yi + 1, state)];
  zmaxmin = zarr[IDX2D(xi + 1, yi, state)];
  zmaxmax = zarr[IDX2D(xi + 1, yi + 1, state)];

  /* get the width and height of the grid cell */
  dx = xmax - xmin;
  dy = ymax - ymin;
  t = (x - xmin)/dx;
  u = (y - ymin)/dy;
  dt = 1./dx; /* partial t / partial x */
  du = 1./dy; /* partial u / partial y */

  zxminmin = state->zx[IDX2D(xi, yi, state)]/dt;
  zxminmax = state->zx[IDX2D(xi, yi + 1, state)]/dt;
  zxmaxmin = state->zx[IDX2D(xi + 1, yi, state)]/dt;
  zxmaxmax = state->zx[IDX2D(xi + 1, yi + 1, state)]/dt;
  zyminmin = state->zy[IDX2D(xi, yi, state)]/du;
  zyminmax = state->zy[IDX2D(xi, yi + 1, state)]/du;
  zymaxmin = state->zy[IDX2D(xi + 1, yi, state)]/du;
  zymaxmax = state->zy[IDX2D(xi + 1, yi + 1, state)]/du;
  zxyminmin = state->zxy[IDX2D(xi, yi, state)]/(dt*du);
  zxyminmax = state->zxy[IDX2D(xi, yi + 1, state)]/(dt*du);
  zxymaxmin = state->zxy[IDX2D(xi + 1, yi, state)]/(dt*du);
  zxymaxmax = state->zxy[IDX2D(xi + 1, yi + 1, state)]/(dt*du);

  t0 = 1;
  t1 = t;
  t2 = t*t;
  t3 = t*t2;
  u0 = 1;
  u1 = u;

  *z_pp = 0;
  v = -3*zminmin + 3*zminmax - 2*zyminmin - zyminmax;
  *z_pp += 2*v*t0*u0;
  v = 2*zminmin-2*zminmax + zyminmin + zyminmax;
  *z_pp += 6*v*t0*u1;
  v = -3*zxminmin + 3*zxminmax - 2*zxyminmin - zxyminmax;
  *z_pp += 2*v*t1*u0;
  v = 2*zxminmin - 2*zxminmax + zxyminmin + zxyminmax;
  *z_pp += 6*v*t1*u1;
  v = 9*zminmin - 9*zmaxmin + 9*zmaxmax - 9*zminmax + 6*zxminmin + 3*zxmaxmin - 3*zxmaxmax - 6*zxminmax + 6*zyminmin - 6*zymaxmin - 3*zymaxmax + 3*zyminmax + 4*zxyminmin + 2*zxymaxmin + zxymaxmax + 2*zxyminmax;
  *z_pp += 2*v*t2*u0;
  v = -6*zminmin + 6*zmaxmin - 6*zmaxmax + 6*zminmax - 4*zxminmin - 2*zxmaxmin + 2*zxmaxmax + 4*zxminmax - 3*zyminmin + 3*zymaxmin + 3*zymaxmax - 3*zyminmax - 2*zxyminmin - zxymaxmin - zxymaxmax - 2*zxyminmax;
  *z_pp += 6*v*t2*u1;
  v = -6*zminmin + 6*zmaxmin - 6*zmaxmax + 6*zminmax - 3*zxminmin - 3*zxmaxmin + 3*zxmaxmax + 3*zxminmax - 4*zyminmin + 4*zymaxmin + 2*zymaxmax - 2*zyminmax - 2*zxyminmin - 2*zxymaxmin - zxymaxmax - zxyminmax;
  *z_pp += 2*v*t3*u0;
  v = 4*zminmin - 4*zmaxmin + 4*zmaxmax - 4*zminmax + 2*zxminmin + 2*zxmaxmin - 2*zxmaxmax - 2*zxminmax + 2*zyminmin - 2*zymaxmin - 2*zymaxmax + 2*zyminmax + zxyminmin + zxymaxmin + zxymaxmax + zxyminmax;
  *z_pp += 6*v*t3*u1;
  *z_pp *= du*du;

  return GSL_SUCCESS;
}

static const gsl_interp2d_type bicubic_type = {
  "bicubic",
  4,
  &bicubic_alloc,
  &bicubic_init,
  &bicubic_eval,
  &bicubic_deriv_x,
  &bicubic_deriv_y,
  &bicubic_deriv_xx,
  &bicubic_deriv_xy,
  &bicubic_deriv_yy,
  &bicubic_free
};

const gsl_interp2d_type * gsl_interp2d_bicubic = &bicubic_type;

#undef IDX2D

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <iostream>
#define ARMA_64BIT_WORD
#include <armadillo>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sf_psi.h>
#include <list>

#include "utils.h"
#include "eval.h"

using namespace std;
using namespace arma;
#include "data.h"

struct model_settings {
    bool verbose;

    string outdir;
    string datadir;

    double a_theta;
    double b_theta;
    double a_beta;
    double b_beta;
    double a_tau;
    double b_tau;
    double a_delta;
    double b_delta;

    bool social_only;
    bool factor_only;
    bool fix_influence;
    bool item_bias;
    bool binary;
    bool directed;

    long   seed;
    int    save_freq;
    int    eval_freq;
    int    conv_freq;
    int    max_iter;
    int    min_iter;
    double likelihood_delta;

    bool svi;
    bool final_pass;
    bool final_pass_test;
    int    sample_size;
    double delay;
    double forget;

    int k;

    void set(bool print, string out, string data, bool use_svi,
             double athe, double bthe, double abet, double bbet,
             double atau, double btau, double adelta, double bdelta,
             bool social, bool factor, bool bias, bool bin, bool dir,
             long rand, int savef, int evalf, int convf,
             int iter_max, int iter_min, double delta,
             bool finalpass, bool finalpasstest,
             int sample, double svi_delay, double svi_forget,
             bool fix, int num_factors) {
        verbose = print;

        outdir = out;
        datadir = data;

        svi = use_svi;

        a_theta = athe;
        b_theta = bthe;
        a_beta  = abet;
        b_beta  = bbet;
        a_tau   = atau;
        b_tau   = btau;
        a_delta = adelta;
        b_delta = bdelta;

        social_only = social;
        factor_only = factor;
        fix_influence = fix;
        item_bias = bias;
        binary = bin;
        directed = dir;

        seed = rand;
        save_freq = savef;
        eval_freq = evalf;
        conv_freq = convf;
        max_iter = iter_max;
        min_iter = iter_min;
        likelihood_delta = delta;

        final_pass = finalpass;
        final_pass_test = finalpasstest;
        sample_size = sample;
        delay = svi_delay;
        forget = svi_forget;

        k = num_factors;
    }

    void set_stochastic_inference(bool setting) {
        svi = setting;
    }

    void set_sample_size(int setting) {
        sample_size = setting;
    }

    void save(string filename) {
        FILE* file = fopen(filename.c_str(), "w");

        fprintf(file, "data directory: %s\n", datadir.c_str());

        fprintf(file, "\nmodel specification:\n");
        if (social_only) {
            fprintf(file, "\tsocial factorization (SF)   [ social factors only ]\n");
        } else if (factor_only) {
            fprintf(file, "\tPoisson factorization (PF)   [ general preference factors only ]\n");
        } else {
            fprintf(file, "\tsocial Poisson factorization (SPF)\n");
        }

        if (fix_influence) {
            fprintf(file,"\tsocial influence parameters fixed to 1\n");
        }

        if (!social_only) {
            fprintf(file, "\tK = %d   (number of latent factors for general preferences)\n", k);
        }

        fprintf(file, "\nshape and rate hyperparameters:\n");
        if (!social_only) {
            fprintf(file, "\ttheta (%f, %f)\n", a_theta, b_theta);
            fprintf(file, "\tbeta  (%f, %f)\n", a_beta, b_beta);
        }
        if (!factor_only) {
            fprintf(file, "\ttau   (%f, %f)\n", a_tau, b_tau);
        }
        if (item_bias) {
            fprintf(file, "\tdelta (%.2f, %.2f)\n", a_delta, b_delta);
        }

        fprintf(file, "\ndata attributes:\n");

        if (binary) {
            fprintf(file, "\tbinary ratings\n");
        } else {
            fprintf(file, "\tinteger ratings\n");
        }

        if (!factor_only) {
            if (directed) {
                fprintf(file, "\tdirected network\n");
            } else {
                fprintf(file, "\tundirected network\n");
            }
        }

        fprintf(file, "\ninference parameters:\n");
        fprintf(file, "\tseed:                                     %d\n", (int)seed);
        fprintf(file, "\tsave frequency:                           %d\n", save_freq);
        fprintf(file, "\tevaluation frequency:                     %d\n", eval_freq);
        fprintf(file, "\tconvergence check frequency:              %d\n", conv_freq);
        fprintf(file, "\tmaximum number of iterations:             %d\n", max_iter);
        fprintf(file, "\tminimum number of iterations:             %d\n", min_iter);
        fprintf(file, "\tchange in log likelihood for convergence: %f\n", likelihood_delta);
        fprintf(file, "\tfinal pass after convergence:             %s\n", final_pass ? "all users" :
            (final_pass_test ? "test users only" : "none"));

        if (svi) {
            fprintf(file, "\nStochastic variational inference parameters\n");
            fprintf(file, "\tsample size:                              %d\n", sample_size);
            fprintf(file, "\tSVI delay (tau):                          %f\n", delay);
            fprintf(file, "\tSVI forgetting rate (kappa):              %f\n", forget);
        } else {
            fprintf(file, "\nusing batch variational inference\n");
        }

        fclose(file);
    }
};

class SPF: protected Model {
    private:
        model_settings* settings;
        Data* data;

        // model parameters
        sp_fmat tau; // user influence
        sp_fmat logtau; // fake "log" user influence
                       // it's really exp(E[log(tau)]) which != E[tau]
        fmat theta;  // user preferences
        fmat beta;   // item attributes
        fmat logtheta;  // log variant of above
        fmat logbeta;   // ditto
        fvec delta;

        // helper parameters
        sp_fmat a_tau;
        sp_fmat b_tau;
        fmat a_theta;
        fmat b_theta;
        fmat a_beta;
        fmat a_beta_user;
        fmat a_beta_old;
        fmat b_beta;
        fvec a_delta;
        float b_delta;
        fvec a_delta_user;
        fvec a_delta_old;

        // random number generator
        gsl_rng* rand_gen;

        void initialize_parameters();
        void reset_helper_params();
        void save_parameters(string label);

        // parameter updates
        void update_shape(int user, int item, int rating);
        double update_tau(int user);
        double update_theta(int user);
        void update_beta(int item);
        void update_delta(int item);

        double get_ave_log_likelihood();
        void log_convergence(int iteration, double ave_ll, double delta_ll);
        void log_time(int iteration, double duration);
        void log_params(int iteration, double tau_change, double theta_change);
        void log_user(FILE* file, int user, int heldout, double rmse,
            double mae, double rank, int first, double crr, double ncrr,
            double ndcg);

        // define how to scale updates (training / sample size) (for SVI)
        double scale;

        // counts of number of times an item has been seen in a sample (for SVI)
        map<int,int> iter_count;

        void evaluate(string label);
        void evaluate(string label, bool write_rankings);

    public:
        SPF(model_settings* model_set, Data* dataset);
        void learn();
        double predict(int user, int item);
        void evaluate();

};

max_stars_repo_name: ajbc/spf

max_issues_repo_name: ajbc/spf

max_forks_repo_name: ajbc/spf


text:
#ifndef _RandomHistObject_H
#define _RandomHistObject_H

#include <gsl/gsl_histogram.h>
#include <QObject>
#include <QVector>
#include "src_config.h"

class RandomNumberGen;

class _HIST_EXPORT_ RandomHistObject : public QObject
{
public:
    static RandomHistObject * getRandomHistObject (QObject * parent=0);
    void GUIRandomProc (QWidget * parent=0, Qt::WindowFlags flags=0);
    void GUIHistogramParams (QWidget * parent=0, Qt::WindowFlags flags=0);
    void calcHist (void);
    void GUIViewHist (QWidget * parent=0, Qt::WindowFlags flags=0);

    enum DistrFuncs
    {
        Uniform=0,
        Gaussian=1,
        Exponential=2
    };

signals:
    void viewResults (QWidget * w);

private:
    //
    // Functions
    //
    RandomHistObject (QObject * parent=0);
    virtual ~RandomHistObject (void);

private:
    //
    // Variables
    //
    static RandomHistObject * self;
    RandomNumberGen * rGenerator;
    QVector<double> resNumb;
    gsl_histogram * hist;
private:
    Q_OBJECT
};

#endif

max_stars_repo_name: YuriyRusinov/Histogram

max_issues_repo_name: YuriyRusinov/Histogram

max_forks_repo_name: YuriyRusinov/Histogram


text:
#pragma once

#include <gsl/span>

#include <cstddef>

namespace imageview {

// Class representing a color in RGBA32 color space.
class RGBA32 {
 public:
  constexpr RGBA32() noexcept = default;

  constexpr RGBA32(unsigned char red, unsigned char green, unsigned char blue, unsigned char alpha) noexcept
      : red(red), green(green), blue(blue), alpha(alpha) {}

  unsigned char red = 0;
  unsigned char green = 0;
  unsigned char blue = 0;
  unsigned char alpha = 0;
};

// Implementation of the PixelFormat concept for RGBA32 pixel format.
// In this pixel format the color is represented via 4 8-bit integers,
// specifying the red, green, blue and alpha channels. When serializing to /
// deserializing from  a byte array, the order of the channels is RGBA (i.e.
// not BGRA).
class PixelFormatRGBA32 {
 public:
  using color_type = RGBA32;
  static constexpr int kBytesPerPixel = 4;

  constexpr color_type read(gsl::span<const std::byte, kBytesPerPixel> data) const;

  constexpr void write(const color_type& color, gsl::span<std::byte, kBytesPerPixel> data) const;
};

constexpr bool operator==(const RGBA32& lhs, const RGBA32& rhs) {
  return lhs.red == rhs.red && lhs.green == rhs.green && lhs.blue == rhs.blue && lhs.alpha == rhs.alpha;
}

constexpr bool operator!=(const RGBA32& lhs, const RGBA32& rhs) { return !(lhs == rhs); }

constexpr PixelFormatRGBA32::color_type PixelFormatRGBA32::read(gsl::span<const std::byte, kBytesPerPixel> data) const {
  return color_type(static_cast<unsigned char>(data[0]), static_cast<unsigned char>(data[1]),
                    static_cast<unsigned char>(data[2]), static_cast<unsigned char>(data[3]));
}

constexpr void PixelFormatRGBA32::write(const color_type& color, gsl::span<std::byte, kBytesPerPixel> data) const {
  data[0] = static_cast<std::byte>(color.red);
  data[1] = static_cast<std::byte>(color.green);
  data[2] = static_cast<std::byte>(color.blue);
  data[3] = static_cast<std::byte>(color.alpha);
}

}  // namespace imageview

max_stars_repo_name: alexanderbelous/imageview

max_issues_repo_name: alexanderbelous/imageview

max_forks_repo_name: alexanderbelous/imageview


text:
#include "matrix.h"
#include <cblas.h>
#include <malloc.h>
#include <stdio.h>
#include <stdlib.h>

Matrix *create_mat(char *path) {
  Matrix *m = malloc(sizeof(Matrix));
  m->row_len = 0;
  m->col_len = 0;
  FILE *myfile = fopen(path, "r");
  int ch = 0;
  while (ch != EOF) {
    ch = fgetc(myfile);
    if (m->row_len == 0 && ch == ' ') {
      m->col_len++;
    }
    if (ch == '\n') {
      m->row_len++;
    }
  };

  m->col_len++;

  m->p_mat = calloc(m->row_len, sizeof(float *));
  for (int i = 0; i < m->row_len; i++) {
    m->p_mat[i] = calloc(m->col_len, sizeof(float));
  }

  rewind(myfile);

  for (int i = 0; i < m->row_len; i++) {
    for (int j = 0; j < m->col_len; j++) {
      if (!fscanf(myfile, "%f", &m->p_mat[i][j])) {
        break;
      }
    }
  }

  fclose(myfile);
  return m;
}

void delete_mat(Matrix *mat) {
  for (int i = 0; i < mat->row_len; i++) {
    free(mat->p_mat[i]);
  }
  free(mat->p_mat);
  free(mat);
}

void copy_mat(Matrix *target, Matrix *source) {
  delete_mat(target);
  target->row_len = source->row_len;
  target->col_len = source->col_len;
  target->p_mat = calloc(target->row_len, sizeof(float *));
  for (int i = 0; i < target->row_len; i++) {
    target->p_mat[i] = calloc(target->col_len, sizeof(float));
  }
  for (int i = 0; i < target->row_len; i++) {
    for (int j = 0; j < target->col_len; j++) {
      target->p_mat[i][j] = source->p_mat[i][j];
    }
  }
}

Matrix *multiple_mat(Matrix *mat1, Matrix *mat2) {
  Matrix *result = malloc(sizeof(Matrix));
  result->row_len = mat1->row_len;
  result->col_len = mat2->col_len;
  result->p_mat = calloc(result->row_len, sizeof(float *));
  for (int i = 0; i < result->row_len; ++i)
    result->p_mat[i] = calloc(result->col_len, sizeof(float));

#pragma omp parallel for
  for (int i = 0; i < result->row_len; i++) {
    for (int k = 0; k < mat1->col_len; k++) {
      float tmp = mat1->p_mat[i][k];
      for (int j = 0; j < result->col_len; j++) {
        result->p_mat[i][j] += tmp * mat2->p_mat[k][j];
      }
    }
  }

  return result;
}
void write_mat(char *path, Matrix *mat) {
  FILE *myfile = fopen(path, "w");
  for (int i = 0; i < mat->row_len; i++) {
    for (int j = 0; j < mat->col_len; j++) {
      if (j == mat->col_len - 1) {
        fprintf(myfile, "%f\n", mat->p_mat[i][j]);
      } else {
        fprintf(myfile, "%f ", mat->p_mat[i][j]);
      }
    }
  }
}

void display_mat(Matrix *m1) {
  for (int i = 0; i < m1->row_len; i++) {
    for (int j = 0; j < m1->col_len; j++) {
      if (j == m1->col_len - 1) {
        printf("%f\n", m1->p_mat[i][j]);
      } else {
        printf("%f ", m1->p_mat[i][j]);
      }
    }
  }
}

Matrix *calc_with_OpenBLAS(Matrix *m1, Matrix *m2) {
  Matrix *result = malloc(sizeof(Matrix));
  result->row_len = m1->row_len;
  result->col_len = m2->col_len;
  result->p_mat = calloc(result->row_len, sizeof(float *));
  for (int i = 0; i < result->row_len; ++i)
    result->p_mat[i] = calloc(result->col_len, sizeof(float));
  const int M = m1->row_len;
  const int N = m2->col_len;
  const int K = m1->col_len;
  int lda = K;
  int ldb = N;
  int ldc = N;
  float *A = calloc(M * K, sizeof(float));
  float *B = calloc(K * N, sizeof(float));
  float *C = calloc(M * N, sizeof(float));
  int k = 0;
  for (int i = 0; i < M; i++) {
    for (int j = 0; j < K; j++) {
      A[k] = m1->p_mat[i][j];
      k += 1;
    }
  }
  k = 0;
  for (int i = 0; i < K; i++) {
    for (int j = 0; j < N; j++) {
      B[k] = m2->p_mat[i][j];
      k += 1;
    }
  }
  cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, M, N, K, 1.0f, A, lda,
              B, ldb, 1.0f, C, ldc);
  k = 0;
  for (int i = 0; i < result->row_len; i++) {
    for (int j = 0; j < result->col_len; j++) {
      result->p_mat[i][j] = C[k];
      k++;
    }
  }
  free(A);
  free(B);
  free(C);
  return result;
}

max_stars_repo_name: Aries-Dawn/Cpp-Program-Design

max_issues_repo_name: Aries-Dawn/Cpp-Program-Design

max_forks_repo_name: Aries-Dawn/Cpp-Program-Design


text:
/* permutation/permutation.c
 *
 * Copyright (C) 2001, 2002 Nicolas Darnis
 *
 * 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.
 */

/* Modified for GSL by Brian Gough.
 * Use in-place algorithms, no need for workspace
 * Use conventions for canonical form given in Knuth (opposite of Sedgewick)
 */

#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_permutation.h>

int
gsl_permutation_linear_to_canonical (gsl_permutation * q,
                                     const gsl_permutation * p)
{
  const size_t n = p->size;
  size_t i, k, s;
  size_t t = n;

  const size_t *const pp = p->data;
  size_t *const qq = q->data;

  if (q->size != p->size)
    {
      GSL_ERROR ("size of q does not match size of p", GSL_EINVAL);
    }

  for (i = 0; i < n; i++)
    {

      k = pp[i];
      s = 1;

      while (k > i)
        {
          k = pp[k];
          s++;
        }

      if (k < i)
        continue;

      /* Now have k == i, i.e the least in its cycle, and s == cycle length */

      t -= s;

      qq[t] = i;

      k = pp[i];
      s = 1;

      while (k > i)
        {
          qq[t + s] = k;
          k = pp[k];
          s++;
        }

      if (t == 0)
        break;
    }

  return GSL_SUCCESS;
}

int
gsl_permutation_canonical_to_linear (gsl_permutation * p,
                                     const gsl_permutation * q)
{
  size_t i, k, kk, first;
  const size_t n = p->size;

  size_t *const pp = p->data;
  const size_t *const qq = q->data;

  if (q->size != p->size)
    {
      GSL_ERROR ("size of q does not match size of p", GSL_EINVAL);
    }

  for (i = 0; i < n; i++)
    {
      pp[i] = i;
    }

  k = qq[0];
  first = pp[k];

  for (i = 1; i < n; i++)
    {
      kk = qq[i];

      if (kk > first)
        {
          pp[k] = pp[kk];
          k = kk;
        }
      else
        {
          pp[k] = first;
          k = kk;
          first = pp[kk];
        }
    }

  pp[k] = first;

  return GSL_SUCCESS;
}

size_t
gsl_permutation_inversions (const gsl_permutation * p)
{
  size_t count = 0;
  size_t i, j;
  const size_t size = p->size;

  for (i = 0; i < size - 1; i++)
    {
      for (j = i + 1; j < size; j++)
        {
          if (p->data[i] > p->data[j])
            {
              count++;
            }
        }
    }

  return count;
}

size_t
gsl_permutation_linear_cycles (const gsl_permutation * p)
{
  size_t i, k;
  size_t count = 0;
  const size_t size = p->size;

  for (i = 0; i < size; i++)
    {

      k = p->data[i];

      while (k > i)
        {
          k = p->data[k];
        }

      if (k < i)
        continue;

      count++;
    }

  return count;
}

size_t
gsl_permutation_canonical_cycles (const gsl_permutation * p)
{
  size_t i;
  size_t count = 1;
  size_t min = p->data[0];

  for (i = 0; i < p->size; i++)
    {
      if (p->data[i] < min)
        {
          min = p->data[i];
          count++;
        }
    }

  return count;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/**
 *
 * @file qwrapper_slaswp.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:44:58 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_slaswp(Quark *quark, Quark_Task_Flags *task_flags,
                       int n, float *A, int lda,
                       int i1,  int i2, const int *ipiv, int inc)
{
    DAG_CORE_LASWP;
    QUARK_Insert_Task(
        quark, CORE_slaswp_quark, task_flags,
        sizeof(int),                      &n,    VALUE,
        sizeof(float)*lda*n,  A,        INOUT | LOCALITY,
        sizeof(int),                      &lda,  VALUE,
        sizeof(int),                      &i1,   VALUE,
        sizeof(int),                      &i2,   VALUE,
        sizeof(int)*n,                     ipiv,     INPUT,
        sizeof(int),                      &inc,  VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slaswp_quark = PCORE_slaswp_quark
#define CORE_slaswp_quark PCORE_slaswp_quark
#endif
void CORE_slaswp_quark(Quark *quark)
{
    int n, lda, i1, i2, inc;
    int *ipiv;
    float *A;

    quark_unpack_args_7(quark, n, A, lda, i1, i2, ipiv, inc);
    LAPACKE_slaswp_work(LAPACK_COL_MAJOR, n, A, lda, i1, i2, ipiv, inc );
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_slaswp_f2(Quark *quark, Quark_Task_Flags *task_flags,
                          int n, float *A, int lda,
                          int i1,  int i2, const int *ipiv, int inc,
                          float *fake1, int szefake1, int flag1,
                          float *fake2, int szefake2, int flag2)
{
    DAG_CORE_LASWP;
    QUARK_Insert_Task(
        quark, CORE_slaswp_f2_quark, task_flags,
        sizeof(int),                        &n,     VALUE,
        sizeof(float)*lda*n,    A,         INOUT | LOCALITY,
        sizeof(int),                        &lda,   VALUE,
        sizeof(int),                        &i1,    VALUE,
        sizeof(int),                        &i2,    VALUE,
        sizeof(int)*n,                       ipiv,      INPUT,
        sizeof(int),                        &inc,   VALUE,
        sizeof(float)*szefake1, fake1,     flag1,
        sizeof(float)*szefake2, fake2,     flag2,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slaswp_f2_quark = PCORE_slaswp_f2_quark
#define CORE_slaswp_f2_quark PCORE_slaswp_f2_quark
#endif
void CORE_slaswp_f2_quark(Quark* quark)
{
    int n, lda, i1, i2, inc;
    int *ipiv;
    float *A;
    void *fake1, *fake2;

    quark_unpack_args_9(quark, n, A, lda, i1, i2, ipiv, inc, fake1, fake2);
    LAPACKE_slaswp_work(LAPACK_COL_MAJOR, n, A, lda, i1, i2, ipiv, inc );
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_slaswp_ontile(Quark *quark, Quark_Task_Flags *task_flags,
                              PLASMA_desc descA, float *Aij,
                              int i1,  int i2, const int *ipiv, int inc, float *fakepanel)
{
    DAG_CORE_LASWP;
    if (fakepanel == Aij) {
        QUARK_Insert_Task(
            quark, CORE_slaswp_ontile_quark, task_flags,
            sizeof(PLASMA_desc),              &descA,     VALUE,
            sizeof(float)*1,      Aij,           INOUT | LOCALITY,
            sizeof(int),                      &i1,        VALUE,
            sizeof(int),                      &i2,        VALUE,
            sizeof(int)*(i2-i1+1)*abs(inc),   ipiv,           INPUT,
            sizeof(int),                      &inc,       VALUE,
            sizeof(float)*1,      fakepanel,     SCRATCH,
            0);
    } else {
        QUARK_Insert_Task(
            quark, CORE_slaswp_ontile_quark, task_flags,
            sizeof(PLASMA_desc),              &descA,     VALUE,
            sizeof(float)*1,      Aij,           INOUT | LOCALITY,
            sizeof(int),                      &i1,        VALUE,
            sizeof(int),                      &i2,        VALUE,
            sizeof(int)*(i2-i1+1)*abs(inc),   ipiv,           INPUT,
            sizeof(int),                      &inc,       VALUE,
            sizeof(float)*1,      fakepanel,     INOUT,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slaswp_ontile_quark = PCORE_slaswp_ontile_quark
#define CORE_slaswp_ontile_quark PCORE_slaswp_ontile_quark
#endif
void CORE_slaswp_ontile_quark(Quark *quark)
{
    int i1, i2, inc;
    int *ipiv;
    float *A, *fake;
    PLASMA_desc descA;

    quark_unpack_args_7(quark, descA, A, i1, i2, ipiv, inc, fake);
    CORE_slaswp_ontile(descA, i1, i2, ipiv, inc);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_slaswp_ontile_f2(Quark *quark, Quark_Task_Flags *task_flags,
                                 PLASMA_desc descA, float *Aij,
                                 int i1,  int i2, const int *ipiv, int inc,
                                 float *fake1, int szefake1, int flag1,
                                 float *fake2, int szefake2, int flag2)
{
    DAG_CORE_LASWP;
    QUARK_Insert_Task(
        quark, CORE_slaswp_ontile_f2_quark, task_flags,
        sizeof(PLASMA_desc),                &descA, VALUE,
        sizeof(float)*1,        Aij,       INOUT | LOCALITY,
        sizeof(int),                        &i1,    VALUE,
        sizeof(int),                        &i2,    VALUE,
        sizeof(int)*(i2-i1+1)*abs(inc),      ipiv,      INPUT,
        sizeof(int),                        &inc,   VALUE,
        sizeof(float)*szefake1, fake1, flag1,
        sizeof(float)*szefake2, fake2, flag2,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slaswp_ontile_f2_quark = PCORE_slaswp_ontile_f2_quark
#define CORE_slaswp_ontile_f2_quark PCORE_slaswp_ontile_f2_quark
#endif
void CORE_slaswp_ontile_f2_quark(Quark *quark)
{
    int i1, i2, inc;
    int *ipiv;
    float *A;
    PLASMA_desc descA;
    void *fake1, *fake2;

    quark_unpack_args_8(quark, descA, A, i1, i2, ipiv, inc, fake1, fake2);
    CORE_slaswp_ontile(descA, i1, i2, ipiv, inc);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_sswptr_ontile(Quark *quark, Quark_Task_Flags *task_flags,
                              PLASMA_desc descA, float *Aij,
                              int i1,  int i2, const int *ipiv, int inc,
                              const float *Akk, int ldak)
{
    DAG_CORE_TRSM;
    QUARK_Insert_Task(
        quark, CORE_sswptr_ontile_quark, task_flags,
        sizeof(PLASMA_desc),              &descA, VALUE,
        sizeof(float)*1,      Aij,       INOUT | LOCALITY,
        sizeof(int),                      &i1,    VALUE,
        sizeof(int),                      &i2,    VALUE,
        sizeof(int)*(i2-i1+1)*abs(inc),    ipiv,      INPUT,
        sizeof(int),                      &inc,   VALUE,
        sizeof(float)*ldak,   Akk,       INPUT,
        sizeof(int),                      &ldak,  VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_sswptr_ontile_quark = PCORE_sswptr_ontile_quark
#define CORE_sswptr_ontile_quark PCORE_sswptr_ontile_quark
#endif
void CORE_sswptr_ontile_quark(Quark *quark)
{
    int i1, i2, inc, ldak;
    int *ipiv;
    float *A, *Akk;
    PLASMA_desc descA;

    quark_unpack_args_8(quark, descA, A, i1, i2, ipiv, inc, Akk, ldak);
    CORE_sswptr_ontile(descA, i1, i2, ipiv, inc, Akk, ldak);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_slaswpc_ontile(Quark *quark, Quark_Task_Flags *task_flags,
                              PLASMA_desc descA, float *Aij,
                              int i1,  int i2, const int *ipiv, int inc, float *fakepanel)
{
    DAG_CORE_LASWP;
    if (fakepanel == Aij) {
        QUARK_Insert_Task(
            quark, CORE_slaswpc_ontile_quark, task_flags,
            sizeof(PLASMA_desc),              &descA,     VALUE,
            sizeof(float)*1,      Aij,           INOUT | LOCALITY,
            sizeof(int),                      &i1,        VALUE,
            sizeof(int),                      &i2,        VALUE,
            sizeof(int)*(i2-i1+1)*abs(inc),   ipiv,           INPUT,
            sizeof(int),                      &inc,       VALUE,
            sizeof(float)*1,      fakepanel,     SCRATCH,
            0);
    } else {
        QUARK_Insert_Task(
            quark, CORE_slaswpc_ontile_quark, task_flags,
            sizeof(PLASMA_desc),              &descA,     VALUE,
            sizeof(float)*1,      Aij,           INOUT | LOCALITY,
            sizeof(int),                      &i1,        VALUE,
            sizeof(int),                      &i2,        VALUE,
            sizeof(int)*(i2-i1+1)*abs(inc),   ipiv,           INPUT,
            sizeof(int),                      &inc,       VALUE,
            sizeof(float)*1,      fakepanel,     INOUT,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slaswpc_ontile_quark = PCORE_slaswpc_ontile_quark
#define CORE_slaswpc_ontile_quark PCORE_slaswpc_ontile_quark
#endif
void CORE_slaswpc_ontile_quark(Quark *quark)
{
    int i1, i2, inc;
    int *ipiv;
    float *A, *fake;
    PLASMA_desc descA;

    quark_unpack_args_7(quark, descA, A, i1, i2, ipiv, inc, fake);
    CORE_slaswpc_ontile(descA, i1, i2, ipiv, inc);
}


text:
//
//  rope.hpp
//  Fans
//
//  Created by tqtifnypmb on 06/12/2017.
//  Copyright © 2017 tqtifnypmb. All rights reserved.
//

#pragma once

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

#include <gsl/gsl>
#include "RopeNode.h"
#include "Range.h"
#include "RopeIter.h"
#include "../types.h"

namespace brick
{

class Rope {
public:
    static constexpr size_t max_leaf_length = 1;

	Rope(Rope&& l, Rope&& r);

	Rope();
    Rope(const Rope&);

	Rope(Rope&& r) = default;
	Rope& operator=(Rope&&) = default;
	Rope& operator=(const Rope&) = delete;

	template <class Converter>
    void insert(gsl::span<const char> bytes, size_t pos);
    void insert(const detail::CodePointList& cp, size_t pos);

    void erase(const Range& range);

    size_t size() const {
        return size_;
    }

    RopeIter begin() const;
    RopeIter end() const;

    std::reverse_iterator<RopeIter> rbegin() const {
        return std::reverse_iterator<RopeIter>(end());
    }

    std::reverse_iterator<RopeIter> rend() const {
        return std::reverse_iterator<RopeIter>(begin());
    }

    RopeIter iterator(size_t index) const;

    bool empty() const {
        return size_ == 0;
    }

    void swap(Rope& other) noexcept {
        std::swap(size_, other.size_);
        root_.swap(other.root_);
    }

    std::string string() const;

    detail::RopeNode* root_test() {
        return root_.get();
    }

    detail::RopeNodePtr nextLeaf_test(detail::RopeNode* current) {
        return nextLeaf(current);
    }

    detail::RopeNodePtr prevLeaf_test(detail::RopeNode* current) {
        return prevLeaf(current);
    }

    void rebalance_test() {
        rebalance();
    }

    std::tuple<detail::RopeNodePtr, size_t> get_test(detail::RopeNode* root, size_t index) {
        return getLeaf(root, index);
    }

    size_t lengthOfWholeRope_test(gsl::not_null<detail::RopeNode*> root) {
        return lengthOfWholeRope(root);
    }

    bool checkHeight();
    bool checkLength();

private:
    friend class RopeIter;

    static const size_t npos = -1;

    Rope(std::vector<std::unique_ptr<detail::RopeNode>>& cplist);

    detail::RopeNodePtr nextLeaf(gsl::not_null<detail::RopeNode*> current) const;
    detail::RopeNodePtr prevLeaf(gsl::not_null<detail::RopeNode*> current) const;
    void removeLeaf(detail::RopeNodePtr node);

    size_t lengthOfWholeRope(gsl::not_null<detail::RopeNode*> root);

    std::vector<std::unique_ptr<detail::RopeNode>> collectLeaves(bool move) const;
	bool needBalance();
	void rebalance();

    std::tuple<detail::RopeNodePtr /* leaf */, size_t /* pos */>
    getLeaf(gsl::not_null<detail::RopeNode*> root, size_t index, detail::RopeNode** lastVisitedNode = nullptr) const;

    void insertSubRope(detail::RopeNodePtr leaf, size_t pos, std::unique_ptr<detail::RopeNode> subRope, size_t len);

    size_t size_ = 0;
	std::unique_ptr<detail::RopeNode> root_;
};

template <class Converter>
void Rope::insert(gsl::span<const char> bytes, size_t pos) {
    insert(Converter::encode(bytes), pos);
}

}   // namespace brick

max_stars_repo_name: tqtifnypmb/brick

max_issues_repo_name: tqtifnypmb/brick

max_forks_repo_name: tqtifnypmb/brick


text:
//
// Autogenerated based on code by L. J. Feldstein
//

#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include "cimple_controller.h"

max_stars_repo_name: shaesaert/TuLiPXML

max_issues_repo_name: shaesaert/TuLiPXML

max_forks_repo_name: shaesaert/TuLiPXML