<commit_before>#include "ghost/config.h"
#include "ghost/omp.h"

#ifdef GHOST_HAVE_MPI
#include <mpi.h> //mpi.h has to be included before stdio.h
#endif
#include <cstdlib>
#include <iostream>
#include <complex>
#include <stdio.h>

#include "ghost/complex.h"
#include "ghost/rand.h"
#include "ghost/util.h"
#include "ghost/densemat_cm.h"
#include "ghost/math.h"
#include "ghost/locality.h"
#include "ghost/log.h"


using namespace std;


template <typename v_t> 
static ghost_error_t ghost_densemat_cm_normalize_tmpl(ghost_densemat_t *vec)
{
    ghost_error_t ret = GHOST_SUCCESS;
    ghost_idx_t v;
    v_t *s = NULL;

    GHOST_CALL_GOTO(ghost_malloc((void **)&s,vec->traits.ncols*sizeof(v_t)),err,ret);
    GHOST_CALL_GOTO(ghost_dot(s,vec,vec),err,ret);

    for (v=0; v<vec->traits.ncols; v++)
    {
        s[v] = (v_t)sqrt(s[v]);
        s[v] = (v_t)(((v_t)1.)/s[v]);
    }
    GHOST_CALL_GOTO(vec->vscale(vec,s),err,ret);

    goto out;
err:

out:
    free(s);

    return ret;

}


template <typename v_t> 
static ghost_error_t ghost_densemat_cm_dotprod_tmpl(ghost_densemat_t *vec, void *res, ghost_densemat_t *vec2)
{ // the parallelization is done manually because reduction does not work with ghost_complex numbers
    if (vec->traits.nrows != vec2->traits.nrows) {
        WARNING_LOG("The input vectors of the dot product have different numbers of rows");
    }
    if (vec->traits.ncols != vec2->traits.ncols) {
        WARNING_LOG("The input vectors of the dot product have different numbers of columns");
    }
    ghost_idx_t i,v,rowidx = 0;

    int nthreads;
#pragma omp parallel
#pragma omp single
    nthreads = ghost_omp_nthread();

    v_t *partsums;
    GHOST_CALL_RETURN(ghost_malloc((void **)&partsums,16*nthreads*sizeof(v_t)));

    for (v=0; v<MIN(vec->traits.ncols,vec2->traits.ncols); v++) {
        v_t sum = 0;
        for (i=0; i<nthreads*16; i++) partsums[i] = (v_t)0.;

#pragma omp parallel 
        {
            int tid = ghost_omp_threadnum();
            ITER_ROWS_BEGIN(vec,i,rowidx)
                partsums[tid*16] += *(v_t *)VECVAL_CM(vec2,vec2->val,v,i) * conjugate((v_t *)(VECVAL_CM(vec,vec->val,v,i)));
            ITER_ROWS_END(rowidx)
        }

        for (i=0; i<nthreads; i++) sum += partsums[i*16];

        ((v_t *)res)[v] = sum;
    }

    free(partsums);
    
    return GHOST_SUCCESS;
}

template <typename v_t> 
static ghost_error_t ghost_densemat_cm_vaxpy_tmpl(ghost_densemat_t *vec, ghost_densemat_t *vec2, void *scale)
{
    if (vec->traits.storage != vec2->traits.storage) {
        ERROR_LOG("Cannot VAXPY densemats with different storage order");
        return GHOST_ERR_INVALID_ARG;
    }
    if (vec->traits.ncols != vec2->traits.ncols) {
        ERROR_LOG("Cannot VAXPY densemats with different number of columns");
        return GHOST_ERR_INVALID_ARG;
    }
    if (vec->traits.nrows != vec2->traits.nrows) {
        ERROR_LOG("Cannot VAXPY densemats with different number of rows");
        return GHOST_ERR_INVALID_ARG;
    }
    ghost_idx_t row1,row2,col,rowidx;
    v_t *s = (v_t *)scale;
    
    ITER2_BEGIN_CM(vec,vec2,col,row1,row2,rowidx)
        *(v_t *)VECVAL_CM(vec,vec->val,col,row1) += *(v_t *)VECVAL_CM(vec2,vec2->val,col,row2) * s[col];
    ITER2_END_CM(rowidx)

    return GHOST_SUCCESS;
}

template <typename v_t> 
static ghost_error_t ghost_densemat_cm_vaxpby_tmpl(ghost_densemat_t *vec, ghost_densemat_t *vec2, void *scale, void *b_)
{
    if (vec->traits.storage != vec2->traits.storage) {
        ERROR_LOG("Cannot VAXPBY densemats with different storage order");
        return GHOST_ERR_INVALID_ARG;
    }
    if (vec->traits.ncols != vec2->traits.ncols) {
        ERROR_LOG("Cannot VAXPBY densemats with different number of columns");
        return GHOST_ERR_INVALID_ARG;
    }
    if (vec->traits.nrows != vec2->traits.nrows) {
        ERROR_LOG("Cannot VAXPBY densemats with different number of rows");
        return GHOST_ERR_INVALID_ARG;
    }
    ghost_idx_t row1,row2,col,rowidx;
    v_t *s = (v_t *)scale;
    v_t *b = (v_t *)b_;
    
    ITER2_BEGIN_CM(vec,vec2,col,row1,row2,rowidx)
        *(v_t *)VECVAL_CM(vec,vec->val,col,row1) = *(v_t *)VECVAL_CM(vec2,vec2->val,col,row2) * s[col] +
            *(v_t *)VECVAL_CM(vec,vec->val,col,row1) * b[col];
    ITER2_END_CM(rowidx)

    return GHOST_SUCCESS;
}
<commit_msg>enhance ghost_densemat_cm_vaxpy_tmpl()<commit_after>#include "ghost/config.h"
#include "ghost/omp.h"

#ifdef GHOST_HAVE_MPI
#include <mpi.h> //mpi.h has to be included before stdio.h
#endif
#include <cstdlib>
#include <iostream>
#include <complex>
#include <stdio.h>

#include "ghost/complex.h"
#include "ghost/rand.h"
#include "ghost/util.h"
#include "ghost/densemat_cm.h"
#include "ghost/math.h"
#include "ghost/locality.h"
#include "ghost/log.h"


using namespace std;


template <typename v_t> 
static ghost_error_t ghost_densemat_cm_normalize_tmpl(ghost_densemat_t *vec)
{
    ghost_error_t ret = GHOST_SUCCESS;
    ghost_idx_t v;
    v_t *s = NULL;

    GHOST_CALL_GOTO(ghost_malloc((void **)&s,vec->traits.ncols*sizeof(v_t)),err,ret);
    GHOST_CALL_GOTO(ghost_dot(s,vec,vec),err,ret);

    for (v=0; v<vec->traits.ncols; v++)
    {
        s[v] = (v_t)sqrt(s[v]);
        s[v] = (v_t)(((v_t)1.)/s[v]);
    }
    GHOST_CALL_GOTO(vec->vscale(vec,s),err,ret);

    goto out;
err:

out:
    free(s);

    return ret;

}


template <typename v_t> 
static ghost_error_t ghost_densemat_cm_dotprod_tmpl(ghost_densemat_t *vec, void *res, ghost_densemat_t *vec2)
{ // the parallelization is done manually because reduction does not work with ghost_complex numbers
    if (vec->traits.nrows != vec2->traits.nrows) {
        WARNING_LOG("The input vectors of the dot product have different numbers of rows");
    }
    if (vec->traits.ncols != vec2->traits.ncols) {
        WARNING_LOG("The input vectors of the dot product have different numbers of columns");
    }
    ghost_idx_t i,v,rowidx = 0;

    int nthreads;
#pragma omp parallel
#pragma omp single
    nthreads = ghost_omp_nthread();

    v_t *partsums;
    GHOST_CALL_RETURN(ghost_malloc((void **)&partsums,16*nthreads*sizeof(v_t)));

    for (v=0; v<MIN(vec->traits.ncols,vec2->traits.ncols); v++) {
        v_t sum = 0;
        for (i=0; i<nthreads*16; i++) partsums[i] = (v_t)0.;

#pragma omp parallel 
        {
            int tid = ghost_omp_threadnum();
            ITER_ROWS_BEGIN(vec,i,rowidx)
                partsums[tid*16] += *(v_t *)VECVAL_CM(vec2,vec2->val,v,i) * conjugate((v_t *)(VECVAL_CM(vec,vec->val,v,i)));
            ITER_ROWS_END(rowidx)
        }

        for (i=0; i<nthreads; i++) sum += partsums[i*16];

        ((v_t *)res)[v] = sum;
    }

    free(partsums);
    
    return GHOST_SUCCESS;
}

template <typename v_t> 
static ghost_error_t ghost_densemat_cm_vaxpy_tmpl(ghost_densemat_t *vec, ghost_densemat_t *vec2, void *scale)
{
    if (vec->traits.storage != vec2->traits.storage) {
        ERROR_LOG("Cannot VAXPY densemats with different storage order");
        return GHOST_ERR_INVALID_ARG;
    }
    if (vec->traits.ncols != vec2->traits.ncols) {
        ERROR_LOG("Cannot VAXPY densemats with different number of columns");
        return GHOST_ERR_INVALID_ARG;
    }
    if (vec->traits.nrows != vec2->traits.nrows) {
        ERROR_LOG("Cannot VAXPY densemats with different number of rows");
        return GHOST_ERR_INVALID_ARG;
    }
    //ghost_idx_t row1,row2,col,rowidx;
    v_t *s = (v_t *)scale;
    
    //ITER2_BEGIN_CM(vec,vec2,col,row1,row2,rowidx)
    //    *(v_t *)VECVAL_CM(vec,vec->val,col,row1) += *(v_t *)VECVAL_CM(vec2,vec2->val,col,row2) * s[col];
    //ITER2_END_CM(rowidx)

    ghost_idx_t i,v,rowidx = 0;

    for (v=0; v<MIN(vec->traits.ncols,vec2->traits.ncols); v++) {
            ITER_ROWS_BEGIN(vec,i,rowidx)
              *(v_t *)VECVAL_CM(vec,vec->val,v,i) += *(v_t *)VECVAL_CM(vec2,vec2->val,v,i) * s[v];
            ITER_ROWS_END(rowidx)
    }
    
    return GHOST_SUCCESS;
}

template <typename v_t> 
static ghost_error_t ghost_densemat_cm_vaxpby_tmpl(ghost_densemat_t *vec, ghost_densemat_t *vec2, void *scale, void *b_)
{
    if (vec->traits.storage != vec2->traits.storage) {
        ERROR_LOG("Cannot VAXPBY densemats with different storage order");
        return GHOST_ERR_INVALID_ARG;
    }
    if (vec->traits.ncols != vec2->traits.ncols) {
        ERROR_LOG("Cannot VAXPBY densemats with different number of columns");
        return GHOST_ERR_INVALID_ARG;
    }
    if (vec->traits.nrows != vec2->traits.nrows) {
        ERROR_LOG("Cannot VAXPBY densemats with different number of rows");
        return GHOST_ERR_INVALID_ARG;
    }
    ghost_idx_t row1,row2,col,rowidx;
    v_t *s = (v_t *)scale;
    v_t *b = (v_t *)b_;
    
    ITER2_BEGIN_CM(vec,vec2,col,row1,row2,rowidx)
        *(v_t *)VECVAL_CM(vec,vec->val,col,row1) = *(v_t *)VECVAL_CM(vec2,vec2->val,col,row2) * s[col] +
            *(v_t *)VECVAL_CM(vec,vec->val,col,row1) * b[col];
    ITER2_END_CM(rowidx)

    return GHOST_SUCCESS;
}

{'fraction_non_alphanumeric': 0.14159785762106672, 'fraction_numerical': 0.012385628207989288, 'mean_word_length': 3.492731829573935, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 81, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1958637', 'n_tokens_mistral': 3667, 'n_tokens_neox': 3520, 'n_words': 657}