Mentions légales du service

Skip to content
Snippets Groups Projects
Commit a47eaafd authored by Mathieu Faverge's avatar Mathieu Faverge
Browse files

MergeDuplicate: update for distributed

parent 87038c1d
No related branches found
No related tags found
1 merge request!36Make spmSymmetrize and spmMergeDuplicate distributed and multi-dofs
Hide whitespace changes
Inline Side-by-side
CMakeLists.txt
+ 1
− 1
View file @ a47eaafd
......@@ -191,7 +191,6 @@ add_custom_target(spm_headers_tgt
set(generated_sources "")
set(SOURCES
src/z_spm.c
src/z_spm_2dense.c
src/z_spm_dof_extend.c
src/z_spm_norm.c
......@@ -206,6 +205,7 @@ set(SOURCES
src/z_spm_integer.c
src/z_spm_laplacian.c
src/z_spm_matrixvector.c
src/z_spm_mergeduplicate.c
src/z_spm_print.c
src/z_spm_sort.c
)
......
src/spm.c
+ 38
− 8
View file @ a47eaafd
......@@ -545,9 +545,10 @@ spmSort( spmatrix_t *spm )
*
* @brief Merge multiple entries in a spm by summing their values together.
*
* The sparse matrix needs to be sorted first (see spmSort()).
* The sparse matrix needs to be sorted first (see z_spmSort()). In distributed,
* only local entries are merged together.
*
* @warning Not implemented for CSR and IJV format.
* @warning Not implemented for IJV format.
*
*******************************************************************************
*
......@@ -565,25 +566,54 @@ spmSort( spmatrix_t *spm )
spm_int_t
spmMergeDuplicate( spmatrix_t *spm )
{
spm_int_t local, global;
switch (spm->flttype) {
case SpmPattern:
return p_spmMergeDuplicate( spm );
local = p_spmMergeDuplicate( spm );
break;
case SpmFloat:
return s_spmMergeDuplicate( spm );
local = s_spmMergeDuplicate( spm );
break;
case SpmDouble:
return d_spmMergeDuplicate( spm );
local = d_spmMergeDuplicate( spm );
break;
case SpmComplex32:
return c_spmMergeDuplicate( spm );
local = c_spmMergeDuplicate( spm );
break;
case SpmComplex64:
return z_spmMergeDuplicate( spm );
local = z_spmMergeDuplicate( spm );
break;
default:
return SPM_ERR_BADPARAMETER;
return (spm_int_t)SPM_ERR_BADPARAMETER;
}
#if defined(SPM_WITH_MPI)
if ( spm->loc2glob ) {
MPI_Allreduce( &local, &global, 1, SPM_MPI_INT, MPI_SUM, spm->comm );
/* Update computed fields */
if( global > 0 ) {
MPI_Allreduce( &(spm->nnz), &(spm->gnnz), 1, SPM_MPI_INT, MPI_SUM, spm->comm );
MPI_Allreduce( &(spm->nnzexp), &(spm->gnnzexp), 1, SPM_MPI_INT, MPI_SUM, spm->comm );
}
}
else
#endif
{
global = local;
if ( global > 0 ) {
spm->gnnz = spm->nnz;
spm->gnnzexp = spm->nnzexp;
}
}
return global;
}
/**
......
src/z_spm.c deleted 100644 → 0
+ 0
− 318
View file @ 87038c1d
/**
*
* @file z_spm.c
*
* SParse Matrix package precision dependent routines.
*
* @copyright 2016-2020 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
* Univ. Bordeaux. All rights reserved.
*
* @version 1.0.0
* @author Mathieu Faverge
* @date 2019-10-29
*
* @precisions normal z -> c d s p
*
**/
#include "common.h"
#include "z_spm.h"
#include <string.h>
/**
*******************************************************************************
*
* @ingroup spm_dev_check
*
* @brief This routine merge the multiple entries in a sparse
* matrix by suming their values together.
*
* The sparse matrix needs to be sorted first (see z_spmSort()).
*
* @warning This function should NOT be called if dof is greater than 1.
*
*******************************************************************************
*
* @param[inout] spm
* On entry, the pointer to the sparse matrix structure.
* On exit, the reduced sparse matrix of multiple entries were present
* in it. The multiple values for a same vertex are sum up together.
*
********************************************************************************
*
* @return The number of vertices that were merged. -1 on error.
*
*******************************************************************************/
spm_int_t
z_spmMergeDuplicate( spmatrix_t *spm )
{
spm_int_t *colptr = spm->colptr;
spm_int_t *oldrow = spm->rowptr;
spm_int_t *newrow = spm->rowptr;
spm_complex64_t *newval = spm->values;
spm_complex64_t *oldval = spm->values;
spm_int_t n = spm->n;
spm_int_t baseval = spm->colptr[0];
spm_int_t dof2 = spm->dof * spm->dof;
spm_int_t idx, i, j, size, savedcolptr;
spm_int_t merge = 0;
#if !defined(PRECISION_p)
spm_int_t d;
#endif
assert( spm->dof >= 1 );
assert( spm->fmttype == SpmCSC );
if ( spm->fmttype == SpmCSC ) {
idx = baseval;
savedcolptr = colptr[0];
for (i=0; i<n; i++, colptr++)
{
size = colptr[1] - savedcolptr;
savedcolptr = colptr[1];
for (j = 0; j < size;
j++, oldrow++, oldval+=dof2, newrow++, newval+=dof2, idx++)
{
if ( newrow != oldrow ) {
newrow[0] = oldrow[0];
#if !defined(PRECISION_p)
memcpy( newval, oldval, dof2 * sizeof(spm_complex64_t) );
#endif
}
while( ((j+1) < size) && (newrow[0] == oldrow[1]) ) {
j++;
oldrow++;
oldval += dof2;
#if !defined(PRECISION_p)
/* Merge the two sets of values */
for (d = 0; d < dof2; d++) {
newval[d] += oldval[d];
}
#endif
merge++;
}
}
assert( ( (merge == 0) && (colptr[1] == idx) ) ||
( (merge != 0) && (colptr[1] > idx) ) );
colptr[1] = idx;
}
assert( ((merge == 0) && (spm->nnz == (idx-baseval))) ||
((merge != 0) && (spm->nnz - merge == (idx-baseval))) );
if (merge > 0) {
spm->nnz = spm->nnz - merge;
spm->gnnz = spm->nnz;
newrow = malloc( spm->nnz * sizeof(spm_int_t) );
memcpy( newrow, spm->rowptr, spm->nnz * sizeof(spm_int_t) );
free(spm->rowptr);
spm->rowptr = newrow;
#if !defined(PRECISION_p)
newval = malloc( spm->nnz * dof2 * sizeof(spm_complex64_t) );
memcpy( newval, spm->values, spm->nnz * dof2 * sizeof(spm_complex64_t) );
free(spm->values);
spm->values = newval;
#endif
}
}
return merge;
}
/**
*******************************************************************************
*
* @ingroup spm_dev_check
*
* @brief This routine corrects the sparse matrix structure if it's
* pattern is not symmetric.
*
* It returns the new symmetric pattern with zeroes on
* the new entries.
*
*******************************************************************************
*
* @param[inout] spm
* On entry, the pointer to the sparse matrix structure.
* On exit, the same sparse matrix with extra entries that makes it
* pattern symmetric.
*
*******************************************************************************
*
* @retval Return the number of entries added to the matrix.
*
*******************************************************************************/
spm_int_t
z_spmSymmetrize( spmatrix_t *spm )
{
spm_complex64_t *oldval, *valtmp, *newval = NULL;
spm_int_t *oldcol, *coltmp, *newcol = NULL;
spm_int_t *oldrow, *rowtmp, *newrow = NULL;
spm_int_t *toaddtab = NULL;
spm_int_t *toaddtmp, toaddcnt, toaddsze;
spm_int_t n = spm->n;
spm_int_t dof2 = spm->dof * spm->dof;
spm_int_t i, j, k, r, size;
spm_int_t baseval;
toaddcnt = 0;
toaddsze = 0;
if ( (spm->fmttype == SpmCSC) || (spm->fmttype == SpmCSR) ) {
if (spm->fmttype == SpmCSC) {
oldcol = spm->colptr;
coltmp = spm->colptr;
oldrow = spm->rowptr;
rowtmp = spm->rowptr;
}
else {
oldcol = spm->rowptr;
coltmp = spm->rowptr;
oldrow = spm->colptr;
rowtmp = spm->colptr;
}
baseval = oldcol[0];
for (i=0; i<n; i++, coltmp++)
{
size = coltmp[1] - coltmp[0];
for (r=0; r<size; r++, rowtmp++ )
{
j = rowtmp[0]-baseval;
if ( i != j ) {
/* Look for the element (j, i) */
spm_int_t frow = oldcol[ j ] - baseval;
spm_int_t lrow = oldcol[ j+1 ] - baseval;
int found = 0;
for (k = frow; k < lrow; k++)
{
if (i == (oldrow[k]-baseval))
{
found = 1;
break;
}
else if ( i < (oldrow[k]-baseval))
{
/* The spm is sorted so we will not find it later */
break;
}
}
if ( !found ) {
if ( newcol == NULL ) {
newcol = malloc( (spm->n+1) * sizeof(spm_int_t) );
for (k=0; k<spm->n; k++) {
newcol[k] = oldcol[k+1] - oldcol[k];
}
/* Let's start with a buffer that will contain 5% of extra elements */
toaddsze = spm_imax( 0.05 * (double)spm->nnz, 1 );
toaddtab = malloc( 2*toaddsze * sizeof(spm_int_t) );
}
if (toaddcnt >= toaddsze) {
toaddsze *= 2;
toaddtab = (spm_int_t*)realloc( toaddtab, 2*toaddsze*sizeof(spm_int_t) );
}
/* Newcol stores the number of elements per column */
newcol[ j ]++;
/* toaddtab stores the couples (j, i) to be added in the final sparse matrix */
toaddtab[ toaddcnt * 2 ] = j;
toaddtab[ toaddcnt * 2 + 1 ] = i;
toaddcnt++;
}
}
}
}
if (toaddcnt > 0) {
spm_int_t newsize, oldsize;
/* Sort the array per column */
spmIntSort2Asc1(toaddtab, toaddcnt);
spm->nnz = spm->nnz + toaddcnt;
spm->gnnz = spm->nnz;
newrow = malloc( spm->nnz * sizeof(spm_int_t) );
#if !defined(PRECISION_p)
newval = malloc( spm->nnz * dof2 * sizeof(spm_complex64_t) );
#endif
coltmp = newcol;
rowtmp = newrow;
valtmp = newval;
oldval = spm->values;
toaddtmp = toaddtab;
newsize = coltmp[0];
coltmp[0] = baseval;
for (i=0; i<n; i++, coltmp++, oldcol++)
{
/* Copy the existing elements */
oldsize = oldcol[1] - oldcol[0];
memcpy( rowtmp, oldrow, oldsize * sizeof(spm_int_t) );
#if !defined(PRECISION_p)
memcpy( valtmp, oldval, oldsize * dof2 * sizeof(spm_complex64_t) );
#endif
oldrow += oldsize;
rowtmp += oldsize;
oldval += oldsize * dof2;
valtmp += oldsize * dof2;
/* Some elements have been added */
assert( newsize >= oldsize );
if ( newsize > oldsize ) {
int added = 0;
int tobeadded = newsize - oldsize;
/* At least one element is equal to i */
assert( toaddtmp[0] == i );
/* Let's add the new vertices */
while( (added < tobeadded) && (toaddtmp[0] == i) )
{
rowtmp[0] = toaddtmp[1] + baseval;
rowtmp++;
toaddtmp+=2;
added++;
}
assert( added == tobeadded );
#if !defined(PRECISION_p)
/* Add 0.s for the associated values */
memset( valtmp, 0, added * dof2 * sizeof(spm_complex64_t) );
valtmp += added * dof2;
#endif
}
/* Use oldsize as temporary variable to update the new colptr */
oldsize = newsize;
newsize = coltmp[1];
coltmp[1] = coltmp[0] + oldsize;
}
assert( coltmp[0]-baseval == spm->nnz );
free( toaddtab );
free( spm->colptr );
free( spm->rowptr );
free( spm->values );
if (spm->fmttype == SpmCSC) {
spm->colptr = newcol;
spm->rowptr = newrow;
}
else {
spm->colptr = newrow;
spm->rowptr = newcol;
}
spm->values = newval;
}
}
return toaddcnt;
}
src/z_spm_mergeduplicate.c 0 → 100644
+ 148
− 0
View file @ a47eaafd
/**
*
* @file z_spm_mergeduplicate.c
*
* SParse Matrix package precision dependent routines.
*
* @copyright 2016-2020 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
* Univ. Bordeaux. All rights reserved.
*
* @version 1.0.0
* @author Mathieu Faverge
* @author Tony Delarue
* @date 2019-10-29
*
* @precisions normal z -> c d s p
*
**/
#include "common.h"
#include "z_spm.h"
#include <string.h>
/**
*******************************************************************************
*
* @ingroup spm_dev_check
*
* @brief This routine merge the multiple entries in a sparse matrix by summing
* their values together.
*
* The sparse matrix needs to be sorted first (see z_spmSort()). In distributed,
* only local entries are merged together.
*
*******************************************************************************
*
* @param[inout] spm
* On entry, the pointer to the sparse matrix structure.
* On exit, the reducton of the input sparse matrix where multiple
* occurences of a same element are summed up together.
*
********************************************************************************
*
* @return The number of vertices that were merged. -1 on error.
*
*******************************************************************************/
spm_int_t
z_spmMergeDuplicate( spmatrix_t *spm )
{
spm_int_t *colptr = (spm->fmttype == SpmCSC) ? spm->colptr : spm->rowptr;
spm_int_t *oldrow = (spm->fmttype == SpmCSC) ? spm->rowptr : spm->colptr;
spm_int_t *newrow = oldrow;
spm_complex64_t *newval = spm->values;
spm_complex64_t *oldval = spm->values;
spm_int_t *loc2glob = spm->loc2glob;
spm_int_t merge = 0;
spm_int_t n = spm->n;
spm_int_t baseval = spmFindBase( spm );
spm_int_t ig, jl, jg, dofi, dofj, dof2;
spm_int_t k, idx, size, valsize, savedcolptr;
#if !defined(PRECISION_p)
spm_int_t d;
#endif
if ( (spm->fmttype != SpmCSC) &&
(spm->fmttype != SpmCSR) )
{
fprintf(stderr, "Error : MergeDuplicate can only be called with SpmCSC or SpmCSR\n");
return SPM_ERR_BADPARAMETER;
}
idx = baseval;
valsize = 0;
savedcolptr = colptr[0];
for (jl=0; jl<n; jl++, colptr++, loc2glob++)
{
jg = (spm->loc2glob == NULL) ? jl : *loc2glob - baseval;
dofj = (spm->dof > 0) ? spm->dof : spm->dofs[jg+1] - spm->dofs[jg];
size = colptr[1] - savedcolptr;
savedcolptr = colptr[1];
for ( k=0; k<size; k++, idx++ )
{
ig = *newrow - baseval;
dofi = (spm->dof > 0) ? spm->dof : spm->dofs[ig+1] - spm->dofs[ig];
dof2 = dofi * dofj;
valsize += dof2;
/*
* A shift has been introduced, we need to first compact the structure
*/
if ( newrow != oldrow ) {
newrow[0] = oldrow[0];
#if !defined(PRECISION_p)
memcpy( newval, oldval, dof2 * sizeof(spm_complex64_t) );
#endif
}
/*
* Let's sum together all identical elements
*/
while( ((k+1) < size) && (newrow[0] == oldrow[1]) ) {
k++;
oldrow++;
oldval += dof2;
#if !defined(PRECISION_p)
/* Merge the two sets of values */
for ( d=0; d<dof2; d++ ) {
newval[d] += oldval[d];
}
#endif
merge++;
}
/* Shift arrays */
oldrow++;
newrow++;
oldval += dof2;
newval += dof2;
}
assert( ( (merge == 0) && (colptr[1] == idx) ) ||
( (merge != 0) && (colptr[1] > idx) ) );
colptr[1] = idx;
}
assert( ((merge == 0) && (spm->nnz == (idx-baseval))) ||
((merge != 0) && (spm->nnz - merge == (idx-baseval))) );
/*
* Realloc the arrays if they have been compacted
*/
if ( merge > 0 ) {
spm->nnz = spm->nnz - merge;
spm->nnzexp = valsize;
if ( spm->fmttype == SpmCSC ) {
spm->rowptr = realloc( spm->rowptr, spm->nnz * sizeof( spm_int_t ) );
}
else {
spm->colptr = realloc( spm->colptr, spm->nnz * sizeof( spm_int_t ) );
}
#if !defined(PRECISION_p)
spm->values = realloc( spm->values, valsize * sizeof( spm_complex64_t ) );
#endif
}
return merge;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment