text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#ifdef _WIN32
    #include "blas.c"
#else
    #include "blas.c"
    // #include <cblas.h>
#endif
#include <moss.h>
#include <objects/string.h>
#include <objects/list.h>
#include <modules/bs.h>
#include <modules/la.h>

double mf_float(mt_object* x, int* e);
extern mt_table* mv_type_array;

void mf_array_delete(mt_array* a){
    if(a->data->refcount==1){
        mf_free(a->data);
    }else{
        a->data->refcount--;
    }
    mf_free(a);
}

void mf_array_dec_refcount(mt_array* a){
    if(a->refcount==1){
        mf_array_delete(a);
    }else{
        a->refcount--;
    }
}

static
mt_array* mf_new_array(long size, int type){
    mt_array* a = mf_malloc(sizeof(mt_array));
    if(type==mv_array_float){
        a->data = mf_malloc(sizeof(mt_array_data)+size*sizeof(double));
    }else{
        abort();
    }
    a->data->refcount = 1;
    a->data->type = type;
    a->base = a->data->a;
    a->refcount = 1;
    a->size = size;
    return a;
}

static
mt_array* mf_array_copy(mt_array* a){
    mt_array* b;
    if(a->n==2){
        if(a->plain){
            goto plain;
        }else{
            unsigned long m=a->shape[0];
            unsigned long n=a->shape[1];
            long si=a->stride[0];
            long sj=a->stride[1];
            b = mf_new_array(m*n,mv_array_float);
            b->n = a->n; b->plain = 1;
            b->shape[0] = m; b->shape[1] = n;
            b->stride[0] = 1; b->stride[1] = m;
            double* pa = (double*)a->base;
            double* pb = (double*)b->base;
            unsigned long i,j,mj;
            long sjj;
            for(j=0; j<n; j++){
                mj=m*j; sjj=sj*j;
                for(i=0; i<m; i++){
                    pb[i+mj]=pa[si*i+sjj];
                }
            }
            return b;
        }
    }else if(a->n==1){
        b = mf_new_array(a->shape[0],mv_array_float);
        b->n = 1;
        b->shape[0] = a->shape[0];
        b->stride[0] = 1;
        cblas_dcopy(a->shape[0],(double*)a->base,a->stride[0],(double*)b->base,1);
        return b;
    }else{
        if(a->plain) goto plain;
        abort();
    }

    plain:
    b = mf_new_array(a->size,mv_array_float);
    b->n = a->n;
    b->plain = 1;
    unsigned long i;
    for(i=0; i<a->n; i++){
        b->shape[i] = a->shape[i];
        b->stride[i] = a->stride[i];
    }
    cblas_dcopy(a->size,(double*)a->base,1,(double*)b->base,1);
    return b;
}

void mf_ensure_plain(mt_array* a){
    if(a->n==2){
        if(a->plain) return;
        mt_array_data* data = mf_malloc(sizeof(mt_array_data)+a->size*sizeof(double));
        data->refcount = 1;
        data->type = mv_array_float;
        double* pa = (double*)a->base;
        double* pb = (double*)data->a;

        unsigned long m = a->shape[0];
        unsigned long n = a->shape[1];
        long si = a->stride[0];
        long sj = a->stride[1];
        unsigned long i,j,mj;
        long sjj;
        for(j=0; j<n; j++){
            mj = m*j; sjj = sj*j;
            for(i=0; i<m; i++){
                pb[i+mj] = pa[si*i+sjj];
            }
        }
        if(a->data->refcount==1){
            mf_free(a->data);
        }else{
            a->data->refcount--;
        }
        a->stride[0] = 1;
        a->stride[1] = m;
        a->data = data;
        a->base = data->a;
        a->plain = 1;
    }else if(a->n==1){
        abort();
    }else{
        if(a->plain) return;
        abort();
    }
}

static
mt_array* mf_array_map_dd(mt_array* a, double(*f)(double)){
    mt_array* b;
    if(a->n==2){
        if(a->plain){
            goto plain;
        }else{
            unsigned long m = a->shape[0];
            unsigned long n = a->shape[1];
            long si = a->stride[0];
            long sj = a->stride[1];
            b = mf_new_array(m*n,mv_array_float);
            b->n = a->n; b->plain = 1;
            b->shape[0] = m; b->shape[1] = n;
            b->stride[0] = 1; b->stride[1] = m;
            double* pa = (double*)a->base;
            double* pb = (double*)b->base;
            unsigned long i,j,mj;
            long sjj;
            for(j=0; j<n; j++){
                mj = m*j; sjj = sj*j;
                for(i=0; i<m; i++){
                    pb[i+mj] = f(pa[si*i+sjj]);
                }
            }
            return b;
        }
    }else if(a->n==1){
        unsigned long n=a->shape[0];
        b = mf_new_array(n,mv_array_float);
        b->n = 1; b->shape[0] = n;
        b->stride[0] = 1;
        long s = a->stride[0];
        unsigned long k;
        double* pa = (double*)a->base;
        double* pb = (double*)b->base;
        for(k=0; k<n; k++){
            pb[k] = f(pa[k*s]);
        }
        return b;
    }else{
        if(a->plain) goto plain;
        abort();
    }

    plain:
    b = mf_new_array(a->size,mv_array_float);
    b->n = a->n;
    b->plain = 1;
    unsigned long i;
    for(i=0; i<a->n; i++){
        b->shape[i] = a->shape[i];
        b->stride[i] = a->stride[i];
    }
    double* pa = (double*)a->base;
    double* pb = (double*)b->base;
    unsigned long size = a->size;
    unsigned long k;
    for(k=0; k<size; k++){
        pb[k] = f(pa[k]);
    }
    return b;
}

long mf_array_size(mt_array* a){
    long p;
    unsigned long k;
    switch(a->n){
    case 1:
        return a->shape[0];
    case 2:
        return a->shape[0]*a->shape[1];
    default:
        p=1;
        for(k=0; k<a->n; k++){
            p = p*a->shape[k];
        }
        return p;
    }
}

mt_array* mf_array(long size, mt_object* v){
    mt_array* a;
    if(size>0 && v[0].type==mv_list){
        long i,j;
        mt_list* list = (mt_list*)v[0].value.p;
        long n = list->size;
        a = mf_new_array(size*n,mv_array_float);
        a->n = 2;
        a->shape[0] = size; a->shape[1] = n;
        a->stride[0] = 1; a->stride[1] = size;
        a->plain = 1;
        double* b = (double*)a->base;
        int e=0;
        for(i=0; i<size; i++){
            if(v[i].type!=mv_list) abort();
            list = (mt_list*)v[i].value.p;
            if(list->size!=n) abort();
            for(j=0; j<n; j++){
                b[i+size*j] = mf_float(list->a+j,&e);
                if(e) abort();
            }
        }
    }else{
        a = mf_new_array(size,mv_array_float);
        a->n = 1;
        a->shape[0] = size;
        a->stride[0] = 1;
        long i;
        double* b = (double*)a->base;
        int e=0;
        for(i=0; i<size; i++){
            b[i] = mf_float(v+i,&e);
        }
    }
    return a;
}

int la_array(mt_object* x, int argc, mt_object* v){
    if(argc!=1){
        mf_argc_error(argc,1,1,"array");
        return 1;
    }
    if(v[1].type!=mv_list){
        mf_type_error1("in array(a): a (type: %s) is not a list.",v+1);
        return 1;
    }
    mt_list* list = (mt_list*)v[1].value.p;
    mt_array* a = mf_array(list->size, list->a);
    if(a==NULL) return 1;
    x->type = mv_array;
    x->value.p = (mt_basic*)a;
    return 0;
}

int la_vector(mt_object* x, int argc, mt_object* v){
    mt_array* a = mf_array(argc,v+1);
    if(a==NULL) return 1;
    x->type = mv_array;
    x->value.p = (mt_basic*)a;
    return 0;
}

int la_matrix(mt_object* x, int argc, mt_object* v){
    mt_array* a = mf_array(argc,v+1);
    if(a==NULL) return 1;
    if(a->n==1){
        a->n = 2;
        unsigned long size = a->shape[0];
        a->shape[0] = size;
        a->shape[1] = 1;
        a->stride[0] = 1;
        a->stride[1] = size;
        a->plain = 0;
    }
    x->type = mv_array;
    x->value.p = (mt_basic*)a;
    return 0;
}

static
void vector_str(mt_bs* bs, mt_array* T){
    char buffer[100];
    unsigned long i;
    double* a = (double*)T->base;
    unsigned long size = T->shape[0];
    long step = T->stride[0];

    mf_bs_push_cstr(bs,"[");
    for(i=0; i<size; i++){
        snprintf(buffer,100,i==0?"%.4g":", %.4g",a[step*i]);
        mf_bs_push_cstr(bs,buffer);
    }
    mf_bs_push_cstr(bs,"]");
}

static
void matrix_str(mt_bs* bs, mt_array* T){
    char buffer[100];
    unsigned long i,j;
    double* a = (double*)T->base;
    unsigned long m,n,is,js;
    m = T->shape[0];
    n = T->shape[1];
    is = T->stride[0];
    js = T->stride[1];
    int max = 0;
    int size;
    for(i=0; i<m; i++){
        for(j=0; j<n; j++){
            size = snprintf(buffer,100,"%.4g",a[i*is+j*js]);
            if(size>max) max=size;
        }
    }
    mf_bs_push_cstr(bs,"[");
    for(i=0; i<m; i++){
        mf_bs_push_cstr(bs,i==0? "[": ",\n [");
        for(j=0; j<n; j++){
            snprintf(buffer,100,j==0?"%*.4g":", %*.4g",max,a[i*is+j*js]);
            mf_bs_push_cstr(bs,buffer);
        }
        mf_bs_push_cstr(bs,"]");
    }
    mf_bs_push_cstr(bs,"]");
}

static
int array_str(mt_object* x, int argc, mt_object* v){
    if(argc!=0){
        mf_argc_error(argc,0,0,"array.type.str");
        return 1;
    }
    if(v[0].type!=mv_array){
        mf_type_error1("in a.str(): a (type: %s) is not an array.",&v[0]);
        return 1;
    }
    mt_array* T = (mt_array*)v[0].value.p;
    mt_bs bs;
    mf_bs_init(&bs);
    if(T->n==1){
        vector_str(&bs,T);
    }else if(T->n==2){
        matrix_str(&bs,T);
    }else{
        abort();
    }

    x->type = mv_string;
    x->value.p = (mt_basic*)mf_str_new(bs.size,(const char*)bs.a);
    mf_free(bs.a);
    return 0;
}

static
int array_list(mt_object* x, int argc, mt_object* v){
    if(argc!=0){
        mf_argc_error(argc,0,0,"array.type.list");
        return 1;
    }
    if(v[0].type!=mv_array){
        mf_type_error1("in a.list(): a (type: %s) is not an array.",&v[0]);
        return 1;
    }
    mt_array* T = (mt_array*)v[0].value.p;
    mt_list* list = mf_raw_list(T->size);
    unsigned long i;
    double* a = (double*)T->base;
    unsigned long size = T->shape[0];
    long step = T->stride[0];
    for(i=0; i<size; i++){
        list->a[i].type = mv_float;
        list->a[i].value.f = a[step*i];
    }

    x->type = mv_list;
    x->value.p = (mt_basic*)list;
    return 0;
}

double addc;
double add(double x){
    return x+addc;
}
double subc;
double sub(double x){
    return x-subc;
}

static
int array_add(mt_object* x, int argc, mt_object* v){
    if(argc!=1){
        mf_argc_error(argc,1,1,"Array.add");
        return 1;
    }
    if(v[0].type!=mv_array){
        mf_type_error1("in a+b: a (type: %s) is not an array.",&v[0]);
        return 1;
    }
    mt_array* a = (mt_array*)v[0].value.p;
    mt_array* c;
    if(v[1].type!=mv_array){
        int e = 0;
        addc = mf_float(&v[1],&e);
        if(e){
            mf_type_error1("in a+r: cannot convert r (type: %s) to float.",&v[1]);
            return 1;
        }
        c = mf_array_map_dd(a,add);
        if(c==NULL) return 1;
        x->type = mv_array;
        x->value.p = (mt_basic*)c;
        return 0;
    }
    mt_array* b = (mt_array*)v[1].value.p;
    if(a->n!=b->n){
        mf_value_error("value error in a+b: a and b have unequal order.");
        return 1;
    }
    if(a->n==1){
        unsigned long size=a->shape[0];
        if(size!=b->shape[0]){
            mf_value_error("Value error in a+b: a and b have unequal shape.");
            return 1;
        }
        c = mf_new_array(size,mv_array_float);
        c->n = 1; c->stride[0] = 1;
        c->shape[0] = size;
        cblas_dcopy(size,(double*)a->base,a->stride[0],(double*)c->base,1);
        cblas_daxpy(size,1,(double*)b->base,b->stride[0],(double*)c->base,1);
    }else{
        unsigned long i;
        for(i=0; i<a->n; i++){
            if(a->shape[i]!=b->shape[i]) abort();
        }
        c = mf_array_copy(a);
        mf_ensure_plain(b);
        cblas_daxpy(c->size,1,(double*)b->base,1,(double*)c->base,1);
    }
    x->type = mv_array;
    x->value.p = (mt_basic*)c;
    return 0;
}

static
int array_radd(mt_object* x, int argc, mt_object* v){
    if(argc!=1){
        mf_argc_error(argc,1,1,"Array.radd");
        return 1;
    }
    if(v[1].type!=mv_array){
        mf_type_error1("in r*a: a (type: %s) is not an array.",&v[1]);
        return 1;
    }
    int e = 0;
    addc = mf_float(&v[0],&e);
    if(e){
        mf_type_error1("in r+a: cannot convert r (type: %s) to float.",&v[0]);
        return 1;
    }
    mt_array* a = (mt_array*)v[1].value.p;
    mt_array* y = mf_array_map_dd(a,add);
    if(y==NULL) return 1;
    x->type = mv_array;
    x->value.p = (mt_basic*)y;
    return 0;
}

static
int array_sub(mt_object* x, int argc, mt_object* v){
    if(argc!=1){
        mf_argc_error(argc,1,1,"Array.sub");
        return 1;
    }
    if(v[0].type!=mv_array){
        mf_type_error1("in a-b: a (type: %s) is not an array.",&v[0]);
        return 1;
    }
    mt_array* a = (mt_array*)v[0].value.p;
    mt_array* c;
    if(v[1].type!=mv_array){
        int e = 0;
        subc = mf_float(&v[1],&e);
        if(e){
            mf_type_error1("in a-r: cannot convert r (type: %s) to float.",&v[1]);
            return 1;
        }
        c = mf_array_map_dd(a,sub);
        if(c==NULL) return 1;
        x->type = mv_array;
        x->value.p = (mt_basic*)c;
        return 0;
    }
    mt_array* b = (mt_array*)v[1].value.p;
    if(a->n==1){
        unsigned long size = a->shape[0];
        if(size!=b->shape[0]){
            mf_value_error("Value error in a-b: a and b have unequal shape.");
            return 1;
        }

        c = mf_new_array(size,mv_array_float);
        c->n = 1;
        c->shape[0] = size;
        c->stride[0] = 1;
        cblas_dcopy(size,(double*)a->base,a->stride[0],(double*)c->base,1);
        cblas_daxpy(size,-1,(double*)b->base,b->stride[0],(double*)c->base,1);
    }else{
        unsigned long i;
        for(i=0; i<a->n; i++){
            if(a->shape[i]!=b->shape[i]) abort();
        }
        c = mf_array_copy(a);
        mf_ensure_plain(b);
        cblas_daxpy(c->size,-1,(double*)b->base,1,(double*)c->base,1);
    }
    x->type = mv_array;
    x->value.p = (mt_basic*)c;
    return 0;
}

static
mt_array* mf_array_scal(double r, mt_array* a){
    mt_array* b;
    if(a->n==1){
        unsigned long size = a->shape[0];
        b = mf_new_array(size,mv_array_float);
        b->n = 1;
        b->shape[0] = size;
        b->stride[0] = 1;
        cblas_dcopy(size,(double*)a->base,a->stride[0],(double*)b->base,1);
    }else{
        b = mf_array_copy(a);
    }
    cblas_dscal(b->size,r,(double*)b->base,1);
    return b;
}

static
mt_array* array_mpy_array(mt_array* a, mt_array* b){
    if(b->n==1){
        unsigned long m = a->shape[0];
        unsigned long n = a->shape[1];
        if(n!=b->shape[0]){
            mf_value_error("Value error in A*x: shape does not match.");
            return NULL;
        }
        if(a->stride[0]!=1){
            mf_ensure_plain(a);
        }
        mt_array* c = mf_new_array(m,mv_array_float);
        c->n = 1;
        c->shape[0] = m;
        c->stride[0] = 1;
        cblas_dgemv(CblasColMajor,CblasNoTrans,m,n,1,
            (double*)a->base, a->stride[1],
            (double*)b->base, b->stride[0],
            0,(double*)c->base,1
        );
        return c;
    }else if(b->n==2){
        if(a->n==2){
            if(a->shape[1]!=b->shape[0]){
                mf_value_error("Value error in A*B: incompatible matrices.");
                return NULL;
            }
            if(a->stride[0]!=1){
                mf_ensure_plain(a);
            }
            if(b->stride[0]!=1){
                mf_ensure_plain(b);
            }

            unsigned long m = a->shape[0];
            unsigned long k = a->shape[1];
            unsigned long n = b->shape[1];
            mt_array* c = mf_new_array(m*n,mv_array_float);
            c->n = 2;
            c->shape[0] = m;
            c->shape[1] = n;
            c->stride[0] = 1;
            c->stride[1] = m;
            c->plain = 1;
            cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n, k, 1,
                (double*)a->base, a->stride[1],
                (double*)b->base, b->stride[1],
                0, (double*)c->base, m
            );
            return c;
        }else{
            abort();
        }
    }else{
        abort();
    }
}

static
mt_array* matrix_identity(unsigned long n){
    mt_array* a = mf_new_array(n*n,mv_array_float);
    a->n = 2;
    a->plain = 1;
    a->shape[0] = n;
    a->shape[1] = n;
    a->stride[0] = 1;
    a->stride[1] = n;
    double* b = (double*)a->base;
    unsigned long i,j,nj;
    for(j=0; j<n; j++){
        nj = n*j;
        for(i=0; i<n; i++){
            b[i+nj] = i==j? 1: 0;
        }
    }
    return a;
}

static
mt_array* square_matrix_mpy(mt_array* a, mt_array* b){
    if(a->stride[0]!=1){
        mf_ensure_plain(a);
    }
    if(b->stride[0]!=1){
        mf_ensure_plain(b);
    }
    unsigned long n = a->shape[0];
    mt_array* c = mf_new_array(n*n,mv_array_float);
    c->n = 2;
    c->shape[0] = n;
    c->shape[1] = n;
    c->stride[0] = 1;
    c->stride[1] = n;
    c->plain = 1;
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, n, n, n, 1,
        (double*)a->base, a->stride[1],
        (double*)b->base, b->stride[1],
        0, (double*)c->base, n
    );
    return c;
}

static
mt_array* square_matrix_pow(mt_array* a, unsigned long n){
    if(n==0){
        return matrix_identity(a->shape[0]);
    }
    mt_array *p,*h;
    p = mf_array_copy(a);
    a = p;
    a->refcount++;
    n--;
    while(n){
        if(n&1){
            h = p;
            p = square_matrix_mpy(p,a);
            mf_array_dec_refcount(h);
        }
        n>>=1;
        if(n==0) break;
        h = a;
        a = square_matrix_mpy(a,a);
        mf_array_dec_refcount(h);
    }
    mf_array_dec_refcount(a);
    return p;
}

static
int array_mul(mt_object* x, int argc, mt_object* v){
    if(argc!=1){
        mf_argc_error(argc,1,1,"Array.mul");
        return 1;
    }
    if(v[0].type!=mv_array){
        mf_type_error1("in a*b: a (type: %s) is not an array.",&v[0]);
        return 1;
    }
    mt_array* a = (mt_array*)v[0].value.p;
    mt_array* b;
    double r;
    if(v[1].type==mv_array){
        b=(mt_array*)v[1].value.p;
        if(a->n==2){
            if(b->n==2 && a->shape[0]==1 && b->shape[1]==1){
                if(a->shape[1]!=b->shape[0]){
                    mf_value_error("Value error in a*b: incompatible matrices.");
                    return 1;
                }
                // todo: strides
                r = cblas_ddot(a->shape[1],
                    (double*)a->base,a->stride[1],
                    (double*)b->base,b->stride[0]
                );
                x->type = mv_float;
                x->value.f = r;
                return 0;
            }
            mt_array* c = array_mpy_array(a,b);
            if(c==NULL) return 1;
            x->type = mv_array;
            x->value.p = (mt_basic*)c;
            return 0;
        }
        if(a->n!=1){
            mf_value_error("Value error in a*b: a and b have incompatible shape.");
            return 1;
        }
        if(b->n==2){
            mt_array* c = array_mpy_array(b,a);
            if(c==NULL) return 1;
            x->type = mv_array;
            x->value.p = (mt_basic*)c;
            return 0;
        }
        if(a->shape[0]!=b->shape[0]){
            mf_value_error("Value error in a*b: a and b have unequal shape.");
            return 1;
        }
        r = cblas_ddot(a->shape[0],
            (double*)a->base,a->stride[0],
            (double*)b->base,b->stride[0]
        );
        x->type = mv_float;
        x->value.f = r;
        return 0;
    }else{
        int e = 0;
        r = mf_float(&v[1],&e);
        if(e){
            mf_type_error1("in a*r: cannot convert r (type: %s) to float.",&v[1]);
            return 1;
        }
        x->type = mv_array;
        x->value.p = (mt_basic*)mf_array_scal(r,a);
        return 0;
    }
}

static
int array_pow(mt_object* x, int argc, mt_object* v){
    if(argc!=1){
        mf_argc_error(argc,1,1,"Array.pow");
        return 1;
    }
    if(v[0].type!=mv_array){
        mf_type_error1("in a^n: a (type: %s) is not an array.",&v[0]);
        return 1;
    }
    mt_array* a = (mt_array*)v[0].value.p;

    double r;
    if(a->n!=1){
        if(a->n==2 && v[1].type==mv_int){
            long n = v[1].value.i;
            if(a->shape[0]!=a->shape[1]){
                mf_value_error("Value error in A^n: A is not a square matrix");
                return 1;
            }
            x->type = mv_array;
            x->value.p = (mt_basic*)square_matrix_pow(a,(unsigned long)n);
            return 0;
        }
        mf_value_error("Value error in a^n: a is not a vector.");
        return 1;
    }
    r = cblas_ddot(a->shape[0],
        (double*)a->base,a->stride[0],
        (double*)a->base,a->stride[0]
    );
    x->type = mv_float;
    if(v[1].type==mv_int && v[1].value.i==2){
        x->value.f = r;
    }else{
        int e = 0;
        double n = mf_float(&v[1],&e);
        if(e){
            mf_type_error1("in a^n: cannot convert n (type: %s) to float.",&v[1]);
            return 1;
        }
        x->value.f = pow(r,n);
    }
    return 0;
}

static
int array_rmul(mt_object* x, int argc, mt_object* v){
    if(argc!=1){
        mf_argc_error(argc,1,1,"Array.rmul");
        return 1;
    }
    if(v[1].type!=mv_array){
        mf_type_error1("in r*a: a (type: %s) is not an array.",&v[1]);
        return 1;
    }
    int e=0;
    double r = mf_float(&v[0],&e);
    if(e){
        mf_type_error1("in r*a: cannot convert r (type: %s) to float.",&v[0]);
        return 1;
    }
    mt_array* a = (mt_array*)v[1].value.p;
    x->type = mv_array;
    x->value.p = (mt_basic*)mf_array_scal(r,a);;
    return 0;
}

static
int array_neg(mt_object* x, int argc, mt_object* v){
    if(argc!=0){
        mf_argc_error(argc,0,0,"Array.neg");
        return 1;
    }
    if(v[0].type!=mv_array){
        mf_type_error1("in -a: a (type: %s) is not an array.",&v[1]);
        return 1;
    }
    mt_array* a = (mt_array*)v[0].value.p;
    x->type = mv_array;
    x->value.p = (mt_basic*)mf_array_scal(-1,a);;
    return 0;
}

static
int array_div(mt_object* x, int argc, mt_object* v){
    if(argc!=1){
        mf_argc_error(argc,1,1,"Array.div");
        return 1;
    }
    if(v[0].type!=mv_array){
        mf_type_error1("in a/r: a (type: %s) is not an array.",&v[1]);
        return 1;
    }
    int e = 0;
    double r = mf_float(&v[1],&e);
    if(e){
        mf_type_error1("in a/r: cannot convert r (type: %s) to float.",&v[0]);
        return 1;
    }
    mt_array* a = (mt_array*)v[0].value.p;
    x->type = mv_array;
    x->value.p = (mt_basic*)mf_array_scal(1/r,a);;
    return 0;
}

static
int array_T(mt_object* x, int argc, mt_object* v){
    if(argc!=0){
        mf_argc_error(argc,0,0,"array.type.T");
        return 1;
    }
    if(v[0].type!=mv_array){
        mf_type_error1("in a.T(): a (type: %s) is not an array.",&v[0]);
        return 1;
    }
    mt_array* a = (mt_array*)v[0].value.p;
    if(a->n==2){
        mt_array* b = mf_malloc(sizeof(mt_array));
        b->refcount = 1;
        b->n = 2;
        a->data->refcount++;
        b->data = a->data;
        b->base = a->base;
        b->size = a->size;
        b->plain = 0;
        b->shape[0] = a->shape[1];
        b->shape[1] = a->shape[0];
        b->stride[0] = a->stride[1];
        b->stride[1] = a->stride[0];
        x->type = mv_array;
        x->value.p = (mt_basic*)b;
    }else{
        a->refcount++;
        x->type = mv_array;
        x->value.p = (mt_basic*)a;
    }
    return 0;
}

static
int array_plain(mt_object* x, int argc, mt_object* v){
    if(v[0].type!=mv_array) abort();
    mt_array* a = (mt_array*)v[0].value.p;
    mf_ensure_plain(a);
    a->refcount++;
    x->type = mv_array;
    x->value.p = (mt_basic*)a;
    return 0;
}

static
int array_copy(mt_object* x, int argc, mt_object* v){
    if(v[0].type!=mv_array) abort();
    mt_array* a = (mt_array*)v[0].value.p;
    x->type = mv_array;
    x->value.p = (mt_basic*)mf_array_copy(a);
    return 0;
}

static
mt_array* vector_slice(mt_array* a, mt_range* r){
    long i,j;
    if(r->a.type!=mv_int || r->b.type!=mv_int){
        mf_type_error("Type error in a[r]: r is not an integer range.");
        return NULL;
    }
    i = r->a.value.i;
    j = r->b.value.i;
    unsigned long n = a->shape[0];
    if(i<0 || j<0){
        mf_value_error("Value error in a[r]: range with negative values.");
        return NULL;
    }
    if((unsigned long)i>=n || (unsigned long)j>=n){
        mf_value_error("Value error in a[r]: range is out of bounds.");
        return NULL;
    }
    mt_array* b = mf_malloc(sizeof(mt_array));
    b->refcount = 1;
    b->size = a->size;
    a->data->refcount++;
    b->data = a->data;
    a->data = b->data;
    b->n = 1;
    b->base = (unsigned char*)((double*)a->base+a->stride[0]*i);
    b->shape[0] = j-i+1;
    b->stride[0] = a->stride[0];
    return b;
}

static
int array_get(mt_object* x, int argc, mt_object* v){
    if(v[0].type!=mv_array){
        mf_type_error1("in a[i]: a (type: %s) is not an array.",&v[0]);
        return 1;
    }
    mt_array* a = (mt_array*)v[0].value.p;
    if(argc==1){
        if(a->n!=1){
            mf_value_error("Value error in a[i]: a is not a vector.");
            return 1;
        }
        if(v[1].type!=mv_int){
            if(v[1].type==mv_range){
                 mt_array* b = vector_slice(a,(mt_range*)v[1].value.p);
                 if(b==NULL) return 1;
                 x->type = mv_array;
                 x->value.p = (mt_basic*)b;
                 return 0;
            }
            mf_type_error1("in a[i]: i (type: %s) is not an integer.",&v[1]);
            return 1;
        }
        long index = v[1].value.i;
        unsigned long n = a->shape[0];
        if(index<0 || (unsigned long)index>=n){
            mf_value_error("Value error in a[i]: i is out of bounds.");
            return 1;
        }
        double* b = (double*)a->base;
        double t = b[index*a->stride[0]];
        x->type = mv_float;
        x->value.f = t;
        return 0;
    }else{
        mf_argc_error(argc,1,1,"Array.get");
        return 1;
    }
}

static
int la_idm(mt_object* x, int argc, mt_object* v){
    if(argc!=1){
        mf_argc_error(argc,1,1,"idm");
        return 1;
    }
    if(v[1].type!=mv_int){
        mf_type_error1("in idm(n): n (type: %s) is not an integer.",v+1);
        return 1;
    }
    long n = v[1].value.i;
    x->type = mv_array;
    x->value.p = (mt_basic*)matrix_identity(n);
    return 0;
}

static
mt_array* diag(mt_array* v){
    if(v->n!=1){
        mf_value_error("Value error in diag(v): v is not a vector.");
        return NULL;
    }
    unsigned long n = v->shape[0];
    mt_array* a = mf_new_array(n*n,mv_array_float);
    a->n = 2;
    a->plain = 1;
    a->shape[0] = n;
    a->shape[1] = n;
    a->stride[0] = 1;
    a->stride[1] = n;
    double* b = (double*)a->base;
    double* vb = (double*)v->base;
    long vs = v->stride[0];
    unsigned long i,j,nj;
    for(j=0; j<n; j++){
        nj = n*j;
        for(i=0; i<n; i++){
            b[i+nj] = i==j? vb[i*vs]: 0;
        }
    }
    return a;
}

static
int la_diag(mt_object* x, int argc, mt_object* v){
    if(argc!=1){
        mf_argc_error(argc,1,1,"diag");
        return 1;
    }
    if(v[1].type!=mv_array){
        mf_type_error1("in diag(a): a is not an array.",v+1);
        return 1;
    }
    mt_array* u = (mt_array*)v[1].value.p;
    mt_array* a = diag(u);
    if(a==NULL) return 1;
    x->type = mv_array;
    x->value.p = (mt_basic*)a;
    return 0;
}

static
mt_array* diag_slice(mt_array* a){
    if(a->n!=2 || a->shape[0]!=a->shape[1]){
        mf_value_error("Value error in a.diag(): a is not an quadratic matrix.");
        return NULL;
    }
    if(a->stride[1]!=1){
        mf_ensure_plain(a);
    }
    mt_array* v = mf_malloc(sizeof(mt_array));
    v->refcount = 1;
    v->n = 1;
    a->data->refcount++;
    v->data = a->data;
    v->base = a->base;
    v->shape[0] = a->shape[0];
    v->stride[0] = (a->shape[0]+1)*a->stride[0];
    return v;
}

static
int array_diag(mt_object* x, int argc, mt_object* v){
    if(argc!=0){
        mf_argc_error(argc,0,0,"diag");
        return 1;
    }
    if(v[0].type!=mv_array){
        mf_type_error1("in a.diag(): a (type: %s) is not an array.",&v[0]);
        return 1;
    }
    mt_array* a = (mt_array*)v[0].value.p;
    mt_array* u = diag_slice(a);
    if(u==NULL) return 1;
    x->type = mv_array;
    x->value.p = (mt_basic*)u;
    return 0;
}

int mf_array_map(mt_array* a, mt_function* f){
    mt_object argv[2];
    argv[0].type = mv_null;
    double* ba = (double*)a->base;
    long size = a->size;
    long i;
    mt_object x;
    int e = 0;
    for(i=0; i<size; i++){
        argv[1].type = mv_float;
        argv[1].value.f = ba[i];
        if(mf_call(f,&x,1,argv)){
            mf_traceback("map");
            return 1;
        }
        ba[i] = mf_float(&x,&e);
        if(e){
            mf_type_error1("in a.map(f): cannot convert f(a[k]) (type: %s) to float.",&x);
            return 1;
        }
    }
    return 0;
}

static
int array_map(mt_object* x, int argc, mt_object* v){
    if(argc!=1){
        mf_argc_error(argc,1,1,"map");
        return 1;
    }
    if(v[0].type!=mv_array){
        mf_type_error1("in a.map(f): a (type: %s) is not an array.",&v[0]);
        return 1;
    }
    if(v[1].type!=mv_function){
        mf_type_error1("in a.map(f): f (type: %s) is not a function.",&v[1]);
        return 1;
    }
    mt_array* a = (mt_array*)v[0].value.p;
    mt_function* f = (mt_function*)v[1].value.p;
    mt_array* b = mf_array_copy(a);
    if(mf_array_map(b,f)) return 1;
    x->type = mv_array;
    x->value.p = (mt_basic*)b;
    return 0;
}

static
int mf_trace(mt_object* x, mt_array* a){
    if(a->n!=2 || a->shape[0]!=a->shape[1]){
        mf_value_error("Value error in trace(a): a is not a quadratic matrix.");
        return 1;
    }
    double s = 0;
    long n = a->shape[0];
    double* ba = (double*)a->base;
    long i,m;
    m = a->stride[0]+a->stride[1];
    for(i=0; i<n; i++){
        s+=ba[i*m];
    }
    x->type = mv_float;
    x->value.f = s;
    return 0;
}

static
int la_trace(mt_object* x, int argc, mt_object* v){
    if(argc!=1){
        mf_argc_error(argc,1,1,"trace");
        return 1;
    }
    if(v[1].type!=mv_array){
        mf_type_error1("in la.trace(a): a (type: %s) is not an array.",&v[1]);
        return 1;
    }
    mt_array* a = (mt_array*)v[1].value.p;
    if(mf_trace(x,a)) return 1;
    return 0;
}

mt_table* mf_la_load(){
    mt_map* m = mv_type_array->m;
    mf_insert_function(m,0,0,"str",array_str);
    mf_insert_function(m,0,0,"list",array_list);
    mf_insert_function(m,1,1,"add",array_add);
    mf_insert_function(m,1,1,"radd",array_radd);
    mf_insert_function(m,1,1,"sub",array_sub);
    mf_insert_function(m,1,1,"mul",array_mul);
    mf_insert_function(m,1,1,"rmul",array_rmul);
    mf_insert_function(m,0,0,"neg",array_neg);
    mf_insert_function(m,1,1,"div",array_div);
    mf_insert_function(m,1,1,"pow",array_pow);
    mf_insert_function(m,1,-1,"get",array_get);
    mf_insert_function(m,0,0,"T",array_T);
    mf_insert_function(m,0,0,"plain",array_plain);
    mf_insert_function(m,0,0,"copy",array_copy);
    mf_insert_function(m,0,0,"diag",array_diag);
    mf_insert_function(m,1,1,"map",array_map);

    mt_table* la = mf_table(NULL);
    la->name = mf_cstr_to_str("module la");
    la->m = mf_empty_map();
    m = la->m;
    mf_insert_function(m,1,1,"array",la_array);
    mf_insert_function(m,1,1,"vector",la_vector);
    mf_insert_function(m,1,1,"matrix",la_matrix);
    mf_insert_function(m,1,1,"idm",la_idm);
    mf_insert_function(m,1,1,"diag",la_diag);
    mf_insert_function(m,1,1,"trace",la_trace);
    m->frozen = 1;
    return la;
}

max_stars_repo_name: JohnBSmith/moss-c

max_issues_repo_name: JohnBSmith/moss-c

max_forks_repo_name: JohnBSmith/moss-c


text:
/**
 * @file batchg_zgemv.c
 *
 * Part of API test for Batched BLAS routines.
 *
 * @author  Samuel  D. Relton
 * @author  Pedro   V. Lara
 * @author  Mawussi Zounon
 * @date 2016-06-01
 *
 * @precisions normal z -> c d s
 *
 **/

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

#define COMPLEX
void batch_zgemv(
	const enum BBLAS_TRANS *trans,
	const int *m, const int *n,
	const BBLAS_Complex64_t *alpha,
	const BBLAS_Complex64_t **arrayA, const int *lda,
	const BBLAS_Complex64_t **arrayx, const int *incx,
	const BBLAS_Complex64_t *beta,
	BBLAS_Complex64_t **arrayy, const int *incy,
	const int batch_count, const enum BBLAS_OPTS batch_opts,
	int* info)
{
	/* Local variables */
	char func_name[15] = "batch_zgemv";

	/* Check input arguments */
	if (batch_count < 0)
    {
		xerbla_batch(func_name, BBLAS_ERR_BATCH_COUNT, -1);
    }

	if (batch_opts == BBLAS_FIXED)
	{
		/* Call fixed size code */
		batchf_zgemv(trans[0],
					 m[0], n[0],
					 alpha[0],
					 arrayA, lda[0],
					 arrayx, incx[0],
					 beta[0],
					 arrayy, incy[0],
					 batch_count, info[0]);
	}
	else if (batch_opts == BBLAS_VARIABLE)
	  {
		/* Call variable size code */
		batchv_zgemv(trans,
					 m, n,
					 alpha,
					 arrayA, lda,
					 arrayx, incx,
					 beta,
					 arrayy, incy,
					 batch_count, info);
	}
	else
	{
		xerbla_batch(func_name, BBLAS_ERR_BATCH_OPTS, -1);
	}
}

ext: c

max_stars_repo_name: mawussi/BBLAS-group


text:
#include <stdio.h>
#include "sqlite3.h"
#include <stdlib.h>
#include <time.h>
#include <sys/time.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

// returned row
static int callback(void *NotUsed, int argc, char **argv, char **azColName){
	int i;
	for(i=0; i<argc; i++){
		printf("%s = %s\t", azColName[i], argv[i] ? argv[i] : "NULL");
	}
	printf("\n");
	return 0;
}

// create test_int table in database
void create_int_tables(sqlite3 *db) {
	int r;
	char *err = 0;

	r = sqlite3_exec(db, "CREATE TABLE IF NOT EXISTS test_int (id INTEGER PRIMARY KEY, uniform INTEGER, normal5 INTEGER, normal20 INTEGER, t20 INTEGER)",
		callback, 0, &err);

	if(r != SQLITE_OK) {
		fprintf(stderr, "SQL error: %s\n", err);
		sqlite3_free(err);
	}
}

// create test table in database
void create_tables(sqlite3 *db) {
	int r;
	char *err = 0;

	r = sqlite3_exec(db, "CREATE TABLE IF NOT EXISTS test (id INTEGER PRIMARY KEY, uniformi INTEGER, normali5 INTEGER, normali20 INTEGER, uniformf FLOAT, normalf5 FLOAT, normalf20 FLOAT)",
		callback, 0, &err);

	if(r != SQLITE_OK) {
		fprintf(stderr, "SQL error: %s\n", err);
		sqlite3_free(err);
	}
}

// insert rows into test with varying random distributions
void insert_data(sqlite3 *db, int rows, gsl_rng *ran) {
	char* err = 0;
	int i;
	int r;

	r = sqlite3_exec(db, "PRAGMA synchronous = 0", NULL, 0, &err);

	if(r != SQLITE_OK) {
		fprintf(stderr, "SQL error: %s\n", err);
		sqlite3_free(err);
	}

	fprintf(stderr, "executed pragma\n");

	r = sqlite3_exec(db, "BEGIN TRANSACTION", NULL, 0, &err);

	if(r != SQLITE_OK) {
		fprintf(stderr, "SQL error: %s\n", err);
		sqlite3_free(err);
	}

	fprintf(stderr, "started transaction\n");

	sqlite3_stmt *stmt;
	sqlite3_prepare_v2(db,
		"INSERT INTO test (id, uniformi, normali5, normali20, uniformf, normalf5, normalf20) VALUES (?,?,?,?,?,?,?)",
		512,
		&stmt,
		0);

	for(i = 0; i < rows; i++) {
		sqlite3_bind_int(stmt, 1, i);
		sqlite3_bind_int(stmt, 2, (int)gsl_ran_flat(ran, -100, 100));
		sqlite3_bind_int(stmt, 3, (int)gsl_ran_gaussian(ran, 5));
		sqlite3_bind_int(stmt, 4, (int)gsl_ran_gaussian(ran, 20));
		sqlite3_bind_int(stmt, 5, (float)gsl_ran_flat(ran, -100, 100));
		sqlite3_bind_int(stmt, 6, (float)gsl_ran_gaussian(ran, 5));
		sqlite3_bind_int(stmt, 7, (float)gsl_ran_gaussian(ran, 20));

		int r = sqlite3_step(stmt);

		if(r != SQLITE_DONE && r != SQLITE_OK) {
			fprintf(stderr, "SQL error: %s\n", sqlite3_errmsg(db));
			sqlite3_free(err);
			return;
		}

		sqlite3_reset(stmt);

		if(i % 10000 == 0)
			printf("%i\n", i);
	}

	sqlite3_finalize(stmt);

	r=sqlite3_exec(db, "CREATE TEMPORARY TABLE tmp AS SELECT * FROM test", 0, 0, &err);

	if(r != SQLITE_OK) {
		fprintf(stderr, "SQL error: %s\n", err);
		sqlite3_free(err);
	}
        r=sqlite3_exec(db, "UPDATE tmp SET id = NULL;", 0, 0, &err);

	if(r != SQLITE_OK) {
		fprintf(stderr, "SQL error: %s\n", err);
		sqlite3_free(err);
	}
        for(int i=0;i<24;i++){
        r=sqlite3_exec(db, "INSERT INTO test SELECT * FROM tmp;", 0, 0, &err);

	if(r != SQLITE_OK) {
		fprintf(stderr, "SQL error: %s\n", err);
		sqlite3_free(err);
	}
}
        r=sqlite3_exec(db, "DROP TABLE tmp;", 0, 0, &err);

	if(r != SQLITE_OK) {
		fprintf(stderr, "SQL error: %s\n", err);
		sqlite3_free(err);
	}
	r=sqlite3_exec(db, "COMMIT", 0, 0, &err);

	if(r != SQLITE_OK) {
		fprintf(stderr, "SQL error: %s\n", err);
		sqlite3_free(err);
	}
}

int main(int argc, char **argv){
	sqlite3 *db;
	char *zErrMsg = 0;
	int rc;

	if(argc != 2) {
		fprintf(stderr, "%s <database name>\n", argv[0]);
		exit(1);
	}

	const gsl_rng_type *type;
	gsl_rng *ran;

	type = gsl_rng_ranlxd2;
	ran = gsl_rng_alloc(type);
	gsl_rng_set(ran, time(0));

	rc = sqlite3_open(argv[1], &db);
	if( rc ){
		fprintf(stderr, "Can't open database: %s\n", sqlite3_errmsg(db));
		sqlite3_close(db);
		exit(1);
	}

	create_tables(db);

	fprintf(stderr, "tables created\n");

	insert_data(db, 100000, ran);

	if(rc != SQLITE_OK) {
		fprintf(stderr, "SQL error: %s\n", zErrMsg);
		sqlite3_free(zErrMsg);
	}

	sqlite3_close(db);
	gsl_rng_free(ran);

	return 0;
}

max_stars_repo_name: mehrdadnikbakhtjam/sphyraena-extension

max_issues_repo_name: mehrdadnikbakhtjam/sphyraena-extension

max_forks_repo_name: mehrdadnikbakhtjam/sphyraena-extension


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

#pragma once

#include <opencv2\core\core.hpp>
#include <opencv2\features2d\features2d.hpp>

#include <utility>
#include <gsl\gsl>
#include <vector>

namespace mage
{
    class CameraCalibration;

    struct Projection
    {
        cv::Point2f Point;
        float Distance;
    };

    void ProjectPoints(const gsl::span<const cv::Point3f>& points3D, const cv::Matx34f& cameraPose, const cv::Matx33f& calibrationMatrix, std::vector<Projection>& points2D);

    Projection ProjectUndistorted(const cv::Matx34f& viewMatrix, const cv::Matx33f& calibration, const cv::Point3f& world);
    Projection ProjectDistorted(const CameraCalibration& cameraCalibration, const cv::Matx34f& viewMatrix, const cv::Point3f& world);

    std::vector<cv::Point2f> ProjectDistorted(
        const mage::CameraCalibration& calibration,
        const cv::Matx34f& viewMatrix,
        gsl::span<const cv::Vec3f> points3D);
}

ext: h

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


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

#include <gsl/gsl_math.h>
#include <gsl/gsl_filter.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>

int
main(void)
{
  const size_t N = 1000;                       /* length of time series */
  const size_t K = 61;                         /* window size */
  const double alpha = 3.0;                    /* Gaussian kernel has +/- 3 standard deviations */
  gsl_vector *x = gsl_vector_alloc(N);         /* input vector */
  gsl_vector *y = gsl_vector_alloc(N);         /* filtered output vector */
  gsl_vector *dy = gsl_vector_alloc(N);        /* first derivative filtered vector */
  gsl_vector *d2y = gsl_vector_alloc(N);       /* second derivative filtered vector */
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
  gsl_filter_gaussian_workspace *gauss_p = gsl_filter_gaussian_alloc(K);
  size_t i;

  /* generate input signal */
  for (i = 0; i < N; ++i)
    {
      double xi = (i > N / 2) ? 0.5 : 0.0;
      double ei = gsl_ran_gaussian(r, 0.1);

      gsl_vector_set(x, i, xi + ei);
    }

  /* apply filters */
  gsl_filter_gaussian(GSL_FILTER_END_PADVALUE, alpha, 0, x, y, gauss_p);
  gsl_filter_gaussian(GSL_FILTER_END_PADVALUE, alpha, 1, x, dy, gauss_p);
  gsl_filter_gaussian(GSL_FILTER_END_PADVALUE, alpha, 2, x, d2y, gauss_p);

  /* print results */
  for (i = 0; i < N; ++i)
    {
      double xi = gsl_vector_get(x, i);
      double yi = gsl_vector_get(y, i);
      double dyi = gsl_vector_get(dy, i);
      double d2yi = gsl_vector_get(d2y, i);
      double dxi;

      /* compute finite difference of x vector */
      if (i == 0)
        dxi = gsl_vector_get(x, i + 1) - xi;
      else if (i == N - 1)
        dxi = gsl_vector_get(x, i) - gsl_vector_get(x, i - 1);
      else
        dxi = 0.5 * (gsl_vector_get(x, i + 1) - gsl_vector_get(x, i - 1));

      printf("%.12e %.12e %.12e %.12e %.12e\n",
             xi,
             yi,
             dxi,
             dyi,
             d2yi);
    }

  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(dy);
  gsl_vector_free(d2y);
  gsl_rng_free(r);
  gsl_filter_gaussian_free(gauss_p);

  return 0;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/* -----------------------------------------------------------------------------
 * Copyright 2021 Jonathan Haigh
 * SPDX-License-Identifier: MIT
 * ---------------------------------------------------------------------------*/

#ifndef SQ_INCLUDE_GUARD_parser_TokenView_h_
#define SQ_INCLUDE_GUARD_parser_TokenView_h_

#include "core/Token.h"
#include "core/typeutil.h"

#include <gsl/gsl>
#include <optional>
#include <range/v3/view/facade.hpp>
#include <string_view>

namespace sq::parser {

/**
 * A ranges::view of the input query split into tokens.
 */
class TokenView : public ranges::view_facade<TokenView> {
public:
  explicit TokenView(std::string_view str) noexcept;

  TokenView() noexcept = default;
  TokenView(const TokenView &) noexcept = default;
  TokenView(TokenView &&) noexcept = default;
  TokenView &operator=(const TokenView &) noexcept = default;
  TokenView &operator=(TokenView &&) noexcept = default;
  ~TokenView() noexcept = default;

  // required for ranges::view_facade
  friend ranges::range_access;

  SQ_ND const Token &read() const;
  SQ_ND bool equal(ranges::default_sentinel_t other) const noexcept;
  void next();

private:
  gsl::index whitespace_length() const;

  std::string_view str_;
  gsl::index pos_ = 0;
  mutable std::optional<Token> cache_ = std::nullopt;
};

} // namespace sq::parser
#endif // SQ_INCLUDE_GUARD_parser_TokenView_h_

max_stars_repo_name: jonathanhaigh/sq

max_issues_repo_name: jonathanhaigh/sq

max_forks_repo_name: jonathanhaigh/sq


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

#ifndef __GSL_RNG_H__
#define __GSL_RNG_H__
#include <stdlib.h>
#include <gsl/gsl_errno.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
  {
    const char *name;
    unsigned long int max;
    unsigned long int min;
    size_t size;
    void (*set) (void *state, unsigned long int seed);
    unsigned long int (*get) (void *state);
    double (*get_double) (void *state);
  }
gsl_rng_type;

typedef struct
  {
    const gsl_rng_type * type;
    void *state;
  }
gsl_rng;

/* These structs also need to appear in default.c so you can select
   them via the environment variable GSL_RNG_TYPE */

extern const gsl_rng_type *gsl_rng_cmrg;
extern const gsl_rng_type *gsl_rng_gfsr4;
extern const gsl_rng_type *gsl_rng_minstd;
extern const gsl_rng_type *gsl_rng_mrg;
extern const gsl_rng_type *gsl_rng_mt19937;
extern const gsl_rng_type *gsl_rng_mt19937_1998;
extern const gsl_rng_type *gsl_rng_r250;
extern const gsl_rng_type *gsl_rng_ran0;
extern const gsl_rng_type *gsl_rng_ran1;
extern const gsl_rng_type *gsl_rng_ran2;
extern const gsl_rng_type *gsl_rng_ran3;
extern const gsl_rng_type *gsl_rng_rand48;
extern const gsl_rng_type *gsl_rng_rand;
extern const gsl_rng_type *gsl_rng_random128_bsd;
extern const gsl_rng_type *gsl_rng_random128_glibc2;
extern const gsl_rng_type *gsl_rng_random128_libc5;
extern const gsl_rng_type *gsl_rng_random256_bsd;
extern const gsl_rng_type *gsl_rng_random256_glibc2;
extern const gsl_rng_type *gsl_rng_random256_libc5;
extern const gsl_rng_type *gsl_rng_random32_bsd;
extern const gsl_rng_type *gsl_rng_random32_glibc2;
extern const gsl_rng_type *gsl_rng_random32_libc5;
extern const gsl_rng_type *gsl_rng_random64_bsd;
extern const gsl_rng_type *gsl_rng_random64_glibc2;
extern const gsl_rng_type *gsl_rng_random64_libc5;
extern const gsl_rng_type *gsl_rng_random8_bsd;
extern const gsl_rng_type *gsl_rng_random8_glibc2;
extern const gsl_rng_type *gsl_rng_random8_libc5;
extern const gsl_rng_type *gsl_rng_random_bsd;
extern const gsl_rng_type *gsl_rng_random_glibc2;
extern const gsl_rng_type *gsl_rng_random_libc5;
extern const gsl_rng_type *gsl_rng_randu;
extern const gsl_rng_type *gsl_rng_ranf;
extern const gsl_rng_type *gsl_rng_ranlux389;
extern const gsl_rng_type *gsl_rng_ranlux;
extern const gsl_rng_type *gsl_rng_ranlxd1;
extern const gsl_rng_type *gsl_rng_ranlxd2;
extern const gsl_rng_type *gsl_rng_ranlxs0;
extern const gsl_rng_type *gsl_rng_ranlxs1;
extern const gsl_rng_type *gsl_rng_ranlxs2;
extern const gsl_rng_type *gsl_rng_ranmar;
extern const gsl_rng_type *gsl_rng_slatec;
extern const gsl_rng_type *gsl_rng_taus;
extern const gsl_rng_type *gsl_rng_transputer;
extern const gsl_rng_type *gsl_rng_tt800;
extern const gsl_rng_type *gsl_rng_uni32;
extern const gsl_rng_type *gsl_rng_uni;
extern const gsl_rng_type *gsl_rng_vax;
extern const gsl_rng_type *gsl_rng_zuf;

const gsl_rng_type ** gsl_rng_types_setup(void);

extern const gsl_rng_type *gsl_rng_default;
extern unsigned long int gsl_rng_default_seed;

gsl_rng *gsl_rng_alloc (const gsl_rng_type * T);
int gsl_rng_memcpy (gsl_rng * dest, const gsl_rng * src);
gsl_rng *gsl_rng_clone (const gsl_rng * r);

void gsl_rng_free (gsl_rng * r);

void gsl_rng_set (const gsl_rng * r, unsigned long int seed);
unsigned long int gsl_rng_max (const gsl_rng * r);
unsigned long int gsl_rng_min (const gsl_rng * r);
const char *gsl_rng_name (const gsl_rng * r);
size_t gsl_rng_size (const gsl_rng * r);
void * gsl_rng_state (const gsl_rng * r);

void gsl_rng_print_state (const gsl_rng * r);

const gsl_rng_type * gsl_rng_env_setup (void);

unsigned long int gsl_rng_get (const gsl_rng * r);
double gsl_rng_uniform (const gsl_rng * r);
double gsl_rng_uniform_pos (const gsl_rng * r);
unsigned long int gsl_rng_uniform_int (const gsl_rng * r, unsigned long int n);

#ifdef HAVE_INLINE
extern inline unsigned long int gsl_rng_get (const gsl_rng * r);

extern inline unsigned long int
gsl_rng_get (const gsl_rng * r)
{
  return (r->type->get) (r->state);
}

extern inline double gsl_rng_uniform (const gsl_rng * r);

extern inline double
gsl_rng_uniform (const gsl_rng * r)
{
  return (r->type->get_double) (r->state);
}

extern inline double gsl_rng_uniform_pos (const gsl_rng * r);

extern inline double
gsl_rng_uniform_pos (const gsl_rng * r)
{
  double x ;
  do
    {
      x = (r->type->get_double) (r->state) ;
    }
  while (x == 0) ;

  return x ;
}

extern inline unsigned long int gsl_rng_uniform_int (const gsl_rng * r, unsigned long int n);

extern inline unsigned long int
gsl_rng_uniform_int (const gsl_rng * r, unsigned long int n)
{
  unsigned long int offset = r->type->min;
  unsigned long int range = r->type->max - offset;
  unsigned long int scale = range / n;
  unsigned long int k;

  if (n > range)
    {
      GSL_ERROR_VAL ("n exceeds maximum value of generator",
			GSL_EINVAL, 0) ;
    }

  do
    {
      k = (((r->type->get) (r->state)) - offset) / scale;
    }
  while (k >= n);

  return k;
}
#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_RNG_H__ */


text:
/**
 * File: specmodel_utils.h
 * Header file for specmodel_utils.c
 *
 */
#ifndef _SPECMODEL_UTILS_H
#define _SPECMODEL_UTILS_H
#include <gsl/gsl_interp.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include "fitsio.h"
#include "aXe_grism.h"

// interpolation type used for the
// spectral models
#define SMODEL_INTERP_TYPE gsl_interp_linear

/*
 * Struct: energy_distrib
 */
typedef struct
{
  int npoints;
  double *wavelength;
  double *flux;

  gsl_interp *interp;
  gsl_interp_accel *accel;

}
energy_distrib;

/*
 * Struct: spectral_models
 */
typedef struct
{
  int n_models;

  energy_distrib **SEDlist;
}
spectral_models;

typedef struct
{
  int dim_x;               // the number of pixels in x-direction
  int dim_y;               // the number of pixels in y-direction
  double xmean;           // mean x-coordinate of the object
  double ymean;           // mean y-coordinate of the object

  gsl_matrix *modimage;   // matrix with the direct image
}
  dirim_emission;

typedef struct
{
  int n_models;              // the number of direct image object

  dirim_emission **obj_list; // the list of direct image objects
}
object_models;

extern object_models *
load_object_models(const char object_models_file[]);

extern dirim_emission **
load_diremission_list(const char object_models_file[], const int n_models);

extern dirim_emission *
load_diremission(const char object_models_file[], const int n_extension);

extern void
free_object_models(object_models *objmodels);

extern void
free_dirim_emission(dirim_emission *diremission);

extern void
print_object_models(const object_models *objmodels);

extern void
print_dirim_emission(const dirim_emission *diremission);

extern dirim_emission *
get_dirim_emission(const object_models *objmodels, const int objspec);

extern int
has_aperture_dirim(const object_models *objmodels, const object *actobject);

extern double
get_diremission_value(const dirim_emission *diremission,
		      const double xpos, const double ypos);

extern double
bilin_interp_matrix(const gsl_matrix *modimage, const double x, const double y);

extern spectral_models *
load_spectral_models(const char spectral_models_file[]);

extern void
print_spectral_models(const spectral_models *smodels);

extern energy_distrib **
load_SED_list(const char spectral_models_file[], const int n_models);

extern energy_distrib *
load_SED_from_fitsext(const char spectral_models_file[], fitsfile *s_models);

extern int
get_num_extensions(const char spectral_models_file[]);

extern void
free_spectral_models(spectral_models *smodels);

extern void
free_specmodels_SEDs(spectral_models *smodels);

extern energy_distrib *
get_model_sed(const spectral_models *spec_mod, const int modspec);

#endif

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


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

/* Author: G. Jungman */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include "tridiag.h"
#include "gsl_linalg.h"

/* for description of method see [Engeln-Mullges + Uhlig, p. 92]
 *
 *     diag[0]  offdiag[0]             0   .....
 *  offdiag[0]     diag[1]    offdiag[1]   .....
 *           0  offdiag[1]       diag[2]
 *           0           0    offdiag[2]   .....
 */
static
int
solve_tridiag(
  const double diag[], size_t d_stride,
  const double offdiag[], size_t o_stride,
  const double b[], size_t b_stride,
  double x[], size_t x_stride,
  size_t N)
{
  int status;
  double *gamma = (double *) malloc (N * sizeof (double));
  double *alpha = (double *) malloc (N * sizeof (double));
  double *c = (double *) malloc (N * sizeof (double));
  double *z = (double *) malloc (N * sizeof (double));

  if (gamma == 0 || alpha == 0 || c == 0 || z == 0)
    {
      status = GSL_ENOMEM;
    }
  else
    {
      size_t i, j;

      /* Cholesky decomposition
         A = L.D.L^t
         lower_diag(L) = gamma
         diag(D) = alpha
       */
      alpha[0] = diag[0];
      gamma[0] = offdiag[0] / alpha[0];

      for (i = 1; i < N - 1; i++)
	{
	  alpha[i] = diag[d_stride * i] - offdiag[o_stride*(i - 1)] * gamma[i - 1];
	  gamma[i] = offdiag[o_stride * i] / alpha[i];
	}

      if (N > 1)
        {
          alpha[N - 1] = diag[d_stride * (N - 1)] - offdiag[o_stride*(N - 2)] * gamma[N - 2];
        }

      /* update RHS */
      z[0] = b[0];
      for (i = 1; i < N; i++)
	{
	  z[i] = b[b_stride * i] - gamma[i - 1] * z[i - 1];
	}
      for (i = 0; i < N; i++)
	{
	  c[i] = z[i] / alpha[i];
	}

      /* backsubstitution */
      x[x_stride * (N - 1)] = c[N - 1];
      if (N >= 2)
	{
	  for (i = N - 2, j = 0; j <= N - 2; j++, i--)
	    {
	      x[x_stride * i] = c[i] - gamma[i] * x[x_stride * (i + 1)];
	    }
	}

      status = GSL_SUCCESS;
    }

  if (z != 0)
    free (z);
  if (c != 0)
    free (c);
  if (alpha != 0)
    free (alpha);
  if (gamma != 0)
    free (gamma);

  return status;
}

/* for description of method see [Engeln-Mullges + Uhlig, p. 96]
 *
 *      diag[0]  offdiag[0]             0   .....  offdiag[N-1]
 *   offdiag[0]     diag[1]    offdiag[1]   .....
 *            0  offdiag[1]       diag[2]
 *            0           0    offdiag[2]   .....
 *          ...         ...
 * offdiag[N-1]         ...
 *
 */
static
int
solve_cyc_tridiag(
  const double diag[], size_t d_stride,
  const double offdiag[], size_t o_stride,
  const double b[], size_t b_stride,
  double x[], size_t x_stride,
  size_t N)
{
  int status;
  double * delta = (double *) malloc (N * sizeof (double));
  double * gamma = (double *) malloc (N * sizeof (double));
  double * alpha = (double *) malloc (N * sizeof (double));
  double * c = (double *) malloc (N * sizeof (double));
  double * z = (double *) malloc (N * sizeof (double));

  if (delta == 0 || gamma == 0 || alpha == 0 || c == 0 || z == 0)
    {
      status = GSL_ENOMEM;
    }
  else
    {
      size_t i, j;
      double sum = 0.0;

      /* factor */

      alpha[0] = diag[0];
      gamma[0] = offdiag[0] / alpha[0];
      delta[0] = offdiag[o_stride * (N-1)] / alpha[0];

      for (i = 1; i < N - 2; i++)
	{
	  alpha[i] = diag[d_stride * i] - offdiag[o_stride * (i-1)] * gamma[i - 1];
	  gamma[i] = offdiag[o_stride * i] / alpha[i];
	  delta[i] = -delta[i - 1] * offdiag[o_stride * (i-1)] / alpha[i];
	}
      for (i = 0; i < N - 2; i++)
	{
	  sum += alpha[i] * delta[i] * delta[i];
	}
      alpha[N - 2] = diag[d_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * gamma[N - 3];
      gamma[N - 2] = (offdiag[o_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * delta[N - 3]) / alpha[N - 2];
      alpha[N - 1] = diag[d_stride * (N - 1)] - sum - offdiag[o_stride * (N - 2)] * gamma[N - 2] * gamma[N - 2];

      /* update */
      z[0] = b[0];
      for (i = 1; i < N - 1; i++)
	{
	  z[i] = b[b_stride * i] - z[i - 1] * gamma[i - 1];
	}
      sum = 0.0;
      for (i = 0; i < N - 2; i++)
	{
	  sum += delta[i] * z[i];
	}
      z[N - 1] = b[b_stride * (N - 1)] - sum - gamma[N - 2] * z[N - 2];
      for (i = 0; i < N; i++)
	{
	  c[i] = z[i] / alpha[i];
	}

      /* backsubstitution */
      x[x_stride * (N - 1)] = c[N - 1];
      x[x_stride * (N - 2)] = c[N - 2] - gamma[N - 2] * x[x_stride * (N - 1)];
      if (N >= 3)
	{
	  for (i = N - 3, j = 0; j <= N - 3; j++, i--)
	    {
	      x[x_stride * i] = c[i] - gamma[i] * x[x_stride * (i + 1)] - delta[i] * x[x_stride * (N - 1)];
	    }
	}

      status = GSL_SUCCESS;
    }

  if (z != 0)
    free (z);
  if (c != 0)
    free (c);
  if (alpha != 0)
    free (alpha);
  if (gamma != 0)
    free (gamma);
  if (delta != 0)
    free (delta);

  return status;
}

int
gsl_linalg_solve_symm_tridiag(
  const gsl_vector * diag,
  const gsl_vector * offdiag,
  const gsl_vector * rhs,
  gsl_vector * solution)
{
  if(diag->size != rhs->size ||
          (offdiag->size != rhs->size && offdiag->size != rhs->size-1) ||
	  (solution->size != rhs->size)
          ) {
    return GSL_EBADLEN;
  }
  else {
    return solve_tridiag(diag->data, diag->stride,
                         offdiag->data, offdiag->stride,
			 rhs->data, rhs->stride,
	                 solution->data, solution->stride,
	                 diag->size);
  }
}

int
gsl_linalg_solve_symm_cyc_tridiag(
  const gsl_vector * diag,
  const gsl_vector * offdiag,
  const gsl_vector * rhs,
  gsl_vector * solution)
{
  if(diag->size != rhs->size ||
          offdiag->size != rhs->size ||
          solution->size != rhs->size
          ) {
    return GSL_EBADLEN;
  }
  else {
    return solve_cyc_tridiag(diag->data, diag->stride,
                             offdiag->data, offdiag->stride,
			     rhs->data, rhs->stride,
	                     solution->data, solution->stride,
	                     diag->size);
  }
}


text:
#include <gsl/gsl_spline.h>

typedef struct FastPMFDInterp{
    size_t size;
    gsl_interp * F;
    gsl_interp * DF;
    gsl_interp * DDF;
    gsl_interp_accel * acc;
} FastPMFDInterp;

void fastpm_fd_interp_init(FastPMFDInterp * FDinterp);

double fastpm_do_fd_interp(FastPMFDInterp * FDinterp, int F_id, double y);

void fastpm_fd_interp_destroy(FastPMFDInterp * FDinterp);

max_stars_repo_name: sbird/FastPMRunner

max_issues_repo_name: sbird/FastPMRunner

max_forks_repo_name: sbird/FastPMRunner


text:
#ifndef _MASSINTEGRAND__H
#define _MASSINTEGRAND__H

#include <vector>
#include <algorithm>
#include <numeric>
#include <functional>
#include <cstdint>

#ifdef HAS_MKL
#include <mkl.h>
#else
#include <cblas.h>
#endif

#ifdef __SSE4_2__
#include <emmintrin.h>
#include <mm_malloc.h>
#endif

#include "SparseAssemblyNative.h"
#include "_det_inv_.h"
#include "_matmul_.h"

#ifndef LL_TYPES
#define LL_TYPES
using Real = double;
using Integer = std::int64_t;
using UInteger = std::uint64_t;
#endif

/*---------------------------------------------------------------------------------------------*/
#ifndef CUSTOM_ALLOCATION_
#define CUSTOM_ALLOCATION_
template<typename T>
FASTOR_INLINE T *allocate(Integer size) {
#if defined(__AVX__)
    T *out = (T*)_mm_malloc(sizeof(T)*size,32);
#elif defined(__SSE__)
    T *out = (T*)_mm_malloc(sizeof(T)*size,16);
#else
    T *out = (T*)malloc(sizeof(T)*size);
#endif
    return out;
}

template<typename T>
FASTOR_INLINE void deallocate(T *a) {
#if defined(__SSE__)
    _mm_free(a);
#else
    free(a);
#endif
}
#endif
/*---------------------------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------------------------*/
#ifndef SPARSE_TRIPLET_FILLER
#define SPARSE_TRIPLET_FILLER
// IJV Filler
FASTOR_INLINE
void fill_triplet(  const Integer *i,
                    const Integer *j,
                    const Real *coeff,
                    int *I,
                    int *J,
                    Real *V,
                    Integer elem,
                    Integer nvar,
                    Integer nodeperelem,
                    const UInteger *elements,
                    Integer i_shape,
                    Integer j_shape
                    ) {

    Integer ndof = nvar*nodeperelem;
    Integer *current_row_column = allocate<Integer>(nvar*nodeperelem);

    Integer const_elem_retriever;
    for (Integer counter=0; counter<nodeperelem; ++counter) {
        const_elem_retriever = nvar*elements[elem*nodeperelem+counter];
        for (Integer ncounter=0; ncounter<nvar; ++ncounter) {
            current_row_column[nvar*counter+ncounter] = const_elem_retriever+ncounter;
        }
    }

    Integer icounter = 0;
    Integer ncounter = ndof*ndof*elem;

    Integer const_I_retriever;
    for (Integer counter=0; counter<ndof; ++counter) {
        const_I_retriever = current_row_column[counter];
        for (Integer iterator=0; iterator<ndof; ++iterator) {
            I[ncounter] = const_I_retriever;
            J[ncounter] = current_row_column[iterator];
            V[ncounter] = coeff[icounter];
            ncounter++;
            icounter++;
        }
    }

    deallocate(current_row_column);
}

FASTOR_INLINE
void fill_global_data(
                    const Integer *i,
                    const Integer *j,
                    const Real *coeff,
                    int *I,
                    int *J,
                    Real *V,
                    Integer elem,
                    Integer nvar,
                    Integer nodeperelem,
                    const UInteger *elements,
                    Integer i_shape,
                    Integer j_shape,
                    int recompute_sparsity_pattern,
                    int squeeze_sparsity_pattern,
                    const int *data_local_indices,
                    const int *data_global_indices,
                    const UInteger *sorted_elements,
                    const Integer *sorter
    )
{

    if (recompute_sparsity_pattern) {
        fill_triplet(i,j,coeff,I,J,V,elem,nvar,nodeperelem,elements,i_shape,j_shape);
    }
    else {
        if (squeeze_sparsity_pattern) {
            SparseAssemblyNativeCSR_RecomputeDataIndex_(
                coeff,
                J,
                I,
                V,
                elem,
                nvar,
                nodeperelem,
                sorted_elements,
                sorter);
        }
        else {
            const int ndof = nvar*nodeperelem;
            const int local_capacity = ndof*ndof;
            SparseAssemblyNativeCSR_(
                coeff,
                data_local_indices,
                data_global_indices,
                elem,
                local_capacity,
                V
                );
        }
    }
}

#endif
/*---------------------------------------------------------------------------------------------*/

inline void _MassIntegrand_Filler_(Real *mass,
    const Real* bases,
    const Real* detJ,
    int ngauss,
    int noderpelem,
    int ndim,
    int nvar,
    Real rho) {

    int local_size = nvar*noderpelem;

    Real *N     = allocate<Real>(nvar*local_size);
    Real *rhoNN = allocate<Real>(local_size*local_size);

    std::fill(N,N+nvar*local_size,0.);

    for (int igauss = 0; igauss < ngauss; ++igauss) {

        // Fill mass integrand
        for (int j=0; j<noderpelem; ++j) {
            const Real bases_j = bases[j*ngauss+igauss];
            for (int ivar=0; ivar<ndim; ++ivar) {
                N[j*nvar*nvar+ivar*nvar+ivar] = bases_j;
            }
        }

        cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
            local_size, local_size, nvar, rho, N, nvar, N, nvar, 0.0, rhoNN, local_size);

        // Multiply mass with detJ
        const Real detJ_igauss = detJ[igauss];
        for (int i=0; i<local_size*local_size; ++i) {
            mass[i] += rhoNN[i]*detJ_igauss;
        }

    }

    deallocate(N);
    deallocate(rhoNN);
}

inline void _ConstantMassIntegrand_Filler_(Real *mass,
    const Real* constant_mass_integrand,
    const Real* detJ,
    int ngauss,
    int local_capacity) {

    for (int igauss = 0; igauss < ngauss; ++igauss) {
        // Multiply mass with detJ
        const Real detJ_igauss = detJ[igauss];
        for (int i=0; i<local_capacity; ++i) {
            mass[i] += constant_mass_integrand[igauss*local_capacity+i]*detJ_igauss;
        }
    }
}

inline void _GenericConstantMassIntegrand_(
    const UInteger* elements,
    const Real* points,
    const Real* Jm,
    const Real* AllGauss,
    const Real* constant_mass_integrand,
    Integer nelem,
    Integer ndim,
    Integer nvar,
    Integer ngauss,
    Integer nodeperelem,
    Integer local_capacity,
    Integer mass_type,
    const Integer* local_rows_mass,
    const Integer* local_cols_mass,
    int *I_mass,
    int *J_mass,
    Real *V_mass,
    Real *mass,
    int recompute_sparsity_pattern,
    int squeeze_sparsity_pattern,
    const int *data_local_indices,
    const int *data_global_indices,
    const UInteger *sorted_elements,
    const Integer *sorter
    ) {

    Integer ndof                    = nodeperelem*nvar;
    Real *LagrangeElemCoords        = allocate<Real>(nodeperelem*ndim);
    Real *MaterialGradient          = allocate<Real>(ndim*nodeperelem);
    // +1 TO AVOID OVERSTEPPONG IN MATMUL UNALIGNED STORE - EXTREMELY DANGEROUSE BUG
    int plus1                       = ndim == 2 ? 0 : 1;
    Real *ParentGradientX           = allocate<Real>(ndim*ndim+plus1);
    Real detJ                       = 0.;
    Real *massel                    = allocate<Real>(local_capacity);
    Real *massel_lumped             = allocate<Real>(ndof);

    // PRE-COMPUTE ISOPARAMETRIC GRADIENTS
    std::vector<std::vector<Real>> current_Jms(ngauss);
    for (int g=0; g<ngauss; ++g) {
        std::vector<Real> current_Jm(ndim*nodeperelem);
        for (int j=0; j<nodeperelem; ++j) {
            for (int k=0; k<ndim; ++k) {
                current_Jm[k*nodeperelem+j] = Jm[k*ngauss*nodeperelem+j*ngauss+g];
            }
        }
        current_Jms[g] = current_Jm;
    }

    // LOOP OVER ELEMETNS
    for (Integer elem=0; elem < nelem; ++elem) {

        // GET THE FIELDS AT THE ELEMENT LEVEL
        for (Integer i=0; i<nodeperelem; ++i) {
            const Integer inode = elements[elem*nodeperelem+i];
            for (Integer j=0; j<ndim; ++j) {
                LagrangeElemCoords[i*ndim+j] = points[inode*ndim+j];
            }
        }

        std::fill_n(massel,local_capacity,0.);

        for (Integer igauss = 0; igauss < ngauss; ++igauss) {

            {
                // USING A STL BASED FILLER REMOVES THE ANNOYING BUG
                std::fill_n(ParentGradientX,ndim*ndim+plus1,0.);
                _matmul_(ndim,ndim,nodeperelem,current_Jms[igauss].data(),LagrangeElemCoords,ParentGradientX);
                const Real detX = _det_(ndim, ParentGradientX);
                detJ = AllGauss[igauss]*std::abs(detX);
            }

            // Multiply constant part of mass with detJ
#ifdef __AVX__
            using V2 = Fastor::SIMDVector<Real>;
            V2 _va, _vb, _vout;
            _vb.set(detJ);
            int Vsize = V2::Size;
            int ROUND_ = ROUND_DOWN(local_capacity,Vsize);
            int i=0;
            for (; i<ROUND_; i+=Vsize) {
                _va.load(&constant_mass_integrand[igauss*local_capacity+i],false);
                _vout.load(&massel[i],false);
#ifdef __FMA__
                _vout = fmadd(_va,_vb,_vout);
#else
                _vout += _va*_vb;
#endif
                _vout.store(&massel[i],false);
            }
            for ( ; i<local_capacity; ++i) {
                massel[i] += constant_mass_integrand[igauss*local_capacity+i]*detJ;
            }
#else
            for (int i=0; i<local_capacity; ++i) {
                massel[i] += constant_mass_integrand[igauss*local_capacity+i]*detJ;
            }
#endif
        }

        // FOR LUMP MASS MATRIX
        if (mass_type == 0) {
            // LUMP MASS
            for (Integer i=0; i<ndof; ++i) {
                massel_lumped[i] = std::accumulate(&massel[i*ndof], &massel[(i+1)*ndof],0.);
            }

            // ASSEMBLE LUMPED MASS
            {
                for (Integer i = 0; i<nodeperelem; ++i) {
                    UInteger T_idx = elements[elem*nodeperelem+i]*nvar;
                    for (Integer iterator = 0; iterator < nvar; ++iterator) {
                        mass[T_idx+iterator] += massel_lumped[i*nvar+iterator];
                    }
                }
            }
        }
        else {
            fill_global_data(   local_rows_mass,
                                local_cols_mass,
                                massel,
                                I_mass,
                                J_mass,
                                V_mass,
                                elem,
                                nvar,
                                nodeperelem,
                                elements,
                                local_capacity,
                                local_capacity,
                                recompute_sparsity_pattern,
                                squeeze_sparsity_pattern,
                                data_local_indices,
                                data_global_indices,
                                sorted_elements,
                                sorter);
        }
    }

    deallocate(ParentGradientX);
    deallocate(MaterialGradient);
    deallocate(LagrangeElemCoords);
    deallocate(massel);
    deallocate(massel_lumped);
}

inline void _SymmetricConstantMassIntegrand_(
    const UInteger* elements,
    const Real* points,
    const Real* Jm,
    const Real* AllGauss,
    const Real* constant_mass_integrand,
    Integer nelem,
    Integer ndim,
    Integer nvar,
    Integer ngauss,
    Integer nodeperelem,
    Integer local_capacity,
    Integer mass_type,
    const Integer* local_rows_mass,
    const Integer* local_cols_mass,
    int *I_mass,
    int *J_mass,
    Real *V_mass,
    Real *mass,
    int recompute_sparsity_pattern,
    int squeeze_sparsity_pattern,
    const int *data_local_indices,
    const int *data_global_indices,
    const UInteger *sorted_elements,
    const Integer *sorter
    ) {

    const Integer ndof                      = nodeperelem*nvar;
    const Integer symmetric_local_capacity  = ndof*(ndof+1)/2;
    Real *LagrangeElemCoords                = allocate<Real>(nodeperelem*ndim);
    Real *MaterialGradient                  = allocate<Real>(ndim*nodeperelem);
    // +1 TO AVOID OVERSTEPPONG IN MATMUL UNALIGNED STORE - EXTREMELY DANGEROUSE BUG
    int plus1                               = ndim == 2 ? 0 : 1;
    Real *ParentGradientX                   = allocate<Real>(ndim*ndim+plus1);
    Real detJ                               = 0.;
    Real *massel                            = allocate<Real>(local_capacity);
    Real *symmetric_massel                  = allocate<Real>(symmetric_local_capacity);
    Real *massel_lumped                     = allocate<Real>(ndof);

    // PRE-COMPUTE SYMMETRIC PART OF CONSTANT MASS WITHOUT ZEROS
    Real *symmetric_constant_mass_integrand = allocate<Real>(ngauss*symmetric_local_capacity);
    Integer counter = 0;
    for (Integer igauss = 0; igauss < ngauss; ++igauss) {
        for (Integer i=0; i<ndof; ++i) {
            for (Integer j=i; j<ndof; ++j) {
                const Real value = constant_mass_integrand[igauss*ndof*ndof+i*ndof+j];
                symmetric_constant_mass_integrand[counter] = value;
                counter++;
            }
        }
    }

    // PRE-COMPUTE ISOPARAMETRIC GRADIENTS
    std::vector<std::vector<Real>> current_Jms(ngauss);
    for (int g=0; g<ngauss; ++g) {
        std::vector<Real> current_Jm(ndim*nodeperelem);
        for (int j=0; j<nodeperelem; ++j) {
            for (int k=0; k<ndim; ++k) {
                current_Jm[k*nodeperelem+j] = Jm[k*ngauss*nodeperelem+j*ngauss+g];
            }
        }
        current_Jms[g] = current_Jm;
    }

    // LOOP OVER ELEMETNS
    for (Integer elem=0; elem < nelem; ++elem) {

        // GET THE FIELDS AT THE ELEMENT LEVEL
        for (Integer i=0; i<nodeperelem; ++i) {
            const Integer inode = elements[elem*nodeperelem+i];
            for (Integer j=0; j<ndim; ++j) {
                LagrangeElemCoords[i*ndim+j] = points[inode*ndim+j];
            }
        }

        std::fill_n(massel,local_capacity,0.);
        std::fill_n(symmetric_massel,symmetric_local_capacity,0.);

        for (Integer igauss = 0; igauss < ngauss; ++igauss) {

            {
                // USING A STL BASED FILLER REMOVES THE ANNOYING BUG
                std::fill_n(ParentGradientX,ndim*ndim+plus1,0.);
                _matmul_(ndim,ndim,nodeperelem,current_Jms[igauss].data(),LagrangeElemCoords,ParentGradientX);
                const Real detX = _det_(ndim, ParentGradientX);
                detJ = AllGauss[igauss]*std::abs(detX);
            }

            // Multiply constant part of mass with detJ
#ifdef __AVX__
            using V2 = Fastor::SIMDVector<Real>;
            V2 _va, _vb, _vout;
            _vb.set(detJ);
            int Vsize = V2::Size;
            int ROUND_ = ROUND_DOWN(symmetric_local_capacity,Vsize);
            int i=0;
            for (; i<ROUND_; i+=Vsize) {
                _va.load(&symmetric_constant_mass_integrand[igauss*symmetric_local_capacity+i],false);
                _vout.load(&symmetric_massel[i],false);
#ifdef __FMA__
                _vout = fmadd(_va,_vb,_vout);
#else
                _vout += _va*_vb;
#endif
                _vout.store(&symmetric_massel[i],false);
            }
            for ( ; i<symmetric_local_capacity; ++i) {
                symmetric_massel[i] += symmetric_constant_mass_integrand[igauss*symmetric_local_capacity+i]*detJ;
            }
#else
            for (int i=0; i<symmetric_local_capacity; ++i) {
                symmetric_massel[i] += symmetric_constant_mass_integrand[igauss*symmetric_local_capacity+i]*detJ;
            }
#endif
        }

        // BUILD TOTAL LOCAL MASS MATRIX
        Integer counter = 0;
        for (Integer i=0; i<ndof; ++i) {
            for (Integer j=i; j<ndof; ++j) {
                const Real value = symmetric_massel[counter];
                massel[i*ndof+j] = value;
                massel[j*ndof+i] = value;
                counter++;
            }
        }

        // FOR LUMP MASS MATRIX
        if (mass_type == 0) {
            // LUMP MASS
            for (Integer i=0; i<ndof; ++i) {
                massel_lumped[i] = std::accumulate(&massel[i*ndof], &massel[(i+1)*ndof],0.);
            }

            // ASSEMBLE LUMPED MASS
            {
                for (Integer i = 0; i<nodeperelem; ++i) {
                    UInteger T_idx = elements[elem*nodeperelem+i]*nvar;
                    for (Integer iterator = 0; iterator < nvar; ++iterator) {
                        mass[T_idx+iterator] += massel_lumped[i*nvar+iterator];
                    }
                }
            }
        }
        else {
            fill_global_data(   local_rows_mass,
                                local_cols_mass,
                                massel,
                                I_mass,
                                J_mass,
                                V_mass,
                                elem,
                                nvar,
                                nodeperelem,
                                elements,
                                local_capacity,
                                local_capacity,
                                recompute_sparsity_pattern,
                                squeeze_sparsity_pattern,
                                data_local_indices,
                                data_global_indices,
                                sorted_elements,
                                sorter);
        }
    }

    deallocate(ParentGradientX);
    deallocate(MaterialGradient);
    deallocate(LagrangeElemCoords);
    deallocate(massel);
    deallocate(symmetric_constant_mass_integrand);
    deallocate(symmetric_massel);
    deallocate(massel_lumped);
}

inline void _SymmetricNonZeroConstantMassIntegrand_(
    const UInteger* elements,
    const Real* points,
    const Real* Jm,
    const Real* AllGauss,
    const Real* constant_mass_integrand,
    Integer nelem,
    Integer ndim,
    Integer nvar,
    Integer ngauss,
    Integer nodeperelem,
    Integer local_capacity,
    Integer mass_type,
    const Integer* local_rows_mass,
    const Integer* local_cols_mass,
    int *I_mass,
    int *J_mass,
    Real *V_mass,
    Real *mass,
    int recompute_sparsity_pattern,
    int squeeze_sparsity_pattern,
    const int *data_local_indices,
    const int *data_global_indices,
    const UInteger *sorted_elements,
    const Integer *sorter
    ) {

    const Integer ndof                      = nodeperelem*nvar;
    Real *LagrangeElemCoords                = allocate<Real>(nodeperelem*ndim);
    Real *MaterialGradient                  = allocate<Real>(ndim*nodeperelem);
    // +1 TO AVOID OVERSTEPPONG IN MATMUL UNALIGNED STORE - EXTREMELY DANGEROUSE BUG
    int plus1                               = ndim == 2 ? 0 : 1;
    Real *ParentGradientX                   = allocate<Real>(ndim*ndim+plus1);
    Real detJ                               = 0.;
    Real *massel                            = allocate<Real>(local_capacity);
    Real *massel_lumped                     = allocate<Real>(ndof);

    // PRE-COMPUTE SPARSITY PATTERN OF THE LOCAL MATRIX -
    // LOCAL MASS MATRICES TYPICALLY HAVE LOTS OF NONZEROS
    constexpr Real tolerance = 1e-12;
    std::vector<Integer> nonzero_i_indices, nonzero_j_indices;
    Integer counter = 0;
    for (Integer i=0; i<ndof; ++i) {
        for (Integer j=i; j<ndof; ++j) {
            if (std::abs(constant_mass_integrand[i*ndof+j]) > tolerance) {
                nonzero_i_indices.push_back(i);
                nonzero_j_indices.push_back(j);
            }
        }
    }
    // PRE-COMPUTE SYMMETRIC PART OF CONSTANT MASS WITHOUT ZEROS
    const Integer symmetric_local_capacity = nonzero_i_indices.size();
    Real *symmetric_massel                  = allocate<Real>(symmetric_local_capacity);
    Real *symmetric_constant_mass_integrand = allocate<Real>(ngauss*symmetric_local_capacity);
    counter = 0;
    for (Integer igauss = 0; igauss < ngauss; ++igauss) {
        for (Integer i=0; i<symmetric_local_capacity; ++i) {
            const Real value = constant_mass_integrand[igauss*ndof*ndof+nonzero_i_indices[i]*ndof+nonzero_j_indices[i]];
            symmetric_constant_mass_integrand[counter] = value;
            counter++;
        }
    }

    // PRE-COMPUTE ISOPARAMETRIC GRADIENTS
    std::vector<std::vector<Real>> current_Jms(ngauss);
    for (int g=0; g<ngauss; ++g) {
        std::vector<Real> current_Jm(ndim*nodeperelem);
        for (int j=0; j<nodeperelem; ++j) {
            for (int k=0; k<ndim; ++k) {
                current_Jm[k*nodeperelem+j] = Jm[k*ngauss*nodeperelem+j*ngauss+g];
            }
        }
        current_Jms[g] = current_Jm;
    }

    // LOOP OVER ELEMETNS
    for (Integer elem=0; elem < nelem; ++elem) {

        // GET THE FIELDS AT THE ELEMENT LEVEL
        for (Integer i=0; i<nodeperelem; ++i) {
            const Integer inode = elements[elem*nodeperelem+i];
            for (Integer j=0; j<ndim; ++j) {
                LagrangeElemCoords[i*ndim+j] = points[inode*ndim+j];
            }
        }

        std::fill_n(massel,local_capacity,0.);
        std::fill_n(symmetric_massel,symmetric_local_capacity,0.);

        for (Integer igauss = 0; igauss < ngauss; ++igauss) {

            {
                // USING A STL BASED FILLER REMOVES THE ANNOYING BUG
                std::fill_n(ParentGradientX,ndim*ndim+plus1,0.);
                _matmul_(ndim,ndim,nodeperelem,current_Jms[igauss].data(),LagrangeElemCoords,ParentGradientX);
                const Real detX = _det_(ndim, ParentGradientX);
                detJ = AllGauss[igauss]*std::abs(detX);
            }

            // Multiply constant part of mass with detJ
#ifdef __AVX__
            using V2 = Fastor::SIMDVector<Real>;
            V2 _va, _vb, _vout;
            _vb.set(detJ);
            int Vsize = V2::Size;
            int ROUND_ = ROUND_DOWN(symmetric_local_capacity,Vsize);
            int i=0;
            for (; i<ROUND_; i+=Vsize) {
                _va.load(&symmetric_constant_mass_integrand[igauss*symmetric_local_capacity+i],false);
                _vout.load(&symmetric_massel[i],false);
#ifdef __FMA__
                _vout = fmadd(_va,_vb,_vout);
#else
                _vout += _va*_vb;
#endif
                _vout.store(&symmetric_massel[i],false);
            }
            for ( ; i<symmetric_local_capacity; ++i) {
                symmetric_massel[i] += symmetric_constant_mass_integrand[igauss*symmetric_local_capacity+i]*detJ;
            }
#else
            for (int i=0; i<symmetric_local_capacity; ++i) {
                symmetric_massel[i] += symmetric_constant_mass_integrand[igauss*symmetric_local_capacity+i]*detJ;
            }
#endif
        }

        // BUILD TOTAL LOCAL MASS MATRIX
        for (Integer i=0; i<symmetric_local_capacity; ++i) {
            const Real value = symmetric_massel[i];
            massel[nonzero_i_indices[i]*ndof+nonzero_j_indices[i]] = value;

        }

        for (Integer i=0; i<ndof; ++i) {
            for (Integer j=i; j<ndof; ++j) {
                massel[j*ndof+i] = massel[i*ndof+j];

            }
        }

        // FOR LUMP MASS MATRIX
        if (mass_type == 0) {
            // LUMP MASS
            for (Integer i=0; i<ndof; ++i) {
                massel_lumped[i] = std::accumulate(&massel[i*ndof], &massel[(i+1)*ndof],0.);
            }

            // ASSEMBLE LUMPED MASS
            {
                for (Integer i = 0; i<nodeperelem; ++i) {
                    UInteger T_idx = elements[elem*nodeperelem+i]*nvar;
                    for (Integer iterator = 0; iterator < nvar; ++iterator) {
                        mass[T_idx+iterator] += massel_lumped[i*nvar+iterator];
                    }
                }
            }
        }
        else {
            fill_global_data(   local_rows_mass,
                                local_cols_mass,
                                massel,
                                I_mass,
                                J_mass,
                                V_mass,
                                elem,
                                nvar,
                                nodeperelem,
                                elements,
                                local_capacity,
                                local_capacity,
                                recompute_sparsity_pattern,
                                squeeze_sparsity_pattern,
                                data_local_indices,
                                data_global_indices,
                                sorted_elements,
                                sorter);
        }
    }

    deallocate(ParentGradientX);
    deallocate(MaterialGradient);
    deallocate(LagrangeElemCoords);
    deallocate(massel);
    deallocate(symmetric_constant_mass_integrand);
    deallocate(symmetric_massel);
    deallocate(massel_lumped);
}

inline void _SymmetricNonZeroOneElementConstantMassIntegrand_(
    const UInteger* elements,
    const Real* points,
    const Real* Jm,
    const Real* AllGauss,
    const Real* constant_mass_integrand,
    Integer nelem,
    Integer ndim,
    Integer nvar,
    Integer ngauss,
    Integer nodeperelem,
    Integer local_capacity,
    Integer mass_type,
    const Integer* local_rows_mass,
    const Integer* local_cols_mass,
    int *I_mass,
    int *J_mass,
    Real *V_mass,
    Real *mass,
    int recompute_sparsity_pattern,
    int squeeze_sparsity_pattern,
    const int *data_local_indices,
    const int *data_global_indices,
    const UInteger *sorted_elements,
    const Integer *sorter
    ) {

    const Integer ndof                      = nodeperelem*nvar;
    Real *LagrangeElemCoords                = allocate<Real>(nodeperelem*ndim);
    Real *MaterialGradient                  = allocate<Real>(ndim*nodeperelem);
    // +1 TO AVOID OVERSTEPPONG IN MATMUL UNALIGNED STORE - EXTREMELY DANGEROUSE BUG
    int plus1                               = ndim == 2 ? 0 : 1;
    Real *ParentGradientX                   = allocate<Real>(ndim*ndim+plus1);
    Real detJ                               = 0.;
    Real *massel                            = allocate<Real>(local_capacity);
    Real *massel_lumped                     = allocate<Real>(ndof);

    // PRE-COMPUTE SPARSITY PATTERN OF THE LOCAL MATRIX -
    // LOCAL MASS MATRICES TYPICALLY HAVE LOTS OF NONZEROS
    constexpr Real tolerance = 1e-12;
    std::vector<Integer> nonzero_i_indices, nonzero_j_indices;
    Integer counter = 0;
    for (Integer i=0; i<ndof; ++i) {
        for (Integer j=i; j<ndof; ++j) {
            if (std::abs(constant_mass_integrand[i*ndof+j]) > tolerance) {
                nonzero_i_indices.push_back(i);
                nonzero_j_indices.push_back(j);
            }
        }
    }
    // PRE-COMPUTE SYMMETRIC PART OF CONSTANT MASS WITHOUT ZEROS
    const Integer symmetric_local_capacity = nonzero_i_indices.size();
    Real *symmetric_massel                  = allocate<Real>(symmetric_local_capacity);
    Real *symmetric_constant_mass_integrand = allocate<Real>(ngauss*symmetric_local_capacity);
    counter = 0;
    for (Integer igauss = 0; igauss < ngauss; ++igauss) {
        for (Integer i=0; i<symmetric_local_capacity; ++i) {
            const Real value = constant_mass_integrand[igauss*ndof*ndof+nonzero_i_indices[i]*ndof+nonzero_j_indices[i]];
            symmetric_constant_mass_integrand[counter] = value;
            counter++;
        }
    }

    // PRE-COMPUTE ISOPARAMETRIC GRADIENTS
    std::vector<std::vector<Real>> current_Jms(ngauss);
    for (int g=0; g<ngauss; ++g) {
        std::vector<Real> current_Jm(ndim*nodeperelem);
        for (int j=0; j<nodeperelem; ++j) {
            for (int k=0; k<ndim; ++k) {
                current_Jm[k*nodeperelem+j] = Jm[k*ngauss*nodeperelem+j*ngauss+g];
            }
        }
        current_Jms[g] = current_Jm;
    }

    // COMPUTE MASS FOR ONE ELEMENT ONLY
    for (Integer elem=0; elem < 1; ++elem) {

        // GET THE FIELDS AT THE ELEMENT LEVEL
        for (Integer i=0; i<nodeperelem; ++i) {
            const Integer inode = elements[elem*nodeperelem+i];
            for (Integer j=0; j<ndim; ++j) {
                LagrangeElemCoords[i*ndim+j] = points[inode*ndim+j];
            }
        }

        std::fill_n(massel,local_capacity,0.);
        std::fill_n(symmetric_massel,symmetric_local_capacity,0.);

        for (Integer igauss = 0; igauss < ngauss; ++igauss) {

            {
                // USING A STL BASED FILLER REMOVES THE ANNOYING BUG
                std::fill_n(ParentGradientX,ndim*ndim+plus1,0.);
                _matmul_(ndim,ndim,nodeperelem,current_Jms[igauss].data(),LagrangeElemCoords,ParentGradientX);
                const Real detX = _det_(ndim, ParentGradientX);
                detJ = AllGauss[igauss]*std::abs(detX);
            }

            // Multiply constant part of mass with detJ
#ifdef __AVX__
            using V2 = Fastor::SIMDVector<Real>;
            V2 _va, _vb, _vout;
            _vb.set(detJ);
            int Vsize = V2::Size;
            int ROUND_ = ROUND_DOWN(symmetric_local_capacity,Vsize);
            int i=0;
            for (; i<ROUND_; i+=Vsize) {
                _va.load(&symmetric_constant_mass_integrand[igauss*symmetric_local_capacity+i],false);
                _vout.load(&symmetric_massel[i],false);
#ifdef __FMA__
                _vout = fmadd(_va,_vb,_vout);
#else
                _vout += _va*_vb;
#endif
                _vout.store(&symmetric_massel[i],false);
            }
            for ( ; i<symmetric_local_capacity; ++i) {
                symmetric_massel[i] += symmetric_constant_mass_integrand[igauss*symmetric_local_capacity+i]*detJ;
            }
#else
            for (int i=0; i<symmetric_local_capacity; ++i) {
                symmetric_massel[i] += symmetric_constant_mass_integrand[igauss*symmetric_local_capacity+i]*detJ;
            }
#endif
        }

        // BUILD TOTAL LOCAL MASS MATRIX
        for (Integer i=0; i<symmetric_local_capacity; ++i) {
            const Real value = symmetric_massel[i];
            massel[nonzero_i_indices[i]*ndof+nonzero_j_indices[i]] = value;
        }

        for (Integer i=0; i<ndof; ++i) {
            for (Integer j=i; j<ndof; ++j) {
                massel[j*ndof+i] = massel[i*ndof+j];
            }
        }

        // FOR LUMP MASS MATRIX
        if (mass_type == 0) {
            // LUMP MASS
            for (Integer i=0; i<ndof; ++i) {
                massel_lumped[i] = std::accumulate(&massel[i*ndof], &massel[(i+1)*ndof],0.);
            }
        }
    }

    // LOOP OVER ELEMETNS
    for (Integer elem=0; elem < nelem; ++elem) {

        // FOR LUMP MASS MATRIX
        if (mass_type == 0) {
            // ASSEMBLE LUMPED MASS
            {
                for (Integer i = 0; i<nodeperelem; ++i) {
                    UInteger T_idx = elements[elem*nodeperelem+i]*nvar;
                    for (Integer iterator = 0; iterator < nvar; ++iterator) {
                        mass[T_idx+iterator] += massel_lumped[i*nvar+iterator];
                    }
                }
            }
        }
        else {
            fill_global_data(   local_rows_mass,
                                local_cols_mass,
                                massel,
                                I_mass,
                                J_mass,
                                V_mass,
                                elem,
                                nvar,
                                nodeperelem,
                                elements,
                                local_capacity,
                                local_capacity,
                                recompute_sparsity_pattern,
                                squeeze_sparsity_pattern,
                                data_local_indices,
                                data_global_indices,
                                sorted_elements,
                                sorter);
        }
    }

    deallocate(ParentGradientX);
    deallocate(MaterialGradient);
    deallocate(LagrangeElemCoords);
    deallocate(massel);
    deallocate(symmetric_constant_mass_integrand);
    deallocate(symmetric_massel);
    deallocate(massel_lumped);
}

#endif

max_stars_repo_name: jdlaubrie/florence

max_issues_repo_name: jdlaubrie/florence

max_forks_repo_name: jdlaubrie/florence


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

#ifndef __GSL_MULTISET_H__
#define __GSL_MULTISET_H__

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

#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

struct gsl_multiset_struct
{
  size_t n;
  size_t k;
  size_t *data;
};

typedef struct gsl_multiset_struct gsl_multiset;

GSL_FUN gsl_multiset *gsl_multiset_alloc (const size_t n, const size_t k);
GSL_FUN gsl_multiset *gsl_multiset_calloc (const size_t n, const size_t k);
GSL_FUN void gsl_multiset_init_first (gsl_multiset * c);
GSL_FUN void gsl_multiset_init_last (gsl_multiset * c);
GSL_FUN void gsl_multiset_free (gsl_multiset * c);
GSL_FUN int gsl_multiset_memcpy (gsl_multiset * dest, const gsl_multiset * src);

GSL_FUN int gsl_multiset_fread (FILE * stream, gsl_multiset * c);
GSL_FUN int gsl_multiset_fwrite (FILE * stream, const gsl_multiset * c);
GSL_FUN int gsl_multiset_fscanf (FILE * stream, gsl_multiset * c);
GSL_FUN int gsl_multiset_fprintf (FILE * stream, const gsl_multiset * c, const char *format);

GSL_FUN size_t gsl_multiset_n (const gsl_multiset * c);
GSL_FUN size_t gsl_multiset_k (const gsl_multiset * c);
GSL_FUN size_t * gsl_multiset_data (const gsl_multiset * c);

GSL_FUN int gsl_multiset_valid (gsl_multiset * c);
GSL_FUN int gsl_multiset_next (gsl_multiset * c);
GSL_FUN int gsl_multiset_prev (gsl_multiset * c);

GSL_FUN INLINE_DECL size_t gsl_multiset_get (const gsl_multiset * c, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
size_t
gsl_multiset_get (const gsl_multiset * c, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= c->k)) /* size_t is unsigned, can't be negative */
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return c->data[i];
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_MULTISET_H__ */

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA