You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
217 lines
6.7 KiB
C
217 lines
6.7 KiB
C
/*
|
|
* Sp3MM_for_AlgebraicMultiGrid
|
|
* (C) Copyright 2021-2022
|
|
* Andrea Di Iorio
|
|
*
|
|
* Redistribution and use in source and binary forms, with or without
|
|
* modification, are permitted provided that the following conditions
|
|
* are met:
|
|
* 1. Redistributions of source code must retain the above copyright
|
|
* notice, this list of conditions and the following disclaimer.
|
|
* 2. Redistributions in binary form must reproduce the above copyright
|
|
* notice, this list of conditions, and the following disclaimer in the
|
|
* documentation and/or other materials provided with the distribution.
|
|
* 3. The name of the Sp3MM_for_AlgebraicMultiGrid or the names of its contributors may
|
|
* not be used to endorse or promote products derived from this
|
|
* software without specific written permission.
|
|
*
|
|
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
|
* ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
|
|
* TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
|
|
* PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE Sp3MM_for_AlgebraicMultiGrid GROUP OR ITS CONTRIBUTORS
|
|
* BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
|
|
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
|
|
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
|
|
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
|
|
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
|
|
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
|
* POSSIBILITY OF SUCH DAMAGE.
|
|
*/
|
|
|
|
//sparse matrix def & aux
|
|
#ifndef SPARSEMATRIX
|
|
#define SPARSEMATRIX
|
|
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
|
|
#include "macros.h"
|
|
#include "config.h"
|
|
|
|
typedef struct {
|
|
idx_t NZ,M,N;
|
|
double *AS;
|
|
idx_t* JA;
|
|
//CSR SPECIFIC
|
|
idx_t* IRP;
|
|
#ifdef ROWLENS
|
|
idx_t* RL; //row lengths
|
|
#endif
|
|
//CUDA SPECIFIC
|
|
idx_t MAX_ROW_NZ;
|
|
|
|
} spmat; //describe a sparse matrix
|
|
|
|
//smart index keeping in a dense map
|
|
typedef struct{
|
|
idx_t len; //num of nnz idx accumulated
|
|
/* nnz index presence packing, implict space enough for all possible indexes*/
|
|
nnz_idxs_flags_t idxsMap;
|
|
uint idxsMapN; //either num of limbs or len of char flag array
|
|
} SPVECT_IDX_DENSE_MAP;
|
|
int initSpVectIdxDenseAcc(idx_t idxMax, SPVECT_IDX_DENSE_MAP* );
|
|
inline void _resetIdxMap(SPVECT_IDX_DENSE_MAP* acc){
|
|
acc->len = 0;
|
|
memset(acc->idxsMap,0,sizeof(*acc->idxsMap)*acc->idxsMapN);
|
|
}
|
|
inline void _freeIdxMap(SPVECT_IDX_DENSE_MAP* acc){
|
|
free(acc->idxsMap);
|
|
free(acc);
|
|
}
|
|
//aux struct for sparse vector-scalar product accumualtion
|
|
typedef struct{
|
|
double* v; //aux accumulating dense vector (sparse)
|
|
idx_t* nnzIdx; //v nnz value's indexes (contiguos)
|
|
idx_t vLen; //size of the aux dense vector //TODO USELESS?
|
|
SPVECT_IDX_DENSE_MAP nnzIdxMap; //
|
|
} ACC_DENSE;
|
|
int allocAccDense(ACC_DENSE* v, ulong size);
|
|
|
|
/*
|
|
**** NNZ IDXS PRESENCE FLAGS ACCESS INTERFACE: ***
|
|
* sp_vect_idx_set(idx,SPVECT_IDX_DENSE_MAP)
|
|
* -> return 0 if idx isn't already setted and in case set it
|
|
*/
|
|
inline int spVect_idx_in(idx_t idx, SPVECT_IDX_DENSE_MAP* idxsMapAcc){
|
|
#if SPVECT_IDX_BITWISE == TRUE
|
|
uint limbID = idx / LIMB_SIZE_BIT; //idx's limb id
|
|
uint limbIdxID = idx % LIMB_SIZE_BIT; //idx's pos in limb
|
|
DEBUGCHECKS assert( limbID < idxsMapAcc->idxsMapN );
|
|
limb_t idxPos = ((limb_t) 1) << limbIdxID;
|
|
if (!(idxsMapAcc->idxsMap[limbID] & idxPos) ){
|
|
idxsMapAcc->idxsMap[limbID] |= idxPos;
|
|
idxsMapAcc->len++;
|
|
return 0;
|
|
}
|
|
#else
|
|
//assert( idx < idxsMapAcc->idxsMapN ); //TODO
|
|
if (!( idxsMapAcc->idxsMap[idx] )){ //TODO usabile primitiva atomica di cmpswap, accorpamento ops a livello HW?
|
|
idxsMapAcc->idxsMap[idx] = 1;
|
|
idxsMapAcc->len++;
|
|
return 0;
|
|
}
|
|
#endif //SPVECT_IDX_BITWISE == TRUE
|
|
return 1;
|
|
}
|
|
inline void _resetAccVect(ACC_DENSE* acc){
|
|
memset(acc->v, 0, acc->vLen * sizeof(*(acc->v)));
|
|
memset(acc->nnzIdx, 0, acc->vLen * sizeof(*(acc->nnzIdx)));
|
|
_resetIdxMap(&acc->nnzIdxMap);
|
|
}
|
|
////Sparse vector accumulator -- corresponding to a matrix portion
|
|
typedef struct{
|
|
//idx_t r; //row index in the corresponding matrix
|
|
//idx_t c; //col index in the corresponding matrix
|
|
idx_t len; //rowLen
|
|
double* AS; //row nnz values
|
|
idx_t* JA; //row nnz colIndexes
|
|
} SPACC;
|
|
|
|
|
|
/*
|
|
* ARRAY BISECTION - RECURSIVE VERSION
|
|
* TODO ASSERT LEN>0 ommitted
|
|
*/
|
|
inline int BISECT_ARRAY(idx_t target, idx_t* arr, idx_t len){
|
|
//if (len == 0) return FALSE;
|
|
if (len <= 1) return *arr == target;
|
|
idx_t middleIdx = (len-1) / 2; //len=5-->2, len=4-->1
|
|
idx_t middle = arr[ middleIdx ];
|
|
if (target == middle) return TRUE;
|
|
else if (target < middle) return BISECT_ARRAY(target,arr,middleIdx);
|
|
else return BISECT_ARRAY(target,arr+middleIdx+1,middleIdx + (len-1)%2);
|
|
}
|
|
|
|
/*
|
|
* return !0 if col @j idx is in row @i of sparse mat @smat
|
|
* bisection used --> O(log_2(ROWLENGHT))
|
|
*/
|
|
inline int IS_NNZ(spmat* smat,idx_t i,idx_t j){
|
|
idx_t rStart = smat->IRP[i];
|
|
idx_t rLen = smat->IRP[i+1] - rStart;
|
|
if (!rLen) return FALSE;
|
|
return BISECT_ARRAY(j,smat->JA + rStart,rLen);
|
|
}
|
|
inline int IS_NNZ_linear(spmat* smat,idx_t i,idx_t j){ //linear -> O(ROWLENGHT)
|
|
int out = 0;
|
|
for (idx_t x=smat->IRP[i]; x<smat->IRP[i+1] && !out; x++){
|
|
out = (j == smat->JA[x]);
|
|
}
|
|
return out;
|
|
}
|
|
////aux functions
|
|
//free sparse matrix
|
|
inline void freeSpmatInternal(spmat* mat){
|
|
if(!mat) return;
|
|
free(mat->AS);
|
|
free(mat->JA);
|
|
free(mat->IRP);
|
|
#ifdef ROWLENS
|
|
free(mat->RL);
|
|
#endif
|
|
}
|
|
|
|
inline void freeSpmat(spmat* mat){
|
|
if (!mat) return;
|
|
freeSpmatInternal(mat);
|
|
free(mat);
|
|
}
|
|
|
|
static inline void zeroSpmat(spmat* m) {
|
|
memset(m->AS, 0, sizeof(*(m->AS)) * m->NZ);
|
|
memset(m->JA, 0, sizeof(*(m->JA)) * m->NZ);
|
|
memset(m->IRP, 0, sizeof(*(m->IRP)) * (m->M + 1));
|
|
}
|
|
|
|
//free max aux structs not NULL pointed
|
|
inline void freeSpAcc(SPACC* r){
|
|
free(r->AS);
|
|
free(r->JA);
|
|
}
|
|
////alloc&init functions
|
|
//alloc&init internal structures only dependent of dimensions @rows,@cols
|
|
inline int allocSpMatrixInternal(idx_t rows, idx_t cols, spmat* mat){
|
|
mat -> M = rows;
|
|
mat -> N = cols;
|
|
if (!(mat->IRP=calloc(mat->M+1,sizeof(*(mat->IRP))))){ //calloc only for 0th
|
|
ERRPRINT("IRP calloc err\n");
|
|
freeSpmatInternal(mat);
|
|
return EXIT_FAILURE;
|
|
}
|
|
#ifdef ROWLENS
|
|
if (!(mat->RL = malloc(mat->M*sizeof(*(mat->RL))))){
|
|
ERRPRINT("RL calloc err\n");
|
|
freeSpmatInternal(mat);
|
|
return EXIT_FAILURE;
|
|
}
|
|
#endif
|
|
return EXIT_SUCCESS;
|
|
}
|
|
|
|
//alloc a sparse matrix of @rows rows and @cols cols
|
|
inline spmat* allocSpMatrix(idx_t rows, idx_t cols){
|
|
|
|
spmat* mat;
|
|
if (!(mat = calloc(1,sizeof(*mat)))) {
|
|
ERRPRINT("mat calloc failed\n");
|
|
return NULL;
|
|
}
|
|
if (allocSpMatrixInternal(rows,cols,mat)){
|
|
free(mat);
|
|
return NULL;
|
|
}
|
|
return mat;
|
|
}
|
|
|
|
#endif
|