text:
/**
 * @file   libhades.c
 * @author Michael Hartmann <michael.hartmann@physik.uni-augsburg.de>
 * @date   February, 2016
 * @brief  library to access low-level LAPACK functions
 */

#include <cblas.h>
#include <errno.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <strings.h>

#include <libhades.h>
#include <libhades/parse_npy_dict.h>

/** \defgroup misc miscellaneous functions
 *  @{
 */

/** @brief malloc wrapper
 *
 * This function uses malloc to allocate size bytes of memory and returns a
 * pointer to the memory. If an error occures, an error message is printed to
 * stderr and the program is aborted.
 *
 * @param [in] size amount of memory to allocate
 * @retval ptr pointer to the allocated memory
 */
void *libhades_malloc(size_t size)
{
    void *ptr = malloc(size);
    if(ptr == NULL)
    {
        const int err = errno;
        fprintf(stderr, "malloc can't allocate %zu bytes of memory: %s (%d)\n", size, strerror(err), err);
        abort();
    }
    return ptr;
}

/** @brief free wrapper
 *
 * This function frees the memory allocated by \ref libhades_malloc or \ref
 * libhades_realloc (or malloc and realloc). If the pointer given is NULL, an error is
 * printed to stderr and the program is aborted.
 *
 * @param [in] ptr pointer to the memory that should be freed
 */
void libhades_free(void *ptr)
{
    if(ptr == NULL)
    {
        fprintf(stderr, "Trying to free NULL pointer\n");
        abort();
    }

    free(ptr);
}

/** @brief realloc wrapper
 *
 * This function changes the size of the memory block pointed to by ptr to size
 * bytes. If ptr is NULL, this function behaves like \ref libhades_malloc. If
 * an error occures, an error is printed to stderr and the program is aboprted.
 *
 * @param [in] ptr pointer to the memory block
 * @param [in] size size of the memory block
 * @retval ptr_new pointer to the new memory block
 */
void *libhades_realloc(void *ptr, size_t size)
{
    void *ptr_new = realloc(ptr, size);

    if(ptr_new == NULL)
    {
        int err = errno;
        fprintf(stderr, "realloc can't allocate %zu bytes of memory: %s (%d)\n", size, strerror(err), err);
        abort();
    }

    return ptr_new;
}

static void *(*malloc_cb)(size_t)          = &libhades_malloc;
static void *(*realloc_cb)(void *, size_t) = &libhades_realloc;
static void  (*free_cb)(void *)            = &libhades_free;

/** macro to create functions argmin, argmax, argabsmin, argabsmax */
#define ARGXXX(FUNCTION_NAME, FUNCTION, RELATION) \
size_t FUNCTION_NAME(double list[], size_t size) \
{ \
    size_t index = 0; \
    for(size_t i = 0; i < size; i++) \
        if(FUNCTION(list[i]) RELATION FUNCTION(list[index])) \
            index = i; \
    return index; \
}

/** @brief Return index of smallest element in list
 *
 * @param [in] list
 * @param [in] size elements in list
 *
 * @retval index
 */
ARGXXX(argmin, +, <)

/** @brief Return index of largest element in list
 *
 * @param [in] list
 * @param [in] size elements in list
 *
 * @retval index
 */
ARGXXX(argmax, +, >)

/** @brief Return index of element with smallest absolute value in list
 *
 * @param [in] list
 * @param [in] size elements in list
 *
 * @retval index
 */
ARGXXX(argabsmin, fabs, <)
/** @brief Return index of element with largest absolute value in list
 *
 * @param [in] list
 * @param [in] size elements in list
 *
 * @retval index
 */
ARGXXX(argabsmax, fabs, >)

/** @}*/

/** \defgroup create Creating, printing and freeing matrices
 *  @{
 */

static inline void _swap_int(int *a, int *b);
static inline void _swap_size_t(size_t *a, size_t *b);

static inline void _swap_int(int *a, int *b)
{
    int c = *a;
    *a = *b;
    *b = c;
}

static inline void _swap_size_t(size_t *a, size_t *b)
{
    size_t c = *a;
    *a = *b;
    *b = c;
}

#define MATRIX_SWAP(FUNCTION_NAME, MATRIX_TYPE, TYPE) \
void FUNCTION_NAME(MATRIX_TYPE *A, MATRIX_TYPE *B) \
{ \
    /* swap pointers */ \
    TYPE *ptr = A->M; \
    A->M = B->M; \
    B->M = ptr; \
\
    /* swap rows */ \
    _swap_int(&A->rows, &B->rows); \
\
    /* swap columns */ \
    _swap_int(&A->columns, &B->columns); \
\
    /* swap min */ \
    _swap_int(&A->min, &B->min); \
\
    /* swap size */ \
    _swap_size_t(&A->size, &B->size); \
\
    /* swap type */ \
    _swap_int(&A->type, &B->type); \
\
    /* swap min */ \
    _swap_int(&A->min, &B->min); \
}

/** @brief Swap matrices A and B
 *
 * This function swaps the matrices A and B. The former content of A will be
 * the content of B and vice versa. No data is copied or moved, but the
 * pointers are swapped.
 *
 * @param [in,out] A matrix A
 * @param [in,out] B matrix B
 */
MATRIX_SWAP(matrix_swap, matrix_t, double);

/** @brief Swap matrices A and B
 *
 * See \ref matrix_swap.
 *
 * @param [in,out] A matrix A
 * @param [in,out] B matrix B
 */
MATRIX_SWAP(matrix_complex_swap, matrix_complex_t, complex_t);

/** macro to create a diagnal matrix out of a vector */
#define MATRIX_DIAG(FUNCTION_NAME, MATRIX_TYPE, ZEROS) \
MATRIX_TYPE *FUNCTION_NAME(MATRIX_TYPE *v) \
{ \
    int dim = v->size; \
    MATRIX_TYPE *A = ZEROS(dim, dim, NULL); \
    if(A == NULL) \
        return NULL; \
\
    for(int i = 0; i < dim; i++) \
        matrix_set(A, i,i, v->M[i]); \
\
    return A; \
}

/** @brief Construct a real diagonal matrix from a row or columns vector
 *
 * @param [in] v row or column vector
 *
 * @retval A A = diag(v) if successfull, NULL otherwise
 */
MATRIX_DIAG(matrix_diag, matrix_t, matrix_zeros)

/** @brief Construct a complex diagonal matrix from a row or columns vector
 *
 * @param [in] v row or column vector
 *
 * @retval A A = diag(v) if successfull, NULL otherwise
 */
MATRIX_DIAG(matrix_complex_diag, matrix_complex_t, matrix_complex_zeros)

/** @brief Set functions to allocate and free memory
 *
 * By default wrappers to malloc and free from <stdlib.h> are used to allocate
 * and free memory. If allocation of memory fails or a NULL pointer is freed,
 * the program will terminate.
 *
 * @param [in] _malloc_cb callback to a malloc-alike function
 * @param [in] _free_cb callback to a free-alike function
 */
void matrix_set_alloc(void *(*_malloc_cb)(size_t), void  (*_free_cb)(void *))
{
    malloc_cb = _malloc_cb;
    free_cb   = _free_cb;
}

/** macro for copying matrices. */
#define MATRIX_COPY(FUNCTION_NAME, MTYPE, TYPE, ALLOC) \
MTYPE *FUNCTION_NAME(MTYPE *A, MTYPE *C) \
{ \
    if(C == NULL) { \
        C = ALLOC(A->rows, A->columns); \
        if(C == NULL) \
            return NULL; \
    } \
\
    C->rows    = A->rows; \
    C->columns = A->columns; \
    C->min     = A->min; \
    C->size    = A->size; \
    C->type    = A->type; \
    C->view    = 0; \
    memcpy(C->M, A->M, C->size*sizeof(TYPE)); \
    return C; \
}

/** @brief Copy real matrix A
 *
 * Copy matrix A into C. If C is NULL, space for the matrix C will be
 * allocated.
 *
 * @param [in]     A real matrix
 * @param [in,out] C real matrix
 *
 * @retval C copy of A
 */
MATRIX_COPY(matrix_copy, matrix_t, double, matrix_alloc)

/** @brief Copy complex matrix A
 *
 * Copy matrix A into C. If C is NULL, space for the matrix C will be
 * allocated.
 *
 * @param [in]     A complex matrix
 * @param [in,out] C complex matrix
 *
 * @retval C copy of A
 */
MATRIX_COPY(matrix_complex_copy, matrix_complex_t, complex_t, matrix_complex_alloc)

/** @brief Copy a real matrix A to a complex matrix C
 *
 * Copy the matrix A to a complex matrix C. The matrix elements of A and C will
 * be identical. If C is NULL, memory for the matrix will be allocated.
 *
 * @param [in] A real matrix
 *
 * @retval C copy of A if successfull, NULL otherwise
 */
matrix_complex_t *matrix_tocomplex(matrix_t *A, matrix_complex_t *C)
{
    if(C == NULL)
    {
        C = matrix_complex_alloc(A->rows, A->columns);
        if(C == NULL)
            return NULL;
    }

    for(size_t i = 0; i < C->size; i++)
        C->M[i] = A->M[i];

    return C;
}

/** @brief Print real matrix M to stream
 *
 * Print the matrix A to the stream given by stream, e.g. stdout or stderr.
 * The format is given by the format string format, the separator of two
 * columns is sep, the separator between lines is given by sep_line. Both sep
 * and sep_line may be NULL.
 *
 * @param [in] stream output stream
 * @param [in] A real matrix
 * @param [in] format output format, e.g. "%lf" or "%g"
 * @param [in] sep separator between columns, e.g. "\t"
 * @param [in] sep_line separator between lines, e.g. "\n"
 */
void matrix_fprintf(FILE *stream, matrix_t *A, const char *format, const char *sep, const char *sep_line)
{
    const int rows = A->rows, columns = A->columns;

    for(int i = 0; i < rows; i++)
    {
        for(int j = 0; j < columns; j++)
        {
            fprintf(stream, format, matrix_get(A, i,j));
            if(sep != NULL)
                fputs(sep, stream);
        }
        if(sep_line != NULL)
            fputs(sep_line, stream);
    }
}

/** @brief Print complex matrix A to stream
 *
 * See matrix_fprintf.
 *
 * @param [in] stream output stream
 * @param [in] A complex matrix
 * @param [in] format output format, e.g. "%+lf%+lfi"
 * @param [in] sep separator between columns, e.g. "\t"
 * @param [in] sep_line separator between lines, e.g. "\n"
 */
void matrix_complex_fprintf(FILE *stream, matrix_complex_t *A, const char *format, const char *sep, const char *sep_line)
{
    const int rows = A->rows, columns = A->columns;

    for(int i = 0; i < rows; i++)
    {
        for(int j = 0; j < columns; j++)
        {
            const complex_t c = matrix_get(A, i,j);
            fprintf(stream, format, CREAL(c), CIMAG(c));
            if(sep != NULL)
                fputs(sep, stream);
        }
        if(sep_line != NULL)
            fputs(sep_line, stream);
    }
}

/** macro for matrix allocations. */
#define MATRIX_ALLOC(FUNCTION_NAME, MTYPE, TYPE) \
MTYPE *FUNCTION_NAME(int rows, int columns) \
{ \
    if(rows <= 0 || columns <= 0) \
        return NULL; \
\
    MTYPE *A = malloc_cb(sizeof(MTYPE)); \
    if(A == NULL) \
        return NULL; \
\
    const size_t size = (size_t)rows*(size_t)columns; \
\
    A->rows    = rows; \
    A->columns = columns; \
    A->min     = MIN(rows, columns); \
    A->size    = size; \
    A->type    = 0; \
    A->view    = 0; \
    A->M       = malloc_cb(size*sizeof(TYPE)); \
    if(A->M == NULL) \
    { \
        free_cb(A); \
        return NULL; \
    } \
\
    return A; \
}

/** @brief Allocate memory for real matrix of m rows and n columns
 *
 * This function will allocate and return a matrix of m lines and n columns.
 * The function given by matrix_set_alloc will be used to allocate memory; by
 * default this is malloc from <stdlib.h>.
 *
 * The matrix elements will be undefined. To allocate a matrix initialized with
 * zeros, see matrix_zeros. To create a unity matrix, see matrix_eye.
 *
 * @param [in] rows rows of matrix M (rows > 0)
 * @param [in] columns columns of matrix M (columns > 0)
 *
 * @retval A real matrix if successful, NULL otherwise
 */
MATRIX_ALLOC(matrix_alloc, matrix_t, double)

/** @brief Allocate memory for complex matrix of m rows and n columns
 *
 * See matrix_alloc.
 *
 * @param [in] rows rows of matrix M (rows > 0)
 * @param [in] columns columns of matrix M (columns > 0)
 *
 * @retval A complex matrix if successful, otherwise NULL
 */
MATRIX_ALLOC(matrix_complex_alloc, matrix_complex_t, complex_t)

/** macro to create zero matrix */
#define MATRIX_ZEROS(FUNCTION_NAME, MTYPE, TYPE, ALLOC, SETALL) \
MTYPE *FUNCTION_NAME(int rows, int columns, MTYPE *A) \
{ \
    if(A == NULL) \
    { \
        A = ALLOC(rows,columns); \
        if(A == NULL) \
            return NULL; \
    } \
\
    SETALL(A, 0); \
\
    return A; \
}

/** @brief Generate a real zero matrix of m lines and n columns
 *
 * Set every element of matrix A to 0. If A is NULL, the matrix will be
 * created. In this case, you have to free the matrix yourself.
 *
 * @param [in] rows rows of matrix M
 * @param [in] columns columns of matrix M
 * @param [in,out] A: matrix
 *
 * @retval A real matrix 0 if successful, otherwise NULL
 */
MATRIX_ZEROS(matrix_zeros, matrix_t, double, matrix_alloc, matrix_setall)

/** @brief Generate a complex zero matrix of m lines and n columns
 *
 * See matrix_zeros.
 *
 * @param [in] rows rows of matrix M
 * @param [in] columns columns of matrix M
 * @param [in,out] A: matrix
 *
 * @retval A complex matrix 0 if successful, otherwise NULL
 */
MATRIX_ZEROS(matrix_complex_zeros, matrix_complex_t, complex_t, matrix_complex_alloc, matrix_complex_setall)

/** macro to create matrix x*Id */
#define MATRIX_SETALL(FUNCTION_NAME, MATRIX_TYPE, TYPE) \
void FUNCTION_NAME(MATRIX_TYPE *A, TYPE x) \
{ \
    const size_t size = A->size; \
    TYPE *M = A->M; \
    for(size_t i = 0; i < size; i++) \
        *M++ = x; \
}

/** @brief Set matrix elements of A to x
 *
 * @param [in,out] A real matrix
 * @param [in] x real number
 */
MATRIX_SETALL(matrix_setall, matrix_t, double)

/** @brief Set matrix elements of A to x
 *
 * @param [in,out] A complex matrix
 * @param [in] x complex number
 */
MATRIX_SETALL(matrix_complex_setall, matrix_complex_t, complex_t)

/** macro to create unity matrix */
#define MATRIX_EYE(FUNCTION_NAME, MTYPE, TYPE, ALLOC, SETALL) \
MTYPE *FUNCTION_NAME(int dim, MTYPE *A) \
{ \
    if(A == NULL) \
    { \
        A = ALLOC(dim,dim); \
        if(A == NULL) \
            return NULL; \
    } \
 \
    const int min = A->min; \
    TYPE *M = A->M; \
    SETALL(A,0); \
    for(int i = 0; i < min; i++) \
        *(M+i*(min+1)) = 1; \
 \
    return A; \
}

/** @brief Create real identity matrix
 *
 * If A == NULL, a identity matrix of dimension dim times dim is created and
 * returned.
 * If A != NULL, the matrix A is set to the identity matrix and the parameter
 * dim is ignored. More specific, for a general (e.g. not square) matrix A, the
 * matrix elements are set to A_ij = Delta_ij.
 *
 * @param [in]     dim dimension of identity matrix (ignored for A == NULL)
 * @param [in,out] A real matrix
 *
 * @retval A identity matrix if successful, NULL otherwise
 */
MATRIX_EYE(matrix_eye, matrix_t, double, matrix_alloc, matrix_setall)

/** @brief Create complex identity matrix
 *
 * See matrix_eye.
 *
 * @param [in]     dim dimension of identity matrix (ignored for A == NULL)
 * @param [in,out] A complex matrix
 *
 * @retval A identity matrix if successful, NULL otherwise
 */
MATRIX_EYE(matrix_complex_eye, matrix_complex_t, complex_t, matrix_complex_alloc, matrix_complex_setall)

/** macro to free matrices */
#define MATRIX_FREE(FUNCTION_NAME, MTYPE) \
void FUNCTION_NAME(MTYPE *A) \
{ \
    if(A != NULL) \
    { \
        if(!A->view && A->M != NULL) \
        { \
            free_cb(A->M); \
            A->M = NULL; \
        } \
        free_cb(A); \
    } \
}

/** @brief Free real matrix
 *
 * This function will free the memory allocated for matrix M. If M is NULL this
 * function will do nothing.
 *
 * @param [in,out] A matrix to free
 */
MATRIX_FREE(matrix_free, matrix_t)

/** @brief Free complex matrix
 *
 * See matrix_free.
 *
 * @param [in,out] A matrix to free
 */
MATRIX_FREE(matrix_complex_free, matrix_complex_t)

/** @}*/

/** \defgroup trdet trace and determinant
 *  @{
 */

/** macro to calculate trace of matrix */
#define MATRIX_TRACE(FUNCTION_NAME, MATRIX_TYPE, TYPE) \
TYPE FUNCTION_NAME(MATRIX_TYPE *A) \
{ \
    const int min = A->min; \
    const TYPE *M = A->M; \
    TYPE trace = 0; \
    for(int i = 0; i < min; i++) \
        trace += *(M+i*(1+min)); \
\
    return trace; \
}

/** @brief Calculate Tr(A) of real matrix A
 *
 * @param [in] A complex matrix
 *
 * @retval x with x=Tr(A)
 */
MATRIX_TRACE(matrix_trace, matrix_t, double)

/** @brief Calculate Tr(A) of complex matrix A
 *
 * @param [in] A complex matrix
 *
 * @retval z with z=Tr(A)
 */
MATRIX_TRACE(matrix_complex_trace, matrix_complex_t, complex_t)

/** @brief Calculate Tr(A*B)
 *
 * This will calculate the trace of A*B: Tr(A*B). Matrix A and B must be square
 * matrices of dimension dim, dim.
 *
 * A and B must not point to the same matrix or the behaviour will be
 * undefined!
 *
 * @param [in] A real matrix
 * @param [in] B real matrix
 *
 * @retval x with x=Tr(A*B)
 */
double matrix_trace_AB(matrix_t *A, matrix_t *B)
{
    const int dim = A->min;
    double sum = 0;
    double *M1 = A->M;
    double *M2 = B->M;

    for(int i = 0; i < dim; i++)
        sum += cblas_ddot(dim, M1+i, dim, M2+i*dim, 1);

    return sum;
}

/** @brief Calculate Tr(A*B) for A,B complex
 *
 * See matrix_trace_AB.
 *
 * @param [in] A complex matrix
 * @param [in] B complex matrix
 *
 * @retval z with z=Tr(A*B)
 */
complex_t matrix_trace_complex_AB(matrix_complex_t *A, matrix_complex_t *B)
{
    const int dim = A->min;
    complex_t sum = 0;

    for(int i = 0; i < dim; i++)
        for(int k = 0; k < dim; k++)
            sum += matrix_get(A, i,k)*matrix_get(B, k,i);

    return sum;
}

/** @brief Calculate Re(Tr(A*B)) for A,B complex
 *
 * See matrix_trace_AB.
 *
 * @param [in] A complex matrix
 * @param [in] B complex matrix
 *
 * @retval x with x=Re(Tr(A*B))
 */
double matrix_trace_complex_AB_real(matrix_complex_t *A, matrix_complex_t *B)
{
    const int dim = A->rows;
    double sum = 0;
    const complex_t *M1 = A->M;
    const complex_t *M2 = B->M;

    for(int i = 0; i < dim; i++)
        for(int k = 0; k < dim; k++)
            sum += CREAL(M1[i*dim+k]*M2[k*dim+i]);

    return sum;
}

/** @}*/

/** \defgroup kron Kronecker product
 *  @{
 */

/** macro to create Kronecker product */
#define MATRIX_KRON(FUNCTION_NAME, MTYPE, TYPE, ALLOC, SETALL) \
MTYPE *FUNCTION_NAME(MTYPE *A, MTYPE *B, MTYPE *C) \
{ \
    const int Am = A->rows, An = A->columns; \
    const int Bm = B->rows, Bn = B->columns; \
    if(C == NULL) \
    { \
        C = ALLOC(Am*Bm, An*Bn); \
        if(C == NULL) \
            return NULL; \
    } \
    SETALL(C, 0); \
\
    for(int m = 0; m < Am; m++) \
        for(int n = 0; n < An; n++) \
        { \
            const TYPE c = matrix_get(A,m,n); \
            if(c != 0) \
            { \
                for(int im = 0; im < Bm; im++) \
                    for(int in = 0; in < Bn; in++) \
                        matrix_set(C, m*Bm+im, n*Bn+in, c*matrix_get(B,im,in)); \
            } \
        } \
\
    return C; \
}

/** @brief Calculate Kronecker product of real matrices A and B
 *
 * Matrix A has dimension Am,An, matrix B has dimension Bm,Bn.
 *
 * If C is not NULL, the Kronecker product will be stored in C. C must have
 * dimension Am+Bm,An+Bn. If C is NULL, memory for the matrix will be allocated
 * and the matrix will be returned. You have to free the memory for the
 * returned matrix yourself.
 *
 * @param [in] A real matrix
 * @param [in] B real matrix
 * @param [in,out] C Kronecker product
 *
 * @retval C Kronecker product of A and B, NULL if no memory could be allocated
 */
MATRIX_KRON(matrix_kron, matrix_t, double, matrix_alloc, matrix_setall)

/** @brief Calculate Kronecker product of complex matrices A and B
 *
 * See matrix_kron.
 *
 * @param [in] A complex matrix
 * @param [in] B complex matrix
 * @param [in,out] C Kronecker product
 *
 * @retval C Kronecker product of A and B
 */
MATRIX_KRON(matrix_complex_kron, matrix_complex_t, complex_t, matrix_complex_alloc, matrix_complex_setall)

/** @}*/

/** \defgroup addsubmult Add, subtract and multiply matrices
 *  @{
 */

/** macro to multiply matrix with scalar factor */
#define MATRIX_MULT_SCALAR(FUNCTION_NAME, MTYPE, TYPE) \
void FUNCTION_NAME(MTYPE *A, TYPE alpha) \
{ \
    const size_t size = A->size; \
    TYPE *M = A->M; \
    for(size_t i = 0; i < size; i++) \
        M[i] *= alpha; \
}

/** @brief Multiply real matrix with a real scalar
 *
 * alpha*A -> A
 *
 * @param [in,out] A real matrix
 * @param [in] alpha real scalar
 * @retval A, A=alpha*A
 */
MATRIX_MULT_SCALAR(matrix_mult_scalar, matrix_t, double)

/** @brief Multiply complex matrix with a complex scalar
 *
 * alpha*A -> A
 *
 * @param [in,out] A complex matrix
 * @param [in] alpha complex scalar
 * @retval A, A=alpha*A
 */
MATRIX_MULT_SCALAR(matrix_complex_mult_scalar, matrix_complex_t, complex_t)

/** @brief Multiply real matrix with a complex scalar
 *
 * alpha*A -> C
 *
 * If C is NULL, memory for the matrix C will be allocated. If C is not NULL, C
 * must have the correct dimension.
 *
 * @param [in] A real matrix
 * @param [in] alpha complex scalar
 * @param [in,out] C complex matrix
 *
 * @retval C, C = alpha*A
 */
matrix_complex_t *matrix_mult_complex_scalar(matrix_t *A, complex_t alpha, matrix_complex_t *C)
{
    const int rows    = A->rows;
    const int columns = A->columns;
    const double *AM  = A->M;
    complex_t *CM;
    if(C == NULL)
    {
        C = matrix_complex_alloc(rows, columns);
        if(C == NULL)
            return NULL;
    }
    CM = C->M;

    for(size_t i = 0; i < C->size; i++)
        CM[i] = alpha*AM[i];

    return C;
}

/* compute alpha*A*B */
#define MATRIX_MULT(FUNCTION_NAME, MATRIX_TYPE, TYPE, ALLOC, BLAS_xGEMM, PTR) \
MATRIX_TYPE *FUNCTION_NAME(MATRIX_TYPE *A, MATRIX_TYPE *B, TYPE alpha, MATRIX_TYPE *C) \
{ \
    TYPE beta = 0; \
\
    if(A->columns != B->rows) \
        return NULL; \
\
    if(C == NULL) \
    { \
        C = ALLOC(A->rows, B->columns); \
        if(C == NULL) \
            return NULL; \
    } \
\
    /* xGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC) */ \
    BLAS_xGEMM(CblasColMajor, /* column order */ \
               CblasNoTrans,  /* don't transpose/conjugate A */ \
               CblasNoTrans,  /* don't transpose/conjugate B */ \
               A->rows,       /* M: rows of A and C */ \
               B->columns,    /* N: columns of B and C */ \
               A->columns,    /* K: columns of A and rows of B */ \
               PTR alpha,     /* alpha: scalar */ \
               A->M,          /* A: matrix A */ \
               A->rows,       /* LDA: leading dimension of A (columns) */ \
               B->M,          /* B: matrix B */ \
               A->columns,    /* LDB: leading dimension of B (columns) */ \
               PTR beta,      /* beta: scalar */ \
               C->M,          /* C: matrix C */ \
               C->rows        /* LDC: leading dimension of C (columns) */ \
    ); \
\
    return C; \
}

/** @brief Multiply real matrices
 *
 * alpha*A*B -> C
 *
 * If C is NULL, memory for the matrix C will be allocated.
 *
 * @param [in] A real matrix
 * @param [in] B real matrix
 * @param [in] alpha real scalar
 * @param [in,out] C real matrix
 *
 * @retval C, C = alpha*A*B
 */
MATRIX_MULT(matrix_mult, matrix_t, double, matrix_alloc, cblas_dgemm, +)

/** @brief Multiply complex matrices
 *
 * alpha*A*B -> C
 *
 * If C is NULL, memory for the matrix C will be allocated.
 *
 * @param [in] A complex matrix
 * @param [in] B complex matrix
 * @param [in] alpha complex scalar
 * @param [in,out] C complex matrix
 *
 * @retval C, C = alpha*A*B
 */
MATRIX_MULT(matrix_complex_mult, matrix_complex_t, complex_t, matrix_complex_alloc, cblas_zgemm, &)

/** macro to add two matrices */
#define MATRIX_ADD(FUNCTION_NAME, TYPE1, MTYPE1, TYPE2, MTYPE2) \
int FUNCTION_NAME(MTYPE1 *A, MTYPE2 *B, TYPE1 alpha, MTYPE1 *C) \
{ \
    const size_t size = A->size; \
    TYPE1 *M3; \
    TYPE1 *M1 = A->M; \
    TYPE2 *M2 = B->M; \
\
    if(C == NULL) \
        M3 = A->M; \
    else \
        M3 = C->M; \
\
    if(A->rows != B->rows || A->columns != B->columns) \
        return LIBHADES_ERROR_SHAPE; \
\
    for(size_t i = 0; i < size; i++) \
        M3[i] = M1[i] + (alpha*M2[i]); \
\
    return 0; \
}

/** @brief Add real matrices A and B
 *
 * Calculate A+alpha*B -> C.
 *
 * The result will be stored in C. If C is NULL, the result is stored in A.
 *
 * @param [in,out] A real matrix
 * @param [in] B real matrix
 * @param [in] alpha real scalar
 * @param [in,out] C real matrix or NULL
 *
 * @retval 0 if successfull
 * @retval LIBHADES_ERROR_SHAPE if matrices have wrong shape
 */
MATRIX_ADD(matrix_add, double, matrix_t, double, matrix_t)

/** @brief Add complex matrices A and B
 *
 * Calculate A+alpha*B -> C.
 *
 * The result will be stored in C. If C is NULL, the result is stored in A.
 *
 * @param [in,out] A complex matrix
 * @param [in] B complex matrix
 * @param [in] alpha complex number
 * @param [in,out] C complex matrix or NULL
 *
 * @retval 0 if successfull
 * @retval LIBHADES_ERROR_SHAPE if matrices have wrong shape
 */
MATRIX_ADD(matrix_complex_add, complex_t, matrix_complex_t, complex_t, matrix_complex_t)

/** @brief Add complex matrix A and real matrix B
 *
 * Calculate A+alpha*B -> C.
 *
 * The result will be stored in C. If C is NULL, the result is stored in A.
 *
 * @param [in,out] A complex matrix
 * @param [in] B real matrix
 * @param [in] alpha complex scalar
 * @param [in,out] C complex matrix or NULL
 *
 * @retval 0 if successfull
 * @retval LIBHADES_ERROR_SHAPE if matrices have wrong shape
 */
MATRIX_ADD(matrix_complex_add_real, complex_t, matrix_complex_t, double, matrix_t)

/** @}*/

/** \defgroup transconj Transpose, conjugate
 *  @{
 */

#define MATRIX_TRANSPOSE(FUNCTION_NAME, MTYPE, TYPE) \
void FUNCTION_NAME(MTYPE *A) \
{ \
    const int rows    = A->rows; \
    const int columns = A->columns; \
    for(int im = 0; im < rows; im++) \
        for(int in = im+1; in < columns; in++) \
        { \
            TYPE temp = matrix_get(A, im, in); \
            matrix_set(A, im, in, matrix_get(A, in, im)); \
            matrix_set(A, in, im, temp); \
        } \
\
    A->rows    = columns; \
    A->columns = rows; \
}

/** @brief Transpose real matrix A
 *
 * The matrix will be transposed.
 *
 * @param [in,out] A real matrix
 *
 * @retval C with C=A^T
 */
MATRIX_TRANSPOSE(matrix_transpose, matrix_t, double)

/** @brief Transpose complex matrix A
 *
 * The matrix will be transposed.
 *
 * @param [in,out] A complex matrix
 *
 * @retval C with C=A^T
 */
MATRIX_TRANSPOSE(matrix_complex_transpose, matrix_complex_t, complex_t)

/** @}*/

/** \defgroup ev Eigenvalue problems
 *  @{
 */

/** @brief Compute eigenvalues and optionally eigenvectors of symmetric matrix A
 *
 * This function computes all eigenvalues and, optionally, eigenvectors of a
 * real symmetric matrix A.
 *
 * See dsyev.
 *
 * @param [in] A real matrix
 * @param [in] JOBZ 'N': only eigenvalues, 'V' eigenvalues and eigenvectors
 * @param [in] UPLO 'U': upper triangle part of A is stored; 'L': lower triangle part of A is stored
 * @param [in] w real matrix of dimension (dim,1) (i.e. a vector); the eigenvalues will be stored in w
 *
 * @retval 0 on success
 */
int eig_sym(matrix_t *A, char *JOBZ, char *UPLO, matrix_t *w)
{
    int info, lwork = -1, N = A->min;
    double workopt;
    double *work;

    dsyev_(JOBZ, UPLO, &N, A->M, &N, w->M, &workopt, &lwork, &info);
    if(info != 0)
        return info;

    lwork = workopt;
    work = malloc_cb(lwork*sizeof(double));
    if(work == NULL)
        return LIBHADES_ERROR_OOM;

    dsyev_(JOBZ, UPLO, &N, A->M, &N, w->M, work, &lwork, &info);

    free_cb(work);

    return info;
}

/** @brief Compute eigenvalues and optionally eigenvectors of Hermitian matrix A
 *
 * This function computes all eigenvalues and, optionally, eigenvectors of a
 * Hermitian symmetric matrix A.
 *
 * See zheev.
 *
 * @param [in] A complex matrix
 * @param [in] JOBZ 'N': only eigenvalues, 'V' eigenvalues and eigenvectors
 * @param [in] UPLO 'U': upper triangle part of A is stored; 'L': lower triangle part of A is stored
 * @param [in] w real matrix of dimension (dim,1) (i.e. a vector); the eigenvalues will be stored in w
 *
 * @retval 0 on success
 */
int eig_herm(matrix_complex_t *A, char *JOBZ, char *UPLO, matrix_t *w)
{
    int info, lwork = -1, N = A->min;
    complex_t workopt;
    complex_t *work = NULL;
    double *rwork;

    rwork = malloc_cb(MAX(1, 3*N-2)*sizeof(double));
    if(rwork == NULL)
        return LIBHADES_ERROR_OOM;

    zheev_(JOBZ, UPLO, &N, A->M, &N, w->M, &workopt, &lwork, rwork, &info);
    if(info != 0)
    {
        free_cb(rwork);
        return info;
    }

    lwork = CREAL(workopt);
    work = malloc_cb(lwork*sizeof(complex_t));
    if(work == NULL)
    {
        free_cb(rwork);
        return LIBHADES_ERROR_OOM;
    }

    zheev_(JOBZ, UPLO, &N, A->M, &N, w->M, work, &lwork, rwork, &info);

    free_cb(rwork);
    free_cb(work);

    return info;
}

/** @brief Compute eigenvalues and optionally eigenvectors of matrix A
 *
 * Compute for an N-by-N complex nonsymmetric matrix A, the eigen- values and,
 * optionally, right eigenvectors
 *
 * See zgeev.
 *
 * @param [in] A real matrix
 * @param [in] w list containing the eigenvalues of A
 * @param [in] vr if vr != NULL, right eigenvectors are computed and will be stored in vr; vr must be a complex matrix of dimension (dim,dim)
 * @param [in] vl if vl != NULL, lefft eigenvectors are computed and will be stored in vl; vl must be a complex matrix of dimension (dim,dim)
 *
 * @retval 0 on success
 */

/* Parameters */
int eig_complex_generic(matrix_complex_t *A, matrix_complex_t *w, matrix_complex_t *vl, matrix_complex_t *vr)
{
    int N = A->min;
    int lwork = -1;
    char jobvr = 'N';
    char jobvl = 'N';
    int info;
    complex_t *evr  = NULL;
    complex_t *evl  = NULL;
    complex_t *work = NULL;
    double *rwork   = NULL;
    complex_t wopt;

    /* if vr is not NULL, calculate right eigenvectors */
    if(vr != NULL)
    {
        jobvr = 'V';
        evr   = vr->M;
    }
    /* if vl is not NULL, calculate left eigenvectors */
    if(vl != NULL)
    {
        jobvl = 'V';
        evl   = vl->M;
    }

    /* get the optimal size for workspace work */
    zgeev_(&jobvl, &jobvr, &N, A->M, &N, w->M, evl, &N, evr, &N, &wopt, &lwork, rwork, &info);

    if(info != 0)
        return info;

    lwork = CREAL(wopt);

    rwork = malloc_cb(2*N*sizeof(double));
    work  = malloc_cb(lwork*sizeof(complex_t));

    if(rwork == NULL || work == NULL)
    {
        if(rwork != NULL)
            free(rwork);
        if(work != NULL)
            free(work);

        return LIBHADES_ERROR_OOM;
    }

    /* SUBROUTINE ZGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO ) */
    zgeev_(
        &jobvl,     /* left eigenvectors of A are not computed */
        &jobvr,     /* calculate/don't calculate right eigenvectors */
        &N,         /* order of matrix A */
        A->M,       /* matrix A */
        &N,         /* LDA - leading dimension of A */
        w->M,       /* eigenvalues */
        evl,        /* left eigenvectors */
        &N,         /* leading dimension of array vl */
        evr,        /* eigenvectors */
        &N,         /* leading dimension of the array VR */
        work,       /* COMPLEX*16 array, dimension (LWORK) */
        &lwork,     /* dimension of the array WORK; LWORK >= max(1,2*N) */
        rwork,      /* (workspace) DOUBLE PRECISION array, dimension (2*N) */
        &info       /* 0 == success */
    );

    free_cb(work);
    free_cb(rwork);

    return info;
}

int eig_complex_vr(matrix_complex_t *M, complex_t lambda, matrix_complex_t *vr)
{
    matrix_complex_t *A = NULL, *b = NULL;
    const int rows    = M->rows;
    const int columns = M->columns;

    /* lapack */
    char trans = 'N';
    int info, lwork;
    complex_t work_size, *work = NULL;
    int m = rows+1, n = columns, nrhs = 1;
    int lda = m, ldb = MAX(m,n);

    if((A = matrix_complex_alloc(rows+1, columns)) == NULL)
        return LIBHADES_ERROR_OOM;

    if((b = matrix_complex_zeros(rows+1,1,NULL)) == NULL)
    {
        matrix_complex_free(A);
        return LIBHADES_ERROR_OOM;
    }
    matrix_set(b, rows,0, 1);

    for(int j = 0; j < columns; j++)
    {
        for(int i = 0; i < rows; i++)
            matrix_set(A, i, j, matrix_get(M, i,j));

        /* - lambda Id */
        matrix_set(A, j,j, matrix_get(A,j,j)-lambda);

        /* set nomalization condition */
        matrix_set(A, rows,j, 1);
    }

    /* determine work size */
    lwork = -1;
    zgels_(&trans, &m, &n, &nrhs, A->M, &lda, b->M, &ldb, &work_size, &lwork, &info);
    lwork = (int)CREAL(work_size);
    if((work = malloc_cb(lwork*sizeof(complex_t))) == NULL)
    {
        matrix_complex_free(A);
        matrix_complex_free(vr);
        return LIBHADES_ERROR_OOM;
    }

    /* calculate eigenvector */
    zgels_(
        &trans, /* find least squares solution of overdetermined system */
        &m,     /* number of rows of matrix A */
        &n,     /* number of columns of matrix A */
        &nrhs,  /* number of columns of matrix B */
        A->M,   /* matrix A (will be overwritten) */
        &lda,   /* leading dimension of A, LDA >= max(1,M) */
        b->M,   /* on entry: right hand side vectors, on exit: solution */
        &ldb,   /* leading dimension of B, LDB >= MAX(1,M,N) */
        work,   /* workspace */
        &lwork, /* dimension of work */
        &info   /* status */
    );

    free_cb(work);
    matrix_complex_free(A);

    for(int i = 0; i < rows; i++)
        matrix_set(vr, i,0, matrix_get(b, i,0));

    matrix_complex_free(b);

    return info;
}

/** @}*/

/** \defgroup exp Calculate matrix exponential
 *  @{
 */

/** macro to calculate matrix norm */
#define MATRIX_NORM(FUNCTION_NAME, MTYPE, LAPACK_FUNC) \
int FUNCTION_NAME(MTYPE *A, char norm_type, double *norm) \
{ \
    double *work = NULL; \
\
    if(norm_type == 'I') \
    { \
        work = malloc_cb(A->rows*sizeof(double)); \
        if(work == NULL) \
            return LIBHADES_ERROR_OOM; \
    } \
\
    *norm = LAPACK_FUNC(&norm_type, &A->rows, &A->columns, A->M, &A->rows, work); \
\
    if(work != NULL) \
        free_cb(work); \
\
    return 0; \
}

/** @brief Compute matrix norm for real matrix
 *
 * See dlange.
 *
 * Returns
 *      max(abs(A(i,j))), norm_type = 'M' or 'm'
 *      norm1(A),         norm_type = '1', 'O' or 'o'
 *      normI(A),         norm_type = 'I' or 'i'
 *      normF(A),         norm_type = 'F', 'f', 'E' or 'e'
 *
 * @param [in]  A real matrix
 * @param [in]  norm_type type of norm, e.g. 'F' for Frobenius norm
 * @param [out] norm value of norm
 *
 * @retval ret 0 if successfull, <0 otherwise
 */
MATRIX_NORM(matrix_norm, matrix_t, dlange_);

/** @brief Compute matrix norm for complex matrix
 *
 * See matrix_norm.
 *
 * @param [in]  A complex matrix
 * @param [in]  norm_type type of norm, e.g. 'F' for Frobenius norm
 * @param [out] norm value of norm
 *
 * @retval ret 0 if successfull, <0 otherwise
 */
MATRIX_NORM(matrix_complex_norm, matrix_complex_t, zlange_);

matrix_complex_t *matrix_complex_exp_taylor(matrix_complex_t *A, int order)
{
    const int rows = A->min;
    matrix_complex_t *C = matrix_complex_alloc(rows,rows);

    /* B = E+A/order */
    matrix_complex_t *B = matrix_complex_copy(A,NULL);
    matrix_complex_mult_scalar(B, 1./order);

    for(int i = 0; i < rows; i++)
        matrix_set(B, i,i, 1+matrix_get(B,i,i));

    for(int k = order-1; k > 0; k--)
    {
        matrix_complex_mult(A, B, 1./order, C);

        /* swap */
        {
            matrix_complex_t *temp;
            temp = C;
            C = B;
            B = temp;
        }

        for(int i = 0; i < rows; i++)
            matrix_set(B, i,i, 1+matrix_get(B,i,i));
    }

    matrix_complex_free(C);

    return B;
}

/** @}*/

/** \defgroup LA LU decomposition, inverting
 *  @{
 */

#define LU_DECOMPOSITION(FUNCTION_NAME, MATRIX_TYPE, XGETRF) \
int FUNCTION_NAME(MATRIX_TYPE *A, int ipiv[]) \
{ \
    int info; \
\
    XGETRF( \
        &A->rows,    /* M number of rows of A */ \
        &A->columns, /* N number of columns of A */ \
        A->M,        /* matrix A to be factored */ \
        &A->columns, /* LDA: leading dimension of A */ \
        ipiv,        /* pivot indices of dimension (min(M,N)) */ \
        &info \
    ); \
\
    return info; \
}

/** @brief Compute LU decomposition of real matrix A
 *
 * See dgetrf.
 *
 * The factorization has the form
 *    A = P * L * U
 * where P is a permutation matrix, L is lower triangular with unit
 * diagonal elements (lower trapezoidal if rows > columns), and U is upper
 * triangular (upper trapezoidal if m < n).
 *
 * @param [in,out] A real matrix
 * @param [out] ipiv pivot indices; array of dimension MIN(rows,columns)
 *
 * @retval INFO
 */
LU_DECOMPOSITION(matrix_lu_decomposition, matrix_t, dgetrf_)

/** @brief Compute LU decomposition of complex matrix A
 *
 * See matrix_lu_decomposition.
 *
 * @param [in,out] A complex matrix
 * @param [out] ipiv pivot indices; array of dimension MIN(rows,columns)
 *
 * @retval INFO
 */
LU_DECOMPOSITION(matrix_complex_lu_decomposition, matrix_complex_t, zgetrf_)

#define MATRIX_INVERT(FUNCTION_NAME, TYPE, MATRIX_TYPE, LU_DECOMPOSITION, XGETRI) \
int FUNCTION_NAME(MATRIX_TYPE *A) \
{ \
    int info, lwork, dim = A->min; \
    int *ipiv = NULL; \
    TYPE *work = NULL; \
    TYPE workopt; \
\
    ipiv = malloc_cb(dim*sizeof(int)); \
    if(ipiv == NULL) \
        return LIBHADES_ERROR_OOM; \
\
    info = LU_DECOMPOSITION(A, ipiv); \
    if(info != 0) \
        goto out; \
\
    lwork = -1; \
    XGETRI(&dim, A->M, &dim, ipiv, &workopt, &lwork, &info); \
    if(info != 0) \
        goto out; \
\
    lwork = (int)workopt; \
    work = malloc_cb(lwork*sizeof(TYPE)); \
    if(work == NULL) \
    { \
        info = LIBHADES_ERROR_OOM; \
        goto out; \
    } \
\
    XGETRI( \
        &dim,   /* order of matrix A */ \
        A->M,   /* factors L and U from LU decomposition */ \
        &dim,   /* LDA: leading dimension of A */ \
        ipiv,   /* pivot indices */ \
        work,   /* workspace of dimension LWORK */ \
        &lwork, /* length of work */ \
        &info \
    ); \
\
out: \
    free_cb(ipiv); \
    if(work != NULL) \
        free_cb(work); \
\
    return info; \
}

/** @brief Invert real matrix A
 *
 * The inverse of A is computed using the LU factorzation of A
 *
 * @param [in] A real matrix
 *
 * @retval INFO
 */
MATRIX_INVERT(matrix_invert, double, matrix_t, matrix_lu_decomposition, dgetri_)

/** @brief Invert complex matrix A
 *
 * See matrix_invert.
 *
 * @param [in] A complex matrix
 *
 * @retval INFO
 */
MATRIX_INVERT(matrix_complex_invert, complex_t, matrix_complex_t, matrix_complex_lu_decomposition, zgetri_)

#define MATRIX_SOLVE(FUNCTION_NAME, MATRIX_TYPE, LU_DECOMPOSITION, XGETRS) \
int FUNCTION_NAME(MATRIX_TYPE *A, MATRIX_TYPE *b) \
{ \
    int N = A->min; \
    char trans = 'N'; \
    int nrhs = b->columns; \
    int info = -1; \
    int *ipiv = malloc_cb(N*sizeof(int)); \
\
    if(ipiv == NULL) \
        return LIBHADES_ERROR_OOM; \
\
    LU_DECOMPOSITION(A, ipiv); \
    XGETRS(&trans, &N, &nrhs, A->M, &N, ipiv, b->M, &N, &info); \
\
    free_cb(ipiv); \
\
    return info; \
}

/** @brief Solve system of linear equations
 *
 * Solve the system of linear equations:
 *      A*x = b
 *
 * @param [in,out] A matrix
 * @param [in,out] b vector/matrix
 *
 * @retval INFO
 */

/** @}*/
MATRIX_SOLVE(matrix_solve, matrix_t, matrix_lu_decomposition, dgetrs_);

/** @brief Solve system of linear equations
 *
 * Solve the system of complex linear equations:
 *      A*x = b
 *
 * @param [in,out] A complex matrix
 * @param [in,out] b complex vector/matrix
 *
 * @retval INFO
 */

/** @}*/
MATRIX_SOLVE(matrix_complex_solve, matrix_complex_t, matrix_complex_lu_decomposition, zgetrs_);

/*
static void _cblas_zaxpy(const int N, const double alpha, const void *X, const int incX, void *Y, const int incY)
{
    complex_t beta = alpha;
    cblas_zaxpy(N, &beta, X, incX, Y, incY);
}
*/

/** @brief Calculate dot product of vectors x and y
 *
 * Calculate dot product of first column of x,y. If x and y have different
 * rows, the minimum is used.
 *
 * @param [in] x vector
 * @param [in] y vector
 * @retval x*y
 */
double vector_dot(matrix_t *x, matrix_t *y)
{
    int incx = 1;
    int incy = 1;
    int N = MIN(x->rows, y->rows);

    return ddot_(&N, x->M, &incx, y->M, &incy);
}

/** @brief Calculate dot product of vectors x and y
 *
 * Calculate dot product of first column of x,y. If x and y have different
 * rows, the minimum is used.
 *
 * @param [in] x vector
 * @param [in] y vector
 * @retval x*y
 */
complex_t vector_complex_dot(matrix_complex_t *x, matrix_complex_t *y)
{
    /*
    int incx = 1;
    int incy = 1;
    int N = MIN(x->rows, y->rows);
    complex_t z = 0;

    zdotc_(&z, &N, x->M, &incx, y->M, &incy);
    return z;
    */

    int N = MIN(x->rows, y->rows);
    complex_t z = 0;
    for(int i = 0; i < N; i++)
        z += matrix_get(x,i,0)*matrix_get(y,i,0);
    return z;
}

#define MATRIX_GET_COLUMN(FUNCTION_NAME, TYPE, MATRIX_TYPE) \
MATRIX_TYPE *FUNCTION_NAME(MATRIX_TYPE *A, int i) \
{ \
    const int rows = A->rows; \
    MATRIX_TYPE *v = malloc_cb(sizeof(MATRIX_TYPE)); \
    if(v == NULL) \
        return NULL; \
\
    v->rows    = rows; \
    v->columns = 1; \
    v->min     = 1; \
    v->size    = rows; \
    v->type    = 0; \
    v->view    = 1; \
    v->M       = &A->M[i*rows]; \
\
    return v; \
}

/** @brief Get i-th column of matrix A
 *
 * @param [in] A matrix
 * @param [in] i column number
 * @retval x*y
 */
MATRIX_GET_COLUMN(matrix_get_column,         double,    matrix_t);

/** @brief Get i-th column of matrix A
 *
 * @param [in] A matrix
 * @param [in] i column number
 * @retval x*y
 */
MATRIX_GET_COLUMN(matrix_complex_get_column, complex_t, matrix_complex_t);

/** \defgroup sparse Functions for sparse matrices
 *  @{
 */

#ifdef SUPPORT_SPARSE

/** @brief Calculate eigenvalues of a sparse complex matrix
 *
 * @param [in] N     number of columns/rows of matrix
 * @param [in] nev   number of eigenvalues to compute
 * @param [in] which LM (largest magnitude), SM (smallest magnitude), LR (largest real part), SR (smallest real part), LI (largest imaginary part), SI (smallest imaginary part)
 * @param [in] Av callback function that implements the matrix-vector operation Av; the input vector is given as in, the vector Av must be written in out
 * @param [in,out] d on exit d contains the Rith approximations (must be of length nev+1)
 * @param [in] mxiter maximum number of Arnoldi update iterations allowed
 * @param [in] tol relative accuracy of the Ritz value
 * @param [in] data pointer that is given to callback function Av
 *
 * @retval 0 if successful
 */
int sparse_complex_eig(int N, int nev, char *which, void (*Av)(int N, complex_t *in, complex_t *out, void *data), complex_t *d, int mxiter, double tol, void *data)
{
    int ret = LIBHADES_ERROR_OOM;
    int info = 0;
    int ido = 0;
    char *bmat = "I";          /* standard eigenproblem */
    int ncv = MIN(2*nev+2, N);

    int ishift = 1;
    int mode = 1;
    int iparam[11] = { ishift, 0, mxiter, 1, 0, 0, mode, 0, 0, 0, 0 };

    int ipntr[14];
    int lworkl = ncv*(3*ncv + 5);
    int rvec = 0;
    char *howmny = "P";

    complex_t *workd = NULL, *workl = NULL, *resid = NULL, *v = NULL, *workev = NULL;
    int *select = NULL;
    double *rwork = NULL;

    /* allocate memory */
    ret = LIBHADES_ERROR_OOM;

    workd = malloc_cb(3*N*sizeof(complex_t));
    if(workd == NULL)
        goto out;
    workl = malloc_cb(lworkl*sizeof(complex_t));
    if(workl == NULL)
        goto out;
    rwork = malloc_cb(ncv*sizeof(double));
    if(rwork == NULL)
        goto out;
    resid  = malloc_cb(N*sizeof(complex_t));
    if(resid == NULL)
        goto out;
    v = malloc_cb(N*ncv*sizeof(complex_t));
    if(v == NULL)
        goto out;
    select = malloc_cb(ncv*sizeof(int));
    if(select == NULL)
        goto out;
    workev = malloc_cb((2*ncv)*sizeof(complex_t));
    if(workev == NULL)
        goto out;

    /* loop */
    while(1)
    {
        /* http://www.caam.rice.edu/software/ARPACK/UG/node138.html */
        znaupd_(
            &ido, bmat, &N, which, &nev, &tol, resid, &ncv, v, &N, iparam,
            ipntr, workd, workl, &lworkl, rwork, &info, strlen(bmat),
            strlen(which)
        );

        if(ido == 1 || ido == -1)
            Av(N, &workd[ipntr[0]-1], &workd[ipntr[1]-1], data);
        else if(ido == 99)
            break;
        else
        {
            ret = ido;
            goto out;
        }
    }

    if(info != 0)
    {
        ret = info;
        goto out;
    }

    /* http://www.mathkeisan.com/usersguide/man/zneupd.html */
    zneupd_(
        &rvec, howmny, select, d, NULL, &N, NULL, workev, bmat, &N, which,
        &nev, &tol, resid, &ncv, v, &N, iparam, ipntr, workd, workl, &lworkl,
        rwork, &info, strlen(howmny), strlen(bmat), strlen(which)
    );

    ret = info;

out:
    if(resid != NULL)
        free_cb(resid);
    if(v != NULL)
        free_cb(v);
    if(workd != NULL)
        free_cb(workd);
    if(workl != NULL)
        free_cb(workl);
    if(rwork != NULL)
        free_cb(rwork);
    if(select != NULL)
        free_cb(select);
    if(workev != NULL)
        free_cb(workev);

    return ret;
}

#endif
/** @}*/

/** \defgroup io Save/load matrices
 *  @{
 */

#define MATRIX_LOAD_FROM_STREAM(FUNCTION_NAME, MATRIX_TYPE, TYPE, ALLOC, TRANSPOSE, IS_COMPLEX) \
MATRIX_TYPE *FUNCTION_NAME(FILE *stream, int *ret) \
{ \
    MATRIX_TYPE *M; \
    uint16_t len; \
    int rows, columns, fortran_order, is_complex; \
    char header[10] = { 0 }; \
    char dict[2048] = { 0 }; \
\
    if(ret != NULL) \
        *ret = 0; \
\
    /* read magic string, major and minor number */ \
    fread(header, 8, 1, stream); \
    if(memcmp(header, "\x93NUMPY\x01\x00", 8) != 0) \
    { \
        if(ret != NULL) \
            *ret = LIBHADES_ERROR_HEADER; \
        return NULL; \
    } \
\
    /* read length of dict */ \
    fread(&len, sizeof(uint16_t), 1, stream); \
\
    if(len >= sizeof(dict)/sizeof(dict[0])) \
    { \
        if(ret != NULL) \
            *ret = LIBHADES_ERROR_INV_LENGTH; \
        return NULL; \
    } \
\
    fread(dict, sizeof(char), len, stream); \
\
    if(npy_dict_get_fortran_order(dict, &fortran_order) != 0) \
    { \
        if(ret != NULL) \
            *ret = LIBHADES_ERROR_ORDER; \
        return NULL; \
    } \
\
    if(npy_dict_get_shape(dict, &rows, &columns) != 0) \
    { \
        if(ret != NULL) \
            *ret = LIBHADES_ERROR_SHAPE; \
        return NULL; \
    } \
\
    if(npy_dict_get_descr(dict, &is_complex) != 0) \
    { \
        if(ret != NULL) \
            *ret = LIBHADES_ERROR_DESCR; \
        return NULL; \
    } \
\
    if(is_complex != IS_COMPLEX) \
    { \
        if(ret != NULL) \
            *ret = LIBHADES_ERROR_FORMAT; \
        return NULL; \
    } \
\
    M = ALLOC(rows,columns); \
    fread(M->M, sizeof(TYPE), rows*columns, stream); \
\
    if(!fortran_order) \
        TRANSPOSE(M); \
\
    M->min  = MIN(rows,columns); \
    M->size = (size_t)rows*(size_t)columns; \
    M->view = 0; \
    M->type = 0; \
\
    return M; \
}

/** @brief Load real matrix from stream
 *
 * Load real matrix A from a stream. This function will also allocate memory
 * for the matrix.
 *
 * If error != NULL, error will be set to:
 *  0                         if successful
 *  LIBHADES_ERROR_HEADER     if magic or major/minor number is invalid
 *  LIBHADES_ERROR_INV_LENGTH if length of dictionary is invalid (too long)
 *  LIBHADES_ERROR_ORDER      if order is invalid (Fortran/C order)
 *  LIBHADES_ERROR_SHAPE      if shape is invalid (rows/columns)
 *  LIBHADES_ERROR_DESCR      if dtype is wrong/not supported
 *  LIBHADES_ERROR_FORMAT     if wrong format (real instead of complex)
 *  0 rows/columns wrong
 *
 * @param [in] stream file handle of a opened file
 * @param [out] error error code
 * @retval A matrix
 * @retval NULL if an error occured
 */
MATRIX_LOAD_FROM_STREAM(matrix_load_from_stream, matrix_t, double, matrix_alloc, matrix_transpose, 0);

/** @brief Load complex matrix from stream
 *
 * Load complex matrix A from a stream. This function will also allocate memory
 * for the matrix.
 *
 * If error != NULL, error will be set to:
 *  0                         if successful
 *  LIBHADES_ERROR_HEADER     if magic or major/minor number is invalid
 *  LIBHADES_ERROR_INV_LENGTH if length of dictionary is invalid (too long)
 *  LIBHADES_ERROR_ORDER      if order is invalid (Fortran/C order)
 *  LIBHADES_ERROR_SHAPE      if shape is invalid (rows/columns)
 *  LIBHADES_ERROR_DESCR      if dtype is wrong/not supported
 *  LIBHADES_ERROR_FORMAT     if wrong format (complex instead of real)
 *  0 rows/columns wrong
 *
 * @param [in] stream file handle of a opened file
 * @param [out] error error code
 * @retval A matrix
 * @retval NULL if an error occured
 */
MATRIX_LOAD_FROM_STREAM(matrix_complex_load_from_stream, matrix_complex_t, complex_t, matrix_complex_alloc, matrix_complex_transpose, 1);

#define MATRIX_LOAD(FUNCTION_NAME, MATRIX_TYPE, LOAD_FUNCTION) \
MATRIX_TYPE *FUNCTION_NAME(const char *filename, int *ret) \
{ \
    FILE *stream; \
    MATRIX_TYPE *M; \
\
    if((stream = fopen(filename, "r")) == NULL) \
    { \
        if(ret != NULL) \
            *ret = LIBHADES_ERROR_IO; \
        return NULL; \
    } \
\
    M = LOAD_FUNCTION(stream, ret); \
\
    fclose(stream); \
\
    return M; \
}

/** @brief Load real matrix from file
 *
 * Load real matrix A from file given by filename. This function will also
 * allocate memory for the matrix. See \ref matrix_load_from_stream for errors.
 *
 * @param [in] filename path to the file
 * @param [out] error error code
 * @retval A matrix
 * @retval NULL if an error occured
 */
MATRIX_LOAD(matrix_load, matrix_t, matrix_load_from_stream);

/** @brief Load complex matrix from file
 *
 * Load complex matrix A from file given by filename. This function will also
 * allocate memory for the matrix. See \ref matrix_complex_load_from_stream for
 * errors.
 *
 * @param [in] filename path to the file
 * @param [out] error error code
 * @retval A matrix
 * @retval NULL if an error occured
 */
MATRIX_LOAD(matrix_complex_load, matrix_complex_t, matrix_complex_load_from_stream);

#define MATRIX_SAVE_TO_STREAM(FUNCTION_NAME, TYPE, MATRIX_TYPE, DTYPE) \
void FUNCTION_NAME(MATRIX_TYPE *M, FILE *stream) \
{ \
    char d_str[512] = { 0 }; \
    uint16_t len = 0; \
    const int rows = M->rows, columns = M->columns; \
\
    /* write magic string, major number and minor number */ \
    fwrite("\x93NUMPY\x01\x00", sizeof(char), 8, stream); \
\
    /* write length of header and header */ \
    snprintf(d_str, sizeof(d_str)/sizeof(d_str[0]), "{'descr': '%s', 'fortran_order': True, 'shape': (%d, %d), }", DTYPE, rows, columns); \
\
    len = strlen(d_str); \
\
    fwrite(&len,  sizeof(len),  1,   stream); \
    fwrite(d_str, sizeof(char), len, stream); \
\
    /* write matrix */ \
    fwrite(M->M, sizeof(TYPE), M->size, stream); \
}

/** @brief Save real matrix A to stream
 *
 * Save real matrix A to a file handle given by stream. The datatype
 * corresponds to Numpy's npy file format.
 *
 * @param [in] A      matrix to be dumped to file
 * @param [in] stream file handle of opened file
 */
MATRIX_SAVE_TO_STREAM(matrix_save_to_stream, double, matrix_t, "<f8");

/** @brief Save complex matrix A to stream
 *
 * Save complex matrix A to a file handle given by stream. The datatype
 * corresponds to Numpy's npy file format.
 *
 * @param [in] A      matrix to be dumped to file
 * @param [in] stream file handle of opened file
 */
MATRIX_SAVE_TO_STREAM(matrix_complex_save_to_stream, complex_t, matrix_complex_t, "<c16");

#define MATRIX_SAVE(FUNCTION_NAME, MATRIX_TYPE, SAVE_FUNCTION) \
int FUNCTION_NAME(MATRIX_TYPE *M, const char *filename) \
{ \
    FILE *stream = fopen(filename, "w"); \
    if(stream == NULL) \
        return LIBHADES_ERROR_IO; \
\
    SAVE_FUNCTION(M, stream); \
\
    fclose(stream); \
\
    return 0; \
};

/** @brief Save real matrix A to file
 *
 * Save real matrix A to file given by filename. The datatype corresponds to
 * Numpy's .npy file format.
 *
 * @param [in] A        matrix to be dumped to file
 * @param [in] filename path to the file
 * @retval 0 if successful
 * @retval LIBHADES_ERROR_IO if file could not be opened
 */
MATRIX_SAVE(matrix_save, matrix_t, matrix_save_to_stream);

/** @brief Save complex matrix A to file
 *
 * Save complex matrix A to file given by filename. The datatype corresponds to
 * Numpy's .npy file format.
 *
 * @param [in] A        matrix to be dumped to file
 * @param [in] filename path to the file
 * @retval 0 if successful
 * @retval LIBHADES_ERROR_IO if file could not be opened
 */
MATRIX_SAVE(matrix_complex_save, matrix_complex_t, matrix_complex_save_to_stream);

/** @}*/

ext: c

max_stars_repo_name: michael-hartmann/libhades

max_issues_repo_name: michael-hartmann/libhades

max_forks_repo_name: michael-hartmann/libhades


text:
#pragma once

#include <algorithm>
#include <array>
#include <cstdint>
#include <initializer_list>
#include <iterator>
#include <sstream>
#include <string>
#include <tuple>

#include <gsl/gsl>

#include "chainerx/axes.h"
#include "chainerx/constant.h"
#include "chainerx/dtype.h"
#include "chainerx/error.h"
#include "chainerx/macro.h"
#include "chainerx/shape.h"

namespace chainerx {

class Strides : public Dims {
    using BaseVector = Dims;

public:
    using const_iterator = BaseVector::const_iterator;
    using const_reverse_iterator = BaseVector::const_reverse_iterator;

    Strides() = default;

    ~Strides() = default;

    // Creates strides for contiguous array.
    Strides(const Shape& shape, Dtype dtype) : Strides{shape, GetItemSize(dtype)} {}
    Strides(const Shape& shape, int64_t item_size);

    // by iterators
    template <typename InputIt>
    Strides(InputIt first, InputIt last) {
        if (std::distance(first, last) > kMaxNdim) {
            throw DimensionError{"too many dimensions: ", std::distance(first, last)};
        }
        insert(begin(), first, last);
    }

    // by gsl:span
    explicit Strides(gsl::span<const int64_t> dims) : Strides{dims.begin(), dims.end()} {}

    // by initializer list
    Strides(std::initializer_list<int64_t> dims) : Strides{dims.begin(), dims.end()} {}

    // copy
    Strides(const Strides&) = default;
    Strides& operator=(const Strides&) = default;

    // move
    Strides(Strides&&) = default;
    Strides& operator=(Strides&&) = default;

    std::string ToString() const;

    int8_t ndim() const noexcept { return gsl::narrow_cast<int8_t>(size()); }

    const int64_t& operator[](int8_t index) const {
        if (!(0 <= index && static_cast<size_t>(index) < size())) {
            throw DimensionError{"Stride index ", index, " out of bounds for strides with ", size(), " size."};
        }
        return this->StackVector::operator[](index);
    }

    int64_t& operator[](int8_t index) {
        if (!(0 <= index && static_cast<size_t>(index) < size())) {
            throw DimensionError{"Stride index ", index, " out of bounds for strides with ", size(), " size."};
        }
        return this->StackVector::operator[](index);
    }

    // span
    gsl::span<const int64_t> span() const { return {*this}; }

    // Rearranges strides in the order specified by the axes.
    //
    // The size of given axes may be fewer than the size of strides.
    // In that case, new strides will be composed by only given axes.
    //
    // It is the caller's responsibility to ensure validity of permutation.
    // If the permutation is invalid, the behavior is undefined.
    Strides Permute(const Axes& axes) const {
        CHAINERX_ASSERT(axes.size() <= size());
        Strides new_strides{};
        for (int8_t axe : axes) {
            new_strides.emplace_back(operator[](axe));
        }
        return new_strides;
    }
};

std::ostream& operator<<(std::ostream& os, const Strides& strides);

void CheckEqual(const Strides& lhs, const Strides& rhs);

// Returns a pair of lower and upper byte offsets to store the data.
// This forumula always holds: lower <= 0 < item_size <= upper
std::tuple<int64_t, int64_t> GetDataRange(const Shape& shape, const Strides& strides, size_t item_size);

}  // namespace chainerx

ext: h

max_stars_repo_name: prabhatnagarajan/chainer

max_issues_repo_name: prabhatnagarajan/chainer

max_forks_repo_name: prabhatnagarajan/chainer


text:
// my_gsl.h is included in bayesMLR.c
// Contains wrapper funtions
// References:
// https://www.gnu.org/software/gsl/manual/html_node/Matrices.html
// https://github.com/Blei-Lab/diln/blob/master/gsl_wrapper.c
#include <stdio.h>             // standard input/output
#include <stdlib.h>            // malloc
#include <math.h>              // fabs, sqrt, etc.
#include <time.h>              // time / set seed
#include <gsl/gsl_rng.h>       // GNU Scientific Library
#include <gsl/gsl_cdf.h>       // GNU Scientific Library
#include <gsl/gsl_randist.h>   // GNU Scientific Library
#include <gsl/gsl_linalg.h>    // GNU Scientific Library
#include <gsl/gsl_matrix.h>    // GNU Scientific Library
#include <gsl/gsl_statistics.h>// GNU Scientific Library
#include <gsl/gsl_blas.h>      // GSL Basic Linear Algebra
#define pi 3.14159265358979323846

// Reading Matrices ////////////////////////////////////////////////
int countFileRows(char* filename) {
  // Count number of rows in file "filename"
  FILE *fp;
  int count = 0;  // Line counter (result)
  char c;  // To store a character read from file

  fp = fopen(filename, "r");
  // Extract characters from file and store in character c
  if (fp) { // check file existence
  for (c = getc(fp); c != EOF; c = getc(fp))
      if (c == '\n') // Increment count if this character is newline
          count = count + 1;
  }
  fclose(fp);

  return count;
}

int countFileCols(char* filename, char delim) {
  // Count number of columns in file "filename" (with delimiter delim)
  FILE *fp = fopen(filename, "r");
  char c;
  int count = 0;

  while( c=fgetc(fp) ) {
    if(c == EOF) {
      break; /* break if end of file */
    } else if (c == delim) {
      count += 1; /* otherwise add one to the count */
    } else if (c == '\n') {
      count += 1;
      break;
    }
  }
  fclose(fp);

  return count;
}

void read_csv(char* filename, gsl_matrix* m) {
  // Read a csv into matrix m
  FILE* f = fopen(filename, "r");
  gsl_matrix_fscanf(f, m);
  fclose(f);
}

// Printing Matrix Info ////////////////////////////
void print_vector(gsl_matrix* v, char* filename) {
  FILE* f;
  if (filename == "") {
    f = stdout;
    printf("\n");
  } else {
    f = fopen(filename, "w");
  }

  gsl_vector_fprintf(f,v,"%g");
}

void print_matrix(gsl_matrix* m, char* filename) {
  // print matrix to filename or stdout if filename == ""
  FILE* f;
  if (filename == "") {
    f = stdout;
    printf("\n");
  } else {
    f = fopen(filename, "w");
  }
  int n = m->size1;
  int k = m->size2;

  for (int i = 0; i < n; i++) {
    for (int j = 0; j < k; j++) {
      fprintf (f, "%g\t", gsl_matrix_get (m, i, j));
    }
    fprintf (f, "\n");
  }

  if (filename == "") {
    printf("\n");
  } else {
    fclose(f);
  }
}

void print_matrix_dims(gsl_matrix* m) {
  // print dimensions of matrix m
  printf("%d%s%d\n",m->size1,"x",m->size2);
}

// Access Matrix Elements: //////////////////////////////////
void mat_sub(gsl_matrix* X, size_t r1, size_t r2, size_t c1, size_t c2, gsl_matrix* X_new) {
  gsl_matrix_view x = gsl_matrix_submatrix(X,r1,c1,r2-r1+1,c2-c1+1);
  gsl_matrix_memcpy(X_new,&x.matrix);
}

// Random Number Functions: /////////////////////////////////
double runif() {
  // Returns a random number b/w 0 & 1
  return (double) rand() / (double) RAND_MAX;
}

void mvrnorm(gsl_vector* m, gsl_matrix* cholS, gsl_rng* r, gsl_vector* out) {
  // vector addition, vector prod
  // m + cholS * e

  int n = m->size;
  double tmp;
  gsl_vector* e = gsl_vector_alloc(n);

  for (int i=0; i<n; i++) {
    tmp = 0;
    for (int j=0; j<n; j++) {
      if (i==0) gsl_vector_set(e,j,gsl_ran_gaussian(r,1.0));
      tmp += gsl_matrix_get(cholS,i,j) * gsl_vector_get(e,j);
    }
    gsl_vector_set(out,i,gsl_vector_get(m,i) + tmp);
  }

  gsl_vector_free(e);
}

// Linear Algebra Wrappers //////////////////////////////////
// mm: matrix-matrix
// mv: matrix-vector
// ms: matrix-scalar
//  d: double real
// ge: general matrix
void mm_prod(gsl_matrix* A, gsl_matrix* B, double a, double b, gsl_matrix* out) {
}

void mm_add(gsl_matrix* A, gsl_matrix* B, double a, double b, gsl_matrix* out) {
}

void mv_prod(gsl_matrix* A, gsl_vector* x, double a, gsl_vector* out) {
  gsl_blas_dgemv(CblasNoTrans,a,A,x,0.0,out);
}

void vv_add(gsl_vector* x, gsl_vector* y, gsl_vector *out) {
  // sum two vectors
  if (x->size != y->size)
    error("Error in vv_add(): input vectors must have the same dimensionality");
  else if (x->size != out->size || y->size != out->size)
    error("Error in vv_add(): output vector has incorrect dimensions");

  double tmp;
  int n = x->size;

  for (int i=0; i<n; i++) {
    tmp = gsl_vector_get(x,i) + gsl_vector_get(y,i);
    gsl_vector_set(out,i,tmp);
  }
}

double vv_prod(gsl_vector* x, gsl_vector* y) {
  // dot product
  if (x->size != y->size)
    error("Error in vv_prod(): input vectors must have the same dimensionality");

  double out = 0;
  int n = x->size;

  for (int i=0; i<n; i++) {
    out += gsl_vector_get(x,i) * gsl_vector_get(y,i);
  }

  return out;
}

void mat_inv(gsl_matrix * X, gsl_matrix * invX) {
// invert the matrix X
    int signum;
    gsl_matrix * LU = gsl_matrix_calloc(X->size1,X->size2);
    gsl_permutation * p = gsl_permutation_calloc(X->size1);
    gsl_matrix_memcpy(LU,X);

    gsl_linalg_LU_decomp (LU,p,&signum);
    gsl_linalg_LU_invert(LU,p,invX);

    gsl_matrix_free(LU);
    gsl_permutation_free(p);
}

void xxi_m(gsl_matrix* X, gsl_matrix* xxi) {
  int n = X->size1;
  int k = X->size2;
  double dot;
  gsl_matrix* xx = gsl_matrix_alloc(k,k);

  for (int i=0; i<k; i++) {
    for (int j=0; j<=i; j++) {
      dot = 0;
      for (int z=0; z<n; z++) {
        dot += gsl_matrix_get(X,z,i) * gsl_matrix_get(X,z,j);
      }
      gsl_matrix_set(xx,j,i,dot);
      gsl_matrix_set(xx,i,j,dot);
    }
  }

  mat_inv(xx,xxi);
  gsl_matrix_free(xx);
}

int my_gsl_linalg_cholesky_decomp (gsl_matrix * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR("cholesky decomposition requires square matrix", GSL_ENOTSQR);
    }
  else
    {
      size_t i,j,k;
      int status = 0;

      /* Do the first 2 rows explicitly.  It is simple, and faster.  And
       * one can return if the matrix has only 1 or 2 rows.
       */

      double A_00 = gsl_matrix_get (A, 0, 0);

      double L_00 = sqrt(A_00);

      if (A_00 <= 0)
        {
          status = GSL_EDOM ;
        }

      gsl_matrix_set (A, 0, 0, L_00);

      if (M > 1)
        {
          double A_10 = gsl_matrix_get (A, 1, 0);
          double A_11 = gsl_matrix_get (A, 1, 1);

          double L_10 = A_10 / L_00;
          double diag = A_11 - L_10 * L_10;
          double L_11 = sqrt(diag);

          if (diag <= 0)
            {
              status = GSL_EDOM;
            }

          gsl_matrix_set (A, 1, 0, L_10);
          gsl_matrix_set (A, 1, 1, L_11);
        }

      for (k = 2; k < M; k++)
        {
          double A_kk = gsl_matrix_get (A, k, k);

          for (i = 0; i < k; i++)
            {
              double sum = 0;

              double A_ki = gsl_matrix_get (A, k, i);
              double A_ii = gsl_matrix_get (A, i, i);

              gsl_vector_view ci = gsl_matrix_row (A, i);
              gsl_vector_view ck = gsl_matrix_row (A, k);

              if (i > 0) {
                gsl_vector_view di = gsl_vector_subvector(&ci.vector, 0, i);
                gsl_vector_view dk = gsl_vector_subvector(&ck.vector, 0, i);

                gsl_blas_ddot (&di.vector, &dk.vector, &sum);
              }

              A_ki = (A_ki - sum) / A_ii;
              gsl_matrix_set (A, k, i, A_ki);
            }

          {
            gsl_vector_view ck = gsl_matrix_row (A, k);
            gsl_vector_view dk = gsl_vector_subvector (&ck.vector, 0, k);

            double sum = gsl_blas_dnrm2 (&dk.vector);
            double diag = A_kk - sum * sum;

            double L_kk = sqrt(diag);

            if (diag <= 0)
              {
                status = GSL_EDOM;
              }

            gsl_matrix_set (A, k, k, L_kk);
          }
        }

      /* Now copy the transposed lower triangle to the upper triangle,
       * the diagonal is common.
       */

      //for (i = 1; i < M; i++)
      //  {
      //    for (j = 0; j < i; j++)
      //      {
      //        double A_ij = 0;//gsl_matrix_get (A, i, j);
      //        gsl_matrix_set (A, j, i, A_ij);
      //      }
      //  }

      if (status == GSL_EDOM)
        {
          GSL_ERROR ("matrix must be positive definite", GSL_EDOM);
        }

      return GSL_SUCCESS;
    }
}

void chol(gsl_matrix* X, gsl_matrix* cholS) {
  int n = X->size1;
  gsl_matrix_memcpy(cholS,X);
  my_gsl_linalg_cholesky_decomp(cholS);

  // lower tri stored in cholS.
  // So: mvrnorm <- function(M,S,n=nrow(S)) M + chol(S) %*% rnorm(n)
}

max_stars_repo_name: luiarthur/padawan

max_issues_repo_name: luiarthur/padawan

max_forks_repo_name: luiarthur/padawan


text:
/**
 *
 * @file testing_zcungesv.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Emmanuel Agullo
 * @date 2010-11-15
 * @precisions mixed zc -> ds
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

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

static int check_solution(int, int, PLASMA_Complex64_t*, int, PLASMA_Complex64_t*, PLASMA_Complex64_t*, int, double);

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

    int N     = atoi(argv[0]);
    int LDA   = atoi(argv[1]);
    int NRHS  = atoi(argv[2]);
    int LDB   = atoi(argv[3]);
    int ITER;

    double eps;
    int info, info_solution;
    int i,j;
    int LDAxN    = LDA*N;
    int LDBxNRHS = LDB*NRHS;

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

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

    eps = LAPACKE_dlamch_work('e');

    /*----------------------------------------------------------
    *  TESTING ZCGELS
    */

    /* Initialize A1 and A2 */
    LAPACKE_zlarnv_work(IONE, ISEED, LDAxN, A1);
    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
            A2[LDA*j+i] = A1[LDA*j+i] ;

    /* Initialize B1 and B2 */
    LAPACKE_zlarnv_work(IONE, ISEED, LDBxNRHS, B1);
    for (i = 0; i < N; i++)
        for (j = 0; j < NRHS; j++)
             B2[LDB*j+i] = B1[LDB*j+i] ;

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

    /* PLASMA ZCUNGESV */
    info = PLASMA_zcungesv(PlasmaNoTrans, N, NRHS, A2, LDA, B1, LDB, B2, LDB, &ITER);
    if (info != PLASMA_SUCCESS ) {
        printf("PLASMA_zcposv is not completed: info = %d\n", info);
        info_solution = 1;
    } else {
        printf(" Solution obtained with %d iterations\n", ITER);

        /* Check the orthogonality, factorization and the solution */
        info_solution = check_solution(N, NRHS, A1, LDA, B1, B2, LDB, eps);
    }

    if (info_solution == 0) {
        printf("***************************************************\n");
        printf(" ---- TESTING ZCUNGESV.................... PASSED !\n");
        printf("***************************************************\n");
    }
    else {
        printf("************************************************\n");
        printf(" - TESTING ZCUNGESV .. FAILED !\n");
        printf("************************************************\n");
    }

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

    return 0;
}

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

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

    alpha = 1.0;
    beta  = -1.0;

    Anorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N,    A1, LDA, work);
    Xnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, NRHS, B2, LDB, work);
    Bnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, NRHS, B1, LDB, work);

    cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, N, NRHS, N, CBLAS_SADDR(alpha), A1, LDA, B2, LDB, CBLAS_SADDR(beta), B1, LDB);
    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, NRHS, B1, LDB, work);

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

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

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

    free(work);

    return info_solution;
}


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

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_trig.h>
#include <gsl/gsl_sf_expint.h>

#include "error.h"

#include "chebyshev.h"
#include "cheb_eval.c"

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

/* based on SLATEC r9sifg.f, W. Fullerton */

/*
 series for f1   on the interval  2.00000e-02 to  6.25000e-02
                                        with weighted error   2.82e-17
                                         log weighted error  16.55
                               significant figures required  15.36
                                    decimal places required  17.20
*/
static double f1_data[20] = {
   -0.1191081969051363610,
   -0.0247823144996236248,
    0.0011910281453357821,
   -0.0000927027714388562,
    0.0000093373141568271,
   -0.0000011058287820557,
    0.0000001464772071460,
   -0.0000000210694496288,
    0.0000000032293492367,
   -0.0000000005206529618,
    0.0000000000874878885,
   -0.0000000000152176187,
    0.0000000000027257192,
   -0.0000000000005007053,
    0.0000000000000940241,
   -0.0000000000000180014,
    0.0000000000000035063,
   -0.0000000000000006935,
    0.0000000000000001391,
   -0.0000000000000000282
};
static cheb_series f1_cs = {
  f1_data,
  19,
  -1, 1,
  10
};

/*

 series for f2   on the interval  0.00000e+00 to  2.00000e-02
                                        with weighted error   4.32e-17
                                         log weighted error  16.36
                               significant figures required  14.75
                                    decimal places required  17.10
*/
static double f2_data[29] = {
   -0.0348409253897013234,
   -0.0166842205677959686,
    0.0006752901241237738,
   -0.0000535066622544701,
    0.0000062693421779007,
   -0.0000009526638801991,
    0.0000001745629224251,
   -0.0000000368795403065,
    0.0000000087202677705,
   -0.0000000022601970392,
    0.0000000006324624977,
   -0.0000000001888911889,
    0.0000000000596774674,
   -0.0000000000198044313,
    0.0000000000068641396,
   -0.0000000000024731020,
    0.0000000000009226360,
   -0.0000000000003552364,
    0.0000000000001407606,
   -0.0000000000000572623,
    0.0000000000000238654,
   -0.0000000000000101714,
    0.0000000000000044259,
   -0.0000000000000019634,
    0.0000000000000008868,
   -0.0000000000000004074,
    0.0000000000000001901,
   -0.0000000000000000900,
    0.0000000000000000432
};
static cheb_series f2_cs = {
  f2_data,
  28,
  -1, 1,
  14
};

/*

 series for g1   on the interval  2.00000e-02 to  6.25000e-02
                                        with weighted error   5.48e-17
                                         log weighted error  16.26
                               significant figures required  15.47
                                    decimal places required  16.92
*/
static double g1_data[21] = {
   -0.3040578798253495954,
   -0.0566890984597120588,
    0.0039046158173275644,
   -0.0003746075959202261,
    0.0000435431556559844,
   -0.0000057417294453025,
    0.0000008282552104503,
   -0.0000001278245892595,
    0.0000000207978352949,
   -0.0000000035313205922,
    0.0000000006210824236,
   -0.0000000001125215474,
    0.0000000000209088918,
   -0.0000000000039715832,
    0.0000000000007690431,
   -0.0000000000001514697,
    0.0000000000000302892,
   -0.0000000000000061400,
    0.0000000000000012601,
   -0.0000000000000002615,
    0.0000000000000000548
};
static cheb_series g1_cs = {
  g1_data,
  20,
  -1, 1,
  13
};

/*

 series for g2   on the interval  0.00000e+00 to  2.00000e-02
                                        with weighted error   5.01e-17
                                         log weighted error  16.30
                               significant figures required  15.12
                                    decimal places required  17.07
*/
static double g2_data[34] = {
   -0.0967329367532432218,
   -0.0452077907957459871,
    0.0028190005352706523,
   -0.0002899167740759160,
    0.0000407444664601121,
   -0.0000071056382192354,
    0.0000014534723163019,
   -0.0000003364116512503,
    0.0000000859774367886,
   -0.0000000238437656302,
    0.0000000070831906340,
   -0.0000000022318068154,
    0.0000000007401087359,
   -0.0000000002567171162,
    0.0000000000926707021,
   -0.0000000000346693311,
    0.0000000000133950573,
   -0.0000000000053290754,
    0.0000000000021775312,
   -0.0000000000009118621,
    0.0000000000003905864,
   -0.0000000000001708459,
    0.0000000000000762015,
   -0.0000000000000346151,
    0.0000000000000159996,
   -0.0000000000000075213,
    0.0000000000000035970,
   -0.0000000000000017530,
    0.0000000000000008738,
   -0.0000000000000004487,
    0.0000000000000002397,
   -0.0000000000000001347,
    0.0000000000000000801,
   -0.0000000000000000501
};
static cheb_series g2_cs = {
  g2_data,
  33,
  -1, 1,
  20
};

/* x >= 4.0 */
static void fg_asymp(const double x, gsl_sf_result * f, gsl_sf_result * g)
{
  /*
      xbig = sqrt (1.0/r1mach(3))
      xmaxf = exp (amin1(-alog(r1mach(1)), alog(r1mach(2))) - 0.01)
      xmaxg = 1.0/sqrt(r1mach(1))
      xbnd = sqrt(50.0)
  */
  const double xbig  = 1.0/GSL_SQRT_DBL_EPSILON;
  const double xmaxf = 1.0/GSL_DBL_MIN;
  const double xmaxg = 1.0/GSL_SQRT_DBL_MIN;
  const double xbnd  = 7.07106781187;

  const double x2 = x*x;

  if(x <= xbnd) {
    gsl_sf_result result_c1;
    gsl_sf_result result_c2;
    cheb_eval_e(&f1_cs, (1.0/x2-0.04125)/0.02125, &result_c1);
    cheb_eval_e(&g1_cs, (1.0/x2-0.04125)/0.02125, &result_c2);
    f->val = (1.0 + result_c1.val)/x;
    g->val = (1.0 + result_c2.val)/x2;
    f->err = result_c1.err/x  + 2.0 * GSL_DBL_EPSILON * fabs(f->val);
    g->err = result_c2.err/x2 + 2.0 * GSL_DBL_EPSILON * fabs(g->val);
  }
  else if(x <= xbig) {
    gsl_sf_result result_c1;
    gsl_sf_result result_c2;
    cheb_eval_e(&f2_cs, 100.0/x2-1.0, &result_c1);
    cheb_eval_e(&g2_cs, 100.0/x2-1.0, &result_c2);
    f->val = (1.0 + result_c1.val)/x;
    g->val = (1.0 + result_c2.val)/x2;
    f->err = result_c1.err/x  + 2.0 * GSL_DBL_EPSILON * fabs(f->val);
    g->err = result_c2.err/x2 + 2.0 * GSL_DBL_EPSILON * fabs(g->val);
  }
  else {
    f->val = (x < xmaxf ? 1.0/x  : 0.0);
    g->val = (x < xmaxg ? 1.0/x2 : 0.0);
    f->err = 2.0 * GSL_DBL_EPSILON * fabs(f->val);
    g->err = 2.0 * GSL_DBL_EPSILON * fabs(g->val);
  }

  return;
}

/* based on SLATEC si.f, W. Fullerton

 series for si   on the interval  0.00000e+00 to  1.60000e+01
                                        with weighted error   1.22e-17
                                         log weighted error  16.91
                               significant figures required  16.37
                                    decimal places required  17.45
*/

static double si_data[12] = {
  -0.1315646598184841929,
  -0.2776578526973601892,
   0.0354414054866659180,
  -0.0025631631447933978,
   0.0001162365390497009,
  -0.0000035904327241606,
   0.0000000802342123706,
  -0.0000000013562997693,
   0.0000000000179440722,
  -0.0000000000001908387,
   0.0000000000000016670,
  -0.0000000000000000122
};

static cheb_series si_cs = {
  si_data,
  11,
  -1, 1,
  9
};

/*
 series for ci   on the interval  0.00000e+00 to  1.60000e+01
                                        with weighted error   1.94e-18
                                         log weighted error  17.71
                               significant figures required  17.74
                                    decimal places required  18.27
*/
static double ci_data[13] = {
   -0.34004281856055363156,
   -1.03302166401177456807,
    0.19388222659917082877,
   -0.01918260436019865894,
    0.00110789252584784967,
   -0.00004157234558247209,
    0.00000109278524300229,
   -0.00000002123285954183,
    0.00000000031733482164,
   -0.00000000000376141548,
    0.00000000000003622653,
   -0.00000000000000028912,
    0.00000000000000000194
};
static cheb_series ci_cs = {
  ci_data,
  12,
  -1, 1,
  9
};

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int gsl_sf_Si_e(const double x, gsl_sf_result * result)
{
  double ax = fabs(x);

  /* CHECK_POINTER(result) */

  if(ax < GSL_SQRT_DBL_EPSILON) {
    result->val = x;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(ax <= 4.0) {
    gsl_sf_result result_c;
    cheb_eval_e(&si_cs, (x*x-8.0)*0.125, &result_c);
    result->val  =  x * (0.75 + result_c.val);
    result->err  = ax * result_c.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    /* Note there is no loss of precision
     * here bcause of the leading constant.
     */
    gsl_sf_result f;
    gsl_sf_result g;
    fg_asymp(ax, &f, &g);
    result->val  = 0.5 * M_PI - f.val*cos(ax) - g.val*sin(ax);
    result->err  = f.err + g.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    if(x < 0.0) result->val = -result->val;
    return GSL_SUCCESS;
  }
}

int gsl_sf_Ci_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x <= 4.0) {
    const double lx = log(x);
    const double y  = (x*x-8.0)*0.125;
    gsl_sf_result result_c;
    cheb_eval_e(&ci_cs, y, &result_c);
    result->val  = lx - 0.5 + result_c.val;
    result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(lx) + 0.5) + result_c.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result sin_result;
    gsl_sf_result cos_result;
    int stat_sin = gsl_sf_sin_e(x, &sin_result);
    int stat_cos = gsl_sf_cos_e(x, &cos_result);
    gsl_sf_result f;
    gsl_sf_result g;
    fg_asymp(x, &f, &g);
    result->val  = f.val*sin_result.val - g.val*cos_result.val;
    result->err  = fabs(f.err*sin_result.val);
    result->err += fabs(g.err*cos_result.val);
    result->err += fabs(f.val*sin_result.err);
    result->err += fabs(g.val*cos_result.err);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_ERROR_SELECT_2(stat_sin, stat_cos);
  }
}

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

#include "eval.h"

double gsl_sf_Si(const double x)
{
  EVAL_RESULT(gsl_sf_Si_e(x, &result));
}

double gsl_sf_Ci(const double x)
{
  EVAL_RESULT(gsl_sf_Ci_e(x, &result));
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#ifndef __GSL_SORT_H__
#define __GSL_SORT_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <gsl/gsl_sort_long_double.h>
#include <gsl/gsl_sort_double.h>
#include <gsl/gsl_sort_float.h>

#include <gsl/gsl_sort_ulong.h>
#include <gsl/gsl_sort_long.h>

#include <gsl/gsl_sort_uint.h>
#include <gsl/gsl_sort_int.h>

#include <gsl/gsl_sort_ushort.h>
#include <gsl/gsl_sort_short.h>

#include <gsl/gsl_sort_uchar.h>
#include <gsl/gsl_sort_char.h>

#endif /* __GSL_SORT_H__ */

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
// run with OPENBLAS_NUM_THREADS=1 and OMP_NUM_THREADS=n
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <cblas.h>
#include <omp.h>
#include <pthread.h>
#define MIN_SIZE 5
#define MAX_SIZE 60
#define NB_SIZE 10

// number of loop for a 1x1 matrix. Lower it if the test is
// too slow on you computer.
#define NLOOP 2e7

typedef struct {
    int matrix_size;
    int n_loop;
    void (* bench_func)();
    void (* blas_func)();
    void * (* create_matrix)(int size);
} BenchParam;

void * s_create_matrix(int size) {
    float * r = malloc(size * sizeof(double));
    int i;
    for(i = 0; i < size; i++)
        r[i] = 1e3 * i / size;
    return r;
}

void * c_create_matrix(int size) {
    float * r = malloc(size * 2 * sizeof(double));
    int i;
    for(i = 0; i < 2 * size; i++)
        r[i] = 1e3 * i / size;
    return r;
}

void * z_create_matrix(int size) {
    double * r = malloc(size * 2 * sizeof(double));
    int i;
    for(i = 0; i < 2 * size; i++)
        r[i] = 1e3 * i / size;
    return r;
}

void * d_create_matrix(int size) {
    double * r = malloc(size * sizeof(double));
    int i;
    for(i = 0; i < size; i++)
        r[i] = 1e3 * i / size;
    return r;
}

void trmv_bench(BenchParam * param)
{
    int i, n;
    int size = param->matrix_size;
    n = param->n_loop / size;
    int one = 1;
    void * A = param->create_matrix(size * size);
    void * y = param->create_matrix(size);
    for(i = 0; i < n; i++) {
        param->blas_func("U", "N", "N", &size, A, &size, y, &one);
    }
    free(A);
    free(y);
}

void gemv_bench(BenchParam * param)
{
    int i, n;
    int size = param->matrix_size;
    n = param->n_loop / size;
    double v = 1.01;
    int one = 1;
    void * A = param->create_matrix(size * size);
    void * y = param->create_matrix(size);
    for(i = 0; i < n; i++) {
        param->blas_func("N", &size, &size, &v, A, &size, y, &one, &v, y, &one);
    }
    free(A);
    free(y);
}

void ger_bench(BenchParam * param) {
    int i, n;
    int size = param->matrix_size;
    n = param->n_loop / size;
    double v = 1.01;
    int one = 1;
    void * A = param->create_matrix(size * size);
    void * y = param->create_matrix(size);
    for(i = 0; i < n; i++) {
        param->blas_func(&size, &size, &v, y, &one, y, &one, A, &size);
    }
    free(A);
    free(y);
}

#ifndef _WIN32
void * pthread_func_wrapper(void * param) {
    ((BenchParam *)param)->bench_func(param);
    pthread_exit(NULL);
}
#endif

#define NB_TESTS 5
void * TESTS[4 * NB_TESTS] = {
    trmv_bench, ztrmv_, z_create_matrix, "ztrmv",
    gemv_bench, dgemv_, d_create_matrix, "dgemv",
    gemv_bench, zgemv_, z_create_matrix, "zgemv",
    ger_bench, dger_, d_create_matrix, "dger",
    ger_bench, zgerc_, z_create_matrix, "zgerc",
};

inline static double delta_time(struct timespec tick) {
    struct timespec tock;
	clock_gettime(CLOCK_MONOTONIC, &tock);
	return (tock.tv_sec - tick.tv_sec) + (tock.tv_nsec - tick.tv_nsec) / 1e9;
}

double pthread_bench(BenchParam * param, int nb_threads)
{
#ifdef _WIN32
    return 0;
#else
    BenchParam threaded_param = *param;
    pthread_t threads[nb_threads];
    int t, rc;
    struct timespec tick;
    threaded_param.n_loop /= nb_threads;
    clock_gettime(CLOCK_MONOTONIC, &tick);
    for(t=0; t<nb_threads; t++){
        rc = pthread_create(&threads[t], NULL, pthread_func_wrapper, &threaded_param);
        if (rc){
            printf("ERROR; return code from pthread_create() is %d\n", rc);
            exit(-1);
        }
    }
    for(t=0; t<nb_threads; t++){
        pthread_join(threads[t], NULL);
    }
	return delta_time(tick);
#endif
}

double seq_bench(BenchParam * param) {
    struct timespec tick;
    clock_gettime(CLOCK_MONOTONIC, &tick);
    param->bench_func(param);
    return delta_time(tick);
}

double omp_bench(BenchParam * param) {
    BenchParam threaded_param = *param;
    struct timespec tick;
    int t;
    int nb_threads = omp_get_max_threads();
    threaded_param.n_loop /= nb_threads;
    clock_gettime(CLOCK_MONOTONIC, &tick);
    #pragma omp parallel for
    for(t = 0; t < nb_threads; t ++){
        param->bench_func(&threaded_param);
    }
    return delta_time(tick);
}

int main(int argc, char * argv[]) {
    double inc_factor = exp(log((double)MAX_SIZE / MIN_SIZE) / NB_SIZE);
    BenchParam param;
    int test_id;
    printf ("Running on %d threads\n", omp_get_max_threads());
    for(test_id = 0; test_id < NB_TESTS; test_id ++) {
        double size = MIN_SIZE;
        param.bench_func = TESTS[test_id * 4];
        param.blas_func = TESTS[test_id * 4 + 1];
        param.create_matrix = TESTS[test_id * 4 + 2];
        printf("\nBenchmark of %s\n", (char*)TESTS[test_id * 4 + 3]);
        param.n_loop = NLOOP;
        while(size <= MAX_SIZE) {
            param.matrix_size = (int)(size + 0.5);
            double seq_time = seq_bench(&param);
            double omp_time = omp_bench(&param);
            double pthread_time = pthread_bench(&param, omp_get_max_threads());
            printf("matrix size %d, sequential %gs, openmp %gs, speedup %g, "
                   "pthread %gs, speedup %g\n",
                   param.matrix_size, seq_time,
                   omp_time, seq_time / omp_time,
                   pthread_time, seq_time / pthread_time);
            size *= inc_factor;
        }
    }
    return(0);
}

max_stars_repo_name: drhpc/OpenBLAS

max_issues_repo_name: drhpc/OpenBLAS

max_forks_repo_name: drhpc/OpenBLAS


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>
#include "cblas.h"

void
cblas_zaxpy (const int N, const void *alpha, const void *X, const int incX,
             void *Y, const int incY)
{
#define BASE double
#include "source_axpy_c.h"
#undef BASE
}

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet