text:
/*
 * Copyright (c) 2010-2020 The University of Tennessee and The University
 *                         of Tennessee Research Foundation.  All rights
 *                         reserved.
 *
 * @precisions normal z -> s d c
 *
 */

#include "dplasma.h"
#include "dplasma/types.h"
#include "dplasmaaux.h"
#include "cores/core_blas.h"
#include "parsec/data_dist/matrix/matrix.h"
#include "parsec/private_mempool.h"

#include <math.h>
#include <stdlib.h>
#include <cblas.h>
#include <string.h>

#include "butterfly_map.h"
#include "zhebut.h"
#include "zgebut.h"
#include "zgebmm.h"
#include <lapacke.h>

#if (PARSEC_zhebut_ADT_IDX_MAX != 1) || (PARSEC_zgebut_ADT_IDX_MAX != 1) || (PARSEC_zgebmm_ADT_IDX_MAX != 1)
#warning Current zhebut can work only if not using named types.
#endif

#define CREATE_N_ENQUEUE 0x0
#define DESTRUCT         0x1

static uint32_t parsec_rbt_rank_of(parsec_data_collection_t *desc, ...){
    int m_seg, n_seg, m_tile, n_tile;
    uintptr_t offset;
    va_list ap;
    parsec_seg_dc_t *segA;
    parsec_data_collection_t *A;

    va_start(ap, desc);
    m_seg = va_arg(ap, int);
    n_seg = va_arg(ap, int);
    va_end(ap);

    segA = (parsec_seg_dc_t *)desc;
    A = (parsec_data_collection_t *)(segA->A_org);

    segment_to_tile(segA, m_seg, n_seg, &m_tile, &n_tile, &offset);

    /* TODO: if not distributed, return 0 */

    return A->rank_of(A, m_tile, n_tile);
}

/*
 * Segments can be handled in two ways:
 * Case 1: The MPI datatype starts from the beginning of the tile (of the original dc) and
 *         uses an offset to get to the beginning of the data of the segment (and a stride).
 * Case 2: The MPI datatype starts from the beginning of the data of the segment (and uses a
 *         stride so it has mb as lda).
 *
 * In case 1, parsec_rbt_data_of() should return a pointer to the beginning of the original tile,
 * i.e., it should return the same thing as data_of() of the original dc for the tile that
 * the given segment falls in.
 * In case 2, parsec_rbt_data_of() should return a pointer to the beginning of the segment,
 * i.e. add the offset to the return value of data_of() of the original dc.
 * The choice between case 1 and case 2 is made in dplasma_datatype_define_subarray(), so
 * these two functions must always correspond.
 *
 * Currently we are using Case 2.
 */
static parsec_data_t *parsec_rbt_data_of(parsec_data_collection_t *desc, ...){
    int m_seg, n_seg, m_tile, n_tile;
    uintptr_t offset, data_start;
    va_list ap;
    parsec_seg_dc_t *segA;
    parsec_data_collection_t *A;

    va_start(ap, desc);
    m_seg = va_arg(ap, int);
    n_seg = va_arg(ap, int);
    va_end(ap);

    segA = (parsec_seg_dc_t *)desc;
    A = &segA->A_org->super;

    segment_to_tile(segA, m_seg, n_seg, &m_tile, &n_tile, &offset);

    data_start = offset*sizeof(dplasma_complex64_t) + (uintptr_t)A->data_of(A, m_tile, n_tile);

    /*
    fprintf(stderr, "Dataof (%d, %d) -> (%d, %d): %p + %llu * %u = %p\n",
            m_seg, n_seg, m_tile, n_tile,
            A->data_of(A, m_tile, n_tile), offset, sizeof(dplasma_complex64_t), data_start);
    */

    return (void *)data_start;
}

/* HE for Hermitian */

/*
 * dplasma_zhebut_New()
 */
parsec_taskpool_t*
dplasma_zhebut_New( parsec_tiled_matrix_t *A, dplasma_complex64_t *U_but_vec, int i_block, int j_block, int level, int *info)
{
    parsec_taskpool_t *parsec_zhebut = NULL;
    parsec_seg_dc_t *seg_descA;
    parsec_memory_pool_t* pool_0;
    dplasma_complex64_t *U_before, *U_after;
    int i, mt, nt, N;

    (void)info;

    seg_descA = (parsec_seg_dc_t *)calloc(1, sizeof(parsec_seg_dc_t));

    /* copy the parsec_tiled_matrix_t part of A into seg_descA */
    memcpy(seg_descA, A, sizeof(parsec_tiled_matrix_t));
    /* overwrite the rank_of() and data_of() */
    seg_descA->super.super.rank_of = parsec_rbt_rank_of;
    seg_descA->super.super.data_of = parsec_rbt_data_of;
    /* store a pointer to A itself */
    seg_descA->A_org = A;
    /* store the level */
    seg_descA->level = level;
    /* store the segment info */
    seg_descA->seg_info = parsec_rbt_calculate_constants(A, level, i_block, j_block);

    N  = A->lm;
    mt = seg_descA->seg_info.tot_seg_cnt_m;
    nt = seg_descA->seg_info.tot_seg_cnt_n;

    pool_0 = (parsec_memory_pool_t*)malloc(sizeof(parsec_memory_pool_t));
    parsec_private_memory_init( pool_0, A->mb * A->nb * sizeof(dplasma_complex64_t) );

    U_before = &U_but_vec[level*N];
    U_after  = &U_but_vec[level*N];

    parsec_zhebut = (parsec_taskpool_t *)parsec_zhebut_new(seg_descA,U_before, U_after, nt, mt, pool_0);

    for(i=0; i<36; i++){
        parsec_arena_datatype_t *adt;
        int type_exists;
        unsigned int m_sz, n_sz;

        type_exists = type_index_to_sizes(&seg_descA->seg_info, i, &m_sz, &n_sz);

        if( type_exists ){
            adt = &((parsec_zhebut_taskpool_t*)parsec_zhebut)->arenas_datatypes[i];
            parsec_add2arena_rect( adt, parsec_datatype_double_complex_t, m_sz, A->nb, A->mb );
        }
    }

    return parsec_zhebut;
}

void
dplasma_zhebut_Destruct( parsec_taskpool_t *tp )
{
    int i;
    parsec_zhebut_taskpool_t *obut = (parsec_zhebut_taskpool_t *)tp;

    for(i=0; i<36; i++){
        parsec_del2arena( &obut->arenas_datatypes[i] );
    }

    parsec_taskpool_free(tp);
}

/* GE for General */

/*
 * dplasma_zgebut_New()
 */
parsec_taskpool_t*
dplasma_zgebut_New( parsec_tiled_matrix_t *A, dplasma_complex64_t *U_but_vec, int i_block, int j_block, int level, int *info)
{
    parsec_taskpool_t *parsec_zgebut = NULL;
    parsec_seg_dc_t *seg_descA;
    parsec_memory_pool_t *pool_0;
    int i, mt, nt, N;
    dplasma_complex64_t *U_before, *U_after;

    (void)info;

    seg_descA = (parsec_seg_dc_t *)calloc(1, sizeof(parsec_seg_dc_t));

    /* copy the parsec_tiled_matrix_t part of A into seg_descA */
    memcpy(seg_descA, A, sizeof(parsec_tiled_matrix_t));
    /* overwrite the rank_of() and data_of() */
    seg_descA->super.super.rank_of = parsec_rbt_rank_of;
    seg_descA->super.super.data_of = parsec_rbt_data_of;
    /* store a pointer to A itself */
    seg_descA->A_org = A;
    /* store the level */
    seg_descA->level = level;
    /* store the segment info */
    seg_descA->seg_info = parsec_rbt_calculate_constants(A, level, i_block, j_block);

    N  = A->lm;
    mt = seg_descA->seg_info.tot_seg_cnt_m;
    nt = seg_descA->seg_info.tot_seg_cnt_n;

    U_before = &U_but_vec[level*N];
    U_after  = &U_but_vec[level*N];

    pool_0 = (parsec_memory_pool_t*)malloc(sizeof(parsec_memory_pool_t));
    parsec_private_memory_init( pool_0, A->mb * A->nb * sizeof(dplasma_complex64_t) );

    parsec_zgebut = (parsec_taskpool_t *)parsec_zgebut_new(seg_descA, U_before, U_after, nt, mt, pool_0);

    for(i=0; i<36; i++){
        parsec_arena_datatype_t *adt;
        int type_exists;
        unsigned int m_sz, n_sz;

        type_exists = type_index_to_sizes(&seg_descA->seg_info, i, &m_sz, &n_sz);

        if( type_exists ){
            adt = &((parsec_zgebut_taskpool_t*)parsec_zgebut)->arenas_datatypes[i];
            parsec_add2arena_rect( adt, parsec_datatype_double_complex_t, m_sz, A->nb, A->mb );
        }
    }

    return parsec_zgebut;
}

void
dplasma_zgebut_Destruct( parsec_taskpool_t *tp )
{
    int i;
    parsec_zgebut_taskpool_t *obut = (parsec_zgebut_taskpool_t *)tp;

    for(i=0; i<36; i++){
        dplasma_matrix_del2arena( &obut->arenas_datatypes[i] );
    }

    parsec_taskpool_free(tp);
}

/*
 * dplasma_zgebmm_New()
 */
parsec_taskpool_t*
dplasma_zgebmm_New( parsec_tiled_matrix_t *A, dplasma_complex64_t *U_but_vec, int i_block, int j_block, int level, dplasma_enum_t trans, int *info)
{
    parsec_taskpool_t *parsec_zgebmm = NULL;
    parsec_seg_dc_t *seg_descA;
    parsec_memory_pool_t *pool_0;
    int i, mt, nt, N;

    (void)info;

    seg_descA = (parsec_seg_dc_t *)calloc(1, sizeof(parsec_seg_dc_t));

    /* copy the parsec_tiled_matrix_t part of A into seg_descA */
    memcpy(seg_descA, A, sizeof(parsec_tiled_matrix_t));
    /* overwrite the rank_of() and data_of() */
    seg_descA->super.super.rank_of = parsec_rbt_rank_of;
    seg_descA->super.super.data_of = parsec_rbt_data_of;
    /* store a pointer to A itself */
    seg_descA->A_org = A;
    /* store the level */
    seg_descA->level = level;
    /* store the segment info */
    seg_descA->seg_info = parsec_rbt_calculate_constants(A, level, i_block, j_block);

    /*
    printf("Apllying zgebmm() in block %d,%d\n",i_block, j_block);
    */

    N  = A->lm;
    U_but_vec = &U_but_vec[level*N];

    mt = seg_descA->seg_info.tot_seg_cnt_m;
    nt = seg_descA->seg_info.tot_seg_cnt_n;

    pool_0 = (parsec_memory_pool_t*)malloc(sizeof(parsec_memory_pool_t));
    parsec_private_memory_init( pool_0, A->mb * A->nb * sizeof(dplasma_complex64_t) );

    parsec_zgebmm = (parsec_taskpool_t *)parsec_zgebmm_new(seg_descA, U_but_vec, nt, mt, trans, pool_0);

    for(i=0; i<36; i++){
        parsec_arena_datatype_t *adt;
        int type_exists;
        unsigned int m_sz, n_sz;

        type_exists = type_index_to_sizes(&seg_descA->seg_info, i, &m_sz, &n_sz);

        if( type_exists ){
            adt = &((parsec_zgebmm_taskpool_t*)parsec_zgebmm)->arenas_datatypes[i];
            parsec_add2arena_rect( adt, parsec_datatype_double_complex_t, m_sz, A->nb, A->mb );
        }
    }

    return parsec_zgebmm;
}

void
dplasma_zgebmm_Destruct( parsec_taskpool_t *tp )
{
    int i;
    parsec_zgebmm_taskpool_t *obmm = (parsec_zgebmm_taskpool_t *)tp;

    for(i=0; i<36; i++){
        dplasma_matrix_del2arena( &obmm->arenas_datatypes[i] );
    }

    parsec_taskpool_free(tp);
}

/*
 * Blocking Interface
 */

static parsec_taskpool_t **iterate_ops(parsec_tiled_matrix_t *A, int tmp_level,
                                    int target_level, int i_block, int j_block,
                                    parsec_taskpool_t **subop,
                                    parsec_context_t *parsec,
                                    dplasma_complex64_t *U_but_vec,
                                    int destroy, int *info)
{
    if(tmp_level == target_level){
        if( (i_block == j_block) ){
            if( destroy ){
                dplasma_zhebut_Destruct(*subop);
            }else{
                *subop = dplasma_zhebut_New(A, U_but_vec, i_block, j_block, target_level, info);
            }
        }else{
            if( destroy ){
                dplasma_zgebut_Destruct(*subop);
            }else{
                *subop = dplasma_zgebut_New(A, U_but_vec, i_block, j_block, target_level, info);
            }
        }
        if( !destroy ){
            parsec_context_add_taskpool(parsec, *subop);
        }
        return subop+1;
    }else{
        if( i_block == j_block ){
            subop = iterate_ops(A, tmp_level+1, target_level, 2*i_block,   2*j_block,   subop, parsec, U_but_vec, destroy, info);
            subop = iterate_ops(A, tmp_level+1, target_level, 2*i_block+1, 2*j_block,   subop, parsec, U_but_vec, destroy, info);
            subop = iterate_ops(A, tmp_level+1, target_level, 2*i_block+1, 2*j_block+1, subop, parsec, U_but_vec, destroy, info);
        }else{
            subop = iterate_ops(A, tmp_level+1, target_level, 2*i_block,   2*j_block,   subop, parsec, U_but_vec, destroy, info);
            subop = iterate_ops(A, tmp_level+1, target_level, 2*i_block+1, 2*j_block,   subop, parsec, U_but_vec, destroy, info);
            subop = iterate_ops(A, tmp_level+1, target_level, 2*i_block,   2*j_block+1, subop, parsec, U_but_vec, destroy, info);
            subop = iterate_ops(A, tmp_level+1, target_level, 2*i_block+1, 2*j_block+1, subop, parsec, U_but_vec, destroy, info);
        }
        return subop;
    }

}

static void RBT_zrandom(int N, dplasma_complex64_t *V)
{
    int i;

    for (i=0; i<N; i++){
        V[i] = (dplasma_complex64_t)exp(((random()/(double)RAND_MAX)-0.5)/10.0);
    }

}

int dplasma_zhebut(parsec_context_t *parsec, parsec_tiled_matrix_t *A, dplasma_complex64_t **U_but_ptr, int levels)
{
    parsec_taskpool_t **subop;
    dplasma_complex64_t *U_but_vec, beta;
    int cur_level, N;
    int info = 0;
    int nbhe = 1<<levels;
    int nbge;
    int final_nt = A->nt/(2*nbhe);
#if defined(DEBUG_BUTTERFLY)
    int i;
#endif
    nbge = (levels>0) ? (1<<(levels-1))*((1<<levels)-1) : 0;

    if( final_nt == 0 ){
        fprintf(stderr,"Too many butterflies. Death by starvation.\n");
        return -1;
    }
    if( A->ln%nbhe != 0 ){
        fprintf(stderr,"Please use a matrix size that is divisible by %d\n - Current Matrix size=%d\n - Number of Hermitian Blocks for this level of RBT=%d\n", 2*nbhe, A->ln, nbhe);
        return -1;
    }

    N = A->lm;

    U_but_vec = (dplasma_complex64_t *)malloc( (levels+1)*N*sizeof(dplasma_complex64_t) );
    *U_but_ptr = U_but_vec;
    srandom(0);
    RBT_zrandom((levels+1)*N, U_but_vec);

    beta = (dplasma_complex64_t)pow(1.0/sqrt(2.0), levels);
    cblas_zscal(levels*N, CBLAS_SADDR(beta), U_but_vec, 1);
#if defined(DEBUG_BUTTERFLY)
    for(i=0; i<levels*N; i++){
        printf("U[%d]: %lf\n",i,creal(U_but_vec[i]));
    }
#endif

    for(cur_level = levels; cur_level >=0; cur_level--){
        nbhe = 1<<cur_level;
        nbge = cur_level > 0 ? (1<<(cur_level-1))*((1<<cur_level)-1) : 0;
        final_nt = A->nt/(2*nbhe);
        if( final_nt == 0 ){
            fprintf(stderr,"Too many butterflies. Death by starvation.\n");
            return -1;
        }
        if( A->ln%nbhe != 0 ){
            fprintf(stderr,"Please use a matrix size that is divisible by %d\n - Current Matrix size=%d\n - Number of Hermitian Blocks for this level of RBT=%d\n", 2*nbhe, A->ln, nbhe);
            return -1;
        }

#if defined(DEBUG_BUTTERFLY)
        printf("\n  =====  Applying Butterfly at level %d\n\n", cur_level);
        fflush(stdout);
#endif

        subop = (parsec_taskpool_t **)malloc((nbhe+nbge) * sizeof(parsec_taskpool_t*));
        (void)iterate_ops(A, 0, cur_level, 0, 0, subop, parsec, U_but_vec, CREATE_N_ENQUEUE, &info);
        dplasma_wait_until_completion(parsec);
        (void)iterate_ops(A, 0, cur_level, 0, 0, subop, parsec, NULL, DESTRUCT, &info);
        free(subop);

#if defined(DEBUG_BUTTERFLY)
        printf("\n\n -+-+-+> Matrix after level %d\n\n", cur_level);
        dplasma_zprint(parsec, dplasmaLower, A);
        printf("\n\n");
#endif

        if( info != 0 ){
            fprintf(stderr,"Terminating the application of butterflies at level: %d\n", cur_level);
            return info;
        }
    }

    return info;
}

max_stars_repo_name: therault/dplasma

max_issues_repo_name: therault/dplasma

max_forks_repo_name: therault/dplasma


text:
/////////////////////////////////////////////////////////////////////
// = NMatrix
//
// A linear algebra library for scientific computation in Ruby.
// NMatrix is part of SciRuby.
//
// NMatrix was originally inspired by and derived from NArray, by
// Masahiro Tanaka: http://narray.rubyforge.org
//
// == Copyright Information
//
// SciRuby is Copyright (c) 2010 - 2012, Ruby Science Foundation
// NMatrix is Copyright (c) 2012, Ruby Science Foundation
//
// Please see LICENSE.txt for additional copyright notices.
//
// == Contributing
//
// By contributing source code to SciRuby, you agree to be bound by
// our Contributor Agreement:
//
// * https://github.com/SciRuby/sciruby/wiki/Contributor-Agreement
//
// == nmatrix.h
//

#ifndef NMATRIX_H
#define NMATRIX_H

#include "nmatrix_config.h"

#include <cblas.h>

#include <math.h>

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ruby.h>
#define RUBY_ZERO INT2FIX(0)

#ifdef BENCHMARK // SOURCE: http://stackoverflow.com/questions/2349776/how-can-i-benchmark-a-c-program-easily
# include <sys/time.h>
# include <sys/resource.h>
#endif

#include "dtypes.h"

#include <stddef.h>
#ifdef HAVE_STDINT_H
# include <stdint.h>
#endif

/*
  Data types used in NArray / NMatrix :
  Please modify these types if your system has any different type.
*/

/* NM_BYTE : unsigned 8-bit integer */
#ifndef HAVE_U_INT8_T
# ifdef HAVE_UINT8_T
typedef uint8_t			u_int8_t;
# else
typedef unsigned char		u_int8_t;
# endif
#endif

//#ifndef HAVE_INT8_T
//typedef char                   int8_t;
//#endif

#ifndef HAVE_INT16_T
# if SIZEOF_SHORT == 2
typedef short                  int16_t;
# else
---->> Please define int16_t manually because sizeof(short) != 2. <<----
# endif
#endif /* HAVE_INT16_T */

#ifndef HAVE_INT32_T
# if SIZEOF_LONG == 4
typedef long                   int32_t;
# else
#  if SIZEOF_INT == 4
typedef int                    int32_t;
#  else
---->> Please define int32_t manually because sizeof(long) != 4. <<----
#  endif
# endif
#endif /* HAVE_INT32_T */

/* unsigned 32-bit integer */
#ifndef HAVE_U_INT32_T
# ifdef HAVE_UINT32_T
typedef uint32_t			u_int32_t;
# else
#  if SIZEOF_LONG == 4
typedef unsigned long                   u_int32_t;
#  else
#   if SIZEOF_INT == 4
typedef unsigned int                    u_int32_t;
#   else
---->> Please define u_int32_t manually because sizeof(long) != 4. <<----
#   endif
#  endif
# endif
#endif /* HAVE_U_INT32_T */

#ifndef HAVE_INT64_T
# if SIZEOF_QUAD == 8
typedef quad                   int64_t;
# else
#  if SIZEOF_LONG == 8
typedef long                   int64_t;
#  else
---->> Please define int64_t manually because sizeof(quad) != 8. <<----
#  endif
# endif
#endif /* HAVE_INT64_T */

/* unsigned 64-bit integer */
#ifndef HAVE_U_INT64_T
# ifdef HAVE_UINT64_T
typedef uint64_t            u_int64_t;
# else
#  if SIZEOF_QUAD == 8
typedef unsigned quad       u_int64_t;
#  else
#   if SIZEOF_LONG == 8
typedef unsigned long       u_int64_t;
#   else
---->> Please define u_int64_t manually because sizeof(quad) != 8. <<----
#   endif
#  endif
# endif
#endif /* HAVE_U_INT64_T */

#ifndef HAVE_SIZE_T /// If you modify this, make sure to modify the definition of y_size_t and Y_SIZE_T!
typedef u_int64_t    size_t;
# define NM_SIZE_T   NM_INT64
#else
# if SIZEOF_SIZE_T == 8
#  define NM_SIZE_T  NM_INT64
# else
#  if SIZEOF_SIZE_T == 4
#   define NM_SIZE_T NM_INT32
#  else
---->> Please define size_t and y_size_t manually because sizeof(size_t) is neither 8 nor 4. <<----
#  endif
# endif
#endif

// for when we need to return array indices.
// This must never be larger than size_t
typedef uint32_t    y_size_t;
#define Y_SIZE_T    NM_INT32

#ifdef HAVE_STDBOOL_H
# include <stdbool.h>
#else
typedef char    bool;
# define true    1;
# define false   0;
#endif

typedef struct { float r,i; } complex64;
typedef struct { double r,i; } complex128;
typedef struct { int16_t n,d; } rational32;
typedef struct { int32_t n,d; } rational64;
typedef struct { int64_t n,d; } rational128;

#if SIZEOF_INT == 8
# define DEFAULT_DTYPE  NM_INT64
#else
# if SIZEOF_INT == 4
#  define DEFAULT_DTYPE NM_INT32
# else
#  define DEFAULT_DTYPE NM_INT16
# endif
#endif

#define YALE_GROWTH_CONSTANT    1.5

enum NMatrix_STypes {
  S_DENSE,
  S_LIST,
  S_YALE,
  S_TYPES
};

// Element-wise operations (see blas/elementwise.template.c)
enum NMatrix_Ops {
  NM_OP_ADD = '+',
  NM_OP_SUB = '-',
  NM_OP_MUL = '*',
  NM_OP_DIV = '/',
  NM_OP_MOD = '%',
  NM_OP_BANG = '!',
  NM_OP_NEG, // unary minus
  NM_OP_EQEQ, // ==
  NM_OP_NEQ, // !=
  NM_OP_GT = '>', // >
  NM_OP_LT = '<', // <
  NM_OP_GTE = ',', // >=
  NM_OP_LTE = '.', // <=
  NM_OP_NOT = '~',
  NM_OP_AND = '&',
  NM_OP_OR = '|',
  NM_OP_XOR = '^',
  NM_OP_LSH, // <<
  NM_OP_RSH // >>
};

/* Singly-linked ordered list
 * - holds keys and values
 * - no duplicate keys
 * - keys are ordered
 * - values may be lists themselves
 */
typedef struct l_node { /* Linked list node */
  size_t key;
  void*  val;
  struct l_node * next; // next
} NODE;

typedef struct l_list {
  NODE* first;
} LIST;

// two vectors and a capacity
typedef struct y_vector {
  void*  ija;
  void*  a;
  size_t capacity;
} VECTOR;

typedef struct common_s { // Common elements found in all _s types.
  int8_t    dtype;
  size_t    rank;
  size_t*   shape;
} STORAGE;

typedef struct list_s {
  int8_t    dtype;
  size_t    rank;
  size_t*   shape;
  void*     default_val;
  LIST*     rows;
} LIST_STORAGE;

typedef struct dense_s {
  int8_t    dtype;
  size_t    rank;
  size_t*   shape;
  void*     elements;
} DENSE_STORAGE;

typedef struct yale_s {
  int8_t    dtype;
  size_t    rank;
  size_t*   shape;
  size_t    ndnz; // strictly non-diagonal non-zero count!
  size_t    capacity;
  int8_t    index_dtype;
  void*     ija;
  void*     a;
} YALE_STORAGE;

typedef struct numeric_matrix {
  int8_t   stype;             /* method of storage (csc, dense, etc) */
  STORAGE* storage;           /* pointer to storage struct */
} NMATRIX;

/* Local */

typedef union {
  u_int8_t b[2];
  int16_t s;
} nm_size16_t;

typedef union {
  u_int8_t b[4];
  int32_t  i;
  float    f;
} nm_size32_t;

typedef union {
  u_int8_t  b[8];
  int64_t   q;
  float     f[2];
  double    d;
  complex64 c;
} nm_size64_t;

typedef union {
  u_int8_t   b[16];
  int64_t    i[2];
  double     d[2];
  float      f[4];
  complex64  c[2];
  complex128 z;
  rational32 r[4];
  rational64 ra[2];
  rational128 rat;
  VALUE      v[2];
} nm_size128_t;

// For calling cblas_gemm functions (see cblas.c)
typedef struct cblas_param_t {
  int M, N, K, lda, ldb, ldc;
  void *A, *B, *C;
  nm_size128_t alpha, beta;
} DENSE_PARAM;

// Formerly in smmp.h:
typedef struct smmp_param_t {
  void *ia, *ja, *a;
  bool diag;
} YALE_PARAM;

// Shouldn't be necessary, as they're defined in nmatrix.h:
// (Oddly, though, these fix the error.)
/*typedef uint8_t   u_int8_t;
typedef uint16_t  u_int16_t;
typedef uint32_t  u_int32_t;
typedef uint64_t  u_int64_t; */

// rational.c
int64_t nmrb_gcd(int64_t x, int64_t y);

// BLAS functions
#define SMMP_MAX_THREE(a,b,c) ((a)>(b) ? ( (a)>(c) ? (a) : (c) ) : ( (b)>(c) ? (b) : (c) ))
#define SMMP_MIN(a,b) ((a)>(b) ? (b) : (a))
#define SMMP_MAX(a,b) ((a)>(b) ? (a) : (b))

void transp(y_size_t n, y_size_t m, void* ia, void* ja, bool diaga, void* a, void* ib, void* jb, void* b, bool move, int8_t itype, int8_t dtype);

void i8_symbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_symbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_symbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_symbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);

void i8_b_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i16_b_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i32_b_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i64_b_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i8_i8_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i16_i8_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i32_i8_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i64_i8_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i8_i16_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i16_i16_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i32_i16_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i64_i16_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i8_i32_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i16_i32_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i32_i32_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i64_i32_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i8_i64_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i16_i64_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i32_i64_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i64_i64_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i8_f32_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i16_f32_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i32_f32_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i64_f32_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i8_f64_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i16_f64_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i32_f64_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i64_f64_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i8_c64_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i16_c64_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i32_c64_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i64_c64_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i8_c128_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i16_c128_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i32_c128_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i64_c128_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i8_r32_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i16_r32_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i32_r32_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i64_r32_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i8_r64_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i16_r64_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i32_r64_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i64_r64_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i8_r128_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i16_r128_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i32_r128_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i64_r128_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i8_v_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i16_v_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i32_v_smmp_sort_columns_(y_size_t, YALE_PARAM);
void i64_v_smmp_sort_columns_(y_size_t, YALE_PARAM);

void i8_b_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_b_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_b_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_b_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_i8_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_i8_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_i8_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_i8_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_i16_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_i16_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_i16_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_i16_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_i32_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_i32_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_i32_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_i32_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_i64_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_i64_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_i64_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_i64_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_f32_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_f32_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_f32_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_f32_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_f64_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_f64_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_f64_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_f64_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_c64_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_c64_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_c64_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_c64_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_c128_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_c128_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_c128_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_c128_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_r32_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_r32_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_r32_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_r32_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_r64_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_r64_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_r64_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_r64_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_r128_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_r128_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_r128_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_r128_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_v_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_v_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_v_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_v_smmp(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);

void i8_b_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_b_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_b_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_b_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_i8_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_i8_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_i8_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_i8_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_i16_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_i16_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_i16_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_i16_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_i32_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_i32_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_i32_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_i32_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_i64_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_i64_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_i64_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_i64_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_f32_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_f32_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_f32_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_f32_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_f64_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_f64_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_f64_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_f64_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_c64_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_c64_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_c64_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_c64_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_c128_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_c128_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_c128_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_c128_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_r32_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_r32_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_r32_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_r32_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_r64_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_r64_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_r64_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_r64_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_r128_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_r128_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_r128_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_r128_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i8_v_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i16_v_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i32_v_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);
void i64_v_numbmm_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, YALE_PARAM);

void i8_b_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i16_b_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i32_b_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i64_b_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i8_i8_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i16_i8_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i32_i8_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i64_i8_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i8_i16_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i16_i16_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i32_i16_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i64_i16_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i8_i32_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i16_i32_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i32_i32_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i64_i32_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i8_i64_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i16_i64_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i32_i64_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i64_i64_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i8_f32_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i16_f32_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i32_f32_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i64_f32_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i8_f64_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i16_f64_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i32_f64_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i64_f64_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i8_c64_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i16_c64_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i32_c64_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i64_c64_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i8_c128_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i16_c128_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i32_c128_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i64_c128_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i8_r32_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i16_r32_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i32_r32_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i64_r32_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i8_r64_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i16_r64_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i32_r64_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i64_r64_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i8_r128_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i16_r128_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i32_r128_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i64_r128_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i8_v_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i16_v_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i32_v_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);
void i64_v_transp_(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool);

// For binary operations involving matrices that need to be casted.
typedef struct storage_pair_t {
  STORAGE* left;
  STORAGE* right;
} STORAGE_PAIR;

#ifndef NMATRIX_C
extern VALUE cNMatrix;

extern const int nm_sizeof[NM_TYPES+1];
#endif

#define NM_MAX_RANK 15

#define UnwrapNMatrix(obj,var)  Data_Get_Struct(obj, struct numeric_matrix, var)
#define IsNMatrix(obj)          (rb_obj_is_kind_of(obj, CNMatrix)==Qtrue)

#define NM_STORAGE(val)         (((struct numeric_matrix*)DATA_PTR(val))->storage)
//#define NM_PTR(a, p)            ((a)->ptr+(p)*nm_sizeof[(a)->type])
#define NM_STRUCT(val)          ((struct numeric_matrix*)DATA_PTR(val))
//#define NM_PTR_TYPE(val,type)   (type)(((struct numeric_matrix*)DATA_PTR(val))->ptr)
#define NM_RANK(val)            (((STORAGE*)(NM_STORAGE(val)))->rank)
#define NM_DTYPE(val)           (((STORAGE*)(NM_STORAGE(val)))->dtype)
#define NM_STYPE(val)           (((struct numeric_matrix*)DATA_PTR(val))->stype)
#define NM_SHAPE(val,i)         (((STORAGE*)(NM_STORAGE(val)))->shape[(i)])
#define NM_SHAPE0(val)          (((struct numeric_matrix*)DATA_PTR(val))->shape[0])
#define NM_SHAPE1(val)          (((struct numeric_matrix*)DATA_PTR(val))->shape[1])
#define NM_SIZEOF_DTYPE(val)    (nm_sizeof[NM_DTYPE(val)])
#define NM_REF(val,coords)      (RefFuncs[NM_STYPE(val)]( NM_STORAGE(val), coords, NM_SIZEOF_DTYPE(val) ))

#define NM_IsNMatrix(obj) (rb_obj_is_kind_of(obj, cNMatrix)==Qtrue)
#define NM_IsArray(obj)   (TYPE(obj)==T_ARRAY || rb_obj_is_kind_of(obj,cNMatrix)==Qtrue)
#define NM_IsROBJ(d) ((d)->dtype==NM_ROBJ)
#define NM_IsINTEGER(a) \
    (NM_DTYPE(a)==NM_BYTE || NM_DTYPE(a)==NM_INT8 || NM_DTYPE(a)==NM_INT16 || NM_DTYPE(a)==NM_INT32 || NM_DTYPE(a)==NM_INT64)
#define NM_IsCOMPLEX(a) \
    (NM_DTYPE(a)==NM_COMPLEX32 || NM_DTYPE(a)==NM_COMPLEX64)
#define NM_MAX(a,b) (((a)>(b))?(a):(b))
#define NM_SWAP(a,b,tmp) {(tmp)=(a);(a)=(b);(b)=(tmp);}

//#define NUM2REAL(v) NUM2DBL( rb_funcall((v),nm_id_real,0) ) // deprecated
#define REAL2DBL(v) NUM2DBL( rb_funcall((v),nm_id_real,0) )
//#define NUM2IMAG(v) NUM2DBL( rb_funcall((v),nm_id_imag,0) ) // deprecated
#define IMAG2DBL(v) NUM2DBL( rb_funcall((v),nm_id_imag,0) )

#define NUM2NUMER(v) NUM2INT( rb_funcall((v), nm_id_numer,0) ) // deprecated
#define NUMER2INT(v) NUM2INT( rb_funcall((v), nm_id_numer,0) )
#define NUM2DENOM(v) NUM2INT( rb_funcall((v), nm_id_denom,0) ) // deprecated
#define DENOM2INT(v) NUM2INT( rb_funcall((v), nm_id_denom,0) )

#define IS_NUMERIC(v)   (FIXNUM_P(v) || TYPE(v) == T_FLOAT || TYPE(v) == T_COMPLEX || TYPE(v) == T_RATIONAL)
#define IS_STRING(v)    (TYPE(v) == T_STRING)

#define CheckNMatrixType(v)   if (TYPE(v) != T_DATA || RDATA(v)->dfree != (RUBY_DATA_FUNC)nm_delete) rb_raise(rb_eTypeError, "expected NMatrix on left-hand side of operation");

//#define YALE_JA_START(sptr)             (((YALE_STORAGE*)(sptr))->shape[0]+1)
#define YALE_IJA(sptr,elem_size,i)          (void*)( (char*)(((YALE_STORAGE*)(sptr))->ija) + i * elem_size )
//#define YALE_JA(sptr,dtype,j)           ((((dtype)*)((YALE_STORAGE*)(sptr))->ija)[(YALE_JA_START(sptr))+j])
#define YALE_ROW_LENGTH(sptr,elem_size,i)   (*(size_t*)YALE_IA((sptr),(elem_size),(i)+1) - *(size_t*)YALE_IJA((sptr),(elem_size),(i)))
#define YALE_A(sptr,elem_size,i)            (void*)((char*)(((YALE_STORAGE*)(sptr))->a) + elem_size * i)
#define YALE_DIAG(sptr, elem_size, i)       ( YALE_A((sptr),(elem_size),(i)) )
//#define YALE_LU(sptr,dtype,i,j)             (((dtype)*)(((YALE_STORAGE*)(sptr))->a)[ YALE_JA_START(sptr) +  ])
#define YALE_MINIMUM(sptr)                  (((YALE_STORAGE*)(sptr))->shape[0]*2 + 1) // arbitrarily defined
#define YALE_SIZE_PTR(sptr,elem_size)       (void*)((char*)((YALE_STORAGE*)(sptr))->ija + ((YALE_STORAGE*)(sptr))->shape[0]*elem_size )
#define YALE_MAX_SIZE(sptr)                 (((YALE_STORAGE*)(sptr))->shape[0] * ((YALE_STORAGE*)(sptr))->shape[1] + 1)
#define YALE_IA_SIZE(sptr)                  ((YALE_STORAGE*)(sptr))->shape[0]

// None of these next three return anything. They set a reference directly.
#define YaleGetIJA(victim,s,i)              (SetFuncs[Y_SIZE_T][(s)->index_dtype](1, &(victim), 0, YALE_IJA((s), nm_sizeof[s->index_dtype], (i)), 0))
#define YaleSetIJA(i,s,from)                (SetFuncs[s->index_dtype][Y_SIZE_T](1, YALE_IJA((s), nm_sizeof[s->index_dtype], (i)), 0, &(from), 0))
#define YaleGetSize(sz,s)                   (SetFuncs[Y_SIZE_T][(s)->index_dtype](1, &sz, 0, (YALE_SIZE_PTR((s), nm_sizeof[(s)->index_dtype])), 0))
//#define YALE_FIRST_NZ_ROW_ENTRY(sptr,elem_size,i)

#if !defined RSTRING_LEN
#define RSTRING_LEN(a) RSTRING(a)->len
#endif
#if !defined RSTRING_PTR
#define RSTRING_PTR(a) RSTRING(a)->ptr
#endif
#if !defined RARRAY_LEN
#define RARRAY_LEN(a) RARRAY(a)->len
#endif
#if !defined RARRAY_PTR
#define RARRAY_PTR(a) RARRAY(a)->ptr
#endif

#define NM_INDEX_TYPES  NM_FLOAT32

typedef void     (*nm_setfunc_t[NM_TYPES][NM_TYPES])(); // copy functions
typedef void     (*nm_incfunc_t[NM_TYPES])();           // increment functions
typedef void*    (*nm_stype_ref_t[S_TYPES])(STORAGE*, size_t*);        // get/ref
typedef VALUE    (*nm_stype_ins_t[S_TYPES])(STORAGE*, size_t*, VALUE); // insert
typedef STORAGE* (*nm_create_storage_t[S_TYPES])();
typedef STORAGE* (*nm_cast_copy_storage_t[S_TYPES])();
typedef STORAGE* (*nm_scast_copy_storage_t[S_TYPES][S_TYPES])();
typedef NMATRIX* (*nm_matrix_multiply_op_t[S_TYPES])();
typedef NMATRIX* (*nm_elementwise_binary_op_casted_t[S_TYPES])();
typedef int      (*nm_d_elementwise_binary_op_t[NM_TYPES])();
typedef int      (*nm_y_elementwise_binary_op_t[NM_TYPES][NM_INDEX_TYPES])();
typedef bool     (*nm_compare_t[S_TYPES])();
typedef void     (*nm_delete_t[S_TYPES])();
typedef void     (*nm_mark_t[S_TYPES])(void*);
typedef void     (*nm_gemm_t[NM_TYPES])();           // general matrix/matrix multiply
typedef void     (*nm_dense_transpose_t[NM_TYPES])(); // dense transpose
typedef void     (*nm_gemv_t[NM_TYPES])();           // general matrix/vector multiply
typedef void     (*nm_smmp_t[NM_TYPES][NM_INDEX_TYPES])(); // sparse (yale) multiply
typedef void     (*nm_smmp_transpose_t[NM_TYPES][NM_INDEX_TYPES])(y_size_t, y_size_t, YALE_PARAM, YALE_PARAM, bool); // sparse (yale) transpose
//typedef void (*nm_setsf_t[S_TYPES][S_TYPES])();
//typedef void (*nm_setdf_t[NM_DTYPES][NM_DTYPES])();

extern nm_setfunc_t SetFuncs;
extern nm_incfunc_t Increment;
extern ID nm_id_real, nm_id_imag;
extern ID nm_id_denom, nm_id_numer;
extern ID nm_id_mult, nm_id_multeq, nm_id_add;

/* blas.c */
int r32gemm(enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const rational32 alpha, const rational32* A, const int lda, const rational32* B, const int ldb, const rational32 beta, rational32* C, const int ldc);
int r32gemv(enum CBLAS_TRANSPOSE Trans, const size_t M, const size_t N, const rational32 alpha, const rational32* A, const size_t lda, const rational32* X, const int incX, const rational32 beta, rational32* Y, const int incY);
int r64gemm(enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const rational64 alpha, const rational64* A, const int lda, const rational64* B, const int ldb, const rational64 beta, rational64* C, const int ldc);
int r64gemv(enum CBLAS_TRANSPOSE Trans, const size_t M, const size_t N, const rational64 alpha, const rational64* A, const size_t lda, const rational64* X, const int incX, const rational64 beta, rational64* Y, const int incY);
int r128gemm(enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const rational128 alpha, const rational128* A, const int lda, const rational128* B, const int ldb, const rational128 beta, rational128* C, const int ldc);
int r128gemv(enum CBLAS_TRANSPOSE Trans, const size_t M, const size_t N, const rational128 alpha, const rational128* A, const size_t lda, const rational128* X, const int incX, const rational128 beta, rational128* Y, const int incY);
int bgemm(enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const u_int8_t alpha, const u_int8_t* A, const int lda, const u_int8_t* B, const int ldb, const u_int8_t beta, u_int8_t* C, const int ldc);
int bgemv(enum CBLAS_TRANSPOSE Trans, const size_t M, const size_t N, const u_int8_t alpha, const u_int8_t* A, const size_t lda, const u_int8_t* X, const int incX, const u_int8_t beta, u_int8_t* Y, const int incY);
int i8gemm(enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const int8_t alpha, const int8_t* A, const int lda, const int8_t* B, const int ldb, const int8_t beta, int8_t* C, const int ldc);
int i8gemv(enum CBLAS_TRANSPOSE Trans, const size_t M, const size_t N, const int8_t alpha, const int8_t* A, const size_t lda, const int8_t* X, const int incX, const int8_t beta, int8_t* Y, const int incY);
int i16gemm(enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const int16_t alpha, const int16_t* A, const int lda, const int16_t* B, const int ldb, const int16_t beta, int16_t* C, const int ldc);
int i16gemv(enum CBLAS_TRANSPOSE Trans, const size_t M, const size_t N, const int16_t alpha, const int16_t* A, const size_t lda, const int16_t* X, const int incX, const int16_t beta, int16_t* Y, const int incY);
int i32gemm(enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const int32_t alpha, const int32_t* A, const int lda, const int32_t* B, const int ldb, const int32_t beta, int32_t* C, const int ldc);
int i32gemv(enum CBLAS_TRANSPOSE Trans, const size_t M, const size_t N, const int32_t alpha, const int32_t* A, const size_t lda, const int32_t* X, const int incX, const int32_t beta, int32_t* Y, const int incY);
int i64gemm(enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const int64_t alpha, const int64_t* A, const int lda, const int64_t* B, const int ldb, const int64_t beta, int64_t* C, const int ldc);
int i64gemv(enum CBLAS_TRANSPOSE Trans, const size_t M, const size_t N, const int64_t alpha, const int64_t* A, const size_t lda, const int64_t* X, const int incX, const int64_t beta, int64_t* Y, const int incY);
int vgemm(enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const VALUE alpha, const VALUE* A, const int lda, const VALUE* B, const int ldb, const VALUE beta, VALUE* C, const int ldc);

// void %%TYPE_ABBREV%%transp(unsigned int M, unsigned int N, const %%TYPE%%* A, int lda, %%TYPE%%* B, int ldb);
void btransp(const unsigned int, const unsigned int, const u_int8_t*, const int, u_int8_t*, const int);
void i8transp(const unsigned int, const unsigned int, const int8_t*, const int, int8_t*, const int);
void i16transp(const unsigned int, const unsigned int, const int16_t*, const int, int16_t*, const int);
void i32transp(const unsigned int, const unsigned int, const int32_t*, const int, int32_t*, const int);
void i64transp(const unsigned int, const unsigned int, const int64_t*, const int, int64_t*, const int);
void f32transp(const unsigned int, const unsigned int, const float*, const int, float*, const int);
void f64transp(const unsigned int, const unsigned int, const double*, const int, double*, const int);
void c64transp(const unsigned int, const unsigned int, const complex64*, const int, complex64*, const int);
void c128transp(const unsigned int, const unsigned int, const complex128*, const int, complex128*, const int);
void r32transp(const unsigned int, const unsigned int, const rational32*, const int, rational32*, const int);
void r64transp(const unsigned int, const unsigned int, const rational64*, const int, rational64*, const int);
void r128transp(const unsigned int, const unsigned int, const rational128*, const int, rational128*, const int);
void vtransp(const unsigned int, const unsigned int, const VALUE*, const int, VALUE*, const int);

int nm_d_b_elementwise(const u_int8_t* A, const u_int8_t* B, u_int8_t* C, size_t n, enum NMatrix_Ops op);
int nm_d_i8_elementwise(const int8_t* A, const int8_t* B, int8_t* C, size_t n, enum NMatrix_Ops op);
int nm_d_i16_elementwise(const int16_t* A, const int16_t* B, int16_t* C, size_t n, enum NMatrix_Ops op);
int nm_d_i32_elementwise(const int32_t* A, const int32_t* B, int32_t* C, size_t n, enum NMatrix_Ops op);
int nm_d_i64_elementwise(const int64_t* A, const int64_t* B, int64_t* C, size_t n, enum NMatrix_Ops op);
int nm_d_f32_elementwise(const float* A, const float* B, float* C, size_t n, enum NMatrix_Ops op);
int nm_d_f64_elementwise(const double* A, const double* B, double* C, size_t n, enum NMatrix_Ops op);
int nm_d_c64_elementwise(const complex64* A, const complex64* B, complex64* C, size_t n, enum NMatrix_Ops op);
int nm_d_c128_elementwise(const complex128* A, const complex128* B, complex128* C, size_t n, enum NMatrix_Ops op);
int nm_d_r32_elementwise(const rational32* A, const rational32* B, rational32* C, size_t n, enum NMatrix_Ops op);
int nm_d_r64_elementwise(const rational64* A, const rational64* B, rational64* C, size_t n, enum NMatrix_Ops op);
int nm_d_r128_elementwise(const rational128* A, const rational128* B, rational128* C, size_t n, enum NMatrix_Ops op);
int nm_d_v_elementwise(const VALUE* A, const VALUE* B, VALUE* C, size_t n, enum NMatrix_Ops op);

// These are in blas.c but are needed by smmp2.c (the smmp template stuff)
rational128 r128_muldiv(int64_t anum, int64_t aden, int64_t bnum, int64_t bden, char k);
rational128 r128_addsub(int64_t anum, int64_t aden, int64_t bnum, int64_t bden, char k);
rational128 r128_mod(int64_t anum, int64_t aden, int64_t bnum, int64_t bden);
rational128 r128_bang(int64_t, int64_t);
rational128 r128_negate(int64_t, int64_t);
rational64 r64_muldiv(int64_t anum, int64_t aden, int64_t bnum, int64_t bden, char k);
rational64 r64_addsub(int64_t anum, int64_t aden, int64_t bnum, int64_t bden, char k);
rational64 r64_mod(int32_t anum, int32_t aden, int32_t bnum, int32_t bden);
rational64 r64_bang(int32_t, int32_t);
rational64 r64_negate(int32_t, int32_t);
rational32 r32_muldiv(int64_t anum, int64_t aden, int64_t bnum, int64_t bden, char k);
rational32 r32_addsub(int64_t anum, int64_t aden, int64_t bnum, int64_t bden, char k);
rational32 r32_mod(int16_t anum, int16_t aden, int16_t bnum, int16_t bden);
rational32 r32_bang(int16_t, int16_t);
rational32 r32_negate(int16_t, int16_t);

rational32 BOOL2R32(bool);
rational64 BOOL2R64(bool);
rational128 BOOL2R128(bool);

/* cblas.c */
DENSE_PARAM init_cblas_params_for_nm_multiply_matrix(int8_t dtype);
void cblas_bgemm_(enum CBLAS_ORDER Order, enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p);
void cblas_bgemv_(enum CBLAS_ORDER Order, enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p);
void cblas_i8gemm_(enum CBLAS_ORDER Order, enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p);
void cblas_i8gemv_(enum CBLAS_ORDER Order, enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p);
void cblas_i16gemm_(enum CBLAS_ORDER Order, enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p);
void cblas_i16gemv_(enum CBLAS_ORDER Order, enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p);
void cblas_i32gemm_(enum CBLAS_ORDER Order, enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p);
void cblas_i32gemv_(enum CBLAS_ORDER Order, enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p);
void cblas_i64gemm_(enum CBLAS_ORDER Order, enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p);
void cblas_i64gemv_(enum CBLAS_ORDER Order, enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p);
void cblas_sgemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p);
void cblas_sgemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p);
void cblas_dgemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p);
void cblas_dgemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p);
void cblas_cgemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p);
void cblas_cgemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p);
void cblas_zgemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p);
void cblas_zgemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p);
void cblas_r32gemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p);
void cblas_r32gemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p);
void cblas_r64gemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p);
void cblas_r64gemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p);
void cblas_r128gemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p);
void cblas_r128gemv_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, DENSE_PARAM p);
void cblas_vgemm_(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, DENSE_PARAM p);

/* smmp2.c */
int i8_b_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int8_t* ija, const u_int8_t* ijb, const u_int8_t* ijc, u_int8_t* a, u_int8_t* b, u_int8_t* c);
int i16_b_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int16_t* ija, const u_int16_t* ijb, const u_int16_t* ijc, u_int8_t* a, u_int8_t* b, u_int8_t* c);
int i32_b_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int32_t* ija, const u_int32_t* ijb, const u_int32_t* ijc, u_int8_t* a, u_int8_t* b, u_int8_t* c);
int i64_b_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int64_t* ija, const u_int64_t* ijb, const u_int64_t* ijc, u_int8_t* a, u_int8_t* b, u_int8_t* c);
int i8_i8_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int8_t* ija, const u_int8_t* ijb, const u_int8_t* ijc, int8_t* a, int8_t* b, int8_t* c);
int i16_i8_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int16_t* ija, const u_int16_t* ijb, const u_int16_t* ijc, int8_t* a, int8_t* b, int8_t* c);
int i32_i8_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int32_t* ija, const u_int32_t* ijb, const u_int32_t* ijc, int8_t* a, int8_t* b, int8_t* c);
int i64_i8_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int64_t* ija, const u_int64_t* ijb, const u_int64_t* ijc, int8_t* a, int8_t* b, int8_t* c);
int i8_i16_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int8_t* ija, const u_int8_t* ijb, const u_int8_t* ijc, int16_t* a, int16_t* b, int16_t* c);
int i16_i16_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int16_t* ija, const u_int16_t* ijb, const u_int16_t* ijc, int16_t* a, int16_t* b, int16_t* c);
int i32_i16_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int32_t* ija, const u_int32_t* ijb, const u_int32_t* ijc, int16_t* a, int16_t* b, int16_t* c);
int i64_i16_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int64_t* ija, const u_int64_t* ijb, const u_int64_t* ijc, int16_t* a, int16_t* b, int16_t* c);
int i8_i32_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int8_t* ija, const u_int8_t* ijb, const u_int8_t* ijc, int32_t* a, int32_t* b, int32_t* c);
int i16_i32_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int16_t* ija, const u_int16_t* ijb, const u_int16_t* ijc, int32_t* a, int32_t* b, int32_t* c);
int i32_i32_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int32_t* ija, const u_int32_t* ijb, const u_int32_t* ijc, int32_t* a, int32_t* b, int32_t* c);
int i64_i32_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int64_t* ija, const u_int64_t* ijb, const u_int64_t* ijc, int32_t* a, int32_t* b, int32_t* c);
int i8_i64_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int8_t* ija, const u_int8_t* ijb, const u_int8_t* ijc, int64_t* a, int64_t* b, int64_t* c);
int i16_i64_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int16_t* ija, const u_int16_t* ijb, const u_int16_t* ijc, int64_t* a, int64_t* b, int64_t* c);
int i32_i64_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int32_t* ija, const u_int32_t* ijb, const u_int32_t* ijc, int64_t* a, int64_t* b, int64_t* c);
int i64_i64_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int64_t* ija, const u_int64_t* ijb, const u_int64_t* ijc, int64_t* a, int64_t* b, int64_t* c);
int i8_f32_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int8_t* ija, const u_int8_t* ijb, const u_int8_t* ijc, float* a, float* b, float* c);
int i16_f32_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int16_t* ija, const u_int16_t* ijb, const u_int16_t* ijc, float* a, float* b, float* c);
int i32_f32_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int32_t* ija, const u_int32_t* ijb, const u_int32_t* ijc, float* a, float* b, float* c);
int i64_f32_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int64_t* ija, const u_int64_t* ijb, const u_int64_t* ijc, float* a, float* b, float* c);
int i8_f64_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int8_t* ija, const u_int8_t* ijb, const u_int8_t* ijc, double* a, double* b, double* c);
int i16_f64_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int16_t* ija, const u_int16_t* ijb, const u_int16_t* ijc, double* a, double* b, double* c);
int i32_f64_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int32_t* ija, const u_int32_t* ijb, const u_int32_t* ijc, double* a, double* b, double* c);
int i64_f64_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int64_t* ija, const u_int64_t* ijb, const u_int64_t* ijc, double* a, double* b, double* c);
int i8_c64_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int8_t* ija, const u_int8_t* ijb, const u_int8_t* ijc, complex64* a, complex64* b, complex64* c);
int i16_c64_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int16_t* ija, const u_int16_t* ijb, const u_int16_t* ijc, complex64* a, complex64* b, complex64* c);
int i32_c64_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int32_t* ija, const u_int32_t* ijb, const u_int32_t* ijc, complex64* a, complex64* b, complex64* c);
int i64_c64_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int64_t* ija, const u_int64_t* ijb, const u_int64_t* ijc, complex64* a, complex64* b, complex64* c);
int i8_c128_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int8_t* ija, const u_int8_t* ijb, const u_int8_t* ijc, complex128* a, complex128* b, complex128* c);
int i16_c128_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int16_t* ija, const u_int16_t* ijb, const u_int16_t* ijc, complex128* a, complex128* b, complex128* c);
int i32_c128_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int32_t* ija, const u_int32_t* ijb, const u_int32_t* ijc, complex128* a, complex128* b, complex128* c);
int i64_c128_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int64_t* ija, const u_int64_t* ijb, const u_int64_t* ijc, complex128* a, complex128* b, complex128* c);
int i8_r32_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int8_t* ija, const u_int8_t* ijb, const u_int8_t* ijc, rational32* a, rational32* b, rational32* c);
int i16_r32_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int16_t* ija, const u_int16_t* ijb, const u_int16_t* ijc, rational32* a, rational32* b, rational32* c);
int i32_r32_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int32_t* ija, const u_int32_t* ijb, const u_int32_t* ijc, rational32* a, rational32* b, rational32* c);
int i64_r32_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int64_t* ija, const u_int64_t* ijb, const u_int64_t* ijc, rational32* a, rational32* b, rational32* c);
int i8_r64_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int8_t* ija, const u_int8_t* ijb, const u_int8_t* ijc, rational64* a, rational64* b, rational64* c);
int i16_r64_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int16_t* ija, const u_int16_t* ijb, const u_int16_t* ijc, rational64* a, rational64* b, rational64* c);
int i32_r64_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int32_t* ija, const u_int32_t* ijb, const u_int32_t* ijc, rational64* a, rational64* b, rational64* c);
int i64_r64_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int64_t* ija, const u_int64_t* ijb, const u_int64_t* ijc, rational64* a, rational64* b, rational64* c);
int i8_r128_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int8_t* ija, const u_int8_t* ijb, const u_int8_t* ijc, rational128* a, rational128* b, rational128* c);
int i16_r128_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int16_t* ija, const u_int16_t* ijb, const u_int16_t* ijc, rational128* a, rational128* b, rational128* c);
int i32_r128_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int32_t* ija, const u_int32_t* ijb, const u_int32_t* ijc, rational128* a, rational128* b, rational128* c);
int i64_r128_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int64_t* ija, const u_int64_t* ijb, const u_int64_t* ijc, rational128* a, rational128* b, rational128* c);
int i8_v_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int8_t* ija, const u_int8_t* ijb, const u_int8_t* ijc, VALUE* a, VALUE* b, VALUE* c);
int i16_v_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int16_t* ija, const u_int16_t* ijb, const u_int16_t* ijc, VALUE* a, VALUE* b, VALUE* c);
int i32_v_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int32_t* ija, const u_int32_t* ijb, const u_int32_t* ijc, VALUE* a, VALUE* b, VALUE* c);
int i64_v_ew(y_size_t n, y_size_t m, enum NMatrix_Ops op, const u_int64_t* ija, const u_int64_t* ijb, const u_int64_t* ijc, VALUE* a, VALUE* b, VALUE* c);

/* dense.c */
DENSE_STORAGE*  create_dense_storage(int8_t dtype, size_t* shape, size_t rank, void* elements, size_t elements_length);
void            delete_dense_storage(DENSE_STORAGE* s);
void            mark_dense_storage(void* s);
DENSE_STORAGE*  cast_copy_dense_storage(DENSE_STORAGE* rhs, int8_t new_dtype);

size_t          count_dense_storage_elements(const DENSE_STORAGE* s);
bool            dense_storage_eqeq(const DENSE_STORAGE*, const DENSE_STORAGE*);

size_t          dense_storage_pos(DENSE_STORAGE* s, size_t* coords);
void*           dense_storage_get(DENSE_STORAGE* s, size_t* coords);
void            dense_storage_set(DENSE_STORAGE* s, size_t* coords, void* val);

/* list.c */
LIST_STORAGE*   create_list_storage(int8_t dtype, size_t* shape, size_t rank, void* init_val);
void            delete_list_storage(LIST_STORAGE* s);
void            mark_list_storage(void* s);
LIST_STORAGE*   cast_copy_list_storage(LIST_STORAGE* rhs, int8_t new_dtype);
size_t          count_storage_max_elements(const STORAGE*);

void*           list_storage_get(LIST_STORAGE* s, size_t* coords);
void*           list_storage_insert(LIST_STORAGE* s, size_t* coords, void* val);
void*           list_storage_remove(LIST_STORAGE* s, size_t* coords);
bool            list_storage_eqeq(const LIST_STORAGE*, const LIST_STORAGE*);

/* yale.c */
void print_vectors(YALE_STORAGE* s);
YALE_STORAGE*   create_yale_storage(int8_t dtype, size_t* shape, size_t rank, size_t init_capacity);
void            init_yale_storage(YALE_STORAGE* s);
void            delete_yale_storage(YALE_STORAGE* s);
void            mark_yale_storage(void* s);
YALE_STORAGE*   cast_copy_yale_storage(YALE_STORAGE* rhs, int8_t new_dtype);
bool            yale_storage_eqeq(const YALE_STORAGE*, const YALE_STORAGE*);

void*           yale_storage_ref(YALE_STORAGE* s, size_t* coords);
char            yale_storage_set(YALE_STORAGE* s, size_t* coords, void* v);

YALE_STORAGE*   create_merged_yale_storage(const YALE_STORAGE*, const YALE_STORAGE*);

size_t          count_list_storage_nd_elements(const LIST_STORAGE*);
size_t          count_list_storage_elements(const LIST_STORAGE*);

/* stype casts */
DENSE_STORAGE* scast_copy_dense_yale(const YALE_STORAGE* rhs, int8_t l_dtype);
DENSE_STORAGE* scast_copy_dense_list(const LIST_STORAGE* rhs, int8_t l_dtype);
YALE_STORAGE* scast_copy_yale_dense(const DENSE_STORAGE* rhs, int8_t l_dtype);
YALE_STORAGE* scast_copy_yale_list(const LIST_STORAGE* rhs, int8_t l_dtype);
LIST_STORAGE* scast_copy_list_yale(const YALE_STORAGE* rhs, int8_t l_dtype);
LIST_STORAGE* scast_copy_list_dense(const DENSE_STORAGE* rhs, int8_t l_dtype);

/* nmatrix.c */
void cast_copy_value_single(void* to, const void* from, int8_t l_dtype, int8_t r_dtype);
int8_t nm_dtypestring_to_dtype(VALUE str);
int8_t nm_dtypesymbol_to_dtype(VALUE sym);
int8_t nm_stypestring_to_stype(VALUE str);
int8_t nm_stypesymbol_to_stype(VALUE sym);
int8_t nm_guess_dtype(VALUE v);
size_t* nm_interpret_shape_arg(VALUE arg, size_t* rank);
NMATRIX* nm_create(int8_t stype, void* storage);
void Init_nmatrix();

#endif

ext: h


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

#include "tasks.h"
#include "utils.h"

void solve_task_par_reconfigure(int nth)
{
    // empty
}

void solve_task_par_finalize(void)
{
    // empty
}

void solve_task_par(void *ptr, int nth, int me)
{
    struct solve_task_arg *arg = (struct solve_task_arg*) ptr;

    int n     = arg->n;
    int nrhs  = arg->nrhs;
    double *L = arg->L;
    int ldL   = arg->ldL;
    double *X = arg->X;
    int ldX   = arg->ldX;

    // Compute nominal block size.
    int blksz = iceil(nrhs, nth);

    // Determine my share of the right-hand sides.
    int my_first_rhs, my_nrhs;
    my_first_rhs = blksz * me;
    my_nrhs = min(blksz, nrhs - my_first_rhs);

    // Solve L * X = B, where B and X share the same memory.
    cblas_dtrsm(CblasColMajor, CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit,
                n, my_nrhs,
                1.0, L, ldL,
                     X + my_first_rhs * ldX, ldX);
}

ext: c

max_stars_repo_name: NLAFET/pcp-runtime

max_issues_repo_name: NLAFET/pcp-runtime

max_forks_repo_name: NLAFET/pcp-runtime


text:
#include <viaio/VImage.h>
#include <viaio/Vlib.h>
#include <viaio/mu.h>
#include <viaio/option.h>

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>

int MeanSigma(double *data,size_t n,double *ave,double *sigma)
{
  size_t i=0;
  double u=0,nx=0,sum1=0,sum2=0,mean=0,var=0;

  for (i=0; i<n; i++) {
    u = data[i];
    sum1 += u;
    sum2 += u*u;
  }
  nx = (double)n;
  mean = sum1/nx;
  var = (sum2 - nx * mean * mean) / (nx - 1.0);
  if (var > TINY) {
    *ave = mean;
    *sigma = sqrt(var);
    return 1;
  }
  else {
    return -1;
  }
}

int NormalizeTimecourse(double *data,size_t n,double *ave,double *sigma)
{
  size_t i=0;
  double u=0,mean=0,sd=0;

  if (MeanSigma(data,n,&mean,&sd) < 0) return -1;

  for (i=0; i<n; i++) {
    u = data[i];
    data[i] = (u-mean)/sd;
  }
  *ave = mean;
  *sigma = sd;
  return 1;
}

void VFreqFilter(VAttrList list,VFloat highpass,VFloat lowpass,VBoolean stop)
{
  int i,j,k,b,r,c;
  double x=0,freq=0;

  /* read functional data */
  int nrows=0,ncols=0,nt=0;
  int nslices = VAttrListNumImages(list);
  VImage *src = VAttrListGetImages(list,nslices);
  VImageDimensions(src,nslices,&nt,&nrows,&ncols);

  /* read header info */
  double *D = (double *)VCalloc(8,sizeof(double));
  for (i=0; i<8; i++) D[i] = 1.0;
  VAttrList geolist = VGetGeoInfo(list);
  if (geolist != NULL) D = VGetGeoPixdim(geolist,D);
  double tr = D[4];

  VFloat xtr=-1;
  if (VGetAttr(VImageAttrList(src[0]),"repetition_time",NULL,VFloatRepn,&xtr) == VAttrFound) {
    tr = (double)xtr;
    if (tr > 100) tr /= 1000.0;  /* convert to seconds */
    if (tr < 0.01) VError("FreqFilter: implausible TR (%f seconds)",tr);
  }
  fprintf(stderr," image dimensions: %d x %d x %d,  nt= %d,  TR= %.3f secs\n",nslices,nrows,ncols,nt,tr);

  /* ini dest image */
  VImage *dst = (VImage *)VCalloc(nslices,sizeof(VImage));
  for (i=0; i<nslices; i++) {
    dst[i] = VCreateImage(nt,nrows,ncols,VFloatRepn);
    VFillImage(dst[i],VAllBands,0);
    VCopyImageAttrs (src[i], dst[i]);
  }
  double smin = VPixelMinValue(dst[0])+TINY;
  double smax = VPixelMaxValue(dst[0])-TINY;

  /* ini fft, zero-padding, add tails */
  size_t tail = 40;
  if (nt < tail) tail = 0;
  size_t ntt = nt + 2*tail;

  size_t nfft = 1;
  while (nfft < ntt) nfft *= 2;
  size_t nc  = (nfft / 2) + 1;
  double alpha = 1.0 / ((double)nfft * tr);

  double *data = (double *)VCalloc(nt,sizeof(double));
  double *data_fft = (double *)VCalloc(nfft,sizeof(double));

  /* bandpass filter function */
  double s=0,xs=1,xf=1;
  double *bandpass = (double *)VCalloc(nc,sizeof(double));
  for (i=1; i<nc; i++) {
    freq = (double)i * alpha;
    s = 1.0/freq;
    xs = xf = 1.0;
    if (lowpass > 0) xf = 1.0 / (1.0 + exp(lowpass-s));
    if (highpass > 0) xs = 1.0 / (1.0 + exp(s-highpass));
    bandpass[i] = xs*xf;
    if (stop == TRUE) bandpass[i] = 1.0-bandpass[i];
  }

  /* main process */
  double ave=0,sigma=0;
  for (b=0; b<nslices; b++) {
    for (r=0; r<nrows; r++) {
      for (c=0; c<ncols; c++) {

	/* get data */
	for (j=0; j<nt; j++) data[j] = VGetPixel(src[b],j,r,c);
	if (NormalizeTimecourse(data,nt,&ave,&sigma) < 0) {
	  for (j=0; j<nt; j++) VSetPixel(src[b],j,r,c,0.0);
	  continue;
	}

	/* zero padding and tail added */
	for (j=0; j<nfft; j++) data_fft[j] = 0;
	k=tail;
	for (j=0; j<tail; j++) {
	  data_fft[j] = data[k--];
	}
	k=0;
	for (j=tail; j<nt+tail; j++) {
	  data_fft[j] = data[k++];
	}
	k=nt;
	for (j=tail+nt; j<nt+2*tail; j++) {
	  data_fft[j] = data[k--];
	}

	/*  fft */
	gsl_fft_real_radix2_transform(data_fft,1,(size_t)nfft);

	/* apply bandpass filter */
	for (i=1; i<nc; i++) {
	  x = bandpass[i];
	  data_fft[i] *= x;
	  data_fft[nfft-i] *= x;
	}

	/* inverse fft */
	gsl_fft_halfcomplex_radix2_inverse(data_fft,1,(size_t)nfft);
	k=0;
	for (j=0; j<nt; j++) {
	  data[k++] = data_fft[j+tail];
	}

	for (j=0; j<nt; j++) {
	  x = data[j]*sigma + ave;
	  if (x < smin) x = smin+TINY;
	  if (x >= smax) x = smax-TINY;
	  VPixel(dst[b],j,r,c,VFloat) = x;
	}
      }
    }
  }
  VFree(data);
  VFree(data_fft);

  /* update attr list */
  VAttrListPosn posn;
  i=0;
  for (VFirstAttr (list, & posn); VAttrExists (& posn); VNextAttr (& posn)) {
    if (VGetAttrRepn (& posn) != VImageRepn) continue;
    VDestroyImage(src[i]);
    if (i >= nslices) break;
    VSetAttrValue (& posn, NULL,VImageRepn,dst[i]);
    i++;
  }
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#include <petsc.h>
/*
#include "libpetsc4py/libpetsc4py.h"
*/
PETSC_EXTERN int import_libpetsc4py(void);
PETSC_EXTERN PetscErrorCode MatPythonGetContext(Mat,void**);
PETSC_EXTERN PetscErrorCode MatPythonSetContext(Mat,void*);
PETSC_EXTERN PetscErrorCode PCPythonGetContext(PC,void**);
PETSC_EXTERN PetscErrorCode PCPythonSetContext(PC,void*);
PETSC_EXTERN PetscErrorCode KSPPythonGetContext(KSP,void**);
PETSC_EXTERN PetscErrorCode KSPPythonSetContext(KSP,void*);
PETSC_EXTERN PetscErrorCode SNESPythonGetContext(SNES,void**);
PETSC_EXTERN PetscErrorCode SNESPythonSetContext(SNES,void*);
PETSC_EXTERN PetscErrorCode TSPythonGetContext(TS,void**);
PETSC_EXTERN PetscErrorCode TSPythonSetContext(TS,void*);
PETSC_EXTERN PetscErrorCode PetscPythonRegisterAll(void);


text:
#include "multi_modelstruct.h"
#include "multivar_support.h"

#include <math.h>
#include <assert.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_errno.h>

/**
 * ccs, 05.06.2012
 * \todo this all needs to be checked carefully against an example that works in R
 * \todo valgrind for mem-leaks
 * \bug no error checking in linalg
 * \bug no error checking in alloc
 * \bug is there error checking in file handles?
 *
 * ccs, 19.06.2012
 * if you allocate a multi_modelstruct m
 * then m->pca_model_array[i]->xmodel == m->xmodel
 * the value in the pca_model_array is just a pointer to m, do we really want this?
 */

/**
 * allocates a multi_modelstruct, like alloc_modelstruct_2,
 * but for multivariate models with t output values at each location
 *
 * @param model_in: (n x d) matrix of the design
 * @param training_matrix: (n x t) matrix of the values of the training values at each of the n locations
 * @param cov_fn_index:  POWEREXPCOVFN, MATERN32, or MATERN52
 * @param regression_order:  0, 1, 2, or 3
 * @param varfrac: the minimum fractional variance that should be retained during the PCA decomp
 *
 * applies a pca decomp to training_matrix to reduce the dimensionality
 *
 */
multi_modelstruct* alloc_multimodelstruct(gsl_matrix *xmodel_in,
																				 gsl_matrix *training_matrix_in,
																				 int cov_fn_index,
																				 int regression_order, double varfrac)
{
	assert(training_matrix_in->size1 == xmodel_in->size1);
	assert(training_matrix_in->size1 > 0);
	assert(training_matrix_in->size2 > 0);
	assert(xmodel_in->size2 > 0);

	int i;
	double mean_temp = 0.0;
	int nt = training_matrix_in->size2;
	int nmodel_points = xmodel_in->size1;
	int nparams = xmodel_in->size2;
	gsl_vector_view col_view;

	/* use default if out of range */
	if (regression_order < 0 || regression_order > 3)
		regression_order = 0;

	/* use a sensible default for the variance fraction */
	if(varfrac < 0 || varfrac > 1)
		varfrac = 0.95;

	/* ntheta is a function of cov_fn_index and nparams */
	int nthetas;
	if ((cov_fn_index == MATERN32) || (cov_fn_index == MATERN52)) {
		nthetas = 3;
	} else if (cov_fn_index == POWEREXPCOVFN) {
		nthetas = nparams + 2;
	} else {
		cov_fn_index = POWEREXPCOVFN;
		nthetas = nparams + 2;
	}

	// this doesn't seem to be allocating correctly, seems to be a broken defn of MallocChecked
	// strangely, code will build in this case...
	//multi_modelstruct * model = (multi_modelstruct*)MallocChecked(sizeof(multi_modelstruct));
	multi_modelstruct * model = (multi_modelstruct*)malloc(sizeof(multi_modelstruct));

	// fill in
	model->nt = nt;
	model->nr = 0; // init at zero
	model->nmodel_points = nmodel_points;
	model->nparams = nparams;
	model->xmodel = xmodel_in;
	model->training_matrix = training_matrix_in;
	model->training_mean = gsl_vector_alloc(nt);
	model->regression_order = regression_order;
	model->cov_fn_index = cov_fn_index;

	/* fill in the mean vector, should probably sum this more carefully... */
	for(i = 0; i < nt; i++){
		col_view = gsl_matrix_column(model->training_matrix, i);
		mean_temp = vector_elt_sum(&col_view.vector, nmodel_points);
		//printf("%lf\n", (mean_temp/((double)nmodel_points)));
		gsl_vector_set(model->training_mean, i, (mean_temp/((double)nmodel_points)) );
	}

	/* carry out the pca decomp on this model, this is defined in multivar_support for now
	 * this will fill in nr, pca_eigenvalues, pca_eigenvectors, pca_evals_r, pca_evecs_r
	 *
	 * this is making a mess if nt = 1
	 */
	gen_pca_decomp(model, varfrac);

	/* fill in pca_model_array */
	gen_pca_model_array(model);

	return model;
}

/**
 * fill in pca_model_array
 * requires:
 * - the pca decomp to have been calculated, so nr and all the pca_... fields are allocated and filled in
 * - m to have been allocated up to: nt, nmodel_points, and xmodel
 *
 * is it possible we could start running out of memory doing all these allocs?
 */
void gen_pca_model_array(multi_modelstruct *m)
{
	int nr = m->nr;
	int i;
	gsl_vector_view col_view;
	gsl_vector* temp_train_vector = gsl_vector_alloc(m->nmodel_points);
	//gsl_matrix* temp_xmodel = gsl_matrix_alloc(m->nmodel_points, m->nparams);

	//gsl_matrix_mempcy(temp_xmodel, m->xmodel);

	// alloc the array of nr model structs
	//m->pca_model_array = (modelstruct**)MallocChecked(sizeof(modelstruct*)*nr);
	m->pca_model_array = (modelstruct**)malloc(sizeof(modelstruct*)*nr);
	// fill in the modelstructs correctly
	for(i = 0; i < nr; i++){
		col_view = gsl_matrix_column(m->pca_zmatrix, i);
		gsl_vector_memcpy(temp_train_vector, &(col_view.vector));

		// this isn't copying in the training vector correctly for somereason
		m->pca_model_array[i] = alloc_modelstruct_2(m->xmodel, temp_train_vector,
																								m->cov_fn_index, m->regression_order);

		// see if brute forcing it will work
		m->pca_model_array[i]->training_vector = gsl_vector_alloc(m->nmodel_points);
		gsl_vector_memcpy(m->pca_model_array[i]->training_vector, temp_train_vector);
	}
	gsl_vector_free(temp_train_vector);
}

/**
 * carries out a pca decomp on m->training_matrix;
 * setting m->nr, m->pca_eigenvalues, m->pca_eigenvectors and
 * initializing m->pca_model_array
 *
 * the pca decomp is pretty simple:
 * let Y_ij = m->y_training
 * let mu_i = (1/nmodel_points) sum_{j=1}^{nmodel_points} Y_{ij}  // the sample mean
 * let ysub_i = Y_i - rep(mu_i,nmodel_points)  //subtract the sample means from each column
 * let sigma_ij = (1/nmodel_points) ( ysub * t(ysub)) // this is the covariance matrix
 *
 * then all we need to do is compute the eigendecomposition
 *
 * sigma_ij = U^{-1} Lambda U
 *
 * where Lambda is a diagonal matrix of t eigenvalues and U is a t x t matrix with the eigenvectors as columns
 *
 * requires:
 * - 0 < vfrac < 1 (0.95 is a good value)
 * - m to have allocated and filled out, nt, nmodel_points, trianing_matrix, (doesn't need xmodel explicitly)
 */
void gen_pca_decomp(multi_modelstruct *m, double vfrac)
{
	FILE *fptr; // for debug output
	int i,j;
	int nt = m->nt;
	int retval;
	double total_variance = 0.0, frac  = 0.0;
	gsl_matrix *pca_zmatrix_temp;
	gsl_matrix *y_sub_mat = gsl_matrix_alloc(m->nmodel_points, nt);
	gsl_matrix *y_temp_mat = gsl_matrix_alloc(m->nmodel_points, nt);

	gsl_matrix *y_cov_mat = gsl_matrix_alloc(nt, nt);

	gsl_vector *evals_temp = gsl_vector_alloc(nt);
	gsl_matrix *evecs_temp = gsl_matrix_alloc(nt,nt);

	gsl_eigen_symmv_workspace *ework = gsl_eigen_symmv_alloc(nt);

	gsl_vector_view col_view;

	// why is this here?
	gsl_matrix_memcpy(y_sub_mat, m->training_matrix);

	// subtract out the mean
	for(i = 0; i < nt; i++){
		//col_view = gsl_matrix_column(y_sub_mat, i);
		printf("# y(%d) mean: %lf\n", i, gsl_vector_get(m->training_mean, i));
		for(j = 0; j < m->nmodel_points; j++){
			gsl_matrix_set(y_sub_mat, j, i, gsl_matrix_get(y_sub_mat, j, i) - gsl_vector_get(m->training_mean, i));
		}
	}

	// compute the sample-variance, by multiplying y_sub_mat, with itself transposed
	gsl_matrix_memcpy(y_temp_mat, y_sub_mat);
	gsl_matrix_set_zero(y_cov_mat);

	/** — Function: int gsl_blas_dgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB, double alpha, const gsl_matrix * A, const gsl_matrix * B, double beta, gsl_matrix * C)  (always forget this one)*/
	/* want C (nt x nt ) so we need to do: (nt x nmodel_points) * (nmodel_points x nt) */
	/**
	 * this is strange, the upper triangle is not obviously correct
	 */
	//retval = gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, y_temp_mat, y_sub_mat, 0.0, y_cov_mat);
	retval = gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, y_sub_mat, y_temp_mat, 0.0, y_cov_mat);
	if(retval){
		printf("# gen_pca_decomp:gsl_blas_dgemm %s\n", gsl_strerror(retval));
		exit(EXIT_FAILURE);
	}

	gsl_matrix_scale(y_cov_mat, (1.0/((double)m->nmodel_points)));

	/**
	 * debug output
	 */
	#ifdef DEBUGPCA
	fptr = fopen("pca-debug.dat","w");
	fprintf(fptr, "# ycov:\n");
	for(j = 0; j < nt; j++){
		for(i = 0; i < nt; i++)
			fprintf(fptr, "%lf ", gsl_matrix_get(y_cov_mat, i, j));
		fprintf(fptr, "\n");
	}
	#endif

	/** now the eigendecomp
	 * y_cov_mat is symmetric and better be real so we can use gsl_eigen_symmv,
	 * note that the ?upper? triangle of y_cov_mat is borked during this process
	 * also: the evals are not sorted by order, but the evectors are in the order of the evalues.
	 * so we need to sort the evalues and the evectors correctly before we can use them
	 * using
	 * — Function: int gsl_eigen_symmv_sort (gsl_vector * eval, gsl_matrix * evec, gsl_eigen_sort_t sort_type)
	 * sort them into descending order (biggest first)
	 */
	gsl_eigen_symmv(y_cov_mat, evals_temp, evecs_temp, ework);
	gsl_eigen_symmv_sort(evals_temp, evecs_temp, GSL_EIGEN_SORT_VAL_DESC);
	/**
	 * eigenvectors are stored in columns of pca_evecs*/
	total_variance = vector_elt_sum(evals_temp, nt);

#ifdef DEBUGPCA
	fprintf(fptr, "# evals:\n");
	for(i = 0; i < nt; i++)
		fprintf(fptr, "%lf %lf\n", gsl_vector_get(evals_temp, i), gsl_vector_get(evals_temp, i) / total_variance);

	fprintf(fptr, "# evecs:\n");
	for(j = 0; j < nt; j++){
		for(i = 0; i < nt; i++)
			fprintf(fptr, "%lf ", gsl_matrix_get(evecs_temp, i, j));
		fprintf(fptr, "\n");
	}

	#endif

	i=0;
	while( frac < vfrac && (i+1) < nt){
		frac = (1.0/total_variance) * vector_elt_sum(evals_temp, i);
		i++;
	}
	m->nr = i;

	if(nt == 1){
		//printf("# 1d case, nr=1\n");
		m->nr = 1;
	}

	m->pca_evals_r = gsl_vector_alloc(m->nr);
	m->pca_evecs_r = gsl_matrix_alloc(m->nt, m->nr);
	// debug...
	fprintf(stderr, "# nr: %d frac: %lf\n", m->nr, frac);

	for(i = 0; i < m->nr; i++){
		gsl_vector_set(m->pca_evals_r, i, gsl_vector_get(evals_temp, i));
		col_view = gsl_matrix_column(evecs_temp, i);
		gsl_matrix_set_col(m->pca_evecs_r, i, &col_view.vector);
	}

	// fill in pca_zmatrix
	m->pca_zmatrix = gsl_matrix_alloc(m->nmodel_points, m->nr);
	pca_zmatrix_temp = gsl_matrix_alloc(m->nmodel_points, m->nr);
	// zmat: (nmodel_points x nr) = (nmodel_points x nt) * ( nt x nr )
	/** — Function: int gsl_blas_dgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB, double alpha, const gsl_matrix * A, const gsl_matrix * B, double beta, gsl_matrix * C)  (always forget this one)*/
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0,  y_sub_mat, m->pca_evecs_r, 0.0, m->pca_zmatrix);

	#ifdef DEBUGPCA
	fprintf(fptr, "# evecs R:\n");
	for(i = 0; i < m->nt; i++){
		for(j = 0; j < m->nr; j++)
			fprintf(fptr, "%lf ", gsl_matrix_get(m->pca_evecs_r, i, j));
		fprintf(fptr, "\n");
	}
	#endif

	gsl_matrix_free(y_temp_mat);
	y_temp_mat = gsl_matrix_alloc(m->nr, m->nr);
	for(i = 0; i < m->nr; i++) // scale the diagonal by the evalue */
		gsl_matrix_set(y_temp_mat, i, i, 1.0/(sqrt(gsl_vector_get(m->pca_evals_r, i))));

	//print_matrix(y_temp_mat, m->nr, m->nr);

	// if nr != nt this won't work!
	//gsl_matrix_memcpy(y_sub_mat, m->pca_zmatrix);
	gsl_matrix_memcpy(pca_zmatrix_temp, m->pca_zmatrix);
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, pca_zmatrix_temp, y_temp_mat, 0.0, m->pca_zmatrix);

	#ifdef DEBUGPCA
	fprintf(fptr, "# zmat:\n");
	for(i = 0; i < 5; i++){
		for(j = 0; j < m->nr; j++)
			fprintf(fptr, "%lf ", gsl_matrix_get(m->pca_zmatrix, i, j));
		fprintf(fptr, "\n");
	}
	fclose(fptr);
	#endif

	gsl_vector_free(evals_temp);
	gsl_matrix_free(evecs_temp);
	gsl_matrix_free(pca_zmatrix_temp);

	gsl_eigen_symmv_free(ework);
	gsl_matrix_free(y_sub_mat);
	gsl_matrix_free(y_temp_mat);
	gsl_matrix_free(y_cov_mat);
}

/**
 * dump the multimodelstruct to fptr, follows from dump_modelstruct_2
 * we first dump the new fields and then iterate through the nr additional models which are dumped as before
 * so we dump a lot of the same info, but this is probably ok, the advantage is that each section defining
 * a model can be pulled out and worked on separately...
 */
void dump_multi_modelstruct(FILE* fptr, multi_modelstruct *m){
	assert(fptr);

	int i,j;
	int nt = m->nt;
	int nr = m->nr;
	int nparams = m->nparams;
	int nmodel_points = m->nmodel_points;
	int cov_fn_index = m->cov_fn_index;
	int regression_order = m->regression_order;

	fprintf(fptr, "%d\n", nt);
	fprintf(fptr, "%d\n", nr);
	fprintf(fptr, "%d\n", nparams);
	fprintf(fptr, "%d\n", nmodel_points);
	fprintf(fptr, "%d\n", cov_fn_index);
	fprintf(fptr, "%d\n", regression_order);

	// multimodel thetas are inside pca_model_array...

	for(i = 0; i < nmodel_points; i++){
		for(j = 0; j < nparams; j++)
			fprintf(fptr, "%.17lf ", gsl_matrix_get(m->xmodel, i, j));
		fprintf(fptr, "\n");
	}

	for(i = 0; i < nmodel_points; i++){
		for(j = 0; j < nt; j++)
			fprintf(fptr, "%.17lf ", gsl_matrix_get(m->training_matrix, i, j));
		fprintf(fptr, "\n");
	}

	// now the rest of the pca information
	for(i = 0; i < nr; i++)
			fprintf(fptr, "%.17lf ", gsl_vector_get(m->pca_evals_r, i));
		fprintf(fptr, "\n");

	for(i = 0; i < nt; i++){
		for(j = 0; j < nr; j++)
			fprintf(fptr, "%.17lf ", gsl_matrix_get(m->pca_evecs_r, i, j));
		fprintf(fptr, "\n");
	}

	for(i = 0; i < nmodel_points; i++){
		for(j = 0; j < nr; j++)
			fprintf(fptr, "%.17lf ", gsl_matrix_get(m->pca_zmatrix, i, j));
		fprintf(fptr, "\n");
	}

	// now the pca_model_array
	// these could go in separate files...
	for(i = 0; i < nr; i++){
		dump_modelstruct_2(fptr, m->pca_model_array[i]);
	}

}

/**
 * loads a multivariate modelstructure from fptr
 */
multi_modelstruct *load_multi_modelstruct(FILE* fptr){
	multi_modelstruct *m = (multi_modelstruct*)malloc(sizeof(multi_modelstruct));
	int i,j;
	int nt, nr;
	int nparams, nmodel_points;
	int cov_fn_index;
	int regression_order;

	double mean_temp;
	gsl_vector_view col_view;

	fscanf(fptr, "%d%*c", & nt);
	fscanf(fptr, "%d%*c", & nr);
	fscanf(fptr, "%d%*c", & nparams);
	fscanf(fptr, "%d%*c", & nmodel_points);
	fscanf(fptr, "%d%*c", & cov_fn_index);
	fscanf(fptr, "%d%*c", & regression_order);

	m->nt = nt;
	m->nr = nr;
	m->nparams = nparams;
	m->nmodel_points = nmodel_points;
	m->cov_fn_index = cov_fn_index;
	m->regression_order = regression_order;

	// now we can allocate everything in m
	m->xmodel = gsl_matrix_alloc(nmodel_points, nparams);
	m->training_matrix = gsl_matrix_alloc(nmodel_points, nt);
	m->training_mean = gsl_vector_alloc(nt); // do we need this? (yes!)
	m->pca_model_array = (modelstruct**)malloc(sizeof(modelstruct*)*nr);
	m->pca_evals_r = gsl_vector_alloc(nr);
	m->pca_evecs_r = gsl_matrix_alloc(nt, nr);
	m->pca_zmatrix = gsl_matrix_alloc(nmodel_points, nr);

	for(i = 0; i < nmodel_points; i++)
		for(j = 0; j < nparams; j++)
			fscanf(fptr, "%lf%*c", gsl_matrix_ptr(m->xmodel, i, j));

	for(i = 0; i < nmodel_points; i++)
		for(j = 0; j < nt; j++)
			fscanf(fptr, "%lf%*c", gsl_matrix_ptr(m->training_matrix, i, j));

	// now the rest of the pca information
	for(i = 0; i < nr; i++)
			fscanf(fptr, "%lf%*c", gsl_vector_ptr(m->pca_evals_r, i));

	for(i = 0; i < nt; i++)
		for(j = 0; j < nr; j++)
			fscanf(fptr, "%lf%*c", gsl_matrix_ptr(m->pca_evecs_r, i, j));

	for(i = 0; i < nmodel_points; i++)
		for(j = 0; j < nr; j++)
			fscanf(fptr, "%lf%*c", gsl_matrix_ptr(m->pca_zmatrix, i, j));

	for(i = 0; i < nr; i++)
		m->pca_model_array[i] = load_modelstruct_2(fptr);

	/* fill in the mean vector */
	for(i = 0; i < nt; i++){
		col_view = gsl_matrix_column(m->training_matrix, i);
		mean_temp = vector_elt_sum(&col_view.vector, nmodel_points);
		gsl_vector_set(m->training_mean, i, (mean_temp/((double)nmodel_points)) );
	}

	return m;
}

/**
 * return the sum of the elements of vec from 0:nstop
 */
double vector_elt_sum(gsl_vector* vec, int nstop)
{
	assert(nstop >= 0);
	assert((unsigned)nstop <= vec->size);
	int i;
	double sum = 0.0;
	for(i = 0; i < nstop; i++){
		sum += gsl_vector_get(vec, i);
	}
	return(sum);
}

/**
 * this free's everything in m
 */
void free_multimodelstruct(multi_modelstruct *m)
{
	int i;
	gsl_vector_free(m->training_mean);
	for(i = 0; i < m->nr; i++){
		free_modelstruct_2(m->pca_model_array[i]);
		//gsl_matrix_free(m->pca_model_array[i]->xmodel);
	}
	free(m->pca_model_array);
	gsl_matrix_free(m->xmodel);
	gsl_matrix_free(m->training_matrix);
	gsl_vector_free(m->pca_evals_r);
	gsl_matrix_free(m->pca_evecs_r);
	gsl_matrix_free(m->pca_zmatrix);
}

max_stars_repo_name: jackdawjackdaw/emulator

max_issues_repo_name: MADAI/MADAIEmulator

max_forks_repo_name: MADAI/MADAIEmulator


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

#include <architecture/ppc/fp_regs.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_errno.h>

int
gsl_ieee_set_mode (int precision, int rounding, int exception_mask)
{
  ppc_fp_scr_t fp_scr = get_fp_scr() ;

  switch (precision)
    {
    case GSL_IEEE_SINGLE_PRECISION:
      GSL_ERROR ("powerpc only supports default precision rounding", GSL_EUNSUP);
      break ;
    case GSL_IEEE_DOUBLE_PRECISION:
      GSL_ERROR ("powerpc only supports default precision rounding", GSL_EUNSUP);
      break ;
    case GSL_IEEE_EXTENDED_PRECISION:
      GSL_ERROR ("powerpc only supports default precision rounding", GSL_EUNSUP);
      break ;
    }

  switch (rounding)
    {
    case GSL_IEEE_ROUND_TO_NEAREST:
      fp_scr.rn = RN_NEAREST ;
      break ;
    case GSL_IEEE_ROUND_DOWN:
      fp_scr.rn = RN_TOWARD_MINUS ;
      break ;
    case GSL_IEEE_ROUND_UP:
      fp_scr.rn = RN_TOWARD_PLUS ;
      break ;
    case GSL_IEEE_ROUND_TO_ZERO:
      fp_scr.rn = RN_TOWARD_ZERO ;
      break ;
    default:
      fp_scr.rn = RN_NEAREST ;
    }

  if (exception_mask & GSL_IEEE_MASK_INVALID)
    fp_scr.ve = 0 ;                             //ve bit:  invalid op exception enable

  if (exception_mask & GSL_IEEE_MASK_DENORMALIZED)
    {
      /* do nothing */
    }
  else
    {
      GSL_ERROR ("powerpc does not support the denormalized operand exception. "
                 "Use 'mask-denormalized' to work around this.", GSL_EUNSUP) ;
    }

  if (exception_mask & GSL_IEEE_MASK_DIVISION_BY_ZERO)
    fp_scr.ze = 0 ;                             //ze bit:  zero divide exception enable

  if (exception_mask & GSL_IEEE_MASK_OVERFLOW)
    fp_scr.oe = 0 ;                             //oe bit:  overflow exception enable

  if (exception_mask & GSL_IEEE_MASK_UNDERFLOW)
    fp_scr.ue  = 0 ;                            //ue bit:  underflow exception enable

  if (exception_mask & GSL_IEEE_TRAP_INEXACT)
    {
      fp_scr.xe = 1 ;                           //xe bit:  inexact exception enable
    }
  else
    {
      fp_scr.xe = 01 ;
    }

  set_fp_scr(fp_scr);

  return GSL_SUCCESS ;

}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
//
// Created by karl_ on 2020-5-3.
//

#include <cblas.h>
#include <imgproc/color-convert.h>

gboolean imgproc_to_grayscale(CoreImage *src, CoreImage **dst) {
    gsize i;
    CoreSize *size;
    gdouble *dst_data_cast, *src_data_cast;
    gpointer src_data, dst_data;
    gsize channel, area;
    CoreColorSpace color_space;
    CorePixelType pixel_type;
    g_return_val_if_fail(src != NULL, FALSE);
    g_return_val_if_fail(dst != NULL, FALSE);

    size = core_image_get_size(src);
    area = core_size_get_area(size);
    color_space = core_image_get_color_space(src);
    pixel_type = core_image_get_pixel_type(src);
    channel = core_image_get_channel(src);
    src_data = core_image_get_data(src);

    if (core_pixel_is_uint8(pixel_type)) {
        src_data_cast = g_malloc(area * sizeof(gdouble) * channel);
        for (i = 0; i < area * channel; ++i) {
            src_data_cast[i] = (gdouble) ((guint8 *) src_data)[i];
        }
    } else {
        /* do not need to cast if pixel type is Dx */
        src_data_cast = src_data;
    }
    dst_data_cast = g_malloc(area * sizeof(gdouble));

    cblas_dcopy(area, src_data_cast, channel, dst_data_cast, 1);
    for (i = 1; i < channel; ++i) {
        cblas_daxpy(area, 1.0, src_data_cast + i, channel, dst_data_cast, 1);
    }
    if (channel > 1) {
        cblas_daxpy(area, 1.0 / channel - 1.0, dst_data_cast, 1, dst_data_cast, 1);
    }

    if (core_pixel_is_uint8(pixel_type)) {
        pixel_type = CORE_PIXEL_U1;
    } else if (core_pixel_is_double(pixel_type)) {
        pixel_type = CORE_PIXEL_D1;
    } else {
        g_return_val_if_fail(FALSE, FALSE);
    }

    /* cast dst_data from double back to uchar if needed */
    if (core_pixel_is_uint8(pixel_type)) {
        for (i = 0; i < area; ++i) {
            ((guint8 *) dst_data_cast)[i] = dst_data_cast[i];
        }
    }
    dst_data = dst_data_cast;
    if (core_pixel_is_uint8(pixel_type)) {
        dst_data = g_realloc(dst_data, sizeof(guint8) * area);
    }

    /* assign result */
    if (*dst == NULL) {
        *dst = core_image_new_with_data(dst_data, CORE_COLOR_SPACE_GRAY_SCALE, pixel_type, size, FALSE);
    } else {
        core_image_assign_data(*dst, dst_data, CORE_COLOR_SPACE_GRAY_SCALE, pixel_type, size, FALSE);
    }

    /* free additional allocated memory */
    if (core_pixel_is_uint8(pixel_type)) {
        g_free(src_data_cast);
    }
    g_object_unref(size);
    return TRUE;
}

gboolean imgproc_to_binary_threshold(CoreImage *src, CoreImage **dst, gdouble threshold, gboolean inverse) {
    CoreImage *gray_src = NULL;
    CoreSize *size;
    gsize i, area;
    gdouble *src_data, *dst_data;
    CoreColorSpace color_space;
    CorePixelType pixel_type;
    g_return_val_if_fail(src != NULL, FALSE);
    g_return_val_if_fail(dst != NULL, FALSE);

    if (core_image_get_channel(src) > 1) {
        imgproc_to_grayscale(src, &gray_src);
    } else {
        gray_src = g_object_ref(src);
    }
    size = core_image_get_size(gray_src);
    area = core_size_get_area(size);
    src_data = core_image_get_data(gray_src);
    dst_data = g_malloc(sizeof(gdouble) * area);
    for (i = 0; i < area; ++i) {
        dst_data[i] = ((unsigned) (src_data[i] > threshold) ^ (unsigned) inverse) * 255.0;
    }
    if (*dst == NULL) {
        core_image_new_with_data(dst_data, 1, CORE_PIXEL_D4, size, FALSE);
    } else {
        core_image_assign_data(*dst, dst_data, 1, CORE_PIXEL_D4, size, FALSE);
    }
    g_object_unref(size);
    g_object_unref(gray_src);
    return TRUE;
}

gboolean imgproc_to_binary_dither() { return FALSE; }

gboolean imgproc_to_binary_ordered_dither() { return FALSE; }

max_stars_repo_name: gtrong/EasyPhotoshop

max_issues_repo_name: gtrong/EasyPhotoshop

max_forks_repo_name: gtrong/EasyPhotoshop