text:
#include "../include/paralleltt.h"

#include <stdio.h>
#include <stdlib.h>
#include <cblas.h>
#include <lapacke.h>
#include <time.h>

flattening_info* flattening_info_init(const MPI_tensor* ten, int flattening, int iscol, int t_v_block)
{
    // Tensor info
    flattening_info* fi = (flattening_info*) malloc(sizeof(flattening_info));
    int t_d = ten->d;

    int* t_nps = (int*) malloc(t_d*sizeof(int));
    int t_Nblocks = 1;
    for (int ii = 0; ii < t_d; ++ii){
        t_nps[ii] = ten->nps[ii];
        t_Nblocks = t_Nblocks * t_nps[ii];
    }

    int* t_t_block = (int*) malloc(t_d*sizeof(int));
    to_tensor_ind(t_t_block, t_v_block, t_nps, t_d);

    long t_N = 1;
    int* t_t_sizes = (int*) malloc(t_d*sizeof(int));
    int* t_t_index = (int*) malloc(t_d*sizeof(int));
    for (int ii = 0; ii < t_d; ++ii){
        int* partition_ii = ten->partitions[ii];
        t_t_index[ii] = partition_ii[t_t_block[ii]];
        t_t_sizes[ii] = partition_ii[t_t_block[ii]+1] - t_t_index[ii];
        t_N = t_N * t_t_sizes[ii];
    }

    // Flattening info
    int f_d = (iscol) ? flattening : t_d-flattening;

    int offset = (iscol) ? 0 : flattening;
    long f_N = 1;
    int f_Nblocks = 1;
    int* f_nps     = (int*) malloc(f_d * sizeof(int));
    int* f_t_block = (int*) malloc(f_d * sizeof(int));
    int* f_t_index = (int*) malloc(f_d * sizeof(int));
    int* f_t_sizes = (int*) malloc(f_d * sizeof(int));
    for (int ii = 0; ii < f_d; ++ii){
        f_nps[ii]     = t_nps[ii+offset];
        f_t_block[ii] = t_t_block[ii+offset];
        f_t_index[ii] = t_t_index[ii+offset];
        f_t_sizes[ii] = t_t_sizes[ii+offset];
        f_N = f_N * f_t_sizes[ii];
        f_Nblocks = f_Nblocks * f_nps[ii];
    }
    int f_v_block = to_vec_ind(f_t_block, f_nps, f_d);

    // Sketching dimensions info
    int s_d = t_d - f_d;
    long s_N = 1;

    offset = (iscol) ? flattening : 0;
    int* s_nps     = (int*) malloc(s_d * sizeof(int));
    int* s_t_index = (int*) malloc(s_d * sizeof(int));
    int* s_t_sizes = (int*) malloc(s_d * sizeof(int));
    for (int ii = 0; ii < s_d; ++ii){
        s_nps[ii]     = t_nps[ii+offset];
        s_t_index[ii] = t_t_index[ii+offset];
        s_t_sizes[ii] = t_t_sizes[ii+offset];
        s_N = s_N * s_t_sizes[ii];
    }

    // Assigning
    fi->t_d       = t_d;
    fi->t_N       = t_N;
    fi->t_Nblocks = t_Nblocks;
    fi->t_nps     = t_nps;
    fi->t_v_block = t_v_block;
    fi->t_t_block = t_t_block;
    fi->t_t_index = t_t_index;
    fi->t_t_sizes = t_t_sizes;

    fi->flattening = flattening;
    fi->iscol      = iscol;
    fi->f_d        = f_d;
    fi->f_N        = f_N;
    fi->f_Nblocks  = f_Nblocks;
    fi->f_nps      = f_nps;
    fi->f_v_block  = f_v_block;
    fi->f_t_block  = f_t_block;
    fi->f_t_index  = f_t_index;
    fi->f_t_sizes  = f_t_sizes;

    fi->s_d       = s_d;
    fi->s_N       = s_N;
    fi->s_nps     = s_nps;
    fi->s_t_index = s_t_index;
    fi->s_t_sizes = s_t_sizes;

    return fi;
}

void flattening_info_free(flattening_info* fi)
{
    free(fi->t_nps);     fi->t_nps = NULL;
    free(fi->t_t_block); fi->t_t_block = NULL;
    free(fi->t_t_index); fi->t_t_index = NULL;
    free(fi->t_t_sizes); fi->t_t_sizes = NULL;
    free(fi->f_nps);     fi->f_nps = NULL;
    free(fi->f_t_block); fi->f_t_block = NULL;
    free(fi->f_t_index); fi->f_t_index = NULL;
    free(fi->f_t_sizes); fi->f_t_sizes = NULL;
    free(fi->s_nps);     fi->s_nps = NULL;
    free(fi->s_t_index); fi->s_t_index = NULL;
    free(fi->s_t_sizes); fi->s_t_sizes = NULL;

    free(fi);
}

void flattening_info_update(flattening_info* fi, const MPI_tensor* ten, int t_v_block)
{
    // Tensor info
    int t_d = ten->d;
    int* t_nps = fi->t_nps;
    int* t_t_block = fi->t_t_block;
    to_tensor_ind(t_t_block, t_v_block, t_nps, t_d);

    long t_N = 1;
    int* t_t_sizes = fi->t_t_sizes;
    int* t_t_index = fi->t_t_index;
    for (int ii = 0; ii < t_d; ++ii){
        int* partition_ii = ten->partitions[ii];
        t_t_index[ii] = partition_ii[t_t_block[ii]];
        t_t_sizes[ii] = partition_ii[t_t_block[ii]+1] - t_t_index[ii];
        t_N = t_N * t_t_sizes[ii];
    }

    // Flattening info
    int flattening = fi->flattening;
    int iscol = fi->iscol;
    int f_d = fi->f_d;

    int offset = (iscol) ? 0 : flattening;
    long f_N = 1;
    int* f_nps     = fi->f_nps;
    int* f_t_block = fi->f_t_block;
    int* f_t_index = fi->f_t_index;
    int* f_t_sizes = fi->f_t_sizes;
    for (int ii = 0; ii < f_d; ++ii){
        f_nps[ii]     = t_nps[ii+offset];
        f_t_block[ii] = t_t_block[ii+offset];
        f_t_index[ii] = t_t_index[ii+offset];
        f_t_sizes[ii] = t_t_sizes[ii+offset];
        f_N = f_N * f_t_sizes[ii];
    }
    int f_v_block = to_vec_ind(f_t_block, f_nps, f_d);

    // Sketching dimensions info
    int s_d = fi->s_d;
    long s_N = 1;

    offset = (iscol) ? flattening : 0;
    int* s_nps     = fi->s_nps;
    int* s_t_index = fi->s_t_index;
    int* s_t_sizes = fi->s_t_sizes;
    for (int ii = 0; ii < s_d; ++ii){
        s_nps[ii]     = t_nps[ii+offset];
        s_t_index[ii] = t_t_index[ii+offset];
        s_t_sizes[ii] = t_t_sizes[ii+offset];
        s_N = s_N * s_t_sizes[ii];
    }

    // Assigning
    fi->t_d = t_d;
    fi->t_N = t_N;
    fi->t_nps = t_nps;
    fi->t_v_block = t_v_block;
    fi->t_t_block = t_t_block;
    fi->t_t_index = t_t_index;
    fi->t_t_sizes = t_t_sizes;

    fi->flattening = flattening;
    fi->iscol = iscol;
    fi->f_d = f_d;
    fi->f_N       = f_N;
    fi->f_nps     = f_nps;
    fi->f_v_block = f_v_block;
    fi->f_t_block = f_t_block;
    fi->f_t_index = f_t_index;
    fi->f_t_sizes = f_t_sizes;

    fi->s_d = s_d;
    fi->s_N = s_N;
    fi->s_nps = s_nps;
    fi->s_t_index = s_t_index;
    fi->s_t_sizes = s_t_sizes;
}

void flattening_info_f_update(flattening_info* fi, const MPI_tensor* ten, int f_v_block)
{
    // tensor info
    int t_d = fi->t_d;

    long t_N = -1;
    int t_v_block = -1;
    int* t_t_block = fi->t_t_block;
    int* t_t_index = fi->t_t_index;
    int* t_t_sizes = fi->t_t_sizes;
    for (int ii = 0; ii < t_d; ++ii){
        t_t_block[ii] = -1;
        t_t_index[ii] = 0;
        t_t_sizes[ii] = 0;
    }
    fi->t_N = t_N;
    fi->t_v_block = t_v_block;

    // flattening info
    int flattening = fi->flattening;
    int iscol = fi->iscol;
    int f_d = fi->f_d;
    int f_Nblocks = fi->f_Nblocks;
    int* f_nps = fi->f_nps;

    int* f_t_block = fi->f_t_block;
    to_tensor_ind(f_t_block, (long) f_v_block, f_nps, f_d);

    long f_N = 1;
    int* f_t_index = fi->f_t_index;
    int* f_t_sizes = fi->f_t_sizes;
    int offset = (iscol) ? 0 : flattening;
    for (int ii = 0; ii < f_d; ++ii){
        int* partition_ii = ten->partitions[ii + offset];
        f_t_index[ii] = partition_ii[f_t_block[ii]];
        f_t_sizes[ii] = partition_ii[f_t_block[ii]+1] - f_t_index[ii];
        f_N = f_N * f_t_sizes[ii];
    }
    fi->f_N = f_N;
    fi->f_v_block = f_v_block;

    // Sketching dimensions info
    int s_d = fi->s_d;
    long s_N = 0;

    int* s_t_index = fi->s_t_index;
    int* s_t_sizes = fi->s_t_sizes;

    for (int ii = 0; ii < s_d; ++ii){
        s_t_index[ii] = 0;
        s_t_sizes[ii] = 0;
    }
    fi->s_N = s_N;
}

void flattening_info_print(flattening_info* fi)
{
    printf("\n~~~~~~~~~~~~~~~ Flattening Info ~~~~~~~~~~~~~~~\n");
    int t_d = fi->t_d;
    printf("t_d = %d\n", t_d);

    printf("t_N = %ld\n", fi->t_N);

    printf("t_Nblocks = %d\n", fi->t_Nblocks);

    if (t_d>0){
        printf("t_nps = [%d", fi->t_nps[0]);
        for (int ii = 1; ii < t_d; ++ii){printf(", %d", fi->t_nps[ii]); }
        printf("]\n");
    }
    else{ printf("t_nps = []\n"); }

    printf("t_v_block = %d\n", fi->t_v_block);

    if (t_d>0){
        printf("t_t_block = [%d", fi->t_t_block[0]);
        for (int ii = 1; ii < t_d; ++ii){printf(", %d", fi->t_t_block[ii]);}
        printf("]\n");
    }
    else{ printf("t_t_block = []\n"); }

    if (t_d>0){
        printf("t_t_index = [%d", fi->t_t_index[0]);
        for (int ii = 1; ii < t_d; ++ii){printf(", %d", fi->t_t_index[ii]);}
        printf("]\n");
    }
    else{ printf("t_t_index = []\n"); }

    if (t_d>0){
        printf("t_t_sizes = [%d", fi->t_t_sizes[0]);
        for (int ii = 1; ii < t_d; ++ii){printf(", %d", fi->t_t_sizes[ii]);}
        printf("]\n");
    }
    else{ printf("t_t_sizes = []\n"); }

    printf("flattening = %d\n", fi->flattening);

    printf("iscol = %d\n", fi->iscol);

    int f_d = fi->f_d;
    printf("f_d = %d\n", f_d);

    printf("f_N = %ld\n", fi->f_N);

    printf("f_Nblocks = %d\n", fi->f_Nblocks);

    if (f_d>0){
        printf("f_nps = [%d", fi->f_nps[0]);
        for (int ii = 1; ii < f_d; ++ii){printf(", %d", fi->f_nps[ii]);}
        printf("]\n");
    }
    else{ printf("f_nps = []\n"); }

    printf("f_v_block = %d\n", fi->f_v_block);

    if (f_d>0){
        printf("f_t_block = [%d", fi->f_t_block[0]);
        for (int ii = 1; ii < f_d; ++ii){printf(", %d", fi->f_t_block[ii]);}
        printf("]\n");
    }
    else{ printf("f_t_block = []\n"); }

    if (f_d>0){
        printf("f_t_index = [%d", fi->f_t_index[0]);
        for (int ii = 1; ii < f_d; ++ii){printf(", %d", fi->f_t_index[ii]);}
        printf("]\n");
    }
    else{ printf("f_t_index = []\n"); }

    if (f_d>0){
        printf("f_t_sizes = [%d", fi->f_t_sizes[0]);
        for (int ii = 1; ii < f_d; ++ii){printf(", %d", fi->f_t_sizes[ii]);}
        printf("]\n");
    }
    else{ printf("f_t_sizes = []\n"); }

    int s_d = fi->s_d;
    printf("s_d = %d\n", s_d);

    printf("s_N = %ld\n", fi->s_N);

    if (s_d>0){
        printf("s_nps = [%d", fi->s_nps[0]);
        for (int ii = 1; ii < s_d; ++ii){printf(", %d", fi->s_nps[ii]);}
        printf("]\n");
    }
    else{ printf("s_nps = []\n"); }

    if (s_d>0){
        printf("s_t_index = [%d", fi->s_t_index[0]);
        for (int ii = 1; ii < s_d; ++ii){printf(", %d", fi->s_t_index[ii]);}
        printf("]\n");
    }
    else{ printf("s_t_index = []\n"); }

    if (s_d>0){
        printf("s_t_sizes = [%d", fi->s_t_sizes[0]);
        for (int ii = 1; ii < s_d; ++ii){printf(", %d", fi->s_t_sizes[ii]);}
        printf("]\n");
    }
    else{ printf("s_t_sizes = []\n"); }
    printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n");
}

int* get_sketch_owners(MPI_tensor* ten, int flattening, int iscol)
{
    int np_flattening = 1;
    int* nps = ten->nps;
    int ii0 = (iscol) ? 0 : flattening;
    int ii1 = (iscol) ? flattening : ten->d;

    for (int ii = ii0; ii < ii1; ++ii){
        np_flattening = np_flattening * nps[ii];
    }

    int* owners = get_partition(ten->comm_size, np_flattening);
    return owners;
}

void get_sketch_height(MPI_tensor* ten, int* owners, int flattening, int iscol, long* X_height_ptr, long* buf_height_ptr)
{
    *X_height_ptr = 0;
    *buf_height_ptr = 0;

    flattening_info* fi = flattening_info_init(ten, flattening, iscol, 0);

    for (int rank = 0; rank < ten->comm_size; ++rank){
        long X_height_tmp = 0;
        for (int jj = owners[rank]; jj < owners[rank + 1]; ++jj){
            flattening_info_f_update(fi, ten, jj);
            X_height_tmp = X_height_tmp + fi->f_N;
            *buf_height_ptr = (fi->f_N > *buf_height_ptr) ? fi->f_N : *buf_height_ptr;
        }

        *X_height_ptr = (X_height_tmp > *X_height_ptr) ? X_height_tmp : *X_height_ptr;
    }

    flattening_info_free(fi);
}

matrix_tt** get_sketch_Omega(const MPI_tensor* ten, int flattening, int r, int buf, int iscol)
{
    // Seed random numbers
    double t = (double) time(NULL);
    MPI_Bcast(&t, 1, MPI_DOUBLE, 0, ten->comm);
    srand((unsigned long) t);

    int ii0 = iscol ? flattening : 0;
    int n_Omega = iscol ? (ten->d) - flattening : flattening;

    matrix_tt** Omegas = (matrix_tt**) malloc(n_Omega * sizeof(matrix_tt*));
    for (int ii = 0; ii < n_Omega; ++ii){
        Omegas[ii] = matrix_tt_init(ten->n[ii+ii0], r+buf);
        matrix_tt_dlarnv(Omegas[ii]);
    }

    return Omegas;
}

void get_KR_info(int *d_KR, int *stride_KR, matrix_tt** KRs, MPI_tensor* ten, int flattening, int r, int buf, int iscol)
{
    int max_m = 1000;

    int assigned = 0;
    *d_KR = 0;
    int N_KR = 1;

    int ii0 = iscol ? flattening : 0;
    int n_Omega = iscol ? (ten->d) - flattening : flattening;

    KRs[2] = matrix_tt_init(1, r+buf);

    for (int ii = 0; ii < n_Omega; ++ii){
        if (assigned == 0){
            int* partition_ii = ten->partitions[ii + ii0];
            int sz_ii = 0;
            for (int jj = 0; jj < ten->nps[ii+ii0]; ++jj){
                int candidate_n = partition_ii[jj+1] - partition_ii[jj];
                sz_ii = (sz_ii > candidate_n) ? sz_ii : candidate_n;
            }
            if (N_KR * sz_ii < max_m){
                *d_KR = ii+1;
                N_KR = N_KR * sz_ii;
            }
            else{
                *stride_KR = max_m / N_KR;
                KRs[0] = matrix_tt_init(N_KR, r+buf);
                KRs[1] = matrix_tt_init(*stride_KR, r+buf);
                N_KR = N_KR * (*stride_KR);
                KRs[3] = matrix_tt_init(N_KR, r+buf);
                assigned = 1;
            }
        }
    }
    if (assigned == 0){
        *stride_KR = 0;
        KRs[0] = matrix_tt_init(N_KR, r+buf);
        KRs[1] = matrix_tt_init(1, r+buf);
        KRs[3] = matrix_tt_init(N_KR, r+buf);
    }
}

sketch* sketch_init(MPI_tensor* ten, int flattening, int r, int buf, int iscol)
{
    sketch* s = (sketch*) malloc(sizeof(sketch));

    s->ten = ten;
    s->flattening = flattening;
    s->r = r;
    s->buf = buf;
    s->iscol = iscol;
    s->owner_partition = get_sketch_owners(ten, flattening, iscol);
    long recv_buf_height = 0;
    get_sketch_height(ten, s->owner_partition, flattening, iscol, &(s->lda), &recv_buf_height);
    s->Omegas = get_sketch_Omega(ten, flattening, r, buf, iscol);
    s->X_size = (s->lda)*(r+buf);
    s->X = (double*) calloc(s->X_size, sizeof(double));
    s->scratch = (double*) calloc(s->X_size, sizeof(double));
    s->recv_buf = (double*) calloc(recv_buf_height * (r+buf), sizeof(double));
    s->fi = flattening_info_init(ten, flattening, iscol, 0);
    s->KRs = (matrix_tt**) calloc(4, sizeof(matrix_tt*));
    get_KR_info(&(s->d_KR), &(s->stride_KR), s->KRs, ten, flattening, r, buf, iscol);

    return s;
}

matrix_tt** copy_sketch_Omega(matrix_tt** Omegas, MPI_tensor* ten, int flattening, int r, int buf, int iscol)
{
    int ii0 = iscol ? flattening : 0;
    int n_Omega = iscol ? (ten->d) - flattening : flattening;

    matrix_tt** Omegas_cp = (matrix_tt**) malloc(n_Omega * sizeof(matrix_tt*));
    for (int ii = 0; ii < n_Omega; ++ii){
        submatrix_update(Omegas[ii], 0, ten->n[ii+ii0], 0, r+buf);
        Omegas_cp[ii] = matrix_tt_copy(Omegas[ii]);
    }

    return Omegas_cp;
}

sketch* sketch_init_with_Omega(MPI_tensor* ten, int flattening, int r, int buf, int iscol, matrix_tt** Omegas)
{
    sketch* s = (sketch*) malloc(sizeof(sketch));
    s->ten = ten;
    s->flattening = flattening;
    s->r = r;
    s->buf = buf;
    s->iscol = iscol;
    s->owner_partition = get_sketch_owners(ten, flattening, iscol);
    long recv_buf_height;
    get_sketch_height(ten, s->owner_partition, flattening, iscol, &(s->lda), &recv_buf_height);
    s->Omegas = copy_sketch_Omega(Omegas, ten, flattening, r, buf, iscol);
    s->X_size = (s->lda)*(r+buf);
    s->X = (double*) calloc(s->X_size, sizeof(double));
    s->scratch = (double*) calloc(s->X_size, sizeof(double));
    s->recv_buf = (double*) calloc(recv_buf_height * (r+buf), sizeof(double));
    s->fi = flattening_info_init(ten, flattening, iscol, 0);
    s->KRs = (matrix_tt**) calloc(4, sizeof(matrix_tt*));
    get_KR_info(&(s->d_KR), &(s->stride_KR), s->KRs, ten, flattening, r, buf, iscol);

    return s;
}

// NOTE: does not free the MPI_tensor. We assume this is shared, and should not be freed this way.
void sketch_free(sketch* s)
{
    MPI_tensor* ten = s->ten; s->ten = NULL;
    int d = ten->d;
    int n_Omega = s->iscol ? d - s->flattening : s->flattening;
    for (int ii = 0; ii < n_Omega; ++ii){
        matrix_tt_free(s->Omegas[ii]); s->Omegas[ii] = NULL;
    }
    for (int ii = 0; ii < 4; ++ii){
        matrix_tt_free(s->KRs[ii]); s->KRs[ii] = NULL;
    }
    free(s->KRs); s->KRs = NULL;

    free(s->owner_partition); s->owner_partition = NULL;
    free(s->X);               s->X = NULL;
    free(s->scratch);         s->scratch = NULL;
    free(s->Omegas);          s->Omegas = NULL;
    free(s->recv_buf);        s->recv_buf = NULL;
    flattening_info_free(s->fi); s->fi = NULL;
    free(s);
}

void sketch_print(sketch* s)
{
    printf("Printing sketch at %p\n", s);

    printf("\nten address: %p\n", s->ten);
    printf("flattening = %d\n",s->flattening);
    printf("r = %d\n", s->r);
    printf("buf = %d\n", s->buf);
    printf("iscol = %d\n", s->iscol);

    int* op = s->owner_partition;
    MPI_tensor* ten = s->ten;
    printf("owner_partition = [%d", op[0]);
    for (int ii = 1; ii < (ten->comm_size) + 1; ++ii){
        printf(", %d", op[ii]);
    }
    printf("]\n");

    printf("lda = %ld\n", s->lda);
    printf("X_size = %ld\n", s->X_size);

    printf("\nX = \n");
    int r = s->r; int buf = s->buf; long X_size = s->X_size; long lda = s->lda;
    matrix_tt* X_mat = matrix_tt_wrap(lda, r+buf, s->X);
    matrix_tt_print(X_mat, 1);
    free(X_mat); X_mat = NULL;

    printf("\nscratch = \n");
    matrix_tt* scratch_mat = matrix_tt_wrap(lda, r+buf, s->scratch);
    matrix_tt_print(scratch_mat, 1);
    free(scratch_mat); scratch_mat = NULL;

    int n_Omega = s->iscol ? (s->ten)->d-s->flattening : s->flattening;
    for (int ii = 0; ii < n_Omega; ++ii){
        printf("\nOmegas[%d] = \n", ii);
        matrix_tt_print(s->Omegas[ii], 1);
    }
}

// Performs C[ii + n_ii * jj,kk] = A[ii, kk] * B[jj, kk]
void submatrix_khatri_rao_outer_product(matrix_tt* A, matrix_tt* B, matrix_tt* C)
{
    for (int kk = 0; kk < A->n; ++kk){
        int A_offset = A->offset + kk*(A->lda);
        int B_offset = B->offset + kk*(B->lda);
        int C_offset = C->offset + kk*(C->lda);

        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                A->m, B->m, 1,
                1.0,
                A->X + A_offset, A->m,
                B->X + B_offset, 1,
                0.0,
                C->X + C_offset, A->m);
    }
}

// General idea of the algorithm:
// Split the Omegas into 3 parts: ii < d_KR, ii == d_KR, ii > d_KR
// For ii <  d_KR, we will pre-multiply the sketch matrices Omega[ii] (Call this KR_1)
// For ii == d_KR, we will loop through blocks of Omega[d_KR]
// For ii > d_KR, we will loop and multiply to get a row vector       (Call this KR_3)

// Then, at each point in the loop, we will perform
//     KR_2 = Omegas[ii] * KR_3 and
//     KR_4 = KR_1 * KR_2,
// where * is the Khatri-Rao product. Finally, we multiply X_mat by KR_4 to update the sketch, and repeat.
void subtensor_khatri_rao(sketch* s, matrix_tt* C, flattening_info* fi, double beta, matrix_tt* X_mat)
{
    int r = s->r + s->buf;
    matrix_tt** Omegas = s->Omegas;

    MPI_tensor* ten = s->ten;
    int d_KR = s->d_KR;
    matrix_tt* KR_1;
    matrix_tt* KR_4;
    if ((d_KR == 0) || (d_KR%2 == 1)){
        KR_1 = s->KRs[0];
        KR_4 = s->KRs[3];
    }
    else{
        KR_1 = s->KRs[3];
        KR_4 = s->KRs[0];
    }
    int s_d = fi->s_d;

    int N_kk = 1;

    int h = 0;

    int Omega_offset_KR;

//    printf("Starting loop to get KR_1\n");
    for (int ii = 0; ii < s_d; ++ii){
        int ii0 = fi->s_t_index[ii];
        int ii1 = ii0 + fi->s_t_sizes[ii];

        submatrix_update(Omegas[ii], ii0, ii1, 0, r);

        if (ii > d_KR){
            N_kk = N_kk * fi->s_t_sizes[ii]; // Number of elements in outer loop
        }
        else if (ii == d_KR){
            Omega_offset_KR = ii0;
        }

        if ((ii == 0) && (d_KR > 0)){
            matrix_tt_reshape(ii1-ii0, r, KR_1);
            matrix_tt_copy_data(KR_1, Omegas[0]);
            h = ii1 - ii0;
        }
        else if (ii < d_KR){ // Multiply inner matrices
            matrix_tt* tmp = KR_1;
            KR_1 = KR_4;
            KR_4 = tmp;

            h = h*(ii1-ii0);
            matrix_tt_reshape(h, r, KR_1);
            submatrix_khatri_rao_outer_product(KR_4, Omegas[ii], KR_1);
        }
    }
//    printf("Finished loop to get KR_1\n");

    if (fi->iscol){
        matrix_tt_wrap_update(X_mat, fi->f_N, fi->s_N, get_X(ten));
    }
    else{
        matrix_tt_wrap_update(X_mat, fi->s_N, fi->f_N, get_X(ten));
    }
    X_mat->transpose = (fi->iscol ? 0 : 1);

    if (d_KR == s_d){ // If we have done everything, just multiply (the sketch is relatively small)
        matrix_tt_dgemm(X_mat, KR_1, C, 1.0, beta);
        return;
    }

    int stride_KR = s->stride_KR;
    matrix_tt* KR_2 = s->KRs[1];
    matrix_tt* KR_3 = s->KRs[2];

    int* t_kk = ten->t_kk;
    int size_KR = fi->s_t_sizes[d_KR];
    int N_ll = 1 + (size_KR - 1) / stride_KR;

    int tensor_offset = 0;

    for (int kk = 0; kk < N_kk; ++kk) {
        to_tensor_ind(t_kk, (long) kk, fi->s_t_sizes + d_KR + 1, s_d - 1 - d_KR);

        for (int ii = d_KR + 1; ii < s_d; ++ii){
            if (ii == d_KR + 1){
                for (int jj = 0; jj < r; ++jj){
                    KR_3->X[jj] = matrix_tt_element(Omegas[ii], t_kk[ii-d_KR-1], jj);
                }
            }
            else{
                for (int jj = 0; jj < r; ++jj){
                    KR_3->X[jj] *= matrix_tt_element(Omegas[ii], t_kk[ii-d_KR-1], jj);
                }
            }
        }

        for (int ll = 0; ll < N_ll; ++ll){
            double bb = ((kk==0) && (ll==0)) ? beta : 1.0;
            int ii0 = Omega_offset_KR + ll * stride_KR;
            int ii1 = ((ll+1)*stride_KR > size_KR) ? size_KR + Omega_offset_KR: (ll+1)*stride_KR + Omega_offset_KR;

            submatrix_update(Omegas[d_KR], ii0, ii1, 0, r);
            if (d_KR + 1 == s_d){ // If we only need KR_1 and Omegas[ii]
                matrix_tt_reshape((KR_1->m) * (ii1-ii0), r, KR_4);
                submatrix_khatri_rao_outer_product(KR_1, Omegas[d_KR], KR_4);
            }
            else if (d_KR == 0){ // If we only need Omegas[ii] and KR_3
                matrix_tt_reshape(Omegas[d_KR]->m, r, KR_4);
                submatrix_khatri_rao_outer_product(Omegas[d_KR], KR_3, KR_4);
            }
            else{ // We need everything...
                matrix_tt_reshape(ii1-ii0, r, KR_2);
                submatrix_khatri_rao_outer_product(Omegas[d_KR], KR_3, KR_2);
                matrix_tt_reshape((ii1-ii0) * (KR_1->m), r, KR_4);
                submatrix_khatri_rao_outer_product(KR_1, KR_2, KR_4);
            }

            if (fi->iscol) {
                submatrix_update(X_mat, 0, fi->f_N, tensor_offset, tensor_offset + KR_4->m);
            }
            else {
                submatrix_update(X_mat, tensor_offset, tensor_offset + KR_4->m, 0, fi->f_N);
            }

            matrix_tt_dgemm(X_mat, KR_4, C, 1.0, bb);
            tensor_offset = tensor_offset + KR_4->m;
        }
    }
}

int get_owner(int block, int* owner_partition, int comm_size)
{
    int color = -1;
    for (int ii = 0; ii < comm_size; ++ii){
        if ((block >= owner_partition[ii]) && (block < owner_partition[ii+1])){
            color = ii;
        }
    }

    return color;
}

int s_get_owner(sketch* s, int f_v_block){
    int* op = s->owner_partition;
    MPI_tensor* ten = s->ten;
    int size = ten->comm_size;

    // Binary search (basically wikipedia)
    int a = 0;
    int b = size-1;
    int c;

    while (a <= b){
        c = (a+b)/2;
        if (f_v_block < op[c]){
            b = c-1;
        }
        else if(f_v_block >= op[c+1]){
            a = c+1;
        }
        else{
            return c;
        }
    }

    return -1;
}

void own_submatrix_update(matrix_tt* mat, sketch* s, int f_v_block, int with_buf){
    if (f_v_block == -1){
        return;
    }

    int sketch_offset = 0;
    int sketch_lda = s->lda;
    MPI_tensor* ten = s->ten;
    int world_rank = ten->rank;
    int* owner_partition = s->owner_partition;

    flattening_info* fi = s->fi;
    if ((f_v_block < owner_partition[world_rank]) || (f_v_block >= owner_partition[world_rank+1])){
        printf("r%d own_submatrix: This core does not own f_v_block = %d\n", world_rank, f_v_block);
    }

    for (int block = s->owner_partition[world_rank]; block < f_v_block; ++block){
        flattening_info_f_update(fi, ten, block);
        sketch_offset = sketch_offset + fi->f_N;
    }
    flattening_info_f_update(fi, ten, f_v_block);

    matrix_tt_wrap_update(mat, sketch_lda, s->r + s->buf, s->X);
    int r = (with_buf) ? s->r + s->buf : s->r;
    submatrix_update(mat, sketch_offset, sketch_offset + fi->f_N, 0, r);
}

matrix_tt* own_submatrix(sketch* s, int f_v_block, int with_buf){
    matrix_tt* mat = (matrix_tt*) malloc(sizeof(matrix_tt));
    own_submatrix_update(mat, s, f_v_block, with_buf);
    return mat;
}

void subtensor_sketch_multiply(sketch* s, flattening_info* fi, matrix_tt** holders)
{
    MPI_tensor* ten = s->ten;
    int* owner_partition = s->owner_partition;
    int world_size = ten->comm_size;
    int world_rank = ten->rank;

    // Actually sketch the thing
    matrix_tt* sketch_mat = holders[0];
    double beta = 0.0;
    if (ten->current_part != -1){
        int current_color = get_owner(fi->f_v_block, owner_partition, world_size);
        if (current_color == world_rank){
            own_submatrix_update(sketch_mat, s, fi->f_v_block, 1);
            beta = 1.0;
        }
        else{
            matrix_tt_wrap_update(sketch_mat, fi->f_N, s->r + s->buf, s->scratch);
            beta = 0.0;
        }

        subtensor_khatri_rao(s, sketch_mat, fi, beta, holders[1]);
    }
}

void subtensor_sketch_communicate(sketch* s, int stream_step, flattening_info* fi, matrix_tt** holders)
{
    MPI_tensor* ten = s->ten;
    int* owner_partition = s->owner_partition;
    int world_size = ten->comm_size;
    int world_rank = ten->rank;
    int* group_ranks = ten->group_ranks;

    flattening_info* fi_tmp = s->fi;
    matrix_tt* sketch_mat = holders[0];

    // For each flattening block
    for (int ii = 0; ii < fi->f_Nblocks; ++ii){
        int kk = 0;
        int owner_ii = get_owner(ii, owner_partition, world_size);
        int group_owner;
        // For each rank in the world
        for (int jj = 0; jj < world_size; ++jj){
            // If the rank owns the block
            if (owner_ii == jj){
                // Add it to the group
                group_ranks[kk] = jj;
                kk = kk+1;
                group_owner = jj;
            }
            else{
                int* schedule_jj = ten->schedule[jj];
                int stream_jj = schedule_jj[stream_step];
                if (stream_jj != -1){
                    flattening_info_update(fi_tmp, ten, schedule_jj[stream_step]);
                    // Else, if the rank just calculated something that adds to the block
                    if (ii == fi_tmp->f_v_block){
                        // Add it to the group
                        group_ranks[kk] = jj;
                        kk = kk+1;
                    }
                }
            }
        }

        if (world_rank == owner_ii){
            matrix_tt* comm_matrix = holders[1];
            own_submatrix_update(comm_matrix, s, ii, 1);
            matrix_tt* buf_mat = holders[2];
            matrix_tt_wrap_update(buf_mat, comm_matrix->m, comm_matrix->n, s->recv_buf);

            matrix_tt_group_reduce(ten->comm, world_rank, comm_matrix, buf_mat, owner_ii, group_ranks, kk);
        }
        else if (sketch_mat->X){
            matrix_tt* buf_mat = holders[2];
            matrix_tt_wrap_update(buf_mat, sketch_mat->m, sketch_mat->n, s->recv_buf);
            matrix_tt_group_reduce(ten->comm, world_rank,sketch_mat, buf_mat, owner_ii, group_ranks, kk);
        }

    }
}

void multi_perform_sketch(sketch** sketches, int n_sketch)
{
    int n_holders = 3;
    matrix_tt** holders = (matrix_tt**) malloc(n_holders*n_sketch*sizeof(matrix_tt*));
    for (int ii = 0; ii < n_holders*n_sketch; ++ii){
        holders[ii] = (matrix_tt*) calloc(1, sizeof(matrix_tt));
    }

    MPI_tensor* ten = sketches[0]->ten; // Tensor
    int rank = ten->rank;
    int* schedule_rank = ten->schedule[rank];

    flattening_info** fis = (flattening_info**) calloc(n_sketch, sizeof(flattening_info*));
    for (int ii = 0; ii < n_sketch; ++ii){
        fis[ii] = flattening_info_init(ten, sketches[ii]->flattening, sketches[ii]->iscol, 0);
    }

    for (int ii = 0; ii < ten->n_schedule; ++ii){
        stream(ten, schedule_rank[ii]);
        for (int jj = 0; jj < n_sketch; ++jj){
            matrix_tt** holders_jj = holders + jj*n_holders;
            flattening_info_update(fis[jj], ten, ten->current_part);
            holders_jj[0]->X = NULL;
            subtensor_sketch_multiply(sketches[jj], fis[jj], holders_jj);
        }

        for (int jj = 0; jj < n_sketch; ++jj){
            matrix_tt** holders_jj = holders + jj*n_holders;
            subtensor_sketch_communicate(sketches[jj], ii, fis[jj], holders_jj);
        }
    }

    for (int ii = 0; ii < n_holders*n_sketch; ++ii){
        free(holders[ii]); holders[ii] = NULL;
    }

    for (int ii = 0; ii < n_sketch; ++ii){
        flattening_info_free(fis[ii]); fis[ii] = NULL;
    }
    free(fis); fis = NULL;
    free(holders); holders = NULL;
}

void sketch_qr(sketch* sketch)
{

    MPI_tensor* ten = sketch->ten;
    int* owner_partition = sketch->owner_partition;
    int flattening = sketch->flattening;
    int iscol = sketch->iscol;
    int r = sketch->r;
    int buf = sketch->buf;

    MPI_Comm comm = ten->comm;
    int size = ten->comm_size;
    int rank = ten->rank;

    // A little pre-processing
    int head = 0;
    flattening_info* fi = flattening_info_init(ten, flattening, iscol, 0);
    int* Ns = (int*) calloc(size, sizeof(int));

    for (int ii = 0; ii < size; ++ii){
        for (int jj = owner_partition[ii]; jj < owner_partition[ii+1]; ++jj){
            flattening_info_f_update(fi, ten, jj);
            Ns[ii] = Ns[ii] + fi->f_N;
        }
    }

    int N_rank = Ns[rank];
    int lda = sketch->lda;
    matrix_tt* Q_big = matrix_tt_wrap(lda, r+buf, sketch->X);
    matrix_tt* Q = submatrix(Q_big, 0, N_rank, 0, r+buf);
    matrix_tt* Q_head = NULL;
    matrix_tt* R = NULL;

    if (rank == head){
        Q_head = matrix_tt_init(size * (r+buf), r+buf);
        R = submatrix(Q_head, 0, r+buf, 0, r+buf);
    }
    else{
        R = matrix_tt_init(r+buf, r+buf);
    }

    matrix_tt_truncated_qr(Q, R, r+buf);

    // Gather the Rs
    for (int jj = 0; jj < r+buf; ++jj){
//        printf("r%d jj%d\n", rank, jj);
        if (rank == head){
            double* col = Q_head->X + jj*(Q_head->lda);
            MPI_Gather(col, r+buf, MPI_DOUBLE, col, r+buf, MPI_DOUBLE, head, comm);
        }
        else{
            MPI_Gather(R->X + jj*(R->lda), r+buf, MPI_DOUBLE, NULL, r+buf, MPI_DOUBLE, head, comm);
        }
    }

    // Take the QR of the Rs
    if (rank == head){
        matrix_tt_truncated_qr(Q_head, NULL, r+buf);
    }

    // Scatter the Q of the preceding step
    for (int jj = 0; jj < r+buf; ++jj){
        if (rank == head){
            double* col = Q_head->X + jj*(Q_head->lda);
            MPI_Scatter(col, r+buf, MPI_DOUBLE, col, r+buf, MPI_DOUBLE, head, comm);
        }
        else{
            MPI_Scatter(NULL, r+buf, MPI_DOUBLE, R->X + jj*(R->lda), r+buf, MPI_DOUBLE, head, comm);
        }
    }

    // Multiply to get the final Q
    matrix_tt* X_big = matrix_tt_wrap(lda, r, sketch->scratch);
    matrix_tt* new_X = submatrix(X_big, 0, N_rank, 0, r);
    matrix_tt* Q_head_sub = submatrix(R, 0, r+buf, 0, r);
    matrix_tt_dgemm(Q, Q_head_sub, new_X, 1.0, 0.0);

    // Switch so the QR lives in X
    double* tmp = sketch->scratch;
    sketch->scratch = sketch->X;
    sketch->X = tmp;

    free(X_big);              X_big = NULL;
    free(new_X);              new_X = NULL;
    free(Q_head_sub);         Q_head_sub = NULL;
    flattening_info_free(fi); fi = NULL;
    free(Ns);                 Ns = NULL;
    free(Q_big);              Q_big = NULL;
    free(Q);                  Q = NULL;
    if (rank == head){
        matrix_tt_free(Q_head); Q_head = NULL;
        free(R);             R = NULL;
    }
    else{
        matrix_tt_free(R); R = NULL;
    }
}

void perform_sketch(sketch* s)
{
    multi_perform_sketch(&s, 1);
}

// Gets the sketch block from the correct owner (stored in fi->f_v_block). Returns a matrix with the correct dimensions
void sendrecv_sketch_block(matrix_tt* mat, sketch* s, flattening_info* fi, int recv_rank, int with_buf)
{
    int f_v_block = fi->f_v_block;

    MPI_tensor* ten = s->ten;
    int rank = ten->rank;
    MPI_Comm comm = ten->comm;

    int owner = s_get_owner(s, f_v_block);
    int r = (with_buf) ? s->r + s->buf : s->r;

    if (rank == recv_rank){
        if (rank == owner){
            own_submatrix_update(mat, s, f_v_block, with_buf);
        }
        else{
            matrix_tt_wrap_update(mat, fi->f_N, r, s->scratch);
            matrix_tt_recv(comm, mat, owner);
        }
    }
    else{
        if(rank == owner){
            own_submatrix_update(mat, s, f_v_block, with_buf);
            matrix_tt_send(comm, mat, recv_rank);
        }
        mat->X = NULL;
    }
}

// Eats s, so be careful!
MPI_tensor* sketch_to_tensor(sketch** s_ptr)
{
    sketch* s = *s_ptr;

    MPI_tensor* sten = (MPI_tensor*) malloc(sizeof(MPI_tensor));
    MPI_tensor* ten = s->ten;
    int rank = ten->rank;
    int size = ten->comm_size;

    flattening_info* fi = flattening_info_init(ten, s->flattening, s->iscol, 0);

    // Get the schedule from the owner_partition
    int n_schedule = 0;
    int* owners = s->owner_partition;
    for (int ii = 0; ii < size; ++ii){
        int sz = owners[ii+1] - owners[ii];
        n_schedule = (n_schedule > sz) ? n_schedule : sz;
    }

    int** schedule = (int**) malloc(size * sizeof(int*));
    int* inverse_schedule = (int*) malloc(fi->f_Nblocks * sizeof(int));
    for (int ii = 0; ii < size; ++ii){
        schedule[ii] = (int*) malloc(n_schedule * sizeof(int));

        int* schedule_ii = schedule[ii];
        int jj0 = owners[ii];
        int sz = owners[ii+1] - jj0;
        for (int jj = 0; jj < n_schedule; ++jj){
            if (jj < sz){
                schedule_ii[jj] = jj+jj0;
                inverse_schedule[jj+jj0] = ii*n_schedule + jj;
            }
            else{
                schedule_ii[jj] = -1;
            }
        }
    }

    int d = fi->f_d + 1;
    int soffset = (s->iscol) ? 0 : 1;
    int offset = (s->iscol) ? 0 : s->flattening;
    int end_ind = (s->iscol) ? d-1 : 0;

    // Inherit the properties of ten
    int* n = (int*) malloc(d * sizeof(int));
    int* nps = (int*) malloc(d * sizeof(int));
    int** partitions = (int**) malloc(d * sizeof(int*));
    for (int ii = 0; ii < d-1; ++ii){
        n[ii + soffset] = ten->n[ii + offset];
        nps[ii + soffset] = ten->nps[ii + offset];
        partitions[ii + soffset] = (int*) malloc( (nps[ii+soffset] + 1) * sizeof(int));

        int* spartition_ii = partitions[ii+soffset];
        int* partition_ii = ten->partitions[ii+offset];
        for (int jj = 0; jj < nps[ii+soffset] + 1; ++jj){
            spartition_ii[jj] = partition_ii[jj];
        }
    }
    // Set the special sketched dimension properties
    n[end_ind] = s->r;
    nps[end_ind] = 1;
    partitions[end_ind] = (int*) malloc( (nps[end_ind] + 1) * sizeof(int));
    int* spartition_end = partitions[end_ind];
    spartition_end[0] = 0; spartition_end[1] = n[end_ind];

    // Get the parameters giving the subtensor locations
    double** subtensors = (double**) malloc(n_schedule * sizeof(double*));
    int* schedule_rank = schedule[rank];
    int with_buf = 0;
    long scratch_offset = 0;
    for (int ii = 0; ii < n_schedule; ++ii){
        if (schedule_rank[ii] != -1){
            matrix_tt* X_submat = own_submatrix(s, schedule_rank[ii], with_buf);
            subtensors[ii] = s->scratch + scratch_offset;
            matrix_tt* scratch_submat = matrix_tt_wrap(X_submat->m, X_submat->n, subtensors[ii]);
            matrix_tt_copy_data(scratch_submat, X_submat);
            scratch_offset = scratch_offset + X_submat->n * X_submat->m;
            free(X_submat);
            free(scratch_submat);
        }
    }
    void* parameters = p_static_init(subtensors);

    // Assigning fields
    sten->d = d;
    sten->n = n;

    sten->comm = ten->comm;
    sten->rank = ten->rank;
    sten->comm_size = ten->comm_size;

    sten->schedule = schedule;
    sten->n_schedule = n_schedule;
    sten->inverse_schedule = inverse_schedule;

    sten->partitions = partitions;
    sten->nps = nps;

    sten->current_part = -1;

    sten->f_ten = NULL;
    sten->parameters = parameters;

    sten->X_size = s->X_size;
    sten->X = s->scratch;

    sten->ind1 = (int*) malloc(d * sizeof(int));
    sten->ind2 = (int*) malloc(d * sizeof(int));
    sten->tensor_part = (int*) malloc(d * sizeof(int));
    sten->group_ranks = (int*) malloc(ten->comm_size * sizeof(int));
    sten->t_kk = (int*) malloc(d * sizeof(int));

    // Freeing sketch things
    int n_Omega = s->iscol ? ten->d - s->flattening : s->flattening;
    for (int ii = 0; ii < n_Omega; ++ii){
        matrix_tt_free(s->Omegas[ii]); s->Omegas[ii] = NULL;
    }
    for (int ii = 0; ii < 4; ++ii){
        matrix_tt_free(s->KRs[ii]); s->KRs[ii] = NULL;
    }
    free(s->KRs); s->KRs = NULL;
    s->ten = NULL;
    free(s->owner_partition); s->owner_partition = NULL;
    free(s->X);               s->X = NULL;
    free(s->Omegas);          s->Omegas = NULL;
    free(s->recv_buf);        s->recv_buf = NULL;
    flattening_info_free(s->fi); s->fi = NULL;
    free(s);                  *s_ptr = NULL;

    flattening_info_free(fi);

    return sten;
}


text:
#include <gsl/gsl_errno.h>
#include <gsl/matrix/gsl_matrix.h>

#define BASE_DOUBLE
#include <gsl/templates_on.h>
#include <gsl/matrix/matrix_source.c>
#include <gsl/templates_off.h>
#undef BASE_DOUBLE

ext: c

max_stars_repo_name: MontyThibault/centre-of-mass-awareness

max_issues_repo_name: MontyThibault/centre-of-mass-awareness

max_forks_repo_name: MontyThibault/centre-of-mass-awareness


text:
/**
  * Copyright 2016 José Manuel Abuín Mosquera <josemanuel.abuin@usc.es>
  *
  * This file is part of Matrix Market Suite.
  *
  * Matrix Market Suite 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.
  *
  * Matrix Market Suite 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 Matrix Market Suite. If not, see <http://www.gnu.org/licenses/>.
  */

#include <lapacke.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>
#include <stdlib.h>

#include "LUDecomposition.h"

void usageLUDecomposition(){

	fprintf(stderr, "\n");
	fprintf(stderr, "Usage: MM-Suite LUDecomposition [options] <input-matrix>\n");
	fprintf(stderr, "\nInput/output options:\n\n");
	fprintf(stderr, "       -o STR        Output file name. Default: stdout\n");
	fprintf(stderr, "       -r            Input format is row per line. Default: False\n");
	fprintf(stderr, "\n");

}

int LUDecomposition(int argc, char *argv[]) {

	int 			ret_code = 1;
	int 			option;

	unsigned long 		*II;

	unsigned long 		*J;
	double 			*values;

	unsigned long 		M;
	unsigned long 		N;
	unsigned long long 	nz;

	int			_M;
	int			_N;

	char			*outputFileName = NULL;

	char			*inputMatrixFile = NULL;

	int			inputFormatRow = 0;

	while ((option = getopt(argc, argv,"ro:")) >= 0) {
		switch (option) {
			case 'o' :
				//free(outputFileName);

				outputFileName = (char *) malloc(sizeof(char)*strlen(optarg)+1);
				strcpy(outputFileName,optarg);

				break;

			case 'r':
				inputFormatRow = 1;
				break;

			default: break;
		}

	}

	if ((optind + 1 > argc) || (optind + 2 <= argc)) {
		usageLUDecomposition();
		return 0;
	}

	if(outputFileName == NULL) {
		outputFileName = (char *) malloc(sizeof(char)*7);
		sprintf(outputFileName,"stdout");
	}

	inputMatrixFile = (char *)malloc(sizeof(char)*strlen(argv[optind])+1);

	strcpy(inputMatrixFile,argv[optind]);

	//Read matrix
	if(inputFormatRow){

		if(!readDenseCoordinateMatrixRowLine(inputMatrixFile,&II,&J,&values,&M,&N,&nz)){
			fprintf(stderr, "[%s] Can not read Matrix\n",__func__);
			return 0;
		}

	}
	else {
		if(!readDenseCoordinateMatrix(inputMatrixFile,&II,&J,&values,&M,&N,&nz)){
			fprintf(stderr, "[%s] Can not read Matrix\n",__func__);
			return 0;
		}
	}

	_M = (int) M;
	_N = (int) N;

	int ipiv[_M];

	/*
	lapack_int LAPACKE_dgetrf( int matrix_layout, lapack_int m, lapack_int n,
                           double* a, lapack_int lda, lapack_int* ipiv );
	*/

        ret_code = LAPACKE_dgetrf(LAPACK_ROW_MAJOR,_M,_N, values,_N, ipiv);

	if(inputFormatRow){
		writeLUCoordinateMatrixRowLine(outputFileName, values,M,N,nz, ipiv);
	}
	else{
		writeLUCoordinateMatrix(outputFileName, values,M,N,nz, ipiv);
	}

	//writeDenseCoordinateMatrix("stdout",values,M,N,nz);
	if(!ret_code)
		return 1;
	else
		return 0;

}


text:
#pragma once
#include <Hypo/Config.h>
#include <Hypo/System/Exports.h>
#include <gsl/span>
#include <Hypo/System/Streams/MemoryStream.h>
#include <Hypo/System/Buffer/Buffer.h>

namespace Hypo
{
	class HYPO_SYSTEM_API BinaryReader
	{
	public:
		BinaryReader(std::istream& stream);

		BinaryReader& operator >>(uInt8& value);
		BinaryReader& operator >>(uInt16& value);
		BinaryReader& operator >>(uInt32& value);
		BinaryReader& operator >>(uInt64& value);
		BinaryReader& operator >>(Int8& value );
		BinaryReader& operator >>(Int16& value);
		BinaryReader& operator >>(Int32& value);
		BinaryReader& operator >>(Int64& value);

		BinaryReader& operator >>(std::string& value);

		template<typename T>
		BinaryReader& operator >>(std::vector<T>& value)
		{
			uInt32 size = 0;
			*this >> size;

			for (int i = 0; i < size; i++)
			{
				T elem;
				*this >> elem;
				value.push_back(elem);
			}
			return *this;
		}

		std::vector<unsigned char> ReadRaw(int size);
		void ReadRaw(void* buffer, int size);

		void ReadBOM();

		void Read7BitEncoded(uInt32& value);
		void Read7BitEncoded(uInt64& value);

		bool Good()const { return m_Stream.good(); }
		bool Bad()const { return m_Stream.bad(); }
		bool Fail() const { return m_Stream.fail(); }
		bool Eof() const { return m_Stream.eof(); }

		std::istream& GetStream() const { return m_Stream; }

	private:
		std::istream& m_Stream;
	};

	template <typename T>
	class BasicMemoryBinaryReader : public BinaryReader
		/// A convenient wrapper for using Buffer and MemoryStream with BinaryReader.
	{
	public:
		BasicMemoryBinaryReader(const Buffer<T>& data) :
			BinaryReader(_istr),
			_data(data),
			_istr(data.begin(), data.Capacity())
		{
		}

		~BasicMemoryBinaryReader()
		{
		}

		const Buffer<T>& data() const
		{
			return _data;
		}

		const MemoryInputStream& stream() const
		{
			return _istr;
		}

		MemoryInputStream& stream()
		{
			return _istr;
		}

	private:
		const Buffer<T>& _data;
		MemoryInputStream _istr;
	};

	typedef BasicMemoryBinaryReader<char> MemoryBinaryReader;

}

ext: h

max_stars_repo_name: TheodorLindberg/Hypo

max_issues_repo_name: TheodorLindberg/Hypo

max_forks_repo_name: TheodorLindberg/Hypo


text:
/**
 * \copyright
 * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
 *            Distributed under a Modified BSD License.
 *              See accompanying file LICENSE.txt or
 *              http://www.opengeosys.org/project/license
 *
 */

#pragma once

/// The LinearSolverLibrarySetup takes care of proper initialization and
/// shutting down of an external linear solver library. The concrete
/// implementation is chosen by the build system.
/// An object of this class must be created at the begining of the scope where
/// it is used. When the scope closes (or the object is destroyed explicitly)
/// library shutting down functions are automatically called.
/// The default implementation is empty providing polymorphic behaviour when
/// using this class.

#include "NumLib/DOF/GlobalMatrixProviders.h"

#if defined(USE_PETSC)
#include <petsc.h>
#include <mpi.h>
namespace ApplicationsLib
{
struct LinearSolverLibrarySetup final
{
    LinearSolverLibrarySetup(int argc, char* argv[])
    {
        MPI_Init(&argc, &argv);
        char help[] = "ogs6 with PETSc \n";
        PetscInitialize(&argc, &argv, nullptr, help);
    }

    ~LinearSolverLibrarySetup()
    {
        NumLib::cleanupGlobalMatrixProviders();
        PetscFinalize();
        MPI_Finalize();
    }
};
}    // ApplicationsLib
#elif defined(USE_LIS)
#include <lis.h>
namespace ApplicationsLib
{
struct LinearSolverLibrarySetup final
{
    LinearSolverLibrarySetup(int argc, char* argv[])
    {
        lis_initialize(&argc, &argv);
    }

    ~LinearSolverLibrarySetup()
    {
        NumLib::cleanupGlobalMatrixProviders();
        lis_finalize();
    }
};
}    // ApplicationsLib
#else
namespace ApplicationsLib
{
struct LinearSolverLibrarySetup final
{
    LinearSolverLibrarySetup(int /*argc*/, char* /*argv*/[]) {}
    ~LinearSolverLibrarySetup()
    {
        NumLib::cleanupGlobalMatrixProviders();
    }
};
}    // ApplicationsLib
#endif

max_stars_repo_name: hosseinsotudeh/ogs

max_issues_repo_name: hosseinsotudeh/ogs

max_forks_repo_name: hosseinsotudeh/ogs


text:
/* GLIB - Library of useful routines for C programming
 * Copyright (C) 1995-1997  Peter Mattis, Spencer Kimball and Josh MacDonald
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

/*
 * Modified by the GLib Team and others 1997-2000.  See the AUTHORS
 * file for a list of people on the GLib Team.  See the ChangeLog
 * files for a list of changes.  These files are distributed with
 * GLib at ftp://ftp.gtk.org/pub/gtk/.
 */

#ifndef __G_LIB_H__
#define __G_LIB_H__

#include <galloca.h>
#include <garray.h>
#include <gasyncqueue.h>
#include <gbacktrace.h>
#include <gcache.h>
#include <gcompat.h>
#include <gcompletion.h>
#include <gconvert.h>
#include <gdataset.h>
#include <gdate.h>
#include <gerror.h>
#include <gfileutils.h>
#include <ghash.h>
#include <ghook.h>
#include <giochannel.h>
#include <glist.h>
#include <gmacros.h>
#include <gmain.h>
#include <gmarkup.h>
#include <gmem.h>
#include <gmessages.h>
#include <gnode.h>
#include <gprimes.h>
#include <gqsort.h>
#include <gquark.h>
#include <gqueue.h>
#include <grand.h>
#include <grel.h>
#include <gscanner.h>
#include <gshell.h>
#include <gslist.h>
#include <gspawn.h>
#include <gstrfuncs.h>
#include <gstring.h>
#include <gthread.h>
#include <gthreadpool.h>
#include <gtimer.h>
#include <gtree.h>
#include <gtypes.h>
#include <gunicode.h>
#include <gutils.h>
#ifdef G_OS_WIN32
#include <gwin32.h>
#endif

#endif /* __G_LIB_H__ */


text:
#include <stdio.h>
#include <gsl/gsl_multiset.h>

void
print_all_multisets ( size_t n, size_t k );

int
main (void)
{
  gsl_multiset * c;
  size_t i;

  printf("all multisets of {0,1,2,3} by size (lex. order)\n") ;
  for(i = 0; i <= 4; i++)
    {
      c = gsl_multiset_calloc (4, i);
      do
        {
          printf("{");
          gsl_multiset_fprintf (stdout, c, " %u");
          printf(" }\n");
        }
      while (gsl_multiset_next(c) == GSL_SUCCESS);
      gsl_multiset_free(c);
    }
  printf("all multisets of {1,2,3,4} by size (reverse lex. order)\n") ;
  for(i = 0; i <= 4; i++)
    {
      c = gsl_multiset_alloc (4, i) ;
      gsl_multiset_init_last(c) ;
      do
        {
          printf("{");
          gsl_multiset_fprintf (stdout, c, " %u");
          printf(" }\n");
        }
      while (gsl_multiset_prev(c) == GSL_SUCCESS);
      gsl_multiset_free(c);
    }
  printf("\n");

  print_all_multisets(5, 3);
  print_all_multisets(5, 0);
  print_all_multisets(5, 5);
  print_all_multisets(1, 1);
  print_all_multisets(3, 1);

  return 0;
}

void
print_all_multisets (size_t n, size_t k)
{
  gsl_multiset * c = gsl_multiset_calloc (n, k);

  printf("multisets %u choose %u (with replacement)\n", n, k);
  do
    {
      gsl_multiset_fprintf (stdout, c, " %u");
      printf("\n");
    }
  while (gsl_multiset_next(c) == GSL_SUCCESS);
  while (gsl_multiset_next(c) == GSL_SUCCESS);
  do
    {
      gsl_multiset_fprintf (stdout, c, " %u");
      printf("\n");
    }
  while (gsl_multiset_prev(c) == GSL_SUCCESS);
  printf("\n");
  gsl_multiset_free(c);
}

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:
#pragma once

#include <gsl/span>

#include "halley/text/halleystring.h"
#include "halley/resources/resource.h"
#include "halley/maths/range.h"
#include "halley/maths/vector4.h"
#include "halley/data_structures/maybe.h"
#include <map>
#include <vector>

#include "halley/utils/type_traits.h"

#if defined(DEV_BUILD)
#define STORE_CONFIG_NODE_PARENTING
#endif

namespace Halley {
	class Serializer;
	class Deserializer;
	class ConfigNode;

	template<typename T>
	struct HasToConfigNode
	{
	private:
		typedef std::true_type yes;
		typedef std::false_type no;
		template<typename U> static auto test(int) -> decltype(std::declval<U>().toConfigNode(), yes());
		template<typename> static no test(...);

	public:
		static constexpr bool value = std::is_same<decltype(test<T>(0)),yes>::value;
	};

	template<typename T>
	struct HasConfigNodeConstructor
	{
	private:
		typedef std::true_type yes;
		typedef std::false_type no;
		template<typename U> static auto test(int) -> decltype(U(std::declval<ConfigNode>()), yes());
		template<typename> static no test(...);

	public:
		static constexpr bool value = std::is_same<decltype(test<T>(0)),yes>::value;
	};

	enum class ConfigNodeType
	{
		Undefined,
		String,
		Sequence,
		Map,
		Int,
		Float,
		Int2,
		Float2,
		Bytes,
		DeltaSequence, // For delta coding
		DeltaMap, // For delta coding
		Noop, // For delta coding
		Idx, // For delta coding
		Del // For delta coding
	};

	template <>
	struct EnumNames<ConfigNodeType> {
		constexpr std::array<const char*, 13> operator()() const {
			return{{
				"undefined",
				"string",
				"sequence",
				"map",
				"int",
				"float",
				"int2",
				"float2",
				"bytes",
				"deltaSequence",
				"deltaMap",
				"noop",
				"idx"
				"del"
			}};
		}
	};

	class ConfigFile;

	class ConfigNode
	{
		friend class ConfigFile;

	public:
		using MapType = std::map<String, ConfigNode, std::less<>>;
		using SequenceType = std::vector<ConfigNode>;

		struct NoopType {};
		struct DelType {};
		struct IdxType {
			int start;
			int len;
			IdxType() = default;
			IdxType(int start, int len) : start(start), len(len) {}
		};

		ConfigNode();
		explicit ConfigNode(const ConfigNode& other);
		ConfigNode(ConfigNode&& other) noexcept;
		ConfigNode(MapType entryMap);
		ConfigNode(SequenceType entryList);
		explicit ConfigNode(String value);
		explicit ConfigNode(const std::string_view& value);
		explicit ConfigNode(bool value);
		explicit ConfigNode(int value);
		explicit ConfigNode(float value);
		explicit ConfigNode(Vector2i value);
		explicit ConfigNode(Vector2f value);
		explicit ConfigNode(Bytes value);
		explicit ConfigNode(NoopType value);
		explicit ConfigNode(DelType value);
		explicit ConfigNode(IdxType value);

		template <typename T>
		explicit ConfigNode(const std::vector<T>& sequence)
		{
			*this = sequence;
		}

		template <typename K, typename V>
		explicit ConfigNode(const std::map<K, V>& values)
		{
			*this = values;
		}

		~ConfigNode();

		ConfigNode& operator=(const ConfigNode& other);
		ConfigNode& operator=(ConfigNode&& other) noexcept;
		ConfigNode& operator=(bool value);
		ConfigNode& operator=(int value);
		ConfigNode& operator=(float value);
		ConfigNode& operator=(Vector2i value);
		ConfigNode& operator=(Vector2f value);

		ConfigNode& operator=(MapType entryMap);
		ConfigNode& operator=(SequenceType entryList);
		ConfigNode& operator=(String value);
		ConfigNode& operator=(Bytes value);
		ConfigNode& operator=(gsl::span<const gsl::byte> bytes);

		ConfigNode& operator=(const char* value);
		ConfigNode& operator=(const std::string_view& value);

		ConfigNode& operator=(NoopType value);
		ConfigNode& operator=(DelType value);
		ConfigNode& operator=(IdxType value);

		template <typename T>
		ConfigNode& operator=(const std::vector<T>& sequence)
		{
			SequenceType seq;
			seq.reserve(sequence.size());
			for (const auto& e: sequence) {
				if constexpr (HasToConfigNode<T>::value) {
					seq.push_back(e.toConfigNode());
				} else {
					seq.push_back(ConfigNode(e));
				}
			}
			return *this = seq;
		}

		template <typename K, typename V>
		ConfigNode& operator=(const std::map<K, V>& values)
		{
			MapType map;
			for (const auto& [k, v]: values) {
				String key = toString(k);

				if constexpr (HasToConfigNode<V>::value) {
					map[std::move(key)] = v.toConfigNode();
				} else {
					map[std::move(key)] = ConfigNode(v);
				}
			}
			return *this = map;
		}

		bool operator==(const ConfigNode& other) const;
		bool operator!=(const ConfigNode& other) const;

		ConfigNodeType getType() const;

		void serialize(Serializer& s) const;
		void deserialize(Deserializer& s);

		int asInt() const;
		float asFloat() const;
		bool asBool() const;
		Vector2i asVector2i() const;
		Vector2f asVector2f() const;
		Vector3i asVector3i() const;
		Vector3f asVector3f() const;
		Vector4i asVector4i() const;
		Vector4f asVector4f() const;
		Range<float> asFloatRange() const;
		String asString() const;
		const Bytes& asBytes() const;

		int asInt(int defaultValue) const;
		float asFloat(float defaultValue) const;
		bool asBool(bool defaultValue) const;
		String asString(const std::string_view& defaultValue) const;
		Vector2i asVector2i(Vector2i defaultValue) const;
		Vector2f asVector2f(Vector2f defaultValue) const;
		Vector3i asVector3i(Vector3i defaultValue) const;
		Vector3f asVector3f(Vector3f defaultValue) const;
		Vector4i asVector4i(Vector4i defaultValue) const;
		Vector4f asVector4f(Vector4f defaultValue) const;

		template <typename T>
		std::vector<T> asVector() const
		{
			if (type == ConfigNodeType::Sequence) {
				std::vector<T> result;
				result.reserve(asSequence().size());
				for (const auto& e : asSequence()) {
					if constexpr (HasConfigNodeConstructor<T>::value) {
						result.emplace_back(T(e));
					} else {
						result.emplace_back(e.convertTo(Tag<T>()));
					}
				}
				return result;
			} else if (type == ConfigNodeType::Undefined) {
				return {};
			} else {
				throw Exception("Can't convert " + getNodeDebugId() + " from " + toString(getType()) + " to std::vector<T>.", HalleyExceptions::Resources);
			}
		}

		template <typename T>
		std::vector<T> asVector(const std::vector<T>& defaultValue) const
		{
			if (type == ConfigNodeType::Sequence) {
				return asVector<T>();
			} else {
				return defaultValue;
			}
		}

		template <typename K, typename V>
		std::map<K, V> asMap() const
		{
			if (type == ConfigNodeType::Map) {
				std::map<K, V> result;
				for (const auto& [k, v] : asMap()) {
					const K key = fromString<K>(k);

					if constexpr (HasConfigNodeConstructor<V>::value) {
						result[std::move(key)] = V(v);
					} else {
						result[std::move(key)] = v.convertTo(Tag<V>());
					}
				}
				return result;
			} else if (type == ConfigNodeType::Undefined) {
				return {};
			} else {
				throw Exception("Can't convert " + getNodeDebugId() + " from " + toString(getType()) + " to std::map<K, V>.", HalleyExceptions::Resources);
			}
		}

		const SequenceType& asSequence() const;
		const MapType& asMap() const;
		SequenceType& asSequence();
		MapType& asMap();

		void ensureType(ConfigNodeType type);

		bool hasKey(const String& key) const;
		void removeKey(const String& key);

		ConfigNode& operator[](const std::string_view& key);
		ConfigNode& operator[](size_t idx);

		const ConfigNode& operator[](const std::string_view& key) const;
		const ConfigNode& operator[](size_t idx) const;

		SequenceType::iterator begin();
		SequenceType::iterator end();

		SequenceType::const_iterator begin() const;
		SequenceType::const_iterator end() const;

		void reset();
		void setOriginalPosition(int line, int column);
		void setParent(const ConfigNode* parent, int idx);
		void propagateParentingInformation(const ConfigFile* parentFile);

		inline void assertValid() const
		{
			Expects(intData != 0xCDCDCDCD);
			Expects(intData != 0xDDDDDDDD);
		}

		struct BreadCrumb {
			const BreadCrumb* prev = nullptr;
			String key;
			OptionalLite<int> idx;
			int depth = 0;

			BreadCrumb() = default;
			BreadCrumb(const BreadCrumb& prev, String key) : prev(&prev), key(std::move(key)), depth(prev.depth + 1) {}
			BreadCrumb(const BreadCrumb& prev, int index) : prev(&prev), idx(index), depth(prev.depth + 1) {}

			bool hasKeyAt(const String& key, int depth) const;
			bool hasIndexAt(int idx, int depth) const;
		};

		class IDeltaCodeHints {
		public:
			virtual ~IDeltaCodeHints() = default;

			virtual std::optional<size_t> getSequenceMatch(const SequenceType& seq, const ConfigNode& newValue, size_t curIdx, const BreadCrumb& breadCrumb) const = 0;
			virtual bool doesSequenceOrderMatter(const BreadCrumb& breadCrumb) const { return true; }
			virtual bool canDeleteKey(const String& key, const BreadCrumb& breadCrumb) const { return true; }
			virtual bool canDeleteAnyKey() const { return true; }
			virtual bool shouldBypass(const BreadCrumb& breadCrumb) const { return false; }
			virtual bool areNullAndEmptyEquivalent(const BreadCrumb& breadCrumb) const { return false; }
		};

		static ConfigNode createDelta(const ConfigNode& from, const ConfigNode& to, const IDeltaCodeHints* hints = nullptr);
		static ConfigNode applyDelta(const ConfigNode& from, const ConfigNode& delta);
		void applyDelta(const ConfigNode& delta);
		void decayDeltaArtifacts();

	private:
		template <typename T>
		class Tag {};

		union {
			String* strData;
			MapType* mapData;
			SequenceType* sequenceData;
			Bytes* bytesData;
			void* rawPtrData;
			int intData;
			float floatData;
			Vector2i vec2iData;
			Vector2f vec2fData;
		};
		ConfigNodeType type = ConfigNodeType::Undefined;
		int auxData = 0; // Used by delta coding

#if defined(STORE_CONFIG_NODE_PARENTING)
		struct ParentingInfo {
			int line = 0;
			int column = 0;
			int idx = 0;
			const ConfigNode* node = nullptr;
			const ConfigFile* file = nullptr;
		};
		std::unique_ptr<ParentingInfo> parent;
#endif

		static ConfigNode undefinedConfigNode;
		static String undefinedConfigNodeName;

		template <typename T> void deserializeContents(Deserializer& s)
		{
			T v;
			s >> v;
			*this = std::move(v);
		}

		String getNodeDebugId() const;
		String backTrackFullNodeName() const;

		int convertTo(Tag<int> tag) const;
		float convertTo(Tag<float> tag) const;
		bool convertTo(Tag<bool> tag) const;
		Vector2i convertTo(Tag<Vector2i> tag) const;
		Vector2f convertTo(Tag<Vector2f> tag) const;
		Vector3i convertTo(Tag<Vector3i> tag) const;
		Vector3f convertTo(Tag<Vector3f> tag) const;
		Vector4i convertTo(Tag<Vector4i> tag) const;
		Vector4f convertTo(Tag<Vector4f> tag) const;
		Range<float> convertTo(Tag<Range<float>> tag) const;
		String convertTo(Tag<String> tag) const;
		const Bytes& convertTo(Tag<Bytes&> tag) const;

		bool isNullOrEmpty() const;

		static ConfigNode doCreateDelta(const ConfigNode& from, const ConfigNode& to, const BreadCrumb& breadCrumb, const IDeltaCodeHints* hints);
		static ConfigNode createMapDelta(const ConfigNode& from, const ConfigNode& to, const BreadCrumb& breadCrumb, const IDeltaCodeHints* hints);
		static ConfigNode createSequenceDelta(const ConfigNode& from, const ConfigNode& to, const BreadCrumb& breadCrumb, const IDeltaCodeHints* hints);
		void applyMapDelta(const ConfigNode& delta);
		void applySequenceDelta(const ConfigNode& delta);

		bool isEquivalent(const ConfigNode& other) const;
		bool isEquivalentStrictOrder(const ConfigNode& other) const;
	};
}

max_stars_repo_name: moonantonio/halley

max_issues_repo_name: moonantonio/halley

max_forks_repo_name: moonantonio/halley


text:
/*
 * Copyright (c) 2015, Aleksas Mazeliauskas and Derek Teaney
 * All rights reserved.
 *
 * rnavier is distributed under MIT license;
 * see the LICENSE file that should be present in the root
 * of the source distribution, or alternately available at:
 * https://github.com/rnavier/rnavier/
 */
#ifndef RN_TEOSLats95p_h
#define RN_TEOSLats95p_h

#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include "TEOS.h"
#include "numeric.h"

// Lattice EoS s95p_v0 by Huovinen and Petreczky. Data files are stored in /s95p, correspondingly
// in *par*.dat (T, \mu_1, \mu_2) and *dens*.dat (e, P, s, n_B, plasma ratio), with natural units,
// e.g., T-GeV, \mu-GeV, e-GeV/fm^3. Each file is
// particularly available in its range.
// file1: T~[0.0159,0.1801], e~[0.001,1.001]
// file2: T~[0.1801,0.2810], e~[1.001,11.00]
// file3: T~[0.2810,0.4243], e~[11.00,61.00]
// file4: T~[0.2423,0.6335], e~[61.00,311.0]
// There are also files avaiable for T<180MeV, with 28 non-zero chemical potential for baryon,
// strangeness, etc. For the time being, and for simplicity we ignore these files. So basically
// we have \mu = 0 and n_B = 0.

class TEOSs95p : public virtual TEOS {

   private:

      //minimum table values
      double fEmin;
      double fSmin;
      double fTmin;
      double fPmin;

      gsl_interp_accel *acc_eVSp;
      gsl_spline *eVSp;

      gsl_interp_accel *acc_tVSx;
      gsl_spline *tVSx;

      gsl_interp_accel *acc_eVSt;
      gsl_spline *eVSt;

      gsl_interp_accel *acc_eVSs;
      gsl_spline *eVSs;

      gsl_interp_accel *acc_sVSe;
      gsl_spline *sVSe;

      void eos(const double &e, double &p, double &cs)  ;
      void st(const double &e, double &s, double &t)  ;

      double fEtaOverS ;    //!< The Shear viscosity/entropy
      double fEtaOverSHad ; //!< The Shear viscosity/entropy
      double fSigmaOverS ;  //!< The Shear viscosity/entropy
      double fKappaTOverS ; //!< The Conductivity/entropy
      double fTPi_EtaST ;  //!< tau_pi / (eta/sT)
      double fL1_EtaTPi ;  //!< lambda_1 /(eta*tau_pi)
      double fL2_EtaTPi ;  //!< lambda_2 /(eta*tau_pi)

   public:

      TEOSs95p ();
      ~TEOSs95p ();

      virtual void read(std::istream &in)  ;

      virtual void write(std::ostream &out) ;

      virtual void eos(const double &e, const double &n,
            double &p, double &cs) ;

      virtual double eofs(const double &s, const double &n) ;

      virtual void stmu(const double &e, const double &n,
            double &s, double &t, double &mu) ;

      virtual void viscosity(const double &e, const double &n,
         double &sigma_overs, double &kappaT_overs, double &eta_overs) ;

      virtual void getBRSSSParams(const double &e, const double &n, double &tpi_etast, double &l1_ntpi, double &l2_ntpi)
      { tpi_etast =  fTPi_EtaST ; l1_ntpi = fL1_EtaTPi ; l2_ntpi = fL2_EtaTPi ; }

      void tmux(const double &e, const double &n, double &t, double &mub, double &mus, double &x) ;

} ;

std::unique_ptr<TEOS> make_eos_eoss95p(TRNavier3DBj *rn, const std::string &icname) ;
#endif

max_stars_repo_name: rnavier/rnavier

max_issues_repo_name: rnavier/rnavier

max_forks_repo_name: rnavier/rnavier


text:
/*
 *
 * Copyright (c) Kresimir Fresl 2002
 *
 * Distributed under the Boost Software License, Version 1.0.
 * (See accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * Author acknowledges the support of the Faculty of Civil Engineering,
 * University of Zagreb, Croatia.
 *
 */

#ifndef BOOST_NUMERIC_BINDINGS_BLAS_DETAIL_CBLAS_H
#define BOOST_NUMERIC_BINDINGS_BLAS_DETAIL_CBLAS_H

#ifdef HAS_OpenBLAS
#define OPENBLAS_CONST_FLOAT_CAST(x) reinterpret_cast<const float*>(x)
#define OPENBLAS_CONST_DOUBLE_CAST(x) reinterpret_cast<const double*>(x)
#define OPENBLAS_FLOAT_CAST(x) reinterpret_cast<float*>(x)
#define OPENBLAS_DOUBLE_CAST(x) reinterpret_cast<double*>(x)
#define OPENBLAS_OPENBLAS_COMPLEX_FLOAT_CAST(x) reinterpret_cast<openblas_complex_float*>(x)
#define OPENBLAS_OPENBLAS_COMPLEX_DOUBLE_CAST(x) reinterpret_cast<openblas_complex_double*>(x)

#else
#define OPENBLAS_CONST_FLOAT_CAST(x) x
#define OPENBLAS_CONST_DOUBLE_CAST(x) x
#define OPENBLAS_FLOAT_CAST(x) x
#define OPENBLAS_DOUBLE_CAST(x) x
#define OPENBLAS_OPENBLAS_COMPLEX_FLOAT_CAST(x) x
#define OPENBLAS_OPENBLAS_COMPLEX_DOUBLE_CAST(x) x
#endif

//
// MKL-specific CBLAS include
//
#if defined BOOST_NUMERIC_BINDINGS_BLAS_MKL

extern "C" {
#include <mkl_cblas.h>
//#include <mkl_service.h>
//
// mkl_types.h defines P4 macro which breaks MPL, undefine it here.
//
#undef P4
}

//
// Default CBLAS include
//
#else

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

#endif
#endif

max_stars_repo_name: ljktest/siconos

max_issues_repo_name: ljktest/siconos

max_forks_repo_name: ljktest/siconos


text:
#pragma once

#ifndef _NOGSL
#include <gsl/gsl_vector.h>
#else
#include "FakeGSL.h"
#endif

#include <fstream>
#include <iostream>
#include <vector>
#include <map>
#include <math.h>
#include "FloatSupport.h"
#include "BooleanDAG.h"
#include "Interface.h"
#include "ConflictGenerator.h"
#include "ActualEvaluators.h"
#include "CommandLineArgs.h"
#include "SmoothEvaluators.h"
#include "NumDebugger.h"

using namespace std;

class NumericalSolver {
protected:
	BooleanDAG* dag;
	Interface* interf;

    int ncontrols;
    map<string, int>& ctrls;
    SmoothEvaluators* smoothEval;
    ActualEvaluators* actualEval;
    OptimizationWrapper* opt;

    SimpleEvaluator* seval;

    set<int> assertConstraints;
    int minimizeNode;

    const vector<vector<int>>& dependentInputs;
    const vector<vector<int>>& dependentCtrls;

    gsl_vector* result;
    LocalState* localState;

	NumDebugger* debugger;

public:
    NumericalSolver(BooleanDAG* _dag, map<string, int>& _ctrls, Interface* _interface, SmoothEvaluators* _smoothEval, ActualEvaluators* _actualEval, OptimizationWrapper* _opt, const vector<vector<int>>& _dependentInputs, const vector<vector<int>>& _dependentCtrls, NumDebugger* _debugger);
    ~NumericalSolver(void);

	// Called by the NumericalSolver
    bool checkSAT(int level, gsl_vector* initState = NULL);

    bool checkFullSAT(gsl_vector* state);

    gsl_vector* getResult();
    LocalState* getLocalState();

};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


text:
 // MIT License

 // Copyright (c) [2017] [Vinay Yuvashankar]

 // Permission is hereby granted, free of charge, to any person obtaining a copy
 // of this software and associated documentation files (the "Software"), to deal
 // in the Software without restriction, including without limitation the rights
 // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 // copies of the Software, and to permit persons to whom the Software is
 // furnished to do so, subject to the following conditions:

 // The above copyright notice and this permission notice shall be included in all
 // copies or substantial portions of the Software.

 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 // SOFTWARE.
/**
	\file "Wavelet.cc"
	\brief This file contains all of the functions that support the ERSP and CWT functions
*/
#include "wavelet.h"
#include <omp.h>
#include <assert.h>
#include <gsl/gsl_statistics.h>

#define TEST 0.00001
#define SETTLING_PERCENTAGE 0.02
#define NORMALIZATION_FACTOR 0.375402913609157562

int Wavelet(double* raw_data, double* scales,
	double sampling_frequency, int n, int J,
	double* result)
{
	//Variable Declarations
	int i, j;
	fftw_complex *data_in, *fft_data;
	fftw_plan plan_forward;

	//Calculate Padding Required
    // const int PADDED_SIZE = CalculatePaddingSize(n, 1);
    const int PADDED_SIZE = n;

    const double dw = (2 * M_PI * sampling_frequency)/(PADDED_SIZE); //NOT IN RAD/SAMPLE in RAD/SEC

    data_in  = (fftw_complex *) fftw_malloc( sizeof( fftw_complex ) * PADDED_SIZE );
	fft_data = (fftw_complex *) fftw_malloc( sizeof( fftw_complex ) * PADDED_SIZE );

	//populate the FFTW data vector
	for (i = 0; i < n; ++i)
    {
    	data_in[i][0] = raw_data[i];
    	// data_in[i + n][0] = raw_data[i];

    	// data_in[i + n][1] = 0.0;
    	data_in[i    ][1] = 0.0;
    }

    // //Force the rest of the data vector to zero just in case
    // for (int i = n; i < PADDED_SIZE; ++i)
    // {
    // 	data_in[i][0] = 0.0;
    // 	data_in[i][1] = 0.0;
    // }

    double *temp = (double*) malloc(n * sizeof(double));
    FILE* debug_file = fopen("debug.log", "w");

	//Calculate the FFT of the data and store it in fft_data
	plan_forward = fftw_plan_dft_1d(PADDED_SIZE, data_in, fft_data,
									FFTW_FORWARD, FFTW_ESTIMATE);
	fftw_execute(plan_forward);

	// #pragma omp parallel num_threads(1) private(i, j) shared (result, sampling_frequency, J, n, scales,  fft_data) default(none)
	// {
		double value;

		fftw_plan plan_backward;
		fftw_complex *filter_convolution, *fftw_result;

		filter_convolution = (fftw_complex *) fftw_malloc( sizeof( fftw_complex )* PADDED_SIZE );
		fftw_result  = 		 (fftw_complex *) fftw_malloc( sizeof( fftw_complex )* PADDED_SIZE );

		// #pragma omp critical (make_plan)
		// {
			//Preapre for the plan backwards
			plan_backward = fftw_plan_dft_1d(PADDED_SIZE, filter_convolution, fftw_result,
				FFTW_BACKWARD, FFTW_ESTIMATE);
		// }

	    // #pragma omp for
		for (i = 0; i < J; ++i)
		{
			//Force the arrays to zero
			memset(filter_convolution, 0.0, sizeof( fftw_complex ) * PADDED_SIZE);
			memset(fftw_result,        0.0, sizeof( fftw_complex ) * PADDED_SIZE);

			//Compute the Fourier Morlet at 0 and N/2
			double norm = sqrt(scales[i]);

			value = CompleteFourierMorlet(0.0, scales[i], norm);

			filter_convolution[0][0] = ( fft_data[0][0] / PADDED_SIZE ) * value;
			filter_convolution[0][1] = ( fft_data[0][1] / PADDED_SIZE ) * value;

			filter_convolution[PADDED_SIZE/2][0] = 0.0;
			filter_convolution[PADDED_SIZE/2][1] = 0.0;

			//Compute the Fourier Morlet Convolution in between
			for (j = 1; j < PADDED_SIZE/2 - 1; ++j)
			{
				value = CompleteFourierMorlet( j * dw , scales[i], norm);

				filter_convolution[j][0] = ( fft_data[j][0] / PADDED_SIZE ) * value;
				filter_convolution[j][1] = ( fft_data[j][1] / PADDED_SIZE ) * value;

				filter_convolution[PADDED_SIZE- j][0] = 0.0;
				filter_convolution[PADDED_SIZE- j][1] = 0.0;
			}

			//Take the inverse FFT.
			fftw_execute(plan_backward);

			//Calculate the power and store it in result
			for (j = 0; j < n; ++j)
			{
				result[i * n + j] = MAGNITUDE(fftw_result[j][0], fftw_result[j][1]) / ( NORMALIZATION_FACTOR * sqrt(scales[i]) );
				temp[j] = result[i * n + j];
			}

			double total = 0;
			for (int j = 0; j < n; ++j)
			{
				total += temp[j];
			}

			double variance = gsl_stats_variance(temp, 1, n);
			fprintf(debug_file, "%.16f\t%.16f\t%.16f\n", SCALE_TO_FREQ(scales[i]), variance, total);

		}

		free(temp);
		fclose(debug_file);

		//FFTW sanitation engineering.
		fftw_destroy_plan(plan_backward);
	    fftw_free(fftw_result);
	    fftw_free(filter_convolution);
	// }
	//Sanitation Engineering
	fftw_destroy_plan(plan_forward);
	fftw_free(fft_data); fftw_free(data_in);
    return(0);
} /*Wavelet */

double* ShortTimeFourierTransform(double * raw_data, double sampling_frequency, int n, int window_size)
{
	//Variable Declarations
	int i, j;
	fftw_complex *data_in, *fft_data;
	fftw_plan plan_forward;
	double* result;

	// FILE * STFT_FILE = fopen("STFT_Result.log", "w");

	int num_windows = ceil((double) n / window_size);
	assert(num_windows > 0);

	// const int PADDED_SIZE = CalculatePaddingSize(window_size, 1);
	const int PADDED_SIZE = window_size;

	result =   (double*)        malloc(num_windows * PADDED_SIZE * sizeof(double));
	data_in  = (fftw_complex *) fftw_malloc( sizeof( fftw_complex ) * PADDED_SIZE );
	fft_data = (fftw_complex *) fftw_malloc( sizeof( fftw_complex ) * PADDED_SIZE );

	//Calculate the FFT of the data and store it in fft_data
	plan_forward = fftw_plan_dft_1d(PADDED_SIZE, data_in, fft_data,
									FFTW_FORWARD, FFTW_ESTIMATE);

	for (i = 0; i < num_windows; ++i)
	{
		memset(data_in,  0.0, sizeof( fftw_complex ) * PADDED_SIZE);
		memset(fft_data, 0.0, sizeof( fftw_complex ) * PADDED_SIZE);

		//Fill Data into FFT Array
		for (j = 0; j < window_size; ++j)
		{
			if (i * window_size + j < n)
			{
				data_in[j][0] = raw_data[i * window_size + j];
				data_in[j][1] = 0.0;
			}
			else
			{
				data_in[j][0] = 0.0;
				data_in[j][1] = 0.0;
			}
		}

		fftw_execute(plan_forward);

		for (j = 0; j < window_size / 2; ++j)
		{
			result[j * num_windows + i] = MAGNITUDE(fft_data[j][0], fft_data[j][1]);
			// fprintf(STFT_FILE, "%d\t%f\t%d\n", j, result[j * num_windows + i], j * num_windows + i);
		}

	}

	// fclose(STFT_FILE);
	fftw_destroy_plan(plan_forward);
	fftw_free(data_in); fftw_free(fft_data);

	return(result);
}

int Find_Peaks(double* array, double* frequency, int n, int J)
{
	FILE*       maximum_file  = fopen("maximum.log", "w");
	ARRAY_DATA *maximum_array = (ARRAY_DATA*) malloc (J * sizeof(ARRAY_DATA));
	int local_maximum_location[J];
	double* temp = (double*) malloc(n * sizeof(double));

	//Find the local maximum at every frequency
	for (int i = 0; i < J; ++i)
	{
		ARRAY_DATA max;
		max.value = array[i * n];
		max.index = i * n;

		for (int j = 0; j < n; ++j)
		{
			if (array[i * n + j] > max.value)
			{
				max.value = array[i * n + j];
				max.index = i * n + j;
			}
		}

		maximum_array[i] = max;
		fprintf(maximum_file, "%f\t%f\n", frequency[i], maximum_array[i].value);
	}

	//Calculate the deravitive of the signal and isolate the peaks
	double sign = (maximum_array[1].value - maximum_array[0].value) / (frequency[1] - frequency[0]);
	int    max_count = 0;
	for (int i = 0; i < J - 1; ++i)
	{
		double slope = (maximum_array[i + 1].value - maximum_array[i].value) / (frequency[i + 1] - frequency[i]);
		if (signbit(slope) != signbit(sign) && sign < 0)
		{
			local_maximum_location[max_count] = i;
			max_count++;
		}

		sign = slope;
	}

	for (int i = 0; i < max_count; ++i)
	{
		int arr_index = local_maximum_location[i];

		//Copy data into memory block
		for (int j = 0; j < n; ++j)
		{
			temp[j] = array[arr_index * n + j];
			// if (i == 1)
			// 	fprintf(maximum_file, "%f\t%.16f\n", (double) j/FS, array[arr_index * n + j]);

		}

		ARRAY_DATA impact_site = Max(temp, n);

		int system_setteled = 0;
		int setteled_index = 0;
		for (int j = impact_site.index; j < n; ++j)
		{
			if (temp[j] < SETTLING_PERCENTAGE * impact_site.value && system_setteled == 0)
			{
				setteled_index = j;
				double setteled_time = (double) (setteled_index - impact_site.index)/FS;
				printf("Frequency[%d]: %f, Settled Time = %f\n", i, frequency[arr_index], setteled_time);
				system_setteled = 1;
			}
		}
	}

	fclose(maximum_file);
	free(maximum_array);
	free(temp);
	return(0);
}

int CalculatePaddingSize(const int array_size, const int pad_flag)
{
	const int pad = ceil(log2(array_size));
	int out = array_size;
	switch(pad_flag)
	{
		case 0: //No Padding what so ever.
			out = array_size;
			break;
		case 1: //Zero - Padding and preforming a Radix-2 Operation
			out = (int) pow(2, pad + 1);
    		break;
    	case 2: //Duplicate array and ramp up and ramp down output
    		out = 2 * array_size;
    		break;

    	case 3: //Repeat the array once.
    		out = 2 * array_size;
    		break;

    	default: //Else return the array size
    		out = array_size;
    		break;
	}
	return(out);
}

void GenerateScalesAndFrequency(const int min_i, const int max_i, const double s_0,
	double* scales, double* frequency)
{
	int count = ( max_i - min_i ) + 1;

	//Populate the scales array
	for (int i = 0; i < count; ++i)
	{
		int counterVariable = min_i + i;
		scales[i] = s_0 * pow(2, counterVariable * D_J);
	}

	IdentifyFrequencies(scales, count, frequency);

}

void IdentifyFrequencies(double* scales, int count, double* frequency)
{
	// double * frequency = (double*) malloc( count * sizeof(double) );

	for (int i = 0; i < count; ++i)
	{
		frequency[i] = (W_0)/(scales[i] * 2 * M_PI);
	}

}

double CompleteFourierMorlet(double w, const double scale, double norm)
{
	// double norm = 1.0/sqrt(scale);
	// double norm = sqrt(scale);

	w = w * scale;
	double out = exp( -0.5 * ( W_0 - w ) * (W_0 - w) )
                                       - K_SIGMA * (exp ( -0.5 * w * w));

	out = norm * C_SIGMA * QUAD_ROOT_PI * out;
	return(out);
}

void TestCases(double *data, const int flag, double freq, double sampling_frequency, int data_size)
{
	// Fit a freq signal at two points
	// double DT = 1./sampling_frequency;
	double fsig = freq/sampling_frequency;
	double dw = 2 * M_PI * fsig;
	double w0 =  0.001; // A SMALL PHASE SHIFT SO ITS NOT ALL INTERGER ALIGNED
	int one_peri = (int)1./fsig;

	double frequency_increment = (MAX_FREQUENCY - MIN_FREQUENCY)/ 3.0; //3.0 seconds.

	int t = 2 * sampling_frequency; //At 2 seconds.

	switch(flag)
	{
		//Impulse at T = 2 seconds
		case 1:
			for (int i = 0; i < data_size; ++i)
			{
				data[i] = 0.0;
			}

			data[t] = 1.0;
			break;

		//Multiple Sines at t = 1500
		case 2:
			for (int i = data_size/2; i < data_size/2 + 2*one_peri; ++i)
			{
				data[i] = sin((i - data_size/2)* dw + w0) + sin((i - data_size/2)* 2* dw + w0);
			}
			break;

		//Multiple Sines at all times
		case 3:
			for (int i = 0; i < data_size; ++i)
			{
				data[i] = sin(i*dw + w0) + sin(i*2*dw + w0);
			}
			break;

		//Single sine at t = 1.0s;
		case 4:
			for (int i = data_size/3; i < data_size/3 + 2 * one_peri; ++i)
			{
				data[i] = sin( (i - data_size/2) * dw + w0 );
			}
			break;

		case 5:
			for (int i = 0; i < data_size; ++i)
			{
				data[i] = cos( i * dw + w0 );

				if (i >= data_size/2 && i <= 2 * (data_size)/3)
				{
					data[i] = 2 * cos(i * dw + w0);
				}
			}
			break;

		case 6:
			for (int i = 0; i < data_size; ++i)
			{
				data[i] = cos(i * dw + w0 );
				if (i >= data_size/3 && i <= data_size/2)
				{
					data[i] = cos(i * (dw - 0.005) + w0);
				}
			}
			break;

		//Frequency Sweep
		case 7:
			for (int i = 0; i < data_size; ++i)
			{
				// fsig = frequency/sampling_frequency;
				// dw = 2*M_PI*fsig;

				data[i] = sin(w0 + 2 * M_PI * (MIN_FREQUENCY + (frequency_increment/2) * pow(i/sampling_frequency, 2)) );

				// frequency += frequency_increment;
			}
			break;

		//Single sine all the way through.
		case 8:
			for (int i = 0; i < data_size; ++i)
			{
				data[i] = cos(i*dw + w0);
			}
			break;
	}
}

int WriteDebug(const double *data, const int length, const int sampling_frequency,
	const char* filename)
{
	FILE* out_file=fopen(filename,"w");
    if (out_file == NULL) return -1;

    double t = 0.0;
    double dt = 1.0/sampling_frequency;

	for (int i = 0; i < length; ++i)
    {
    	// double value = (double) i/length;
    	fprintf(out_file, "%f\t%.16e\n", t, data[i]);
    	t += dt;
    }

    fclose(out_file);
    return 0;
}

void FillData(double * data)
{
	// Fit a FREQ signal at two points
	// double DT = 1./FS;
	double fsig = FREQ/FS;
	double dw = 2*M_PI*fsig;
	double w0 =  0.01; // A SMALL PHASE SHIFT SO ITS NOT ALL INTERGER ALIGNED
	int one_peri = (int)1./fsig;
	// printf("FS  %.2d   Pitch %.f   Discrete Period = %d \n",FS,FREQ,one_peri);

	for (int i = 0; i < DATA_SIZE; ++i)
	{
		data[i] = 0.0;
	}

	// //Impulse Sample
	// data[2000] = 1.0;
	int i;
	// double t=0;
	for(i=0;i<DATA_SIZE;i++){
		// data[i] = sin(i*dw) + sin(i*dw*4);
		data[i]=0.;
		// if((i>200)&(i<400))data[i]=sin( (i-200)*dw+w0);
		if((i>0.25*DATA_SIZE)&(i<0.25*DATA_SIZE+one_peri)) data[i]=sin( (i-200)*dw+w0);
		if((i>0.5*DATA_SIZE)&(i<0.5*DATA_SIZE+2*one_peri))data[i]=sin( (i-1000)*dw+w0);
		if((i>0.75*DATA_SIZE)&(i<0.75*DATA_SIZE+3*one_peri))data[i]=sin( (i-2000)*dw+w0);
	}
}

int GetFileSize(char filename[])
{
	FILE* signalFile = fopen(filename, "r");
	assert(signalFile != NULL);
	// obtain file size:
	fseek (signalFile , 0 , SEEK_END);
	long lSize = ftell (signalFile);
	rewind (signalFile);

	char * buffer = (char*) malloc(sizeof(char)*lSize);
	assert(buffer != NULL);

	int result = fread (buffer, 1, lSize, signalFile);
	assert(result == lSize);

	char * token = strtok(buffer, "\n");

    //Get input from text.
	int counterVariable = 0;
	while (token !=NULL)
    {
    	// data[counterVariable] = atof(token);
    	counterVariable++;
        token = strtok (NULL, "\n");

    }
    fclose(signalFile);

    return (counterVariable);
}

int ReadFile(double data[], char filename[])
{
	FILE* signalFile = fopen(filename, "r");
	assert(signalFile != NULL);
	// obtain file size:
	fseek (signalFile , 0 , SEEK_END);
	long lSize = ftell (signalFile);
	rewind (signalFile);

	char * buffer = (char*) malloc(sizeof(char)*lSize);
	assert(buffer != NULL);

	int result = fread (buffer, 1, lSize, signalFile);
	assert(result == lSize);
	// puts(buffer);

	char * token = strtok(buffer, "\n");

    //Get input from text.
	int counterVariable = 0;
	while (token !=NULL)
    {
    	data[counterVariable] = atof(token);
    	counterVariable++;
        token = strtok (NULL, "\n");

    }
    fclose(signalFile);

    return (counterVariable);
}

double CWT_Cosine_Real(double time, double scale)
{
    double norm = sqrt(scale);
    double w_o = FREQ * 2 * M_PI;
    w_o *= scale;
    double mor = exp( - 0.5 * (W_0 - w_o) * (W_0 - w_o) ) + K_SIGMA * exp( -0.5 * w_o * w_o );
    double cosine = cos(w_o * time);

    double out = 0.5 * C_SIGMA * QUAD_ROOT_PI * norm * mor * cosine;
    return(out);
}

double CWT_Cosine_Complex(double time, double scale)
{
    double norm = sqrt(scale);
    double w_o = FREQ * 2 * M_PI;
    w_o *= scale;
    double mor = exp( - 0.5 * (W_0 - w_o) * (W_0 - w_o) ) + K_SIGMA * exp( -0.5 * w_o * w_o );
    double sine = sin(w_o * time);

    double out = 0.5 * C_SIGMA * QUAD_ROOT_PI * norm * mor * sine;
    return(out);
}

double CWT_Dirac_Real(double time, double scale)
{
    double impulse = 2.0;
    time = (impulse - time )/scale;
    double norm = 1.0/sqrt(scale);
    double out = exp( - 0.5 * time * time ) * ( cos( W_0 * time ) - K_SIGMA );

    out = norm * C_SIGMA * QUAD_ROOT_PI * out;
    return(out);
}

double CWT_Dirac_Complex(double time, double scale)
{
    double impulse = 2.0;
    time = (impulse - time )/scale;
    double norm = 1.0/sqrt(scale);
    double out = exp( - 0.5 * time * time ) * ( sin( W_0 * time ) - K_SIGMA );
    out = norm * C_SIGMA * QUAD_ROOT_PI * out;
    return(out);
}

ARRAY_DATA Max(double * array, int size)
{
	double max = array[0];
	int array_index = 0;

	ARRAY_DATA out;

	for (int i = 0; i < size; ++i)
	{
		if (array[i] > max)
		{
			max = array[i];
			// if (max != max)
			// {
			// 	printf("Naan Alert!\n");
			// }
			array_index = i;
		}
	}

	out.value = max;
	out.index = array_index;

	// printf("Max: Array[%d] = %.17f\n", array_index, array[array_index]);
	return(out);
}

/**
    \fn int OpenFile(const char* fileName, struct edf_hdr_struct *header)

    \brief Openes a .BDF file and allocates it to an edf_hdr_struct.

    \param fileName The name and location of the file to be opened
    \param header The pointer to the edf header structure

    \return 0 if file is opened successfully
    \return -1 if there is an error

*/
int OpenFile(const char* fileName, struct edf_hdr_struct *header)
{
    if(edfopen_file_readonly(fileName, header, EDFLIB_READ_ALL_ANNOTATIONS))
    {
        switch(header->filetype)
        {
            case EDFLIB_MALLOC_ERROR                : printf("\nmalloc error\n\n");
                                                    break;
            case EDFLIB_NO_SUCH_FILE_OR_DIRECTORY   : printf("\ncan not open file, no such file or directory\n\n");
                                                    break;
            case EDFLIB_FILE_CONTAINS_FORMAT_ERRORS : printf("\nthe file is not EDF(+) or BDF(+) compliant\n"
                                                           "(it contains format errors)\n\n");
                                                    break;
            case EDFLIB_MAXFILES_REACHED            : printf("\nto many files opened\n\n");
                                                    break;
            case EDFLIB_FILE_READ_ERROR             : printf("\na read error occurred\n\n");
                                                    break;
            case EDFLIB_FILE_ALREADY_OPENED         : printf("\nfile has already been opened\n\n");
                                                    break;
            default                                 : printf("\nunknown error\n\n");
                                                    break;
        }

    return(-1);
    }
    return(0);
}

/**
    \fn int64_t FindTriggers(const int * statusInput, const int64_t numberOfRecords,
                        int64_t * outputBuffer)

    \brief This function should take an array input and return the rising and falling edges of the triggers.

    \param statusInput: The Status Channel Input from the BDF or EDF flie. use the edfread_digital_samples
    \param numberOfRecords: The size of statusInput
    \param outputBuffer: a 1 x 2 * MAXIMUM_TRIGGERS int64_t array with the odd entries being the
                rising edge and the even entries being the falling edges.

    \return counterVariable The number of triggers that were found.
*/

int FindTriggers(const int * statusInput, const int64_t numberOfRecords,
                        int64_t * outputBuffer)
{
    int64_t counterVariable = 0;
    int counterVariable_int = 0;
    //int i needs to be int64_t because we are recording it into a int64_t arrray.
    int edge = 0;
    for (int64_t i = 0; i < numberOfRecords; ++i)
    {
        //Bit and the lower 16 bits to see if any of the triggers have been triggered to on.
        if ( ((statusInput[i] & 0x0000FFFF) > 0) && (edge == 0) && (statusInput[i-1] != statusInput[i]) ) //Rising Edge Detected.
        {
            outputBuffer[counterVariable] = i;
            // printf("outputBuffer[%" PRId64 "] = %" PRId64"\n", counterVariable, outputBuffer[counterVariable]);

            counterVariable++;
            counterVariable_int++;
            edge = 1;
        }

        if (statusInput[i-1] != statusInput[i] && edge == 1) //Falling Edge Detected.
        {
            edge = 0;
        }

        if (counterVariable > MAXIMUM_TRIGGERS)
            return -1;
    }

    return(counterVariable_int);
}

/**
    \fn int FilterTriggers(const int code, const int button, const int numberOfRecords,
            const int64_t * triggerList,
            const int * readBuffer,
            int * outputBuffer)

    \brief Filteres the triggers coming in, and finds the specified events

    \param code The code of the trigger list that is needed
                Possible Inputs: 1, or 2
    \param button The code of the button that is needed
                Possible Inputs 1, or 2
    \param triggerList The list of all of the possible triggers
    \param readBuffer The Status Channel input from the file

    \param outputBuffer The buffer that FilterTriggers will populate with the location of the location of the triggers that we're looking for

    \return counterVariable The number of triggers found.
*/

int FilterTriggers(const int code,
    const int button,
    const int numberOfRecords,
    const int64_t * triggerList,
    const int * readBuffer,
    int * outputBuffer)
{
    int readCode;
    int buttonCode;
    int counterVariable = 0;
    for (int i = 0; i < numberOfRecords; ++i)
    {
        readCode = readBuffer[i] & 255;
        buttonCode = (readBuffer[i] >> 8) & 3;

        if ( (readCode == code) && (buttonCode == button) )
        {
            outputBuffer[counterVariable] = triggerList[i];
            counterVariable++;
        }
    }
    return counterVariable;
}

void CleanData(double * data, double n)
{
    double mean = gsl_stats_mean(data, 1, n);
    double sDeviation = gsl_stats_sd_m(data, 1, n, mean);
    // printf("Mean: %f, SD: %f\n", mean, sDeviation);

    //Compute the Z-Score or Standard Score
    for (int i = 0; i < n; ++i)
    {
        data[i] = (data[i] - mean)/sDeviation;
    }
}

int  WriteFile(const double *data, const double *frequency, const int x, const int y, int sampling_frequency,
    const char filename[])
{

    FILE* out_file=fopen(filename,"w");
    if (out_file == NULL) return -1;

    //Xticks
    fprintf(out_file, "%d\t", x);
    for (int i = 0; i < y; ++i)
    {
        fprintf(out_file, "%f\t", (double) i/sampling_frequency);
    }
    fprintf(out_file, "\n");

    // double small_eps = 0.00001; //Add a small eps so that logs of zero don't happen.
    for (int i = 0; i < x; ++i)
    {
        //Feed Frequency
        fprintf(out_file, "%f\t", frequency[i]);

        //Feed Data
        for (int j = 0; j < y; ++j)
        {
            // value = Magnitude(result[i*n + j], result[i*n + j]);
            fprintf(out_file, "%.16e\t", data[i*y + j]);
        }
        //Ready for the next line.
        fprintf(out_file, "\n");

    }

    fclose(out_file);
    return(0);
}

int WriteGnuplotScript(const char graph_title[], const char filename[])
{
    FILE* gnuplot_file = fopen("script.gplot", "w");
    if (gnuplot_file == NULL) return -1;

    // fprintf(gnuplot_file, "set term x11\n");

    // fprintf(gnuplot_file, "set logscale z 10\n");

    fprintf(gnuplot_file, "set term pngcairo enhanced font 'arial,12'\n");
    fprintf(gnuplot_file, "%s%s%s","set output '", graph_title,  ".png' \n");
    fprintf(gnuplot_file, "set pm3d map\n");
    fprintf(gnuplot_file, "set logscale y 2\n");
    fprintf(gnuplot_file, "set ticslevel 0\n");
    fprintf(gnuplot_file, "set xlabel \"time (s)\"\n");
    fprintf(gnuplot_file, "set ylabel \"Frequency (Hz)\"\n");

    //Input Graph Title
    fprintf(gnuplot_file, "%s%s%s\n", "set title \"", graph_title, "\"");
    //Plot filename
    // plot "DATA.log" matrix nonuniform with pm3d t ''
    fprintf(gnuplot_file, "%s%s%s\n", "splot \"", filename, "\" matrix nonuniform with pm3d t ''");
    // fprintf(gnuplot_file, "pause -1 \"Hit Return to continue\"");

    return(0);
}

max_stars_repo_name: yuvashankar/Research

max_issues_repo_name: yuvashankar/Research

max_forks_repo_name: yuvashankar/Research


text:
#include <errno.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <time.h>
#include <clapack.h>
#include <cblas.h>
#include <omp.h>

#define DIM 16
#define NB 256

//#define DIM 96
//#define NB 128

extern void spotrf_(char *, int *, float *, int *, int *);

/*#pragma css task input(NB) inout(A[NB][NB]) highpriority */
void smpSs_spotrf_tile(float *A)
{

unsigned char LO='L';
int INFO;
int nn=NB;
   spotrf_(&LO,
          &nn,
          A,&nn,
          &INFO);
}

/*#pragma css task input(A[NB][NB], B[NB][NB], NB) inout(C[NB][NB])*/
void smpSs_sgemm_tile(float  *A, float *B, float *C)
{
unsigned char TR='T', NT='N';
float DONE=1.0, DMONE=-1.0;
    cblas_sgemm(
        CblasColMajor,
        CblasNoTrans, CblasTrans,
        NB, NB, NB,
        -1.0, A, NB,
              B, NB,
         1.0, C, NB);

}

/*#pragma css task input(T[NB][NB], NB) inout(B[NB][NB])*/
void smpSs_strsm_tile(float *T, float *B)
{
unsigned char LO='L', TR='T', NU='N', RI='R';
float DONE=1.0;

    cblas_strsm(
        CblasColMajor,
        CblasRight, CblasLower, CblasTrans, CblasNonUnit,
        NB, NB,
        1.0, T, NB,
             B, NB);

}

/*#pragma css task input(A[NB][NB], NB) inout(C[NB][NB])*/
void smpSs_ssyrk_tile( float *A, float *C)
{
unsigned char LO='L', NT='N';
float DONE=1.0, DMONE=-1.0;

    cblas_ssyrk(
        CblasColMajor,
        CblasLower,CblasNoTrans,
        NB, NB,
        -1.0, A, NB,
         1.0, C, NB);

}

void compute(struct timeval *start, struct timeval *stop, float *A[DIM][DIM])
{
#pragma omp parallel num_threads(48)
#pragma omp single
{
 gettimeofday(start,NULL);
  double t1 = omp_get_wtime();

  long j,k,i;
  for (j = 0; j < DIM; j++)
  {
    for (k= 0; k< j; k++)
    {
      for (i = j+1; i < DIM; i++)
      {
      #pragma omp task shared(A) firstprivate(i,k,j)  depend(in : A[i][k], A[j][k]) depend(inout:A[i][j])
        smpSs_sgemm_tile(  A[i][k], A[j][k], A[i][j]);
      }
    }
    for (i = 0; i < j; i++)
    {
    #pragma omp task shared(A) firstprivate(i,j) depend (in : A[j][i]) depend(inout : A[j][j])
      smpSs_ssyrk_tile( A[j][i], A[j][j]);
    }

   #pragma omp task shared(A) firstprivate(j) depend(inout : A[j][j])
    smpSs_spotrf_tile( A[j][j]);

    for (i = j+1; i < DIM; i++)
    {
  #pragma omp task shared(A) firstprivate(i,j) depend(in : A[j][j]) depend (inout : A[i][j])
      smpSs_strsm_tile( A[j][j], A[i][j]);
    }

  }
#pragma omp taskwait

  double t2 = omp_get_wtime();
  fprintf(stderr,"Time: %f\n",t2-t1);
}
 gettimeofday(stop,NULL);
 exit(0);
}

static void init(int argc, char **argv,  long *N_p);

float **A;
float * Alin;

long N;

int
main(int argc, char *argv[])
{

unsigned char LO='L';
int  INFO;

  struct timeval start;
  struct timeval stop;
  unsigned long elapsed;

  init(argc, argv, &N);
  fprintf(stderr,"Computing cholesky...\n");

  compute(&start, &stop, (void *)A);

  int nn=N;

  elapsed = 1000000 * (stop.tv_sec - start.tv_sec);
  elapsed += stop.tv_usec - start.tv_usec;

  printf ("%lu;\t", elapsed);
  printf("%d\n", (int)((0.33*N*N*N+0.5*N*N+0.17*N)/elapsed));

  printf("par_sec_time_us:%lu\n",elapsed);
  return 0;
}

static void convert_to_blocks(long N, float *Alin, float *A[DIM][DIM])
{
  long i,j;
  for (i = 0; i < N; i++)
  {
    for (j = 0; j < N; j++)
    {
      A[j/NB][i/NB][(i%NB)*NB+j%NB] = Alin[i*N+j];
    }
  }

}

void fill_random(float *Alin, int NN)
{
  int i;
  for (i = 0; i < NN; i++)
  {
    Alin[i]=((float)rand())/((float)RAND_MAX);
  }
}

static void init(int argc, char **argv,  long *N_p)
{
  long ISEED[4] = {0,0,0,1};
  long IONE=1;

  long N = NB*DIM;
  long NN = N * N;

  *N_p = N;

   Alin = (float *) malloc(NN * sizeof(float));

  fill_random(Alin,NN);
  long i;
  for(i=0; i<N; i++)
  {
    Alin[i*N + i] += N;
  }

  A = (float **) malloc(DIM*DIM*sizeof(float *));
  for ( i = 0; i < DIM*DIM; i++)
    {
//     #pragma omp memory
     A[i] = (float *) malloc(NB*NB*sizeof(float));
      int z;
      for (z=0;z<DIM*DIM;z++)
	   { float *zz = (float *) A[i];
	     zz[z] = Alin[(DIM*DIM)*i + z]; }
    }

  convert_to_blocks(N, Alin, (void *)A);

}

max_stars_repo_name: podobas/BLYSK

max_issues_repo_name: podobas/BLYSK

max_forks_repo_name: podobas/BLYSK