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

#include <config.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_errno.h>

/*
 * This module contains routines related to the Cholesky decomposition
 * of a complex Hermitian positive definite matrix.
 */

static void cholesky_complex_conj_vector(gsl_vector_complex *v);

/*
gsl_linalg_complex_cholesky_decomp()
  Perform the Cholesky decomposition on a Hermitian positive definite
matrix. See Golub & Van Loan, "Matrix Computations" (3rd ed),
algorithm 4.2.2.

Inputs: A - (input/output) complex postive definite matrix

Return: success or error

The lower triangle of A is overwritten with the Cholesky decomposition
*/

int
gsl_linalg_complex_cholesky_decomp(gsl_matrix_complex *A)
{
  const size_t N = A->size1;

  if (N != A->size2)
    {
      GSL_ERROR("cholesky decomposition requires square matrix", GSL_ENOTSQR);
    }
  else
    {
      size_t i, j;
      gsl_complex z;
      double ajj;

      for (j = 0; j < N; ++j)
        {
          z = gsl_matrix_complex_get(A, j, j);
          ajj = GSL_REAL(z);

          if (j > 0)
            {
              gsl_vector_complex_const_view aj =
                gsl_matrix_complex_const_subrow(A, j, 0, j);

              gsl_blas_zdotc(&aj.vector, &aj.vector, &z);
              ajj -= GSL_REAL(z);
            }

          if (ajj <= 0.0)
            {
              GSL_ERROR("matrix is not positive definite", GSL_EDOM);
            }

          ajj = sqrt(ajj);
          GSL_SET_COMPLEX(&z, ajj, 0.0);
          gsl_matrix_complex_set(A, j, j, z);

          if (j < N - 1)
            {
              gsl_vector_complex_view av =
                gsl_matrix_complex_subcolumn(A, j, j + 1, N - j - 1);

              if (j > 0)
                {
                  gsl_vector_complex_view aj =
                    gsl_matrix_complex_subrow(A, j, 0, j);
                  gsl_matrix_complex_view am =
                    gsl_matrix_complex_submatrix(A, j + 1, 0, N - j - 1, j);

                  cholesky_complex_conj_vector(&aj.vector);

                  gsl_blas_zgemv(CblasNoTrans,
                                 GSL_COMPLEX_NEGONE,
                                 &am.matrix,
                                 &aj.vector,
                                 GSL_COMPLEX_ONE,
                                 &av.vector);

                  cholesky_complex_conj_vector(&aj.vector);
                }

              gsl_blas_zdscal(1.0 / ajj, &av.vector);
            }
        }

      /* Now store L^H in upper triangle */
      for (i = 1; i < N; ++i)
        {
          for (j = 0; j < i; ++j)
            {
              z = gsl_matrix_complex_get(A, i, j);
              gsl_matrix_complex_set(A, j, i, gsl_complex_conjugate(z));
            }
        }

      return GSL_SUCCESS;
    }
} /* gsl_linalg_complex_cholesky_decomp() */

/*
gsl_linalg_complex_cholesky_solve()
  Solve A x = b where A is in cholesky form
*/

int
gsl_linalg_complex_cholesky_solve (const gsl_matrix_complex * cholesky,
                                   const gsl_vector_complex * b,
                                   gsl_vector_complex * x)
{
  if (cholesky->size1 != cholesky->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (cholesky->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (cholesky->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      gsl_vector_complex_memcpy (x, b);

      /* solve for y using forward-substitution, L y = b */

      gsl_blas_ztrsv (CblasLower, CblasNoTrans, CblasNonUnit, cholesky, x);

      /* perform back-substitution, L^H x = y */

      gsl_blas_ztrsv (CblasLower, CblasConjTrans, CblasNonUnit, cholesky, x);

      return GSL_SUCCESS;
    }
} /* gsl_linalg_complex_cholesky_solve() */

/*
gsl_linalg_complex_cholesky_svx()
  Solve A x = b in place where A is in cholesky form
*/

int
gsl_linalg_complex_cholesky_svx (const gsl_matrix_complex * cholesky,
                                 gsl_vector_complex * x)
{
  if (cholesky->size1 != cholesky->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (cholesky->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* solve for y using forward-substitution, L y = b */

      gsl_blas_ztrsv (CblasLower, CblasNoTrans, CblasNonUnit, cholesky, x);

      /* perform back-substitution, L^H x = y */

      gsl_blas_ztrsv (CblasLower, CblasConjTrans, CblasNonUnit, cholesky, x);

      return GSL_SUCCESS;
    }
} /* gsl_linalg_complex_cholesky_svx() */

/********************************************
 *           INTERNAL ROUTINES              *
 ********************************************/

static void
cholesky_complex_conj_vector(gsl_vector_complex *v)
{
  size_t i;

  for (i = 0; i < v->size; ++i)
    {
      gsl_complex z = gsl_vector_complex_get(v, i);
      gsl_vector_complex_set(v, i, gsl_complex_conjugate(z));
    }
} /* cholesky_complex_conj_vector() */


text:
#include <jni.h>
#include <omp.h>
#ifdef __INTEL_COMPILER
#include <mkl.h>
#include <mkl_trans.h>
#else
#include <cblas.h>
#endif

JNIEXPORT jdouble JNICALL Java_edu_berkeley_bid_CBLAS_ddot
(JNIEnv * env, jobject calling_obj, jint N, jdoubleArray jX, jint incX, jdoubleArray jY, jint incY){
	jdouble * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jdouble * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);
	jdouble returnValue;

	returnValue = cblas_ddot(N, X, incX, Y, incY);

	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
	return returnValue;
}

JNIEXPORT jdouble JNICALL Java_edu_berkeley_bid_CBLAS_ddotxx
(JNIEnv * env, jobject calling_obj, jint N, jdoubleArray jX, jint startX, jdoubleArray jY, jint startY){
	jdouble * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jdouble * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);
	jdouble returnValue;

	returnValue = cblas_ddot(N, X+startX, 1, Y+startY, 1);

	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
	return returnValue;
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_ddotm
(JNIEnv * env, jobject calling_obj, jint nrows, jint ncols, jdoubleArray jX, jint ldx, jdoubleArray jY, jint ldy, jdoubleArray jZ){
	jdouble * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jdouble * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);
	jdouble * Z = (*env)->GetPrimitiveArrayCritical(env, jZ, JNI_FALSE);
        int i;

    for (i = 0; i < ncols; i++) {
      Z[i] = cblas_ddot(nrows, X+i*ldx, 1, Y+i*ldy, 1);
    }

	(*env)->ReleasePrimitiveArrayCritical(env, jZ, Z, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_ddotr
(JNIEnv * env, jobject calling_obj, jint nrows, jint ncols, jdoubleArray jX, jint ldx, jdoubleArray jY, jint ldy, jdoubleArray jZ){
	jdouble * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jdouble * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);
	jdouble * Z = (*env)->GetPrimitiveArrayCritical(env, jZ, JNI_FALSE);
    int i, j;

    for (i = 0; i < ncols; i++) {
#pragma omp parallel for
      for (j = 0; j < nrows; j++) {
        Z[j] += X[j + i*ldx] * Y[j + i*ldy];
      }
    }

	(*env)->ReleasePrimitiveArrayCritical(env, jZ, Z, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_daxpy
(JNIEnv * env, jobject calling_obj, jint N, jdouble a, jdoubleArray jX, jint incX, jdoubleArray jY, jint incY){
	jdouble * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jdouble * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);

        cblas_daxpy(N, a, X, incX, Y, incY);

	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_daxpyxx
(JNIEnv * env, jobject calling_obj, jint N, jdouble a, jdoubleArray jX, jint startX, jdoubleArray jY, jint startY){
	jdouble * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jdouble * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);

        cblas_daxpy(N, a, X+startX, 1, Y+startY, 1);

	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_dgemv
(JNIEnv * env, jobject calling_obj, jint order, jint transA, jint M, jint N, jdouble alpha,
 jdoubleArray jA, jint lda, jdoubleArray jX, jint incX, jdouble beta, jdoubleArray jY, jint incY){
	jdouble * A = (*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE);
	jdouble * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jdouble * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);

	cblas_dgemv((CBLAS_ORDER)order, (CBLAS_TRANSPOSE)transA, M, N, alpha, A, lda, X, incX, beta, Y, incY);

	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_dgemm
(JNIEnv * env, jobject calling_obj, jint order, jint transA, jint transB, jint M, jint N, jint K,
 jdouble alpha, jdoubleArray jA, jint lda, jdoubleArray jB, jint ldb, jdouble beta, jdoubleArray jC, jint ldc){
	jdouble * A = (*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE);
	jdouble * B = (*env)->GetPrimitiveArrayCritical(env, jB, JNI_FALSE);
	jdouble * C = (*env)->GetPrimitiveArrayCritical(env, jC, JNI_FALSE);

	cblas_dgemm((CBLAS_ORDER)order, (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB, M, N, K,
                    alpha, A, lda, B, ldb, beta, C, ldc);

	(*env)->ReleasePrimitiveArrayCritical(env, jC, C, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jB, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_domatcopy
(JNIEnv * env, jobject calling_obj, jstring j_order, jstring j_transA, jint M, jint N,
 jdouble alpha, jdoubleArray j_A, jint lda, jdoubleArray j_B, jint ldb) {
	char * order = (char *)(*env)->GetStringUTFChars(env, j_order, 0);
	char * transA = (char *)(*env)->GetStringUTFChars(env, j_transA, 0);
	jdouble * A = (*env)->GetPrimitiveArrayCritical(env, j_A, JNI_FALSE);
	jdouble * B = (*env)->GetPrimitiveArrayCritical(env, j_B, JNI_FALSE);

#ifdef __INTEL_COMPILER
	mkl_domatcopy(order[0], transA[0], M, N, alpha, A, lda, B, ldb);
#else
	int corder = (order[0] == 'r' || order[0] == 'R') ? CblasRowMajor : CblasColMajor;
	int ctrans = (transA[0] == 't' || transA[0] == 'T') ? CblasTrans : CblasNoTrans;
	cblas_domatcopy(corder, ctrans, M, N, alpha, A, lda, B, ldb);
#endif

	(*env)->ReleasePrimitiveArrayCritical(env, j_B, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_A, A, 0);
	(*env)->ReleaseStringUTFChars(env, j_transA, transA);
	(*env)->ReleaseStringUTFChars(env, j_order, order);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_dmcscm
(JNIEnv * env, jobject calling_obj, jint M, jint N, jdoubleArray j_A, jint lda,
 jdoubleArray j_B, jintArray j_ir, jintArray j_jc, jdoubleArray j_C, jint ldc){
	jdouble * A = (*env)->GetPrimitiveArrayCritical(env, j_A, JNI_FALSE);
	jdouble * B = (*env)->GetPrimitiveArrayCritical(env, j_B, JNI_FALSE);
	jint * ir = (*env)->GetPrimitiveArrayCritical(env, j_ir, JNI_FALSE);
	jint * jc = (*env)->GetPrimitiveArrayCritical(env, j_jc, JNI_FALSE);
	jdouble * C = (*env)->GetPrimitiveArrayCritical(env, j_C, JNI_FALSE);

        int ioff = jc[0];
        int i;
#pragma omp parallel for
        for (i = 0; i < N; i++) {
          int j, ir0;
          for (j = jc[i]-ioff; j < jc[i+1]-ioff; j++) {
            ir0 = ir[j]-ioff;
            cblas_daxpy(M, B[j], A+(ir0*lda), 1, C+(i*ldc), 1);
          }
        }

	(*env)->ReleasePrimitiveArrayCritical(env, j_C, C, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_jc, jc, 0);
    (*env)->ReleasePrimitiveArrayCritical(env, j_ir, ir, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_B, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_A, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_dmcsrm
(JNIEnv * env, jobject calling_obj, jint M, jint N, jdoubleArray j_A, jint lda,
 jdoubleArray j_B, jintArray j_ir, jintArray j_jc, jdoubleArray j_C, jint ldc){
	jdouble * A = (*env)->GetPrimitiveArrayCritical(env, j_A, JNI_FALSE);
	jdouble * B = (*env)->GetPrimitiveArrayCritical(env, j_B, JNI_FALSE);
	jint * ir = (*env)->GetPrimitiveArrayCritical(env, j_ir, JNI_FALSE);
	jint * jc = (*env)->GetPrimitiveArrayCritical(env, j_jc, JNI_FALSE);
	jdouble * C = (*env)->GetPrimitiveArrayCritical(env, j_C, JNI_FALSE);

    int ioff = jc[0];
    int i;
#pragma omp parallel for
    for (i = 0; i < N; i++) {
      int j, k;
      for (j = jc[i]-ioff; j < jc[i+1]-ioff; j++) {
        k = ir[j]-ioff;
        cblas_daxpy(M, B[j], A+(i*lda), 1, C+(k*ldc), 1);
      }
    }

	(*env)->ReleasePrimitiveArrayCritical(env, j_C, C, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_jc, jc, 0);
    (*env)->ReleasePrimitiveArrayCritical(env, j_ir, ir, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_B, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_A, A, 0);
}

JNIEXPORT jfloat JNICALL Java_edu_berkeley_bid_CBLAS_sdot
(JNIEnv * env, jobject calling_obj, jint N, jfloatArray jX, jint incX, jfloatArray jY, jint incY){
	jfloat * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jfloat * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);
	jfloat returnValue;

	returnValue = cblas_sdot(N, X, incX, Y, incY);

	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
	return returnValue;
}

JNIEXPORT jfloat JNICALL Java_edu_berkeley_bid_CBLAS_sdotxx
(JNIEnv * env, jobject calling_obj, jint N, jfloatArray jX, jint startX, jfloatArray jY, jint startY){
	jfloat * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jfloat * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);
	jfloat returnValue;

	returnValue = cblas_sdot(N, X+startX, 1, Y+startY, 1);

	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
	return returnValue;
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_sdotm
(JNIEnv * env, jobject calling_obj, jint nrows, jint ncols, jfloatArray jX, jint ldx, jfloatArray jY, jint ldy, jfloatArray jZ){
	jfloat * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jfloat * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);
	jfloat * Z = (*env)->GetPrimitiveArrayCritical(env, jZ, JNI_FALSE);
    int i, j;

    for (i = 0; i < ncols; i++) {
      Z[i] = cblas_sdot(nrows, X+i*ldx, 1, Y+i*ldy, 1);
    }

	(*env)->ReleasePrimitiveArrayCritical(env, jZ, Z, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_sdotr
(JNIEnv * env, jobject calling_obj, jint nrows, jint ncols, jfloatArray jX, jint ldx, jfloatArray jY, jint ldy, jfloatArray jZ){
	jfloat * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jfloat * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);
	jfloat * Z = (*env)->GetPrimitiveArrayCritical(env, jZ, JNI_FALSE);
    int i, j;

    for (i = 0; i < ncols; i++) {
#pragma omp parallel for
      for (j = 0; j < nrows; j++) {
        Z[j] += X[j + i*ldx] * Y[j + i*ldy];
      }
    }

	(*env)->ReleasePrimitiveArrayCritical(env, jZ, Z, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_sgemv
(JNIEnv * env, jobject calling_obj, jint order, jint transA, jint M, jint N, jfloat alpha,
jfloatArray jA, jint lda, jfloatArray jX, jint incX, jfloat beta, jfloatArray jY, jint incY){
	jfloat * A = (*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE);
	jfloat * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jfloat * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);

	cblas_sgemv((CBLAS_ORDER)order, (CBLAS_TRANSPOSE)transA, M, N, alpha, A, lda, X, incX, beta, Y, incY);

	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_sgemm
(JNIEnv * env, jobject calling_obj, jint order, jint transA, jint transB, jint M, jint N, jint K,
 jfloat alpha, jfloatArray jA, jint lda, jfloatArray jB, jint ldb, jfloat beta, jfloatArray jC, jint ldc){
	jfloat * A = (*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE);
	jfloat * B = (*env)->GetPrimitiveArrayCritical(env, jB, JNI_FALSE);
	jfloat * C = (*env)->GetPrimitiveArrayCritical(env, jC, JNI_FALSE);

	cblas_sgemm((CBLAS_ORDER)order, (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB, M, N, K,
                    alpha, A, lda, B, ldb, beta, C, ldc);

	(*env)->ReleasePrimitiveArrayCritical(env, jC, C, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jB, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_sgemmx
(JNIEnv * env, jobject calling_obj, jint order, jint transA, jint transB, jint M, jint N, jint K,
 jfloat alpha, jfloatArray jA, jint Aoff, jint lda, jfloatArray jB, jint Boff, jint ldb, jfloat beta, jfloatArray jC, jint Coff, jint ldc){
	jfloat * A = (*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE);
	jfloat * B = (*env)->GetPrimitiveArrayCritical(env, jB, JNI_FALSE);
	jfloat * C = (*env)->GetPrimitiveArrayCritical(env, jC, JNI_FALSE);

	cblas_sgemm((CBLAS_ORDER)order, (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB, M, N, K,
                    alpha, A+Aoff, lda, B+Boff, ldb, beta, C+Coff, ldc);

	(*env)->ReleasePrimitiveArrayCritical(env, jC, C, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jB, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_somatcopy
(JNIEnv * env, jobject calling_obj, jstring j_order, jstring j_transA, jint M, jint N,
 jfloat alpha, jfloatArray j_A, jint lda, jfloatArray j_B, jint ldb) {
	char * order = (char *)(*env)->GetStringUTFChars(env, j_order, 0);
	char * transA = (char *)(*env)->GetStringUTFChars(env, j_transA, 0);
	jfloat * A = (*env)->GetPrimitiveArrayCritical(env, j_A, JNI_FALSE);
	jfloat * B = (*env)->GetPrimitiveArrayCritical(env, j_B, JNI_FALSE);

#ifdef __INTEL_COMPILER
	mkl_somatcopy(order[0], transA[0], M, N, alpha, A, lda, B, ldb);
#else
	int corder = (order[0] == 'r' || order[0] == 'R') ? CblasRowMajor : CblasColMajor;
	int ctrans = (transA[0] == 't' || transA[0] == 'T') ? CblasTrans : CblasNoTrans;
	cblas_somatcopy(corder, ctrans, M, N, alpha, A, lda, B, ldb);
#endif

	(*env)->ReleasePrimitiveArrayCritical(env, j_B, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_A, A, 0);
	(*env)->ReleaseStringUTFChars(env, j_transA, transA);
	(*env)->ReleaseStringUTFChars(env, j_order, order);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_iomatcopy
(JNIEnv * env, jobject calling_obj, jstring j_order, jstring j_transA, jint M, jint N,
 jintArray j_A, jint lda, jintArray j_B, jint ldb) {
	char * order = (char *)(*env)->GetStringUTFChars(env, j_order, 0);
	char * transA = (char *)(*env)->GetStringUTFChars(env, j_transA, 0);
	jfloat * A = (*env)->GetPrimitiveArrayCritical(env, j_A, JNI_FALSE);
	jfloat * B = (*env)->GetPrimitiveArrayCritical(env, j_B, JNI_FALSE);

#ifdef __INTEL_COMPILER
	mkl_somatcopy(order[0], transA[0], M, N, 1.0f, A, lda, B, ldb);
#else
	int corder = (order[0] == 'r' || order[0] == 'R') ? CblasRowMajor : CblasColMajor;
	int ctrans = (transA[0] == 't' || transA[0] == 'T') ? CblasTrans : CblasNoTrans;
	cblas_somatcopy(corder, ctrans, M, N, 1.0f, A, lda, B, ldb);
#endif
	(*env)->ReleasePrimitiveArrayCritical(env, j_B, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_A, A, 0);
	(*env)->ReleaseStringUTFChars(env, j_transA, transA);
	(*env)->ReleaseStringUTFChars(env, j_order, order);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_lomatcopy
(JNIEnv * env, jobject calling_obj, jstring j_order, jstring j_transA, jint M, jint N,
 jlongArray j_A, jint lda, jlongArray j_B, jint ldb) {
	char * order = (char *)(*env)->GetStringUTFChars(env, j_order, 0);
	char * transA = (char *)(*env)->GetStringUTFChars(env, j_transA, 0);
	jdouble * A = (*env)->GetPrimitiveArrayCritical(env, j_A, JNI_FALSE);
	jdouble * B = (*env)->GetPrimitiveArrayCritical(env, j_B, JNI_FALSE);

#ifdef __INTEL_COMPILER
	mkl_domatcopy(order[0], transA[0], M, N, 1.0, A, lda, B, ldb);
#else
	int corder = (order[0] == 'r' || order[0] == 'R') ? CblasRowMajor : CblasColMajor;
	int ctrans = (transA[0] == 't' || transA[0] == 'T') ? CblasTrans : CblasNoTrans;
	cblas_domatcopy(corder, ctrans, M, N, 1.0, A, lda, B, ldb);
#endif

	(*env)->ReleasePrimitiveArrayCritical(env, j_B, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_A, A, 0);
	(*env)->ReleaseStringUTFChars(env, j_transA, transA);
	(*env)->ReleaseStringUTFChars(env, j_order, order);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_spermute
(JNIEnv * env, jobject calling_obj, jint M, jint N, jint K, jfloatArray j_A, jfloatArray j_B) {
        int i, offset, step = M*N;
	jfloat * A = (*env)->GetPrimitiveArrayCritical(env, j_A, JNI_FALSE);
	jfloat * B = (*env)->GetPrimitiveArrayCritical(env, j_B, JNI_FALSE);

#ifdef __INTEL_COMPILER
        for (i = 0, offset = 0; i < K; i++, offset += step) {
          mkl_somatcopy('C', 'T', M, N, 1.0f, A+offset, M, B+offset, N);
        }
#else
        for (i = 0, offset = 0; i < K; i++, offset += step) {
          cblas_somatcopy(CblasColMajor, CblasTrans, M, N, 1.0f, A+offset, M, B+offset, N);
        }
#endif
	(*env)->ReleasePrimitiveArrayCritical(env, j_B, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_A, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_ipermute
(JNIEnv * env, jobject calling_obj, jint M, jint N, jint K, jintArray j_A, jintArray j_B) {
        int i, offset, step = M*N;
	jfloat * A = (jfloat *)((*env)->GetPrimitiveArrayCritical(env, j_A, JNI_FALSE));
	jfloat * B = (jfloat *)((*env)->GetPrimitiveArrayCritical(env, j_B, JNI_FALSE));

#ifdef __INTEL_COMPILER
        for (i = 0, offset = 0; i < K; i++, offset += step) {
          mkl_somatcopy('C', 'T', M, N, 1.0f, A+offset, M, B+offset, N);
        }
#else
        for (i = 0, offset = 0; i < K; i++, offset += step) {
          cblas_somatcopy(CblasColMajor, CblasTrans, M, N, 1.0f, A+offset, M, B+offset, N);
        }
#endif
	(*env)->ReleasePrimitiveArrayCritical(env, j_B, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_A, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_dpermute
(JNIEnv * env, jobject calling_obj, jint M, jint N, jint K, jdoubleArray j_A, jdoubleArray j_B) {
        int i, offset, step = M*N;
	jdouble * A = (*env)->GetPrimitiveArrayCritical(env, j_A, JNI_FALSE);
	jdouble * B = (*env)->GetPrimitiveArrayCritical(env, j_B, JNI_FALSE);

#ifdef __INTEL_COMPILER
        for (i = 0, offset = 0; i < K; i++, offset += step) {
          mkl_domatcopy('C', 'T', M, N, 1.0, A+offset, M, B+offset, N);
        }
#else
        for (i = 0, offset = 0; i < K; i++, offset += step) {
          cblas_domatcopy(CblasColMajor, CblasTrans, M, N, 1.0, A+offset, M, B+offset, N);
        }
#endif

	(*env)->ReleasePrimitiveArrayCritical(env, j_B, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_A, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_lpermute
(JNIEnv * env, jobject calling_obj, jint M, jint N, jint K, jlongArray j_A, jlongArray j_B) {
        int i, offset, step = M*N;
	jdouble * A = (jdouble *)((*env)->GetPrimitiveArrayCritical(env, j_A, JNI_FALSE));
	jdouble * B = (jdouble *)((*env)->GetPrimitiveArrayCritical(env, j_B, JNI_FALSE));

#ifdef __INTEL_COMPILER
        for (i = 0, offset = 0; i < K; i++, offset += step) {
          mkl_domatcopy('C', 'T', M, N, 1.0, A+offset, M, B+offset, N);
        }
#else
        for (i = 0, offset = 0; i < K; i++, offset += step) {
          cblas_domatcopy(CblasColMajor, CblasTrans, M, N, 1.0, A+offset, M, B+offset, N);
        }
#endif

	(*env)->ReleasePrimitiveArrayCritical(env, j_B, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_A, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_saxpy
(JNIEnv * env, jobject calling_obj, jint N, jfloat a, jfloatArray jX, jint incX, jfloatArray jY, jint incY){
	jfloat * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jfloat * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);

        cblas_saxpy(N, a, X, incX, Y, incY);

	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_saxpyxx
(JNIEnv * env, jobject calling_obj, jint N, jfloat a, jfloatArray jX, jint startX, jfloatArray jY, jint startY){
	jfloat * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jfloat * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);

        cblas_saxpy(N, a, X+startX, 1, Y+startY, 1);

	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_smcscm
(JNIEnv * env, jobject calling_obj, jint M, jint N, jfloatArray j_A, jint lda,
 jfloatArray j_B, jintArray j_ir, jintArray j_jc, jfloatArray j_C, jint ldc){
	jfloat * A = (*env)->GetPrimitiveArrayCritical(env, j_A, JNI_FALSE);
	jfloat * B = (*env)->GetPrimitiveArrayCritical(env, j_B, JNI_FALSE);
	jint * ir = (*env)->GetPrimitiveArrayCritical(env, j_ir, JNI_FALSE);
	jint * jc = (*env)->GetPrimitiveArrayCritical(env, j_jc, JNI_FALSE);
	jfloat * C = (*env)->GetPrimitiveArrayCritical(env, j_C, JNI_FALSE);

        int ioff = jc[0];
        int i;
#pragma omp parallel for
        for (i = 0; i < N; i++) {
          int j, ir0;
          for (j = jc[i]-ioff; j < jc[i+1]-ioff; j++) {
            ir0 = ir[j]-ioff;
            cblas_saxpy(M, B[j], A+(ir0*lda), 1, C+(i*ldc), 1);
          }
        }

	(*env)->ReleasePrimitiveArrayCritical(env, j_C, C, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_jc, jc, 0);
        (*env)->ReleasePrimitiveArrayCritical(env, j_ir, ir, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_B, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_A, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_smcsrm
(JNIEnv * env, jobject calling_obj, jint M, jint N, jfloatArray j_A, jint lda,
 jfloatArray j_B, jintArray j_ir, jintArray j_jc, jfloatArray j_C, jint ldc){
	jfloat * A = (*env)->GetPrimitiveArrayCritical(env, j_A, JNI_FALSE);
	jfloat * B = (*env)->GetPrimitiveArrayCritical(env, j_B, JNI_FALSE);
	jint * ir = (*env)->GetPrimitiveArrayCritical(env, j_ir, JNI_FALSE);
	jint * jc = (*env)->GetPrimitiveArrayCritical(env, j_jc, JNI_FALSE);
	jfloat * C = (*env)->GetPrimitiveArrayCritical(env, j_C, JNI_FALSE);

        int ioff = jc[0];
        int i;
#pragma omp parallel for
        for (i = 0; i < N; i++) {
          int j, jj, k;
          for (j = jc[i]-ioff; j < jc[i+1]-ioff; j++) {
            jj = ir[j]-ioff;
            if (M == 1) {
              C[jj*ldc] += B[j] * A[i*lda];
            } else if (M > 10) {
              cblas_saxpy(M, B[j], A+(i*lda), 1, C+(jj*ldc), 1);
            } else {
              int iia = i*lda;
              int jjc = jj*ldc;
              float Bj = B[j];
              for (k = 0; k < M; k++) {
                C[jjc+k] += Bj * A[iia+k];
              }
            }
          }
        }

	(*env)->ReleasePrimitiveArrayCritical(env, j_C, C, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_jc, jc, 0);
        (*env)->ReleasePrimitiveArrayCritical(env, j_ir, ir, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_B, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, j_A, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_cdot
(JNIEnv * env, jobject calling_obj, jint N, jfloatArray jX, jint incX, jfloatArray jY, jint incY, jfloatArray jZ){
	jfloat * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jfloat * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);
	jfloat * Z = (*env)->GetPrimitiveArrayCritical(env, jZ, JNI_FALSE);

    cblas_cdotu_sub(N, X, incX, Y, incY, Z);

	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jZ, Z, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_cdotxx
(JNIEnv * env, jobject calling_obj, jint N, jfloatArray jX, jint startX, jfloatArray jY, jint startY, jfloatArray jZ){
	jfloat * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jfloat * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);
	jfloat * Z = (*env)->GetPrimitiveArrayCritical(env, jZ, JNI_FALSE);

    cblas_cdotu_sub(N, X+startX, 1, Y+startY, 1, Z);

	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jZ, Z, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_cdotm
(JNIEnv * env, jobject calling_obj, jint nrows, jint ncols, jfloatArray jX, jint ldx, jfloatArray jY, jint ldy, jfloatArray jZ){
	jfloat * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jfloat * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);
	jfloat * Z = (*env)->GetPrimitiveArrayCritical(env, jZ, JNI_FALSE);
        int i;

    for (i=0; i<2*ncols; i+=2) {
      cblas_cdotu_sub(nrows, X+i*ldx, 1, Y+i*ldy, 1, Z+i);
    }

	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jZ, Z, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_cdotr
(JNIEnv * env, jobject calling_obj, jint nrows, jint ncols, jfloatArray jX, jint ldx, jfloatArray jY, jint ldy, jfloatArray jZ){
	jfloat * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jfloat * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);
	jfloat * Z = (*env)->GetPrimitiveArrayCritical(env, jZ, JNI_FALSE);
    int i, j;

    for (i = 0; i < ncols; i++) {
#pragma omp parallel for
      for (j = 0; j < nrows; j++) {
        int ix, iy;
        ix = 2*(j + i*ldx);
        iy = 2*(j + i*ldy);
        Z[2*j] += X[ix] * Y[ix] - X[ix+1] * Y[ix+1];
        Z[2*j+1] += X[ix] * Y[ix+1] + X[ix+1] * Y[ix];
      }
    }

	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jZ, Z, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_cgemv
(JNIEnv * env, jobject calling_obj, jint order, jint transA, jint M, jint N, jfloatArray jAlpha,
jfloatArray jA, jint lda, jfloatArray jX, jint incX, jfloatArray jBeta, jfloatArray jY, jint incY){
	jfloat * A = (*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE);
	jfloat * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jfloat * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);
	jfloat * alpha = (*env)->GetPrimitiveArrayCritical(env, jAlpha, JNI_FALSE);
	jfloat * beta = (*env)->GetPrimitiveArrayCritical(env, jBeta, JNI_FALSE);

	cblas_cgemv((CBLAS_ORDER)order, (CBLAS_TRANSPOSE)transA, M, N, alpha, A, lda, X, incX, beta, Y, incY);

	(*env)->ReleasePrimitiveArrayCritical(env, jBeta, beta, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jAlpha, alpha, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_cgemm
(JNIEnv * env, jobject calling_obj, jint order, jint transA, jint transB, jint M, jint N, jint K,
jfloatArray jAlpha, jfloatArray jA, jint lda, jfloatArray jB, jint ldb, jfloatArray jBeta, jfloatArray jC, jint ldc){
	jfloat * A = (*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE);
	jfloat * B = (*env)->GetPrimitiveArrayCritical(env, jB, JNI_FALSE);
	jfloat * C = (*env)->GetPrimitiveArrayCritical(env, jC, JNI_FALSE);
	jfloat * alpha = (*env)->GetPrimitiveArrayCritical(env, jAlpha, JNI_FALSE);
	jfloat * beta = (*env)->GetPrimitiveArrayCritical(env, jBeta, JNI_FALSE);

	cblas_cgemm((CBLAS_ORDER)order, (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB, M, N, K,
                    alpha, A, lda, B, ldb, beta, C, ldc);

	(*env)->ReleasePrimitiveArrayCritical(env, jC, C, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jB, B, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_caxpy
(JNIEnv * env, jobject calling_obj, jint N, jfloatArray jA, jfloatArray jX, jint incX,
 jfloatArray jY, jint incY){
	jfloat * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jfloat * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);
	jfloat * a = (*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE);

     cblas_caxpy(N, a, X, incX, Y, incY);

	(*env)->ReleasePrimitiveArrayCritical(env, jA, a, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_caxpyxx
(JNIEnv * env, jobject calling_obj, jint N, jfloatArray jA, jfloatArray jX, jint startX, jfloatArray jY, jint startY){
	jfloat * X = (*env)->GetPrimitiveArrayCritical(env, jX, JNI_FALSE);
	jfloat * Y = (*env)->GetPrimitiveArrayCritical(env, jY, JNI_FALSE);
	jfloat * a = (*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE);

    cblas_caxpy(N, a, X+startX, 1, Y+startY, 1);

	(*env)->ReleasePrimitiveArrayCritical(env, jA, a, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jY, Y, 0);
	(*env)->ReleasePrimitiveArrayCritical(env, jX, X, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_blockSgemm
(JNIEnv *env, jobject obj, jint transA, jint transB, jint nr, jint nc, jint kk, jfloat alpha, jfloatArray jA, jint aoff, jint lda, jint astep,
 jfloatArray jB, jint boff, jint ldb, jint bstep, jfloat beta, jfloatArray jC, jint coff, jint ldc, jint cstep, jint reps)
{
  int i, at, bt;
  jfloat *A = (*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE);
  jfloat *B = (*env)->GetPrimitiveArrayCritical(env, jB, JNI_FALSE);
  jfloat *C = (*env)->GetPrimitiveArrayCritical(env, jC, JNI_FALSE);

  at = (transA) ? CblasTrans : CblasNoTrans;
  bt = (transB) ? CblasTrans : CblasNoTrans;
  A += aoff;
  B += boff;
  C += coff;
  for (i = 0; i < reps; i++) {
    cblas_sgemm(CblasColMajor, at, bt, nr, nc, kk, alpha, A, lda, B, ldb, beta, C, ldc);
    A += astep;
    B += bstep;
    C += cstep;
  }
  (*env)->ReleasePrimitiveArrayCritical(env, jC, C, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jB, B, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_blockDgemm
(JNIEnv *env, jobject obj, jint transA, jint transB, jint nr, jint nc, jint kk, jdouble alpha, jdoubleArray jA, jint aoff, jint lda, jint astep,
 jdoubleArray jB, jint boff, jint ldb, jint bstep, jdouble beta, jdoubleArray jC, jint coff, jint ldc, jint cstep, jint reps)
{
  int i, at, bt;
  jdouble *A = (*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE);
  jdouble *B = (*env)->GetPrimitiveArrayCritical(env, jB, JNI_FALSE);
  jdouble *C = (*env)->GetPrimitiveArrayCritical(env, jC, JNI_FALSE);

  at = (transA) ? CblasTrans : CblasNoTrans;
  bt = (transB) ? CblasTrans : CblasNoTrans;
  A += aoff;
  B += boff;
  C += coff;
  for (i = 0; i < reps; i++) {
    cblas_dgemm(CblasColMajor, at, bt, nr, nc, kk, alpha, A, lda, B, ldb, beta, C, ldc);
    A += astep;
    B += bstep;
    C += cstep;
  }
  (*env)->ReleasePrimitiveArrayCritical(env, jC, C, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jB, B, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_blockSgemm4D
(JNIEnv *env, jobject obj, jint transA, jint transB, jint nr, jint nc, jint kk, jfloat alpha,
 jfloatArray jA, jint aoff, jint lda, jint astep1, jint astep2,
 jfloatArray jB, jint boff, jint ldb, jint bstep1, jint bstep2, jfloat beta,
 jfloatArray jC, jint coff, jint ldc, jint cstep1, jint cstep2, jint reps1, jint reps2)
{
  int i, j, at, bt;
  float *A = (float *)((*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE));
  float *B = (float *)((*env)->GetPrimitiveArrayCritical(env, jB, JNI_FALSE));
  float *C = (float *)((*env)->GetPrimitiveArrayCritical(env, jC, JNI_FALSE));
  float *pA, *pB, *pC;

  at = (transA) ? CblasTrans : CblasNoTrans;
  bt = (transB) ? CblasTrans : CblasNoTrans;
  A += aoff;
  B += boff;
  C += coff;
  for (i = 0; i < reps2; i++) {
    for (j = 0; j < reps1; j++) {
      pA = A + (j * astep1 + i * astep2);
      pB = B + (j * bstep1 + i * bstep2);
      pC = C + (j * cstep1 + i * cstep2);
      cblas_sgemm(CblasColMajor, at, bt, nr, nc, kk, alpha, pA, lda, pB, ldb, beta, pC, ldc);
    }
  }
  (*env)->ReleasePrimitiveArrayCritical(env, jC, C, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jB, B, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_blockDgemm4D
(JNIEnv *env, jobject obj, jint transA, jint transB, jint nr, jint nc, jint kk, jdouble alpha,
 jdoubleArray jA, jint aoff, jint lda, jint astep1, jint astep2,
 jdoubleArray jB, jint boff, jint ldb, jint bstep1, jint bstep2, jdouble beta,
 jdoubleArray jC, jint coff, jint ldc, jint cstep1, jint cstep2, jint reps1, jint reps2)
{
  int i, j, at, bt;
  double *A = (double *)((*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE));
  double *B = (double *)((*env)->GetPrimitiveArrayCritical(env, jB, JNI_FALSE));
  double *C = (double *)((*env)->GetPrimitiveArrayCritical(env, jC, JNI_FALSE));
  double *pA, *pB, *pC;

  at = (transA) ? CblasTrans : CblasNoTrans;
  bt = (transB) ? CblasTrans : CblasNoTrans;
  A += aoff;
  B += boff;
  C += coff;
  for (i = 0; i < reps2; i++) {
    for (j = 0; j < reps1; j++) {
      pA = A + (j * astep1 + i * astep2);
      pB = B + (j * bstep1 + i * bstep2);
      pC = C + (j * cstep1 + i * cstep2);
      cblas_dgemm(CblasColMajor, at, bt, nr, nc, kk, alpha, pA, lda, pB, ldb, beta, pC, ldc);
    }
  }
  (*env)->ReleasePrimitiveArrayCritical(env, jC, C, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jB, B, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_word2vecFwd
(JNIEnv *env, jobject obj, jint nrows, jint ncols, const jint nwa, const jint nwb, jintArray jWA, jintArray jWB,
 jfloatArray jA, jfloatArray jB, jfloatArray jC)
{
  jint * WA = (jint *)((*env)->GetPrimitiveArrayCritical(env, jWA, JNI_FALSE));
  jint * WB = (jint *)((*env)->GetPrimitiveArrayCritical(env, jWB, JNI_FALSE));
  jfloat * A = (jfloat *)((*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE));
  jfloat * B = (jfloat *)((*env)->GetPrimitiveArrayCritical(env, jB, JNI_FALSE));
  jfloat * C = (jfloat *)((*env)->GetPrimitiveArrayCritical(env, jC, JNI_FALSE));
  int i;
#pragma omp parallel for
  for (i = 0; i < ncols; i++) {
    int j, k, c, ia, ib, coff;
    float sum;
    for (j = 0; j < nwa; j++) {
      ia = nrows*WA[j+i*nwa];
      for (k = 0; k < nwb; k++) {
        ib = nrows*WB[k+i*nwb];
        sum = 0;
        for (c = 0; c < nrows; c++) {
          sum += A[c + ia] * B[c + ib];
        }
        coff = nwa * (k + nwb * i);
        C[j + coff] = sum;
      }
    }
  }

  (*env)->ReleasePrimitiveArrayCritical(env, jC, C, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jB, B, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jWB, WB, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jWA, WA, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_word2vecBwd
(JNIEnv *env, jobject obj, jint nrows, jint ncols, jint nwa, jint nwb, jintArray jWA, jintArray jWB,
 jfloatArray jA, jfloatArray jB, jfloatArray jDA, jfloatArray jDB, jfloatArray jC, jfloat lrate)
{
  jint * WA = (jint *)((*env)->GetPrimitiveArrayCritical(env, jWA, JNI_FALSE));
  jint * WB = (jint *)((*env)->GetPrimitiveArrayCritical(env, jWB, JNI_FALSE));
  jfloat * A = (jfloat *)((*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE));
  jfloat * B = (jfloat *)((*env)->GetPrimitiveArrayCritical(env, jB, JNI_FALSE));
  jfloat * DA = (jfloat *)((*env)->GetPrimitiveArrayCritical(env, jDA, JNI_FALSE));
  jfloat * DB = (jfloat *)((*env)->GetPrimitiveArrayCritical(env, jDB, JNI_FALSE));
  jfloat * C = (jfloat *)((*env)->GetPrimitiveArrayCritical(env, jC, JNI_FALSE));
  int i;

#pragma omp parallel for
  for (i = 0; i < ncols; i++) {
    int j, k, c;
    float cv;
    int ia, ib;
    for (j = 0; j < nwa; j++) {
      ia = nrows*WA[j+i*nwa];
      for (k = 0; k < nwb; k++) {
        ib = nrows*WB[k+i*nwb];
        cv = lrate * C[j + nwa * (k + nwb * i)];
        for (c = 0; c < nrows; c++) {
          A[c + ia] += cv * DB[c + ib];
        }
      }
    }
    for (k = 0; k < nwb; k++) {
      ib = nrows*WB[k+i*nwb];
      for (j = 0; j < nwa; j++) {
        ia = nrows*WA[j+i*nwa];
        cv = lrate * C[j + nwa * (k + nwb * i)];
        for (c = 0; c < nrows; c++) {
          B[c + ib] += cv * DA[c + ia];
        }
      }
    }
  }
  (*env)->ReleasePrimitiveArrayCritical(env, jC, C, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jDB, DB, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jDA, DA, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jB, B, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jWB, WB, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jWA, WA, 0);
}

#define mymax(A, B) ((A>B) ? A : B)

#define mymin(A, B) ((A<B) ? A : B)

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_reduceTensorFloat
(JNIEnv *env, jobject obj, jfloatArray jA, jfloatArray jB, jint m, jint n, jint p, jint op)
{
  jfloat * A = (jfloat *)((*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE));
  jfloat * B = (jfloat *)((*env)->GetPrimitiveArrayCritical(env, jB, JNI_FALSE));
  int i, j;

  for (i = 0; i < p; i++) {
    int k;
    for (k = 0; k < m; k++) {
      B[k + m * i] = A[k + m * n * i];
    }
  }
  for (j = 1; j < n; j++) {
    switch (op) {
    case 0 : {
#pragma omp parallel for
      for (i = 0; i < p; i++) {
	int k;
	for (k = 0; k < m; k++) {
	  B[k + m * i] = B[k + m * i] + A[k + m * (j + n * i)];
	}
      }
    }
      break;
    case 1 : {
#pragma omp parallel for
      for (i = 0; i < p; i++) {
	int k;
	for (k = 0; k < m; k++) {
	  B[k + m * i] = B[k + m * i] * A[k + m * (j + n * i)];
	}
      }
    }
      break;
    case 2 : {
#pragma omp parallel for
      for (i = 0; i < p; i++) {
	int k;
	for (k = 0; k < m; k++) {
	  B[k + m * i] = mymax(B[k + m * i], A[k + m * (j + n * i)]);
	}
      }
    }
      break;
    case 3 : {
#pragma omp parallel for
      for (i = 0; i < p; i++) {
	int k;
	for (k = 0; k < m; k++) {
	  B[k + m * i] = mymin(B[k + m * i], A[k + m * (j + n * i)]);
	}
      }
    }
      break;
    default: {}
    }
  }
  (*env)->ReleasePrimitiveArrayCritical(env, jB, B, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_reduceTensorDouble
  (JNIEnv *env, jobject obj, jdoubleArray jA, jdoubleArray jB, jint m, jint n, jint p, jint op)
{
  jdouble * A = (jdouble *)((*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE));
  jdouble * B = (jdouble *)((*env)->GetPrimitiveArrayCritical(env, jB, JNI_FALSE));
  int i, j;

  for (i = 0; i < p; i++) {
    int k;
    for (k = 0; k < m; k++) {
      B[k + m * i] = A[k + m * n * i];
    }
  }
  for (j = 1; j < n; j++) {
    switch (op) {
    case 0 : {
#pragma omp parallel for
      for (i = 0; i < p; i++) {
	int k;
	for (k = 0; k < m; k++) {
	  B[k + m * i] += A[k + m * (j + n * i)];
	}
      }
    }
      break;
    case 1 : {
#pragma omp parallel for
      for (i = 0; i < p; i++) {
	int k;
	for (k = 0; k < m; k++) {
	  B[k + m * i] *= A[k + m * (j + n * i)];
	}
      }
    }
      break;
    case 2 : {
#pragma omp parallel for
      for (i = 0; i < p; i++) {
	int k;
	for (k = 0; k < m; k++) {
	  B[k + m * i] = mymax(B[k + m * i], A[k + m * (j + n * i)]);
	}
      }
    }
      break;
    case 3 : {
#pragma omp parallel for
      for (i = 0; i < p; i++) {
	int k;
	for (k = 0; k < m; k++) {
	  B[k + m * i] = mymin(B[k + m * i], A[k + m * (j + n * i)]);
	}
      }
    }
      break;
    default: {}
    }
  }
  (*env)->ReleasePrimitiveArrayCritical(env, jB, B, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
}

JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_reduceTensorInt
(JNIEnv *env, jobject obj, jintArray jA, jintArray jB, jint m, jint n, jint p, jint op)
{
  jint * A = (jint *)((*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE));
  jint * B = (jint *)((*env)->GetPrimitiveArrayCritical(env, jB, JNI_FALSE));
  int i, j;

  for (i = 0; i < p; i++) {
    int k;
    for (k = 0; k < m; k++) {
      B[k + m * i] = A[k + m * n * i];
    }
  }
  for (j = 1; j < n; j++) {
    switch (op) {
    case 0 : {
#pragma omp parallel for
      for (i = 0; i < p; i++) {
	int k;
	for (k = 0; k < m; k++) {
	  B[k + m * i] += A[k + m * (j + n * i)];
	}
      }
    }
      break;
    case 1 : {
#pragma omp parallel for
      for (i = 0; i < p; i++) {
	int k;
	for (k = 0; k < m; k++) {
	  B[k + m * i] *= A[k + m * (j + n * i)];
	}
      }
    }
      break;
    case 2 : {
#pragma omp parallel for
      for (i = 0; i < p; i++) {
	int k;
	for (k = 0; k < m; k++) {
	  B[k + m * i] = mymax(B[k + m * i], A[k + m * (j + n * i)]);
	}
      }
    }
      break;
    case 3 : {
#pragma omp parallel for
      for (i = 0; i < p; i++) {
	int k;
	for (k = 0; k < m; k++) {
	  B[k + m * i] = mymin(B[k + m * i], A[k + m * (j + n * i)]);
	}
      }
    }
      break;
    default: {}
    }
  }
  (*env)->ReleasePrimitiveArrayCritical(env, jB, B, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
}

  JNIEXPORT void JNICALL Java_edu_berkeley_bid_CBLAS_reduceTensorLong
(JNIEnv *env, jobject obj, jlongArray jA, jlongArray jB, jint m, jint n, jint p, jint op)
{
  jlong * A = (jlong *)((*env)->GetPrimitiveArrayCritical(env, jA, JNI_FALSE));
  jlong * B = (jlong *)((*env)->GetPrimitiveArrayCritical(env, jB, JNI_FALSE));
  int i, j;

  for (i = 0; i < p; i++) {
    int k;
    for (k = 0; k < m; k++) {
      B[k + m * i] = A[k + m * n * i];
    }
  }
  for (j = 1; j < n; j++) {
    switch (op) {
    case 0 : {
#pragma omp parallel for
      for (i = 0; i < p; i++) {
	int k;
	for (k = 0; k < m; k++) {
	  B[k + m * i] += A[k + m * (j + n * i)];
	}
      }
    }
      break;
    case 1 : {
#pragma omp parallel for
      for (i = 0; i < p; i++) {
	int k;
	for (k = 0; k < m; k++) {
	  B[k + m * i] *= A[k + m * (j + n * i)];
	}
      }
    }
      break;
    case 2 : {
#pragma omp parallel for
      for (i = 0; i < p; i++) {
	int k;
	for (k = 0; k < m; k++) {
	  B[k + m * i] = mymax(B[k + m * i], A[k + m * (j + n * i)]);
	}
      }
    }
      break;
    case 3 : {
#pragma omp parallel for
      for (i = 0; i < p; i++) {
	int k;
	for (k = 0; k < m; k++) {
	  B[k + m * i] = mymin(B[k + m * i], A[k + m * (j + n * i)]);
	}
      }
    }
      break;
    default: {}
    }
  }
  (*env)->ReleasePrimitiveArrayCritical(env, jB, B, 0);
  (*env)->ReleasePrimitiveArrayCritical(env, jA, A, 0);
}

ext: c

max_stars_repo_name: BIDData/BIDMat

max_issues_repo_name: BIDData/BIDMat

max_forks_repo_name: BIDData/BIDMat


text:
#pragma once

#include "string.h"
#include <math.h>
#include "flow.h"

extern "C" {
#include <cblas.h>
}

void OpenblasGemm(const CBLAS_TRANSPOSE TransA,
	const CBLAS_TRANSPOSE TransB, const int M, const int N, const int K,
	const float alpha, const float* A, const float* B, const float beta,
	float* C);

void OpenblasGemv(const CBLAS_TRANSPOSE TransA, const int M,
	const int N, const float alpha, const float* A, const float* x,
	const float beta, float* y);

void InitWeightBias(Flow &weight_, Flow &bias_);

//Y[i] = X[i]
void CopyData(const int N, const float* X, float* Y);

//y[i] = exp(a[i])
void ExpData(const int n, const float* a, float* y);

//y[i] = a[i]\b[i]
void DivData(const int n, const float* a, const float* b, float* y);

//X = alpha*X
void ScalData(const int N, const float alpha, float *X);

//Y=alpha*X+Y
void AxpyData(const int N, const float alpha, const float* X, float* Y);

//Y= alpha*X+beta*Y
void AxpbyData(const int N, const float alpha, const float* X, const float beta, float* Y);

ext: h

max_stars_repo_name: samylee/LazyNet


text:
#include "ccl_utils.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "ccl_params.h"
#include "ccl_error.h"
#include <gsl/gsl_errno.h>

/* ------- ROUTINE: ccl_linear spacing ------
INPUTS: [xmin,xmax] of the interval to be divided in N bins
OUTPUT: bin edges in range [xmin,xmax]
*/

double * ccl_linear_spacing(double xmin, double xmax, int N)
{
  double dx = (xmax-xmin)/(N -1.);

  double * x = malloc(sizeof(double)*N);
  if (x==NULL) {
    fprintf(stderr, "ERROR: Could not allocate memory for linear-spaced array (N=%d)\n", N);
    return x;
  }

  for (int i=0; i<N; i++) {
    x[i] = xmin + dx*i;
  }
  x[0]=xmin; //Make sure roundoff errors don't spoil edges
  x[N-1]=xmax; //Make sure roundoff errors don't spoil edges

  return x;
}

/* ------- ROUTINE: ccl_log spacing ------
INPUTS: [xmin,xmax] of the interval to be divided logarithmically in N bins
TASK: divide an interval in N logarithmic bins
OUTPUT: bin edges in range [xmin,xmax]
*/

double * ccl_log_spacing(double xmin, double xmax, int N)
{
  if (N<2) {
    fprintf(stderr, "ERROR: Cannot make log-spaced array with %d points - need at least 2\n", N);
    return NULL;
  }

  if (!(xmin>0 && xmax>0)) {
    fprintf(stderr, "ERROR: Cannot make log-spaced array xmax or xmax non-positive (had %le, %le)\n", xmin, xmax);
    return NULL;
  }

  double log_xmax = log(xmax);
  double log_xmin = log(xmin);
  double dlog_x = (log_xmax - log_xmin) /  (N-1.);

  double * x = malloc(sizeof(double)*N);
  if (x==NULL) {
    fprintf(stderr, "ERROR: Could not allocate memory for log-spaced array (N=%d)\n", N);
    return x;
  }

  for (int i=0; i<N; i++) {
    x[i] = exp(log_xmin + dlog_x*i);
  }
  x[0]=xmin; //Make sure roundoff errors don't spoil edges
  x[N-1]=xmax; //Make sure roundoff errors don't spoil edges

  return x;
}

//Spline creator
//n     -> number of points
//x     -> x-axis
//y     -> f(x)-axis
//y0,yf -> values of f(x) to use beyond the interpolation range
SplPar *ccl_spline_init(int n,double *x,double *y,double y0,double yf)
{
  SplPar *spl=malloc(sizeof(SplPar));
  if(spl==NULL)
    return NULL;

  spl->intacc=gsl_interp_accel_alloc();
  spl->spline=gsl_spline_alloc(gsl_interp_cspline,n);
  int parstatus=gsl_spline_init(spl->spline,x,y,n);
  if(parstatus) {
    gsl_interp_accel_free(spl->intacc);
    gsl_spline_free(spl->spline);
    return NULL;
  }

  spl->x0=x[0];
  spl->xf=x[n-1];
  spl->y0=y0;
  spl->yf=yf;

  return spl;
}

//Evaluates spline at x checking for bound errors
double ccl_spline_eval(double x,SplPar *spl)
{
  if(x<=spl->x0)
    return spl->y0;
  else if(x>=spl->xf)
    return spl->yf;
  else {
    double y;
    int stat=gsl_spline_eval_e(spl->spline,x,spl->intacc,&y);
    if (stat!=GSL_SUCCESS) {
      ccl_raise_exception(stat,"ccl_utils.c: ccl_splin_eval(): gsl error\n");
      return NAN;
    }
    return y;
  }
}

//Spline destructor
void ccl_spline_free(SplPar *spl)
{
  gsl_spline_free(spl->spline);
  gsl_interp_accel_free(spl->intacc);
  free(spl);
}


text:
/*
 * The MIT License (MIT)
 *
 * Copyright (c) 2014 Boone "Bea" Adkins
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */
/**
 * @class pomdp::MultinomialPolicy
 *
 * CompoundPolicy whose action is determined by sampling a <s>multinomial</s> categorical distribution across
 * sub-policies.
 *
 * @date Aug 10, 2013
 * @author Bea Adkins
 */

#ifndef MULTINOMIAL_POLICY_H
#define MULTINOMIAL_POLICY_H

#include <boost/static_assert.hpp>
#include <boost/type_traits.hpp>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include "assert_rng.h"
#include "compound_policy.h"

namespace pomdp
{
template <class Base_Policy>
class MultinomialPolicy : public CompoundPolicySpec<Base_Policy>
{
public:
  MultinomialPolicy() : val_table_(NULL)
  {
    rng_ = gsl_rng_alloc(gsl_rng_taus2);

    // Randomly seed the RNG.
    srand(time(NULL));
    gsl_rng_set (rng_, rand());
  }

  virtual ~MultinomialPolicy()
  {
    gsl_rng_free(rng_);
    if(CompoundPolicySpec<Base_Policy>::isInit())
    {
      gsl_ran_discrete_free(val_table_);
    }
  }

  virtual bool loadFromFiles(const std::vector<std::string>& paths)
  {
    if(!CompoundPolicySpec<Base_Policy>::loadFromFiles(paths))
      return false;
    CompoundPolicySpec<Base_Policy>::is_init_ = false; // Not yet initialized.

    // Now that model number is known, initialize random number variables.
    freshenDistribution();

    CompoundPolicySpec<Base_Policy>::is_init_ = true;
    return true;
  }

  virtual int applyPolicy(const FeatureValue& current) const
  {
    // Sample a policy index.
    int i_sampled_policy = sample();
    if(i_sampled_policy == -1)
      return -1;
    ROS_ASSERT(0 <= i_sampled_policy && i_sampled_policy < CompoundPolicySpec<Base_Policy>::getNumModels());

    // Apply that policy.
    ROS_DEBUG_STREAM("Appying policy #" << i_sampled_policy << ":");
    return CompoundPolicySpec<Base_Policy>::models_[i_sampled_policy]->applyPolicy(current);
  }

  virtual bool setModelLikelihoods(const std::vector<double>& model_likelihoods)
  {
    if(!CompoundPolicySpec<Base_Policy>::setModelLikelihoods(model_likelihoods))
      return false;

    freshenDistribution(); // With new model likelihoods
    return true;
  }
protected:
  /**
   * Random number generator.
   */
  gsl_rng* rng_;

  /**
   * Table for quickly sampling values.
   */
  gsl_ran_discrete_t* val_table_;

  /**
   * @return int within [0, numModels] or -1 on error.
   */
  int sample() const
  {
    if(!CompoundPolicySpec<Base_Policy>::isInit() || CompoundPolicySpec<Base_Policy>::getNumModels() == 0)
      return -1;

    // Sample a value.
    return gsl_ran_discrete(rng_, val_table_);
  }

  /**
   * Call this whenever model_likelihoods_ change, in value or length.
   */
  void freshenDistribution()
  {
    ROS_ASSERT(CompoundPolicySpec<Base_Policy>::model_likelihoods_.size() > 0);

    if(CompoundPolicySpec<Base_Policy>::isInit())
      gsl_ran_discrete_free(val_table_);
    val_table_ = gsl_ran_discrete_preproc(CompoundPolicySpec<Base_Policy>::model_likelihoods_.size(),
                                          &CompoundPolicySpec<Base_Policy>::model_likelihoods_[0]);
  }
};

typedef MultinomialPolicy<MDPPolicy> MultinomialMDPPolicy;
typedef MultinomialPolicy<POMDPPolicy> MultinomialPOMDPPolicy;

} /* namespace pomdp */
#endif /* MULTINOMIAL_POLICY_H */


text:
#ifndef CACHED_IMAGE_H
#define CACHED_IMAGE_H

#ifndef solaris
#  include <stdint.h>
#endif
#include "asf_meta.h"
#include "float_image.h"
#include <stdio.h>
#include <sys/types.h>

#include <glib.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_histogram.h>

//---------------------------------------------------------------------------
// Adding a new supported data type:
//  cache.c:
//    data_size()
//    cached_image_get_pixel()
//    cached_image_get_rgb()
//  stats.c:
//    generate_thumbnail_data()
//  big_image.c:
//    update_pixel_info()
//  read_X.c: (for clients that will handle the data type)
//    open_X_data()
//
typedef enum {
    UNDEFINED = 0,
    GREYSCALE_FLOAT = 1,
    GREYSCALE_BYTE = 2,
    RGB_BYTE = 3,
    RGB_FLOAT = 4
} ssv_data_type_t;

// Stats structure -- this one is used by both greyscale and RGB
typedef struct {
    double map_min, map_max;
    double avg, stddev;
    double act_min, act_max; // absolute min/max of all values
    int hist[256];           // histogram
    double no_data_value;    // value indicating "no data"
    int have_no_data;        // TRUE if "no_data_value" is present
    double no_data_min;
    double no_data_max;
    int have_no_data_range;
    int truncate;
} ImageStats;

// Stats structure -- this one is used only for RGB images, keeps track
//                    of stats for each channel
typedef struct {
    double map_min, map_max;
    double avg, stddev;
    double act_min, act_max; // absolute min/max of all values
    double no_data_value;
    int have_no_data;
    double no_data_max;
    double no_data_min;
    int have_no_data_range;
    int truncate;
} ImageStatsRGB;

// NOTE: At the moment, the Pixbuf that contains the displayed data
// is plain RGB, not RGB-A, so adding support for transparency would
// also involve changing big_image.c/make_big_image(), and you'd have
// to add a method to the cache interface that returned the alpha
// channel in addition to the rgb values.

//---------------------------------------------------------------------------
// This is the client interface -- the set of function pointers, etc,
// that encapsulate how the image cache gets data from the file.
// See: -read.c/read_file() for how the cache is hooked up to the
//          client depending on what the file type is
//      -read_template.c for an example of how to add another client
typedef int ReadClientFn(int row_start, int n_rows_to_get,
                         void *dest, void *read_client_info,
                         meta_parameters *meta, int data_type);
typedef int ThumbFn(int thumb_size_x, int thumb_size_y, meta_parameters *meta,
                    void *read_client_info, void *dest, int data_type);
typedef void FreeFn(void *read_client_info);

typedef struct {
    ReadClientFn *read_fn;
    ThumbFn *thumb_fn;
    FreeFn *free_fn;
    void *read_client_info;
    ssv_data_type_t data_type;
    int require_full_load;
} ClientInterface;

//---------------------------------------------------------------------------
// Here is the ImageCache stuff.  The global ImageCache that holds the
// loaded image is "data_ci".  This is all private data.
typedef struct {
  int nl, ns;               // Image dimensions.
  ClientInterface *client;  // pointers to data read implementations
  int n_tiles;              // Number of tiles in memory
  int reached_max_tiles;    // Have we loaded as many tiles as we can?
  int rows_per_tile;        // Number of rows in each tile
  int entire_image_fits;    // TRUE if we can load the entire image
  int *rowstarts;           // Row numbers starting each tile
  unsigned char **cache;    // Cached values (floats, unsigned chars ...)
  int *access_counts;       // Updated when a tile is accessed
  int n_access;             // used to find oldest tile
  ssv_data_type_t data_type;// type of data we have
  meta_parameters *meta;    // metadata -- don't own this pointer
  ImageStats *stats;        // not owned by us, not populated by us
  ImageStatsRGB *stats_r;   // not owned by us, not populated by us
  ImageStatsRGB *stats_g;   // not owned by us, not populated by us
  ImageStatsRGB *stats_b;   // not owned by us, not populated by us
} CachedImage;

CachedImage * cached_image_new_from_file(
    const char *file, meta_parameters *meta, ClientInterface *client,
    ImageStats *stats, ImageStatsRGB *stats_r, ImageStatsRGB *stats_g,
    ImageStatsRGB *stats_b);

float cached_image_get_pixel (CachedImage *self, int line, int samp);
void cached_image_get_rgb(CachedImage *self, int line, int samp,
                          unsigned char *r, unsigned char *g,
                          unsigned char *b);
void cached_image_get_rgb_float(CachedImage *self, int line, int samp,
                                float *r, float *g, float *b);

void load_thumbnail_data(CachedImage *self, int thumb_size_x, int thumb_size_y,
                         void *dest);

void cached_image_free (CachedImage *self);

#endif

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
#pragma once
#include <gsl/gsl>
#include <memory>
class BassEnhance;
class Eq_5band;
class Reverb;

struct Fx_Parameter {
    const char *name;
    int default_value;
    enum Type { Integer, Boolean };
    Type type;
    enum Flag { HasSeparator = 1 };
    unsigned flags;
};

class Synth_Fx {
public:
    Synth_Fx();
    ~Synth_Fx();

    void init(float sample_rate);
    void clear();
    void compute(float data[], unsigned nframes);

    int get_parameter(size_t index) const;
    void set_parameter(size_t index, int value);

    static gsl::span<const Fx_Parameter> parameters();

    enum {
        P_Bass_Enhance,
        P_Bass_Amount,
        P_Bass_Tone,
        P_Eq_Enable,
        P_Eq_Low,
        P_Eq_Mid_Low,
        P_Eq_Mid,
        P_Eq_Mid_High,
        P_Eq_High,
        P_Reverb_Enable,
        P_Reverb_Amount,
        P_Reverb_Size,
        Parameter_Count,
    };

private:
    std::unique_ptr<BassEnhance> be_;
    std::unique_ptr<Eq_5band> eq_;
    std::unique_ptr<Reverb> rev_;
    bool be_enable_ = false;
    bool eq_enable_ = false;
    bool rev_enable_ = false;
};

max_stars_repo_name: jpcima/smf-dsp

max_issues_repo_name: jpcima/smf-dsp

max_forks_repo_name: jpcima/smf-dsp


text:
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include <math.h>
typedef struct{
  double g;//gravity
  double M;//mass of the car
  double m;//mass of the load
  double r;//length of the rope
  double C;//damping coefficient
  double k;//spring constang

  double F;//input force
  double T;//tension to the rope
  //  double FF;//Fictitious force

  double x,x0;//x position of the load
  double y,y0;//y position of the load
  double X;//X position of the car
  double a;//angle (theta)

  double dx,dx0;
  double dy,dy0;
  double dX,dX0;
  double dr;
  double da,da0;
  double dXmax;
  double ddXmax;
  double Fmax;

  double ddx;
  double ddy;
  double ddX;
  double ddr;
  double dda;
  double h;
  int nh;
  double _h;
  int _dim;
  double *_y;
  double *_y_err;
  double *_dydt_in;
  double *_dydt_out;
  const gsl_odeiv_step_type *_T;
  gsl_odeiv_step *_s;
  gsl_odeiv_system _sys;
  double _t;
} CRANE;

#define dim_crane 4
int cranefunc (double t, const double y[], double f[], void *params)
{
  CRANE *c = (CRANE *)params;

#ifdef ICONIP2014 //wrong dynamics
  f[0] = y[1];
  f[1] = -c->k*(y[0]-y[2])/c->m - c->C*(y[1]-y[3])/c->m;//iconip2014
  f[2] = y[3];
  f[3] = (c->k*(y[0]-y[2]) + c->F)/c->M;//iconip2014
//  f[0] = y[1];
//  f[1] = -c->k*(y[0]-y[2])/c->m - c->C*(y[1]-y[3])/c->m;
//  f[2] = y[3];
//  f[3] = (c->k*(y[0]-y[2]) + c->C*(y[1]-y[3] + c->F)/c->M;
#else //corrected dynamics
  f[0] = y[1];
  f[1] = (-c->k*(y[0]-y[2]) - c->C*(y[1]-y[3])       )/c->m;//iconip2014
  f[2] = y[3];
  f[3] = ( c->k*(y[0]-y[2]) + c->C*(y[1]-y[3]) + c->F)/c->M;//iconip2014
#endif

//  f[0] = y[1];
//  f[1] = -c->k*y[0]/c->m - c->C*y[1]/c->m - c->ddX;
//  f[2] = y[3];
//  f[3] = (c->k*y[0] + c->F)/c->M;
//
//  f[0] = y[1];
//  f[1] = (-(2*c->dr+c->C)*y[1] -c->ddX*cos(y[0]) -c->g*sin(y[0]))/c->r;
//  f[2] = y[3];
//  f[3] = (c->F+c->T*sin(y[0]))/c->M;
  return GSL_SUCCESS;
}
////////////////
#define square(x) ((x)*(x))
CRANE crane;
int initialize()
{
  crane.g=9.8;
  crane.M=100;
  crane.m=20;//10;//20
  crane.r=5;//5;
  crane.C=10;//5;
  crane.k=15;//
  crane.dr=crane.ddr=0;
  crane.x0=0;
  crane.y0=crane.r;
  crane.dx0=crane.dy0=crane.da0=0;
  crane.dXmax=1.;//??ddXmax=(dXmax-dX)/dt;//??
  crane.ddXmax=0.2;//??ddXmax=(dXmax-dX)/dt;//??
  crane.Fmax=30.0;//check//  crane.Fmax=30;//check
#ifdef CRANESUB
  crane.h=AP_tS;
  crane.m=_crane_m;
  crane.r=_crane_r;
  crane.C=_crane_C;
  crane.k=_crane_k;//
  crane.dXmax=_crane_dXmax;
  if(_AP_umax>0) AP_u_max=crane.ddXmax=_AP_umax;
  else AP_u_max=crane.ddXmax;
  AP_u_min=-AP_u_max;
  //  starttime=0;
  //  totaltime=50;
  rr=AP_r=_AP_r;//10
  rr_kyoyou=_rr_kyoyou;
  //  p=(double *)malloc(buffsize*sizeof(double));
  C_MODE=11;
  //  iteration=_iteration;
#else
  crane.h=0.01;
  //  crane.h=0.0001;
#endif
  crane.nh=10;
  crane._h=crane.h/crane.nh;//0.001
  crane._dim=4;
  if(crane._y==NULL){
    crane._y=(double*)malloc(crane._dim*sizeof(double));
    crane._y_err=(double*)malloc(crane._dim*sizeof(double));
    crane._dydt_in=(double*)malloc(crane._dim*sizeof(double));
    crane._dydt_out=(double*)malloc(crane._dim*sizeof(double));
  }
  crane._T= gsl_odeiv_step_rk4;
  crane._s= gsl_odeiv_step_alloc (crane._T, crane._dim);
  crane._sys= (gsl_odeiv_system){cranefunc, NULL, crane._dim, &crane};
  //  crane._sys[0]=cranefunc;
  //  crane._sys[1]=NULL;
  //  crane._sys[2]=crane._dim;
  //  crane._sys[3]=&crane;
  //initialize variables
  crane._t=crane._y[0]=crane._y[1]=crane._y[2]=crane._y[3]=0;
  //  crane.FF=0;
  int i;for(i=0;i<crane._dim;i++) crane._dydt_in[i]=0;
  GSL_ODEIV_FN_EVAL(&(crane._sys), crane._t, crane._y, crane._dydt_in);
//  crane.a  =crane._y[0];
//  crane.da =crane._y[1];
//  crane.X  =crane._y[2];
//  crane.dX =crane._y[3];
//  crane.x  =crane.X+crane.r * sin(crane.a); //output y
//  crane.y  =crane.r*cos(crane.a);
//  crane.dXmax=1.;//??ddXmax=(dXmax-dX)/dt;//??
//  crane.ddXmax=0.2;//??ddXmax=(dXmax-dX)/dt;//??
  return(0);
}
#ifndef CRANESUB
char *fn="crane3io.dat";
FILE *fp;
#endif//#ifndef CRANESUB

//double plant(double uu)
double plant(double uu)
{
  if(uu>crane.ddXmax) uu=crane.ddXmax;
  else if(uu<-crane.ddXmax) uu=-crane.ddXmax;

  int n;
  for(n=0;n<crane.nh;n++){
    crane.ddX=uu;//crane.ddX=(uu<crane.ddXmax)?uu:crane.ddXmax;
    if(crane.dX>=crane.dXmax && crane.ddX>0) crane.ddX=0;
    else if(crane.dX<=-crane.dXmax && crane.ddX<0) crane.ddX=0;

    int status = gsl_odeiv_step_apply (crane._s, crane._t, crane._h,
				       crane._y, crane._y_err,
				       crane._dydt_in,
				       crane._dydt_out,
				       &crane._sys);
    if (status != GSL_SUCCESS) {fprintf(stderr,"### failure of plant calc.\n"); break;}

    int i;for(i=0;i<dim_crane;i++) crane._dydt_in[i]=crane._dydt_out[i];
    crane._t+=crane._h;
    crane.x  =crane._y[0];
    crane.dx =crane._y[1];
    crane.X  =crane._y[2];
    crane.dX =crane._y[3];
    //    crane.x  =crane.X+crane.a; //output y
    //    crane.dx=(crane.x-crane.x0)/crane._h;
    //crane.x  =crane.X+crane.r * sin(crane.a); //output y
    //crane.y  =crane.r*cos(crane.a);
    //    crane.dx=(crane.x-crane.x0)/crane._h;
    //    crane.dy=(crane.y-crane.y0)/crane._h;
    crane.ddx=(crane.dx-crane.dx0)/crane._h;
    //    crane.ddX=(crane.dX-crane.dX0)/crane._h;
    //    crane.ddy=(crane.dy-crane.dy0)/crane._h;
    //    crane.T=crane.m*sqrt(square(crane.ddx)+square(crane.ddy-crane.g));
    crane.F=crane.M*crane.ddX-crane.k*(crane.x-crane.X)-crane.C*(crane.dx-crane.dX);;
    //    crane.F=crane.M*crane.ddX-crane.k*(crane.x-crane.X);
    //    crane.F=crane.M*crane.ddX-crane.T*sin(crane.a);
    //    crane.dda=(crane.da-crane.da0)/crane._h;
    crane.x0=crane.x;//?
    //    crane.y0=crane.y;
    crane.dx0=crane.dx;//?
    //    crane.dy0=crane.dy;
    //    crane.da0=crane.da;
    //    crane.FF=crane.ddX;
    crane.dX0=crane.dX;
  }
#ifndef CRANESUB
  fprintf(fp,"%.7e %.7e %.7e", crane._t,crane._y[0],crane._y[1]);//crane.x,crane.dx
  fprintf(fp," %.7e %.7e %.7e" ,crane.X,crane.dX,crane.ddX);
  fprintf(fp," %.7e %.7e %.7e" ,crane.x,crane.dx,crane.ddx);
  fprintf(fp," %.7e %.7e %.7e" ,crane.y,crane.dy,crane.ddy);
  fprintf(fp," %.7e %.7e" ,crane.T,crane.F);
  fprintf(fp,"\n");
#endif //#ifndef CRANESUB
  return(crane.x);//crane.x in [0,AP_y=10?] ==> rt in [0,1]
}
#ifndef CRANESUB
int main(int argc,char **argv)
{
  ////////////////////////////////////////////////////
  /// method 1                                     ///
  /// input ddX                                    ///
  /// output a,x,y,F                               ///
  ////////////////////////////////////////////////////
  //  double h = 0.001;//,hh=0.1; h = 0.0001;//,hh=0.1;
  double T=5;
  double t0=0;     //stationary
  double t1=T+t0;//accelerate(speed up)
  double t2=T+t1;//free run
  double t3=T+t2;//decelerate(slow down)
  double t4=t0+50;   //free run

  initialize();
  int M=(t4/crane.h)+1;
  double *ddX=(double*)malloc(sizeof(double)*M);
  //  double *F  =(double*)malloc(sizeof(double)*M);
  double Vmax=1;//dX/dt=1[m/s]
  //C:10;K:15;M:100;m:10;TS:0.5;/*数式処理の場合、この行を実行しない.m=10,40,70,100にして行う*/
  /*  crane.M=100;
  crane.C=10;
  crane.k=15;
  crane.m=10;*/
  int i;
  for(i=1;i<argc;i++){
    if(strncmp(argv[i],"m:",2)==0){
      sscanf(&argv[i][2],"%lf",&crane.m);
      //      fprintf(stderr,"#m=%g",crane.m);
    }
  }
  fprintf(stderr,"#C:%g K:%g M:%g m:%g\n",crane.C,crane.k,crane.M,crane.m);

  double TS=0.5;
  double umax=10;
  int n;
  double t;
  ////////////////////////////////////////////////////
  ///   set ddX                                    ///
  ////////////////////////////////////////////////////
  int n0 =t0/crane.h+0.5, n4=t4/crane.h+0.5;
  for(t=0;t<t4;t+=crane.h){
    double u;
    n=t/crane.h;
    if(t<=t0) u=0; //zero
    else if(t<=t1) u=umax*0.8; //accelerate
    else if(t<=t2) u=0;            //free move
    else if(t<=t3) u=-umax*0.8;//decelerate
    else u=0;                      //free move
    ddX[n]=u;//??
    //    ddX[n]=u*crane.h;//??
  }
  ////////////////////////////////////////////////////
  ///   Solve by the Runge Kutta Method of GSL     ///
  ////////////////////////////////////////////////////
  fp=fopen(fn,"w");
  fprintf(fp,"#%.7e %d %d %.7e %.7e %.7e %.7e #h,n0,n4,M,n,r\n",crane.h,n0,n4,crane.M,crane.m,crane.r,crane.C);
  //  double dx0=0,dy0=0,da0=0;
  //  for(n=0;n<n4;){
  for(t=0;t<t4;t+=crane.h){
    n=t/crane.h;
    plant(ddX[n]);
  }
  fclose(fp);
  fprintf(stdout,"Results are stored in '%s'.\n",fn);

  gsl_odeiv_step_free (crane._s);
  return 0;
}
#endif //#ifndef SUB


text:
/**
 * \file FFT.h
 */

#ifndef ATK_UTILITY_FFT_H
#define ATK_UTILITY_FFT_H

#include <complex>
#include <cstdint>
#include <memory>
#include <vector>

#include <gsl/gsl>

#include <ATK/Utility/config.h>

namespace ATK
{
  /// An FFT class
  template<class DataType_>
  class ATK_UTILITY_EXPORT FFT
  {
  public:
    /// Builds the wrapper class
    FFT();
    /// Destructor
    ~FFT();

    FFT(const FFT&) = delete;
    FFT& operator=(const FFT&) = delete;

    /*!
     * @brief Sets a new size for the FFT and creates an associated plan for it
     * @ size is the new size of the plan
     */
    void set_size(gsl::index size);
    gsl::index get_size() const;
    /*!
     * @brief Processes a FFT from a real input array and stores it internally
     * @param input is the real input array of data to process
     * @param input_size is the size of the input array
     */
    void process(const DataType_* input, gsl::index input_size) const;
    /*!
    * @brief Processes a FFT from a real input to a complex output
    * @param input is the real input array of data to process
    * @param output is the complex output array
    * @param input_size is the size of both input and output arrays
    */
    void process_forward(const DataType_* input, std::complex<DataType_>* output, gsl::index input_size) const;
    /*!
    * @brief Processes an inverse FFT from a complex input to a real output
    * @param input is the complex input array of data to process
    * @param output is the real output array
    * @param input_size is the size of both input and output arrays
    */
    void process_backward(const std::complex<DataType_>* input, DataType_* output, gsl::index input_size) const;
    /*!
    * @brief Processes a FFT from a complex input array and stores it internally
    * @param input is the real input array of data to process
    * @param input_size is the size of the input array
    */
    void process(const std::complex<DataType_>* input, gsl::index input_size) const;
    /*!
    * @brief Processes a FFT from a complex input to a complex output
    * @param input is the real input array of data to process
    * @param output is the complex output array
    * @param input_size is the size of both input and output arrays
    */
    void process_forward(const std::complex<DataType_>* input, std::complex<DataType_>* output, gsl::index input_size) const;
    /*!
    * @brief Processes an inverse FFT from a complex input to a complex output
    * @param input is the complex input array of data to process
    * @param output is the real output array
    * @param input_size is the size of both input and output arrays
    */
    void process_backward(const std::complex<DataType_>* input, std::complex<DataType_>* output, gsl::index input_size) const;
    /*!
    * @brief Computes the amplitude of the resulting spectrum
    * @param amp is the output angle container
    */
    void get_amp(std::vector<DataType_>& amp) const;
    /*!
    * @brief Computes the angle of the resulting spectrum
    * @param angle is the output angle container
    */
    void get_angle(std::vector<DataType_>& angle) const;

  protected:
	  gsl::index size;

    class FFTImpl;
    std::unique_ptr<FFTImpl> impl;
  };
}

#endif

max_stars_repo_name: D-J-Roberts/AudioTK

max_issues_repo_name: D-J-Roberts/AudioTK

max_forks_repo_name: D-J-Roberts/AudioTK


text:
/*
 * Copyright (c) 2011-2018 The University of Tennessee and The University
 *                         of Tennessee Research Foundation.  All rights
 *                         reserved.
 *
 * @precisions normal z -> s d c
 *
 */
#include <math.h>
#include <stdlib.h>
#include <core_blas.h>
#include <cblas.h>
#include "dplasma.h"

/*
#define DEBUG_BUTTERFLY
*/

/* Forward declaration of kernels for the butterfly transformation */
void BFT_zQTL( int mb, int nb, int lda, int i_seg, int j_seg, int lvl, int N,
          PLASMA_Complex64_t *tl, PLASMA_Complex64_t *bl,
          PLASMA_Complex64_t *tr, PLASMA_Complex64_t *br,
          PLASMA_Complex64_t *C, PLASMA_Complex64_t *U_before, PLASMA_Complex64_t *U_after, int bl_is_tr_trans, int is_diagonal );

void BFT_zQBL( int mb, int nb, int lda, int i_seg, int j_seg, int lvl, int N,
          PLASMA_Complex64_t *tl, PLASMA_Complex64_t *bl,
          PLASMA_Complex64_t *tr, PLASMA_Complex64_t *br,
          PLASMA_Complex64_t *C, PLASMA_Complex64_t *U_before, PLASMA_Complex64_t *U_after, int bl_is_tr_trans, int is_diagonal );
void BFT_zQTR_trans( int mb, int nb, int lda, int i_seg, int j_seg, int lvl, int N,
          PLASMA_Complex64_t *tl, PLASMA_Complex64_t *bl,
          PLASMA_Complex64_t *tr, PLASMA_Complex64_t *br,
          PLASMA_Complex64_t *C, PLASMA_Complex64_t *U_before, PLASMA_Complex64_t *U_after, int bl_is_tr_trans );
void BFT_zQTR( int mb, int nb, int lda, int i_seg, int j_seg, int lvl, int N,
          PLASMA_Complex64_t *tl, PLASMA_Complex64_t *bl,
          PLASMA_Complex64_t *tr, PLASMA_Complex64_t *br,
          PLASMA_Complex64_t *C, PLASMA_Complex64_t *U_before, PLASMA_Complex64_t *U_after, int bl_is_tr_trans);
void BFT_zQBR( int mb, int nb, int lda, int i_seg, int j_seg, int lvl, int N,
          PLASMA_Complex64_t *tl, PLASMA_Complex64_t *bl,
          PLASMA_Complex64_t *tr, PLASMA_Complex64_t *br,
          PLASMA_Complex64_t *C, PLASMA_Complex64_t *U_before, PLASMA_Complex64_t *U_after, int bl_is_tr_trans, int is_diagonal );

void RBMM_zTOP( int mb, int nb, int lda, int i_seg, int lvl, int N, int trans,
          PLASMA_Complex64_t *top, PLASMA_Complex64_t *btm,
          PLASMA_Complex64_t *C, PLASMA_Complex64_t *U_but_vec);
void RBMM_zBTM( int mb, int nb, int lda, int i_seg, int lvl, int N, int trans,
          PLASMA_Complex64_t *top, PLASMA_Complex64_t *btm,
          PLASMA_Complex64_t *C, PLASMA_Complex64_t *U_but_vec);

/* Bodies of kernels for the butterfly transformation */

void BFT_zQTL( int mb, int nb, int lda, int i_seg, int j_seg, int lvl, int N,
          PLASMA_Complex64_t *tl, PLASMA_Complex64_t *bl,
          PLASMA_Complex64_t *tr, PLASMA_Complex64_t *br,
          PLASMA_Complex64_t *C, PLASMA_Complex64_t *U_before, PLASMA_Complex64_t *U_after, int bl_is_tr_trans, int is_diagonal )
{
    int i, j;

    (void)lvl;
    (void)N;

#if defined(DEBUG_BUTTERFLY)
    printf ("BFT_zQTL(mb:%d, nb:%d, lda:%d, i_seg:%d, j_seg:%d, lvl:%d, N:%d, bl_is_tr_trans:%d, is_diagonal:%d)\n",
            mb, nb, lda, i_seg, j_seg, lvl, N, bl_is_tr_trans, is_diagonal);
#endif

    if( bl_is_tr_trans ){
        for (j=0; j<nb; j++) {
            int start = is_diagonal ? j : 0;
            for (i=start; i<mb; i++) {
                PLASMA_Complex64_t ri = U_before[i_seg+i];
                PLASMA_Complex64_t rj = U_after[j_seg+j];
#if defined(DEBUG_BUTTERFLY)
                printf ("HE A[%d][%d] = U_before[%d]*((tl[%d]+bl[%d]) + (tr[%d]+br[%d])) * U_after[%d]\n", i_seg+i, j_seg+j, i_seg+i, j*lda+i, j*lda+i, i*lda+j, j*lda+i, j_seg+j);
#endif
                C[j*lda+i] = ri * ((tl[j*lda+i]+bl[j*lda+i]) + (tr[i*lda+j]+br[j*lda+i])) * rj;
#if defined(DEBUG_BUTTERFLY)
                printf ("HE %lf %lf %lf %lf %lf %lf %lf\n", creal(C[j*lda+i]), creal(ri), creal(tl[j*lda+i]), creal(bl[j*lda+i]), creal(tr[i*lda+j]), creal(br[j*lda+i]), creal(rj));
#endif
            }
        }
    }else{
        for (j=0; j<nb; j++) {
            int start = is_diagonal ? j : 0;
            for (i=start; i<mb; i++) {
                PLASMA_Complex64_t ri = U_before[i_seg+i];
                PLASMA_Complex64_t rj = U_after[j_seg+j];
#if defined(DEBUG_BUTTERFLY)
                printf ("GE A[%d][%d] = U_before[%d]*((tl[%d]+bl[%d]) + (tr[%d]+br[%d])) * U_after[%d]\n", i_seg+i, j_seg+j, i_seg+i, j*lda+i, j*lda+i, j*lda+i, j*lda+i, j_seg+j);
#endif
                C[j*lda+i] = ri * ((tl[j*lda+i]+bl[j*lda+i]) + (tr[j*lda+i]+br[j*lda+i])) * rj;
#if defined(DEBUG_BUTTERFLY)
                printf ("GE %lf %lf %lf %lf %lf %lf %lf\n", creal(C[j*lda+i]), creal(ri), creal(tl[j*lda+i]), creal(bl[j*lda+i]), creal(tr[j*lda+i]), creal(br[j*lda+i]), creal(rj));
#endif
            }
        }
    }
    return;
}

void BFT_zQBL( int mb, int nb, int lda, int i_seg, int j_seg, int lvl, int N,
          PLASMA_Complex64_t *tl, PLASMA_Complex64_t *bl,
          PLASMA_Complex64_t *tr, PLASMA_Complex64_t *br,
          PLASMA_Complex64_t *C, PLASMA_Complex64_t *U_before, PLASMA_Complex64_t *U_after, int bl_is_tr_trans, int is_diagonal )
{
    int i, j;
    int r_to_s = N/(1<<(lvl+1));

#if defined(DEBUG_BUTTERFLY)
    if( is_diagonal )
        printf ("BFT_zQBL(diag)\n");
    else
        printf ("BFT_zQBL(lower)\n");
#endif

    if( bl_is_tr_trans ){
        for (j=0; j<nb; j++) {
            for (i=0; i<mb; i++) {
                PLASMA_Complex64_t si = U_before[i_seg+r_to_s+i];
                PLASMA_Complex64_t rj = U_after[j_seg+j];
#if defined(DEBUG_BUTTERFLY)
                printf ("HE A[%d][%d] = U_before[%d]*((tl[%d]-bl[%d]) + (tr[%d]-br[%d])) * U_after[%d]\n", i_seg+i+r_to_s, j_seg+j, i_seg+r_to_s+i, j*lda+i, j*lda+i, i*lda+j, j*lda+i, j_seg+j);
#endif
                if( is_diagonal && j>i ){ /* If the tile of the TL quarter is on the diagonal, then we need to take the transpose element for the tiles that came from TL and BR */
                    C[j*lda+i] = si * ((tl[i*lda+j]-bl[j*lda+i]) + (tr[i*lda+j]-br[i*lda+j])) * rj;
                } else{
                    C[j*lda+i] = si * ((tl[j*lda+i]-bl[j*lda+i]) + (tr[i*lda+j]-br[j*lda+i])) * rj;
                }
#if defined(DEBUG_BUTTERFLY)
                printf ("HE %lf %lf %lf %lf %lf %lf %lf\n",creal(C[j*lda+i]), creal(si), creal(tl[j*lda+i]), creal(bl[j*lda+i]), creal(tr[i*lda+j]), creal(br[j*lda+i]), creal(rj));
#endif
            }
        }
    }else{
        for (j=0; j<nb; j++) {
            for (i=0; i<mb; i++) {
                PLASMA_Complex64_t si = U_before[i_seg+r_to_s+i];
                PLASMA_Complex64_t rj = U_after[j_seg+j];
#if defined(DEBUG_BUTTERFLY)
                printf ("GE A[%d][%d] = U_before[%d]*((tl[%d]-bl[%d]) + (tr[%d]-br[%d])) * U_after[%d]\n", i_seg+i+r_to_s, j_seg+j, i_seg+r_to_s+i, j*lda+i, j*lda+i, j*lda+i, j*lda+i, j_seg+j);
#endif
                C[j*lda+i] = si * ((tl[j*lda+i]-bl[j*lda+i]) + (tr[j*lda+i]-br[j*lda+i])) * rj;
#if defined(DEBUG_BUTTERFLY)
                printf ("GE %lf %lf %lf %lf %lf %lf %lf\n", creal(C[j*lda+i]), creal(si), creal(tl[j*lda+i]), creal(bl[j*lda+i]), creal(tr[j*lda+i]), creal(br[j*lda+i]), creal(rj));
#endif
            }
        }
    }
    return;
}

/* This function writes into a transposed tile, so C is always transposed. */
void BFT_zQTR_trans( int mb, int nb, int lda, int i_seg, int j_seg, int lvl, int N,
          PLASMA_Complex64_t *tl, PLASMA_Complex64_t *bl,
          PLASMA_Complex64_t *tr, PLASMA_Complex64_t *br,
          PLASMA_Complex64_t *C, PLASMA_Complex64_t *U_before, PLASMA_Complex64_t *U_after, int bl_is_tr_trans )
{
    int i,j;
    int r_to_s = N/(1<<(lvl+1));
#if defined(DEBUG_BUTTERFLY)
    printf ("BFT_zQTR_trans()\n");
#endif
    if( bl_is_tr_trans ){
        for (j=0; j<nb; j++) {
            for (i=0; i<mb; i++) {
                PLASMA_Complex64_t ri = U_before[i_seg+i];
                PLASMA_Complex64_t sj = U_after[j_seg+r_to_s+j];
#if defined(DEBUG_BUTTERFLY)
                printf ("HE A[%d][%d] = U_before[%d]*((tl[%d]+bl[%d]) - (tr[%d]+br[%d])) * U_after[%d]\n", i_seg+j, j_seg+i+r_to_s, i_seg+i, j*lda+i, j*lda+i, i*lda+j, j*lda+i, j_seg+r_to_s+j);
#endif
                C[i*lda+j] = ri * ((tl[j*lda+i]+bl[j*lda+i]) - (tr[i*lda+j]+br[j*lda+i])) * sj;
#if defined(DEBUG_BUTTERFLY)
                printf ("HE %lf %lf %lf %lf %lf %lf %lf\n", creal(C[i*lda+j]), creal(ri), creal(tl[j*lda+i]), creal(bl[j*lda+i]), creal(tr[i*lda+j]), creal(br[j*lda+i]), creal(sj));
#endif
            }
        }
    }else{
        assert(0); // This should never happen.
    }
    return;
}

void BFT_zQTR( int mb, int nb, int lda, int i_seg, int j_seg, int lvl, int N,
          PLASMA_Complex64_t *tl, PLASMA_Complex64_t *bl,
          PLASMA_Complex64_t *tr, PLASMA_Complex64_t *br,
          PLASMA_Complex64_t *C, PLASMA_Complex64_t *U_before, PLASMA_Complex64_t *U_after, int bl_is_tr_trans )
{
    int i, j;
    int r_to_s = N/(1<<(lvl+1));
#if defined(DEBUG_BUTTERFLY)
    printf ("BFT_zQTR()\n");
#endif
    if( bl_is_tr_trans ){
        for (j=0; j<nb; j++) {
            for (i=0; i<mb; i++) {
                PLASMA_Complex64_t ri = U_before[i_seg+i];
                PLASMA_Complex64_t sj = U_after[j_seg+r_to_s+j];
#if defined(DEBUG_BUTTERFLY)
                printf ("HE A[%d][%d] = U_before[%d]*((tl[%d]+bl[%d]) - (tr[%d]+br[%d])) * U_after[%d]\n", i_seg+i, j_seg+j+r_to_s, i_seg+i, j*lda+i, j*lda+i, i*lda+j, j*lda+i, j_seg+r_to_s+j);
#endif
                C[j*lda+i] = ri * ((tl[j*lda+i]+bl[j*lda+i]) - (tr[i*lda+j]+br[j*lda+i])) * sj;
#if defined(DEBUG_BUTTERFLY)
                printf ("HE %lf %lf %lf %lf %lf %lf %lf\n", creal(C[j*lda+i]), creal(ri), creal(tl[j*lda+i]), creal(bl[j*lda+i]), creal(tr[i*lda+j]), creal(br[j*lda+i]), creal(sj));
#endif
            }
        }
    }else{
        for (j=0; j<nb; j++) {
            for (i=0; i<mb; i++) {
                PLASMA_Complex64_t ri = U_before[i_seg+i];
                PLASMA_Complex64_t sj = U_after[j_seg+r_to_s+j];
#if defined(DEBUG_BUTTERFLY)
                printf ("GE A[%d][%d] = U_before[%d]*((tl[%d]+bl[%d]) - (tr[%d]+br[%d])) * U_after[%d]\n", i_seg+i, j_seg+j+r_to_s, i_seg+i, j*lda+i, j*lda+i, j*lda+i, j*lda+i, j_seg+r_to_s+j);
#endif
                C[j*lda+i] = ri * ((tl[j*lda+i]+bl[j*lda+i]) - (tr[j*lda+i]+br[j*lda+i])) * sj;
#if defined(DEBUG_BUTTERFLY)
                printf ("GE %lf %lf %lf %lf %lf %lf %lf\n", creal(C[j*lda+i]), creal(ri), creal(tl[j*lda+i]), creal(bl[j*lda+i]), creal(tr[j*lda+i]), creal(br[j*lda+i]), creal(sj));
#endif
            }
        }
    }
    return;
}

void BFT_zQBR( int mb, int nb, int lda, int i_seg, int j_seg, int lvl, int N,
          PLASMA_Complex64_t *tl, PLASMA_Complex64_t *bl,
          PLASMA_Complex64_t *tr, PLASMA_Complex64_t *br,
          PLASMA_Complex64_t *C, PLASMA_Complex64_t *U_before, PLASMA_Complex64_t *U_after, int bl_is_tr_trans, int is_diagonal )
{
    int i, j;
    int r_to_s = N/(1<<(lvl+1));

#if defined(DEBUG_BUTTERFLY)
    if( is_diagonal )
        printf ("BFT_zQBR(diag)\n");
    else
        printf ("BFT_zQBR(lower)\n");
#endif

    if( bl_is_tr_trans ){
        for (j=0; j<nb; j++) {
            int start = is_diagonal ? j : 0;
            for (i=start; i<mb; i++) {
                PLASMA_Complex64_t si = U_before[i_seg+r_to_s+i];
                PLASMA_Complex64_t sj = U_after[j_seg+r_to_s+j];
#if defined(DEBUG_BUTTERFLY)
                printf ("HE A[%d][%d] = U_before[%d]*((tl[%d]-bl[%d]) - (tr[%d]-br[%d])) * U_after[%d]\n", i_seg+i+r_to_s, j_seg+j+r_to_s, i_seg+r_to_s+i, j*lda+i, j*lda+i, i*lda+j, j*lda+i, j_seg+r_to_s+j);
#endif
                C[j*lda+i] = si * ((tl[j*lda+i]-bl[j*lda+i]) - (tr[i*lda+j]-br[j*lda+i])) * sj;
#if defined(DEBUG_BUTTERFLY)
                printf ("HE %lf %lf %lf %lf %lf %lf %lf\n", creal(C[j*lda+i]), creal(si), creal(tl[j*lda+i]), creal(bl[j*lda+i]), creal(tr[i*lda+j]), creal(br[j*lda+i]), creal(sj));
#endif
            }
        }
    }else{
        for (j=0; j<nb; j++) {
            int start = is_diagonal ? j : 0;
            for (i=start; i<mb; i++) {
                PLASMA_Complex64_t si = U_before[i_seg+r_to_s+i];
                PLASMA_Complex64_t sj = U_after[j_seg+r_to_s+j];
#if defined(DEBUG_BUTTERFLY)
                printf ("GE A[%d][%d] = U_before[%d]*((tl[%d]-bl[%d]) - (tr[%d]-br[%d])) * U_after[%d]\n", i_seg+i+r_to_s, j_seg+j+r_to_s, i_seg+r_to_s+i, j*lda+i, j*lda+i, j*lda+i, j*lda+i, j_seg+r_to_s+j);
#endif
                C[j*lda+i] = si * ((tl[j*lda+i]-bl[j*lda+i]) - (tr[j*lda+i]-br[j*lda+i])) * sj;
#if defined(DEBUG_BUTTERFLY)
                printf ("GE %lf %lf %lf %lf %lf %lf %lf\n", creal(C[j*lda+i]), creal(si), creal(tl[j*lda+i]), creal(bl[j*lda+i]), creal(tr[j*lda+i]), creal(br[j*lda+i]), creal(sj));
#endif
            }
        }
    }
    return;
}

void RBMM_zTOP( int mb, int nb, int lda, int i_seg, int lvl, int N, int trans,
          PLASMA_Complex64_t *top, PLASMA_Complex64_t *btm,
          PLASMA_Complex64_t *C, PLASMA_Complex64_t *U_but_vec)
{
    int i, j;
    int r_to_s = N/(1<<(lvl+1));
#if defined(DEBUG_BUTTERFLY)
    printf ("RBMM_zTOP()\n");
#endif
    for (j=0; j<nb; j++) {
        for (i=0; i<mb; i++) {
            PLASMA_Complex64_t r = U_but_vec[i_seg+i];
            if( PlasmaConjTrans == trans ){
#if defined(DEBUG_BUTTERFLY)
                printf ("C[%d] = U_but_vec[%d]*(top[%d]+btm[%d])\n", j*lda+i, i_seg+i, j*lda+i, j*lda+i);
#endif
                C[j*lda+i] = r * (top[j*lda+i] + btm[j*lda+i]);
#if defined(DEBUG_BUTTERFLY)
                printf ("%lf %lf %lf %lf\n",creal(C[j*lda+i]), creal(r), creal(top[j*lda+i]), creal(btm[j*lda+i]));
#endif
            }else{
                PLASMA_Complex64_t s = U_but_vec[i_seg+r_to_s+i];
#if defined(DEBUG_BUTTERFLY)
                printf ("C[%d] = U_but_vec[%d]*top[%d] + U_but_vec[%d]*btm[%d]\n", j*lda+i, i_seg+i, j*lda+i, i_seg+r_to_s+i, j*lda+i);
#endif
                C[j*lda+i] =  r*top[j*lda+i] + s*btm[j*lda+i];
#if defined(DEBUG_BUTTERFLY)
                printf ("%lf %lf %lf %lf %lf\n",creal(C[j*lda+i]), creal(r), creal(top[j*lda+i]), creal(s), creal(btm[j*lda+i]));
#endif
            }
        }
    }
    return;
}

void RBMM_zBTM( int mb, int nb, int lda, int i_seg, int lvl, int N, int trans,
          PLASMA_Complex64_t *top, PLASMA_Complex64_t *btm,
          PLASMA_Complex64_t *C, PLASMA_Complex64_t *U_but_vec)
{
    int i, j;
    int r_to_s = N/(1<<(lvl+1));
#if defined(DEBUG_BUTTERFLY)
    printf ("RBMM_zBTM()\n");
#endif
    for (j=0; j<nb; j++) {
        for (i=0; i<mb; i++) {
            PLASMA_Complex64_t s = U_but_vec[i_seg+r_to_s+i];
            if( PlasmaConjTrans == trans ){
#if defined(DEBUG_BUTTERFLY)
                printf ("C[%d] = U_but_vec[%d]*(top[%d]-btm[%d])\n", j*lda+i, i_seg+r_to_s+i, j*lda+i, j*lda+i);
#endif
                C[j*lda+i] = s * (top[j*lda+i] - btm[j*lda+i]);
#if defined(DEBUG_BUTTERFLY)
                printf ("%lf %lf %lf %lf\n",creal(C[j*lda+i]), creal(s), creal(top[j*lda+i]), creal(btm[j*lda+i]));
#endif
            }else{
                PLASMA_Complex64_t r = U_but_vec[i_seg+i];
#if defined(DEBUG_BUTTERFLY)
                printf ("C[%d] = U_but_vec[%d]*top[%d] - U_but_vec[%d]*btm[%d]\n", j*lda+i, i_seg+i, j*lda+i, i_seg+r_to_s+i, j*lda+i);
#endif
                C[j*lda+i] = r*top[j*lda+i] - s*btm[j*lda+i];
#if defined(DEBUG_BUTTERFLY)
                printf ("%lf %lf %lf %lf %lf\n", creal(C[j*lda+i]), creal(r), creal(top[j*lda+i]), creal(s), creal(btm[j*lda+i]));
#endif
            }
        }
    }
    return;
}

max_stars_repo_name: NLAFET/ABFT

max_issues_repo_name: NLAFET/ABFT

max_forks_repo_name: NLAFET/ABFT


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics.h>

/* number of data points to fit */
#define N        200

/* number of fit coefficients */
#define NCOEFFS  12

/* nbreak = ncoeffs + 2 - k = ncoeffs - 2 since k = 4 */
#define NBREAK   (NCOEFFS - 2)

int
main (void)
{
  const size_t n = N;
  const size_t ncoeffs = NCOEFFS;
  const size_t nbreak = NBREAK;
  size_t i, j;
  gsl_bspline_workspace *bw;
  gsl_vector *B;
  double dy;
  gsl_rng *r;
  gsl_vector *c, *w;
  gsl_vector *x, *y;
  gsl_matrix *X, *cov;
  gsl_multifit_linear_workspace *mw;
  double chisq, Rsq, dof, tss;

  gsl_rng_env_setup();
  r = gsl_rng_alloc(gsl_rng_default);

  /* allocate a cubic bspline workspace (k = 4) */
  bw = gsl_bspline_alloc(4, nbreak);
  B = gsl_vector_alloc(ncoeffs);

  x = gsl_vector_alloc(n);
  y = gsl_vector_alloc(n);
  X = gsl_matrix_alloc(n, ncoeffs);
  c = gsl_vector_alloc(ncoeffs);
  w = gsl_vector_alloc(n);
  cov = gsl_matrix_alloc(ncoeffs, ncoeffs);
  mw = gsl_multifit_linear_alloc(n, ncoeffs);

  /* this is the data to be fitted */
  for (i = 0; i < n; ++i)
    {
      double sigma;
      double xi = (15.0 / (N - 1)) * i;
      double yi = cos(xi) * exp(-0.1 * xi);

      sigma = 0.1 * yi;
      dy = gsl_ran_gaussian(r, sigma);
      yi += dy;

      gsl_vector_set(x, i, xi);
      gsl_vector_set(y, i, yi);
      gsl_vector_set(w, i, 1.0 / (sigma * sigma));

      printf("%f %f\n", xi, yi);
    }

  /* use uniform breakpoints on [0, 15] */
  gsl_bspline_knots_uniform(0.0, 15.0, bw);

  /* construct the fit matrix X */
  for (i = 0; i < n; ++i)
    {
      double xi = gsl_vector_get(x, i);

      /* compute B_j(xi) for all j */
      gsl_bspline_eval(xi, B, bw);

      /* fill in row i of X */
      for (j = 0; j < ncoeffs; ++j)
        {
          double Bj = gsl_vector_get(B, j);
          gsl_matrix_set(X, i, j, Bj);
        }
    }

  /* do the fit */
  gsl_multifit_wlinear(X, w, y, c, cov, &chisq, mw);

  dof = n - ncoeffs;
  tss = gsl_stats_wtss(w->data, 1, y->data, 1, y->size);
  Rsq = 1.0 - chisq / tss;

  fprintf(stderr, "chisq/dof = %e, Rsq = %f\n",
                   chisq / dof, Rsq);

  printf("\n\n");

  /* output the smoothed curve */
  {
    double xi, yi, yerr;

    for (xi = 0.0; xi < 15.0; xi += 0.1)
      {
        gsl_bspline_eval(xi, B, bw);
        gsl_multifit_linear_est(B, c, cov, &yi, &yerr);
        printf("%f %f\n", xi, yi);
      }
  }

  gsl_rng_free(r);
  gsl_bspline_free(bw);
  gsl_vector_free(B);
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_matrix_free(X);
  gsl_vector_free(c);
  gsl_vector_free(w);
  gsl_matrix_free(cov);
  gsl_multifit_linear_free(mw);

  return 0;
} /* main() */

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE