diff --git a/CMakeLists.txt b/CMakeLists.txt
index 46098c79fce213f7c4ed876a63389b6d7acb9dc1..66834ebfca0c93922eabd7a02a566af101523956 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -207,6 +207,7 @@ set(SOURCES
src/z_spm_laplacian.c
src/z_spm_matrixvector.c
src/z_spm_print.c
+ src/z_spm_sort.c
)
precisions_rules_py(generated_sources
diff --git a/src/common.h b/src/common.h
index 814be961faa2ff87026dd56ebe3be0b8f78ccb9d..1210be7be3f7811e23147be9f8d75ad9a4214028 100644
--- a/src/common.h
+++ b/src/common.h
@@ -55,6 +55,7 @@ spm_get_datatype( const spmatrix_t *spm )
spm_int_t *spm_get_glob2loc( spmatrix_t *spm, spm_int_t baseval );
int spm_get_distribution( const spmatrix_t *spm );
+spm_int_t *spm_create_asc_values( const spmatrix_t *spm );
/********************************************************************
* Conjuguate/Id functions
diff --git a/src/spm.c b/src/spm.c
index a27458ce4b4e37df9144f9fff994f24871c2e548..82a7fff45721e5b5dddfdf1e778f209874bd7337 100644
--- a/src/spm.c
+++ b/src/spm.c
@@ -558,44 +558,24 @@ spmNorm( spm_normtype_t ntype,
int
spmSort( spmatrix_t *spm )
{
- if ( (spm->dof != 1) && (spm->flttype != SpmPattern) ) {
- switch (spm->flttype) {
- case SpmFloat:
- s_spmSortMultidof( spm );
- break;
- case SpmDouble:
- d_spmSortMultidof( spm );
- break;
- case SpmComplex32:
- c_spmSortMultidof( spm );
- break;
- case SpmComplex64:
- z_spmSortMultidof( spm );
- break;
- default:
- return SPM_ERR_BADPARAMETER;
- }
- }
- else {
- switch (spm->flttype) {
- case SpmPattern:
- p_spmSort( spm );
- break;
- case SpmFloat:
- s_spmSort( spm );
- break;
- case SpmDouble:
- d_spmSort( spm );
- break;
- case SpmComplex32:
- c_spmSort( spm );
- break;
- case SpmComplex64:
- z_spmSort( spm );
- break;
- default:
- return SPM_ERR_BADPARAMETER;
- }
+ switch (spm->flttype) {
+ case SpmPattern:
+ p_spmSort( spm );
+ break;
+ case SpmFloat:
+ s_spmSort( spm );
+ break;
+ case SpmDouble:
+ d_spmSort( spm );
+ break;
+ case SpmComplex32:
+ c_spmSort( spm );
+ break;
+ case SpmComplex64:
+ z_spmSort( spm );
+ break;
+ default:
+ return SPM_ERR_BADPARAMETER;
}
return SPM_SUCCESS;
}
@@ -1699,3 +1679,82 @@ spm_get_distribution( const spmatrix_t *spm )
assert(distribution > 0);
return distribution;
}
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup spm_dev_check
+ *
+ * @brief Create an nnz array that represents the shift of the original
+ * multidof value array.
+ *
+ *******************************************************************************
+ *
+ * @param[in] spm
+ * The sparse matrix structure.
+ *
+ ********************************************************************************
+ *
+ * @return An nnz array which stores the multidof shift of the original
+ * values aray
+ *
+ *******************************************************************************/
+spm_int_t *
+spm_create_asc_values( const spmatrix_t *spm )
+{
+ spm_int_t i, j, ig, jg;
+ spm_int_t baseval, n;
+ spm_int_t dof, dofi, dofj;
+ spm_int_t *colptr = spm->colptr;
+ spm_int_t *rowptr = spm->rowptr;
+ spm_int_t *dofs = spm->dofs;
+ spm_int_t *loc2glob = spm->loc2glob;
+ spm_int_t *values = malloc( (spm->nnz + 1) * sizeof(spm_int_t));
+ spm_int_t *valtmp = values;
+
+ values[0] = 0;
+ baseval = spmFindBase(spm);
+ dof = spm->dof;
+ switch (spm->fmttype)
+ {
+ case SpmCSR:
+ colptr = spm->rowptr;
+ rowptr = spm->colptr;
+
+ spm_attr_fallthrough;
+
+ case SpmCSC:
+ n = spm->n;
+ loc2glob = spm->loc2glob;
+ for ( j = 0; j < n; j++, colptr++, loc2glob++ )
+ {
+ jg = (spm->loc2glob == NULL) ? j : *loc2glob - baseval;
+ dofj = (dof > 0) ? dof : dofs[jg+1] - dofs[jg];
+ for ( i = colptr[0]; i < colptr[1]; i++, rowptr++, valtmp++ )
+ {
+ ig = *rowptr - baseval;
+ dofi = (dof > 0) ? dof : dofs[ig+1] - dofs[ig];
+
+ valtmp[1] = valtmp[0] + (dofj*dofi);
+ }
+ }
+ break;
+
+ case SpmIJV:
+ n = spm->nnz;
+ for ( j = 0; j < n; j++, colptr++, rowptr++, valtmp++ )
+ {
+ jg = *colptr - baseval;
+ dofj = (dof > 0) ? dof : dofs[jg+1] - dofs[jg];
+ ig = *rowptr - baseval;
+ dofi = (dof > 0) ? dof : dofs[ig+1] - dofs[ig];
+
+ valtmp[1] = valtmp[0] + (dofj*dofi);
+ }
+ break;
+ }
+ assert((valtmp - values) == spm->nnz);
+ values = realloc( values, spm->nnz * sizeof(spm_int_t) );
+
+ return values;
+}
diff --git a/src/z_spm.c b/src/z_spm.c
index 366a84c528d0e7d600c79cad7464a78abd86ff96..273593dbd204562b12dd4fd34496626f2b2b2fe6 100644
--- a/src/z_spm.c
+++ b/src/z_spm.c
@@ -18,372 +18,6 @@
#include "z_spm.h"
#include <string.h>
-/**
- *******************************************************************************
- *
- * @ingroup spm_dev_check
- *
- * @brief This routine sorts the spm matrix.
- *
- * For the CSC and CSR formats, the subarray of edges for each vertex are sorted.
- * For the IJV format, the edges are storted first by column indexes, and then
- * by row indexes. To perform a sort first by row, second by column, please swap
- * the colptr and rowptr of the structure before calling the subroutine.
- *
- * @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 same sparse matrix with subarrays of edges sorted by
- * ascending order.
- *
- *******************************************************************************/
-void
-z_spmSort( spmatrix_t *spm )
-{
- spm_int_t *colptr = spm->colptr;
- spm_int_t *rowptr = spm->rowptr;
- spm_complex64_t *values = spm->values;
- void *sortptr[3];
- spm_int_t n = spm->n;
- spm_int_t i, size;
- (void)sortptr;
-
-#if !defined(PRECISION_p)
- assert( spm->dof == 1 );
-#endif
-
- /* Sort in place each subset */
- if ( spm->fmttype == SpmCSC ) {
- for (i=0; i<n; i++, colptr++)
- {
- size = colptr[1] - colptr[0];
-
-#if defined(PRECISION_p)
- spmIntSort1Asc1( rowptr, size );
-#else
- sortptr[0] = rowptr;
- sortptr[1] = values;
- z_spmIntFltSortAsc( sortptr, size );
-#endif
- rowptr += size;
- values += size;
- }
- }
- else if ( spm->fmttype == SpmCSR ) {
- for (i=0; i<n; i++, rowptr++)
- {
- size = rowptr[1] - rowptr[0];
-
-#if defined(PRECISION_p)
- spmIntSort1Asc1( colptr, size );
-#else
- sortptr[0] = colptr;
- sortptr[1] = values;
- z_spmIntFltSortAsc( sortptr, size );
-#endif
- colptr += size;
- values += size;
- }
- }
- else if ( spm->fmttype == SpmIJV ) {
- size = spm->nnz;
-
- sortptr[0] = colptr;
- sortptr[1] = rowptr;
-
-#if defined(PRECISION_p)
- spmIntMSortIntAsc( sortptr, size );
-#else
- sortptr[2] = values;
- z_spmIntIntFltSortAsc( sortptr, size );
-#endif
- }
-}
-
-/**
- * @brief Sort subarrays of rowptr and values
- *
- * @param[in] rowptr
- * The original rowptr.
- *
- * @param[in] rowtmp
- * The sorted copy of the rowptr.
- *
- * @param[in] values
- * The original values array.
- *
- * @param[inout] valtmp
- * The copy of the valptr to sort.
- *
- * @param[in] dofs
- * The pointer to the dofs array.
- *
- * @param[in] dof
- * SPM dof value.
- *
- * @param[in] baseval
- * SPM baseval.
- *
- * @param[in] dofj
- * Current colum dof.
- *
- * @param[in] size
- * Size of the current subarray.
- *
- * @return number of value sorted in the value array
- */
-static inline spm_int_t
-z_spm_sort_multidof_csx_values( const spm_int_t *rowptr,
- const spm_int_t *rowtmp,
- const spm_complex64_t *values,
- spm_complex64_t *valtmp,
- const spm_int_t *dofs,
- spm_int_t dof,
- spm_int_t baseval,
- spm_int_t dofj,
- spm_int_t size )
-{
- spm_int_t i, ig, dofi;
- spm_int_t k = 0;
- spm_int_t memory, count, added = 0;
-
- while (k < size)
- {
- memory = 0;
- while ( (k < (size - 1)) && (rowtmp[k] == rowtmp[k+1]) )
- {
- memory++;
- k++;
- }
-
- count = 0;
- for ( i = 0; i < size; i++)
- {
- ig = rowptr[i] - baseval;
- dofi = (dof > 0) ? dof : dofs[ig+1] - dofs[ig];
- if ( rowtmp[k] != rowptr[i] ) {
- count += dofj * dofi;
- continue;
- }
- /*
- * The matrix isn't merged.
- * We have to make sure that we don't copy the same information.
- */
- memcpy( valtmp + added,
- values + count,
- dofi * dofj * sizeof(spm_complex64_t) );
- added += dofi * dofj;
-
- if ( memory > 0 ) {
- memory--;
- continue;
- }
-
- k++;
- break;
- }
- }
- return added;
-}
-
-/**
- * @brief Sort a IJV matrix.
- *
- * @param[in] spm
- * Pointer to the spm structure.
- *
- * @param[inout] newcol
- * The sorted copy of the colptr.
- *
- * @param[inout] newrow
- * The sorted copy of the rowptr.
- *
- * @param[inout] newval
- * The copy of the valptr to sort.
- */
-static inline void
-z_spm_sort_multidof_ijv_values( const spmatrix_t *spm,
- spm_int_t *newcol,
- spm_int_t *newrow,
- spm_complex64_t *newval )
-{
- spm_int_t *colptr;
- spm_int_t *rowptr;
- spm_complex64_t *values;
- spm_int_t *dofs;
- spm_int_t i, ig, jg, dofi, dofj, dof2;
- spm_int_t size, baseval;
- spm_int_t k = 0;
- spm_int_t count, memory = 0;
-
- values = spm->values;
- dofs = spm->dofs;
- size = spm->nnz;
- baseval = spmFindBase(spm);
- while (k < size)
- {
- while ( (newcol[0] == newcol[1]) && (newrow[0] == newrow[1]) )
- {
- newcol++;
- newrow++;
- memory++;
- k++;
- }
-
- jg = *newcol - baseval;
- dofj = (spm->dof > 0) ? spm->dof : dofs[jg+1] - dofs[jg];
- ig = *newrow - baseval;
- dofi = (spm->dof > 0) ? spm->dof : dofs[ig+1] - dofs[ig];
- dof2 = dofi * dofj;
-
- count = 0;
- colptr = spm->colptr;
- rowptr = spm->rowptr;
- for ( i = 0; i < size; i++, colptr++, rowptr++ )
- {
- jg = *colptr - baseval;
- dofj = (spm->dof > 0) ? spm->dof : dofs[jg+1] - dofs[jg];
- ig = *rowptr - baseval;
- dofi = (spm->dof > 0) ? spm->dof : dofs[ig+1] - dofs[ig];
-
- if ( ((*newcol) != (*colptr)) || ((*newrow) != (*rowptr)) ) {
- count += dofj * dofi;
- continue;
- }
-
- memcpy( newval,
- values + count,
- dof2 * sizeof(spm_complex64_t) );
- newval += dof2;
-
- if( memory > 0 ) {
- memory--;
- continue;
- }
-
- newcol++;
- newrow++;
- k++;
- break;
- }
- assert(memory == 0);
- }
-}
-
-/**
- *******************************************************************************
- *
- * @ingroup spm_dev_check
- *
- * @brief This routine sorts the spm matrix.
- *
- * For the CSC and CSR formats, the subarray of edges for each vertex are sorted.
- * For the IJV format, the edges are sorted first by column indexes, and then
- * by row indexes. To perform a sort first by row, second by column, please swap
- * the colptr and rowptr of the structure before calling the subroutine.
- * This routine is used for multidof matrices. It's way less efficient than the
- * single dof one.
- *
- *******************************************************************************
- *
- * @param[inout] spm
- * On entry, the pointer to the sparse matrix structure.
- * On exit, the same sparse matrix with subarrays of edges sorted by
- * ascending order.
- *
- *******************************************************************************/
-void
-z_spmSortMultidof( spmatrix_t *spm )
-{
- spm_int_t *colptr, *newcol, *coltmp;
- spm_int_t *rowptr, *newrow, *rowtmp;
- spm_complex64_t *values, *newval, *valtmp;
- spm_int_t size, n = spm->n;
-
- newrow = malloc( spm->nnz * sizeof(spm_int_t) );
- newval = malloc( spm->nnzexp * sizeof(spm_complex64_t) );
- values = spm->values;
-
- if ( spm->fmttype != SpmIJV ) {
- spm_int_t *loc2glob = spm->loc2glob;
- spm_int_t *dofs = spm->dofs;
- spm_int_t j, jg, dofj, baseval;
- spm_int_t added;
-
- baseval = spmFindBase(spm);
- rowtmp = newrow;
- valtmp = newval;
- colptr = (spm->fmttype == SpmCSC) ? spm->colptr : spm->rowptr;
- rowptr = (spm->fmttype == SpmCSC) ? spm->rowptr : spm->colptr;
-
- memcpy( newrow, rowptr, spm->nnz * sizeof(spm_int_t) );
- for (j=0; j<n; j++, colptr++, loc2glob++)
- {
- size = colptr[1] - colptr[0];
- jg = (spm->loc2glob == NULL) ? j : *loc2glob - baseval;
- dofj = (spm->dof > 0) ? spm->dof : dofs[jg+1] - dofs[jg];
-
- /* Sort rowptr */
- spmIntSort1Asc1( rowtmp, size );
-
- /* Sort values */
- added = z_spm_sort_multidof_csx_values( rowptr, rowtmp, values, valtmp, dofs,
- spm->dof, baseval, dofj, size );
-
- rowptr += size;
- rowtmp += size;
- values += added;
- valtmp += added;
- }
-
- if(spm->fmttype == SpmCSC) {
- memcpy( spm->rowptr, newrow, spm->nnz * sizeof( spm_int_t ) );
- }
- else {
- memcpy( spm->colptr, newrow, spm->nnz * sizeof( spm_int_t ) );
- }
- }
-
- else {
- void *sortptr[2];
-
- colptr = spm->colptr;
- rowptr = spm->rowptr;
- size = spm->nnz;
- newcol = malloc( size * sizeof(spm_int_t) );
-
- memcpy( newcol, colptr, size * sizeof(spm_int_t) );
- memcpy( newrow, rowptr, size * sizeof(spm_int_t) );
-
- sortptr[0] = newcol;
- sortptr[1] = newrow;
-
- /* Sort the colptr and the rowptr */
- spmIntMSortIntAsc( sortptr, size );
-
- coltmp = newcol;
- rowtmp = newrow;
- valtmp = newval;
-
- /* Sort values */
- z_spm_sort_multidof_ijv_values( spm, coltmp, rowtmp, valtmp );
-
- memcpy(spm->colptr, newcol, spm->nnz * sizeof( spm_int_t ));
- memcpy(spm->rowptr, newrow, spm->nnz * sizeof( spm_int_t ));
-
- free(newcol);
- }
-
- memcpy( spm->values, newval, spm->nnzexp * sizeof( spm_complex64_t ) );
-
- free(newrow);
- free(newval);
-}
-
-
/**
*******************************************************************************
*
diff --git a/src/z_spm.h b/src/z_spm.h
index f2758d4e90a218d803d5a60c0b2d9ee646309498..7cd4c301a664bf653e7ff2f300ed1b898143e94e 100644
--- a/src/z_spm.h
+++ b/src/z_spm.h
@@ -70,7 +70,6 @@ double z_spmNorm( spm_normtype_t ntype, const spmatrix_t *spm );
* Extra routines
*/
void z_spmSort( spmatrix_t *spm );
-void z_spmSortMultidof( spmatrix_t *spm );
spm_int_t z_spmMergeDuplicate( spmatrix_t *spm );
spm_int_t z_spmSymmetrize( spmatrix_t *spm );
diff --git a/src/z_spm_sort.c b/src/z_spm_sort.c
new file mode 100644
index 0000000000000000000000000000000000000000..725a45acd181fe60b2d102f2b3ff37eafab8aeb7
--- /dev/null
+++ b/src/z_spm_sort.c
@@ -0,0 +1,300 @@
+/**
+ *
+ * @file z_spm_sort.c
+ *
+ * SParse Matrix package precision dependent sort 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 2020-09-24
+ *
+ * @precisions normal z -> c d s p
+ *
+ **/
+#include "common.h"
+#include "z_spm.h"
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup spm_dev_check
+ *
+ * @brief This routine sorts the single dof spm matrix.
+ *
+ * For the CSC and CSR formats, the subarray of edges for each vertex are sorted.
+ * For the IJV format, the edges are storted first by column indexes, and then
+ * by row indexes. To perform a sort first by row, second by column, please swap
+ * the colptr and rowptr of the structure before calling the subroutine.
+ *
+ * @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 same sparse matrix with subarrays of edges sorted by
+ * ascending order.
+ *
+ *******************************************************************************/
+static inline void
+z_spmSortNoDof( spmatrix_t *spm )
+{
+ spm_int_t *colptr = spm->colptr;
+ spm_int_t *rowptr = spm->rowptr;
+ spm_complex64_t *values = spm->values;
+ void *sortptr[3];
+ spm_int_t n = spm->n;
+ spm_int_t i, size;
+ (void)sortptr;
+
+ /* Sort in place each subset */
+ if ( spm->fmttype == SpmCSC ) {
+ for (i=0; i<n; i++, colptr++)
+ {
+ size = colptr[1] - colptr[0];
+
+#if defined(PRECISION_p)
+ spmIntSort1Asc1( rowptr, size );
+#else
+ sortptr[0] = rowptr;
+ sortptr[1] = values;
+ z_spmIntFltSortAsc( sortptr, size );
+#endif
+ rowptr += size;
+ values += size;
+ }
+ }
+ else if ( spm->fmttype == SpmCSR ) {
+ for (i=0; i<n; i++, rowptr++)
+ {
+ size = rowptr[1] - rowptr[0];
+
+#if defined(PRECISION_p)
+ spmIntSort1Asc1( colptr, size );
+#else
+ sortptr[0] = colptr;
+ sortptr[1] = values;
+ z_spmIntFltSortAsc( sortptr, size );
+#endif
+ colptr += size;
+ values += size;
+ }
+ }
+ else if ( spm->fmttype == SpmIJV ) {
+ size = spm->nnz;
+
+ sortptr[0] = colptr;
+ sortptr[1] = rowptr;
+
+#if defined(PRECISION_p)
+ spmIntMSortIntAsc( sortptr, size );
+#else
+ sortptr[2] = values;
+ z_spmIntIntFltSortAsc( sortptr, size );
+#endif
+ }
+}
+
+/**
+ * @brief Apply the permutations on the values array for a CSX spm
+ * The values array holds the permutation indexes of the new values array.
+ *
+ * @param[in] spm
+ * The pointer to the spm.
+ *
+ * @param[in] values
+ * The original values array.
+ *
+ * @param[inout] newval
+ * The new values array with the correct permutations.
+ * Must be allocated before this routine.
+ */
+static inline void
+z_spm_sort_multidof_csx_values( const spmatrix_t *spm,
+ const spm_complex64_t *values,
+ spm_complex64_t *newval )
+{
+ spm_int_t i, j, ig, jg, index;
+ spm_int_t size, baseval, dof;
+ spm_int_t dofi, dofj, dof2;
+
+ spm_int_t *colptr = (spm->fmttype == SpmCSC) ? spm->colptr: spm->rowptr;
+ spm_int_t *rowptr = (spm->fmttype == SpmCSC) ? spm->rowptr: spm->colptr;
+ spm_int_t *indexes = spm->values;
+ spm_int_t *dofs = spm->dofs;
+ spm_int_t *loc2glob = spm->loc2glob;
+
+ spm_complex64_t *valtmp = newval;
+
+ size = spm->n;
+ baseval = spmFindBase(spm);
+ dof = spm->dof;
+ for ( j = 0; j < size; j++, colptr++, loc2glob++ )
+ {
+ jg = (spm->loc2glob == NULL) ? j : *loc2glob - baseval;
+ dofj = (dof > 0) ? dof : dofs[jg+1] - dofs[jg];
+
+ for ( i = colptr[0]; i < colptr[1]; i++, rowptr++, indexes++ )
+ {
+ ig = *rowptr - baseval;
+ dofi = (dof > 0) ? dof : dofs[ig+1] - dofs[ig];
+ dof2 = dofi * dofj;
+
+ index = *indexes;
+ memcpy( valtmp, values + index, dof2 * sizeof(spm_complex64_t) );
+ valtmp += dof2;
+ }
+ }
+}
+
+/**
+ * @brief Apply the permutations on the values array fon an IJV spm.
+ * The values array holds the permutation indexes of the new values array.
+ *
+ * @param[in] spm
+ * The pointer to the spm.
+ *
+ * @param[in] values
+ * The original values array.
+ *
+ * @param[inout] newval
+ * The new values array with the correct permutations.
+ * Must be allocated before this routine.
+ */
+static inline void
+z_spm_sort_multidof_ijv_values( const spmatrix_t *spm,
+ const spm_complex64_t *values,
+ spm_complex64_t *newval )
+{
+ spm_int_t i, ig, jg, index;
+ spm_int_t size, baseval, dof;
+ spm_int_t dofi, dofj, dof2;
+
+ spm_int_t *colptr = spm->colptr;
+ spm_int_t *rowptr = spm->rowptr;
+ spm_int_t *indexes = spm->values;
+ spm_int_t *dofs;
+
+ spm_complex64_t *valtmp = newval;
+
+ size = spm->nnz;
+ baseval = spmFindBase(spm);
+ dof = spm->dof;
+ dofs = spm->dofs - baseval;
+ for ( i = 0; i < size; i++, colptr++, rowptr++, indexes++ )
+ {
+ jg = *colptr;
+ dofj = (dof > 0) ? dof : dofs[jg+1] - dofs[jg];
+ ig = *rowptr;
+ dofi = (dof > 0) ? dof : dofs[ig+1] - dofs[ig];
+ dof2 = dofi * dofj;
+
+ index = *indexes;
+
+ memcpy( valtmp, values + index, dof2 * sizeof(spm_complex64_t) );
+ valtmp += dof2;
+ }
+}
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup spm_dev_check
+ *
+ * @brief This routine sorts the multiple dof spm matrix.
+ *
+ * For the CSC and CSR formats, the subarray of edges for each vertex are sorted.
+ * For the IJV format, the edges are sorted first by column indexes, and then
+ * by row indexes. To perform a sort first by row, second by column, please swap
+ * the colptr and rowptr of the structure before calling the subroutine.
+ * This routine is used for multidof matrices. It's way less efficient than the
+ * single dof one.
+ *
+ *******************************************************************************
+ *
+ * @param[inout] spm
+ * On entry, the pointer to the sparse matrix structure.
+ * On exit, the same sparse matrix with subarrays of edges sorted by
+ * ascending order.
+ *
+ *******************************************************************************/
+static inline void
+z_spmSortMultidof( spmatrix_t *spm )
+{
+ spm_int_t dof;
+ spm_coeftype_t flttype;
+ spm_complex64_t *values = spm->values;
+ /* This array will store the solution */
+ spm_complex64_t *newval = malloc( spm->nnzexp * sizeof(spm_complex64_t) );
+
+ /* Create a tmp array composed by the multidof indexes of the valptr */
+ spm_int_t *indexes = spm_create_asc_values(spm);
+
+ /*
+ * Sort the spm as a single dof matrix.
+ * The value array will represent the permutations.
+ */
+ dof = spm->dof;
+ flttype = spm->flttype;
+
+ spm->values = indexes;
+ spm->dof = 1;
+
+ if ( sizeof(spm_int_t) == sizeof(spm_fixdbl_t) ) {
+ spm->flttype = 3; // SpmDouble
+ }
+ else {
+ spm->flttype = 2; // SpmFloat
+ }
+ spmSort( spm );
+
+ spm->dof = dof;
+ spm->flttype = flttype;
+
+ /* Apply the permutations and copy datas in the newval */
+ if ( spm->fmttype != SpmIJV ) {
+ z_spm_sort_multidof_csx_values( spm, values, newval );
+ }
+ else {
+ z_spm_sort_multidof_ijv_values( spm, values, newval );
+ }
+ free(indexes);
+ free(values);
+
+ spm->values = newval;
+}
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup spm_dev_check
+ *
+ * @brief This routine sorts the spm matrix.
+ *
+ * For the CSC and CSR formats, the subarray of edges for each vertex are sorted.
+ * For the IJV format, the edges are storted first by column indexes, and then
+ * by row indexes. To perform a sort first by row, second by column, please swap
+ * the colptr and rowptr of the structure before calling the subroutine.
+ *
+ *******************************************************************************
+ *
+ * @param[inout] spm
+ * On entry, the pointer to the sparse matrix structure.
+ * On exit, the same sparse matrix with subarrays of edges sorted by
+ * ascending order.
+ *
+ *******************************************************************************/
+void
+z_spmSort( spmatrix_t *spm )
+{
+ if ( (spm->dof != 1) && (spm->flttype != SpmPattern) ) {
+ z_spmSortMultidof( spm );
+ }
+ else {
+ z_spmSortNoDof( spm );
+ }
+}
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 7257f221c866401a13394bedf1b904a00130821a..09b01b87e8f5b5c21b2d71e95dea094ff107fe4a 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -17,6 +17,7 @@ set(LIB_SOURCES
z_spm_tests.c
core_zgeadd.c
core_zplrnt.c
+ z_spm_sort_tests.c
)
## reset variables
@@ -55,6 +56,7 @@ if ( SPM_WITH_MPI )
spm_dist_norm_tests.c
spm_dist_genrhs_tests.c
spm_dist_matvec_tests.c
+ spm_dist_sort_tests.c
)
endif()
@@ -66,14 +68,22 @@ endforeach()
## CTest execution
set( SPM_TESTS
- spm_convert_tests spm_norm_tests spm_matvec_tests )
+ spm_convert_tests
+ spm_norm_tests
+ spm_matvec_tests
+ )
set( SPM_DOF_TESTS
- spm_dof_expand_tests spm_dof_norm_tests spm_dof_matvec_tests spm_dof_sort_tests)
+ spm_dof_expand_tests
+ spm_dof_norm_tests
+ spm_dof_matvec_tests
+ spm_dof_sort_tests
+ )
set( SPM_MPI_TESTS
spm_scatter_gather_tests
spm_dist_norm_tests
spm_dist_genrhs_tests
spm_dist_matvec_tests
+ spm_dist_sort_tests
)
# List of run types
diff --git a/tests/spm_dist_sort_tests.c b/tests/spm_dist_sort_tests.c
new file mode 100644
index 0000000000000000000000000000000000000000..f6de343602ed19b525cd34aee38c5f2dd43b8ccb
--- /dev/null
+++ b/tests/spm_dist_sort_tests.c
@@ -0,0 +1,192 @@
+/**
+ *
+ * @file spm_dist_sort_tests.c
+ *
+ * Tests and validate the spm_sort routines with a distributed spm.
+ *
+ * @copyright 2015-2020 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+ * Univ. Bordeaux. All rights reserved.
+ *
+ * @version 1.0.0
+ * @author Mathieu Faverge
+ * @author Delarue Tony
+ * @date 2020-09-07
+ *
+ **/
+#include <stdint.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <assert.h>
+#include <time.h>
+#include <spm_tests.h>
+
+static inline int
+spm_dist_sort_check( spmatrix_t *spm, spmatrix_t *spmdist )
+{
+ spmatrix_t expand1, expand2;
+ spmatrix_t *gathered;
+ int rc;
+
+ spmSort( spm );
+ spmSort( spmdist );
+
+ gathered = spmGather( spmdist, -1 );
+
+ spmExpand( spm, &expand1 );
+ spmExpand( gathered, &expand2 );
+
+ rc = spmCompare( &expand1, &expand2 );
+
+ MPI_Allreduce( MPI_IN_PLACE, &rc, 1, MPI_INT, MPI_SUM, spmdist->comm );
+
+ spmExit( &expand1 );
+ spmExit( &expand2 );
+ spmExit( gathered );
+ free( gathered );
+
+ return rc;
+}
+
+int main (int argc, char **argv)
+{
+ char *filename;
+ spmatrix_t original, *spm, *spmdist;
+ spm_driver_t driver;
+ int clustnbr = 1;
+ int clustnum = 0;
+ spm_mtxtype_t mtxtype;
+ spm_fmttype_t fmttype;
+ int baseval, distbycol;
+ int rc = SPM_SUCCESS;
+ int err = 0;
+ int dof, dofmax = 4;
+ int to_free = 0;
+
+ MPI_Init( &argc, &argv );
+
+ /**
+ * Get options from command line
+ */
+ spmGetOptions( argc, argv,
+ &driver, &filename );
+
+ rc = spmReadDriver( driver, filename, &original );
+ free(filename);
+
+ if ( rc != SPM_SUCCESS ) {
+ fprintf(stderr, "ERROR: Could not read the file, stop the test !!!\n");
+ return EXIT_FAILURE;
+ }
+
+ MPI_Comm_size( MPI_COMM_WORLD, &clustnbr );
+ MPI_Comm_rank( MPI_COMM_WORLD, &clustnum );
+
+ if ( original.flttype == SpmPattern ) {
+ spmGenFakeValues( &original );
+ }
+
+ spmPrintInfo( &original, stdout );
+
+ if ( clustnum == 0 ) {
+ printf(" -- SPM Sort Test --\n");
+ }
+
+ for( fmttype=SpmCSC; fmttype<=SpmIJV; fmttype++ ) {
+
+ distbycol = (fmttype == SpmCSR) ? 0 : 1;
+ if ( spmConvert( fmttype, &original ) != SPM_SUCCESS ) {
+ fprintf( stderr, "Issue to convert to %s format\n", fmtnames[fmttype] );
+ continue;
+ }
+
+ for( dof=-1; dof<2; dof++ )
+ {
+ if ( dof >= 0 ) {
+ spm = spmDofExtend( &original, dof, dofmax );
+ to_free = 1;
+ }
+ else {
+ spm = malloc(sizeof(spmatrix_t));
+ memcpy( spm, &original, sizeof(spmatrix_t) );
+ to_free = 0;
+ }
+
+ if ( spm == NULL ) {
+ fprintf( stderr, "Issue to extend the matrix\n" );
+ continue;
+ }
+
+ spmdist = spmScatter( spm, -1, NULL, distbycol, -1, MPI_COMM_WORLD );
+ if ( spmdist == NULL ) {
+ fprintf( stderr, "Failed to scatter the spm\n" );
+ err++;
+ continue;
+ }
+
+ for( baseval=0; baseval<2; baseval++ )
+ {
+ spmBase( spmdist, baseval );
+ spmBase( spm, baseval );
+
+ for( mtxtype=SpmGeneral; mtxtype<=SpmHermitian; mtxtype++ )
+ {
+ if ( (mtxtype == SpmHermitian) &&
+ ( ((original.flttype != SpmComplex64) &&
+ (original.flttype != SpmComplex32)) ||
+ (original.mtxtype != SpmHermitian) ) )
+ {
+ continue;
+ }
+
+ if ( (mtxtype != SpmGeneral) &&
+ (original.mtxtype == SpmGeneral) )
+ {
+ continue;
+ }
+
+ if ( spm ) {
+ spm->mtxtype = mtxtype;
+ }
+ spmdist->mtxtype = mtxtype;
+
+ if ( clustnum == 0 ) {
+ printf( " Case: %s / %s / %d / %s / %d\n",
+ fltnames[spmdist->flttype],
+ fmtnames[spmdist->fmttype], baseval,
+ mtxnames[mtxtype - SpmGeneral], dof );
+ }
+
+ rc = spm_dist_sort_check( spm, spmdist );
+
+ err = (rc != 0) ? err+1 : err;
+ }
+ }
+ spmExit( spmdist );
+ free( spmdist );
+
+ if ( spm != &original ) {
+ if( to_free ){
+ spmExit( spm );
+ }
+ free( spm );
+ }
+ }
+ }
+
+ spmExit( &original );
+ MPI_Finalize();
+
+ if( err == 0 ) {
+ if(clustnum == 0) {
+ printf(" -- All tests PASSED --\n");
+ }
+ return EXIT_SUCCESS;
+ }
+ else
+ {
+ printf(" -- %d tests FAILED --\n", err);
+ return EXIT_FAILURE;
+ }
+}
diff --git a/tests/spm_dof_sort_tests.c b/tests/spm_dof_sort_tests.c
index 1921bae614cacdd1493b0880460b158600479ac7..146be200d0be9cb612421e0cfe95d156e5520427 100644
--- a/tests/spm_dof_sort_tests.c
+++ b/tests/spm_dof_sort_tests.c
@@ -2,12 +2,12 @@
*
* @file spm_dof_sort_tests.c
*
- * Tests and validate the spm_sort routines when the spm_tests.hold constant and/or variadic dofs.
+ * Tests and validate the spm_sort routines when the spm contains constant and/or variadic dofs.
*
- * @copyright 2015-2017 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+ * @copyright 2015-2020 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
* Univ. Bordeaux. All rights reserved.
*
- * @version 6.0.0
+ * @version 1.0.0
* @author Mathieu Faverge
* @author Delarue Tony
* @date 2020-09-07
@@ -25,31 +25,189 @@
#define PRINT_RES(_ret_) \
if(_ret_) { \
printf("FAILED(%d)\n", _ret_); \
- err++; \
} \
else { \
printf("SUCCESS\n"); \
}
+/**
+ *******************************************************************************
+ *
+ * @brief This routine unsorts the spm matrix to check our sort routine.
+ * It will only change the pattern, the value array doesn't follow it.
+ *
+ *******************************************************************************
+ *
+ * @param[inout] spm
+ * On entry, the pointer to the sparse matrix structure.
+ * On exit, the same sparse matrix with subarrays of edges unsorted
+ *
+ *******************************************************************************/
+static inline void
+spm_unsort( spmatrix_t *spm )
+{
+ spm_int_t i, j, size;
+ spm_int_t index1, index2, count;
+ spm_int_t baseval;
+ spm_int_t coltmp, rowtmp;
+ spm_int_t *colptr = spm->colptr;
+ spm_int_t *rowptr = spm->rowptr;
+
+ if( spm->dof < 0 ) {
+ return;
+ }
+
+ baseval = spmFindBase(spm);
+ switch (spm->fmttype)
+ {
+ case SpmCSR:
+ /* Swap pointers to call CSC */
+ colptr = spm->rowptr;
+ rowptr = spm->colptr;
+
+ spm_attr_fallthrough;
+
+ case SpmCSC:
+ size = spm->n;
+ for ( j = 0; j < size; j++, colptr++ )
+ {
+ count = colptr[1] - colptr[0];
+ for ( i = 0; i < count; i++ )
+ {
+ index1 = ( rand() % count ) + colptr[0] - baseval;
+ index2 = ( rand() % count ) + colptr[0] - baseval;
+
+ rowtmp = rowptr[index1];
+ rowptr[index1] = rowptr[index2];
+ rowptr[index2] = rowtmp;
+ break;
+ }
+ }
+ break;
+
+ case SpmIJV:
+ size = spm->nnz;
+ for ( i = 0; i < size; i++ )
+ {
+ index1 = ( rand() % size );
+ index2 = ( rand() % size );
+
+ coltmp = colptr[index1];
+ rowtmp = rowptr[index1];
+
+ colptr[index1] = colptr[index2];
+ rowptr[index1] = rowptr[index2];
+
+ colptr[index2] = coltmp;
+ rowptr[index2] = rowtmp;
+ }
+ break;
+ }
+}
+
static inline int
-spm_sort_check( spmatrix_t *spm )
+spm_sort_check_csx( const spmatrix_t *spm )
{
- spmatrix_t expand1, expand2;
- int rc;
+ spm_int_t i, j, max;
+ spm_int_t n = spm->n;
+ spm_int_t *colptr = (spm->fmttype == SpmCSC) ? spm->colptr : spm->rowptr;
+ spm_int_t *rowptr = (spm->fmttype == SpmCSC) ? spm->rowptr : spm->colptr;
- spmExpand( spm, &expand1 );
+ for ( j = 0; j < n; j++, colptr++ )
+ {
+ max = (colptr[1] - 1);
+ for ( i = colptr[0]; i < max; i++, rowptr++ )
+ {
+ if( rowptr[0] > rowptr[1] ) {
+ return 1;
+ }
+ }
+ rowptr++;
+ }
+
+ return 0;
+}
+
+static inline int
+spm_sort_check_ijv( const spmatrix_t *spm )
+{
+ spm_int_t k;
+ spm_int_t nnz = spm->nnz - 1;
+ spm_int_t *colptr = spm->colptr;
+ spm_int_t *rowptr = spm->rowptr;
+
+ k = 0;
+ while ( k < nnz )
+ {
+ while ( colptr[0] == colptr[1] )
+ {
+ if( rowptr[0] > rowptr[1] ) {
+ return 1;
+ }
+ k++;
+ colptr++;
+ rowptr++;
+ if (k == nnz) {
+ return 0;
+ }
+ }
+ if( colptr[0] > colptr[1] ) {
+ return 1;
+ }
+ k++;
+ colptr++;
+ rowptr++;
+ }
+ return 0;
+}
+
+static inline int
+spm_sort_check( spmatrix_t *spm)
+{
+ spmatrix_t *spmcpy;
+ int rc1, rc2;
+
+ spm_unsort(spm);
+
+ spmcpy = spmCopy(spm);
+ spmSort( spmcpy );
- spmSort( spm );
- spmSort( &expand1 );
+ /* Check that the matrix pattern is well sorted */
+ if ( spm->fmttype != SpmIJV ) {
+ rc1 = spm_sort_check_csx( spmcpy );
+ }
+ else {
+ rc1 = spm_sort_check_ijv( spmcpy );
+ }
+
+ /* Check that the matrix values follows the original pattern */
+ switch (spm->flttype)
+ {
+ case SpmFloat:
+ rc2 = s_spm_sort_check_values(spm, spmcpy);
+ break;
- spmExpand( spm, &expand2 );
+ case SpmDouble:
+ rc2 = d_spm_sort_check_values(spm, spmcpy);
+ break;
- rc = spmCompare( &expand1, &expand2 );
+ case SpmComplex32:
+ rc2 = c_spm_sort_check_values(spm, spmcpy);
+ break;
+
+ case SpmComplex64:
+ rc2 = z_spm_sort_check_values(spm, spmcpy);
+ break;
+
+ default:
+ rc2 = 0;
+ break;
+ }
- spmExit( &expand1 );
- spmExit( &expand2 );
+ spmExit(spmcpy);
+ free(spmcpy);
- return rc;
+ return rc1 + rc2;
}
int main (int argc, char **argv)
diff --git a/tests/spm_tests.h b/tests/spm_tests.h
index e747b28013d5597aae97571eff05efb22602f971..4a63693ffe9de1af1c92feb017b9a8b63dd32992 100644
--- a/tests/spm_tests.h
+++ b/tests/spm_tests.h
@@ -17,6 +17,7 @@
#define _spm_tests_h_
#include <stdlib.h>
+#include <string.h>
#include <assert.h>
#include <stdio.h>
#include <spm.h>
@@ -86,6 +87,7 @@ int z_spm_dist_norm_check( const spmatrix_t *spm, const spmatrix_t *spmdist );
int z_spm_dist_genrhs_check( const spmatrix_t *spm, spm_int_t nrhs,
const spm_complex64_t *bloc, const spm_complex64_t *bdst, int root );
int z_spm_dist_matvec_check( spm_int_t baseval, spm_trans_t trans, const spmatrix_t *spm );
+int z_spm_sort_check_values( const spmatrix_t *spm1, const spmatrix_t *spm2 );
void c_spm_print_check( char *filename, const spmatrix_t *spm );
int c_spm_matvec_check( spm_trans_t trans, const spmatrix_t *spm );
@@ -94,6 +96,7 @@ int c_spm_dist_norm_check( const spmatrix_t *spm, const spmatrix_t *spmdist );
int c_spm_dist_genrhs_check( const spmatrix_t *spm, spm_int_t nrhs,
const spm_complex32_t *bloc, const spm_complex32_t *bdst, int root );
int c_spm_dist_matvec_check( spm_int_t baseval, spm_trans_t trans, const spmatrix_t *spm );
+int c_spm_sort_check_values( const spmatrix_t *spm1, const spmatrix_t *spm2 );
void d_spm_print_check( char *filename, const spmatrix_t *spm );
int d_spm_matvec_check( spm_trans_t trans, const spmatrix_t *spm );
@@ -102,6 +105,7 @@ int d_spm_dist_norm_check( const spmatrix_t *spm, const spmatrix_t *spmdist );
int d_spm_dist_genrhs_check( const spmatrix_t *spm, spm_int_t nrhs,
const double *bloc, const double *bdst, int root );
int d_spm_dist_matvec_check( spm_int_t baseval, spm_trans_t trans, const spmatrix_t *spm );
+int d_spm_sort_check_values( const spmatrix_t *spm1, const spmatrix_t *spm2 );
void s_spm_print_check( char *filename, const spmatrix_t *spm );
int s_spm_matvec_check( spm_trans_t trans, const spmatrix_t *spm );
@@ -110,6 +114,7 @@ int s_spm_dist_norm_check( const spmatrix_t *spm, const spmatrix_t *spmdist );
int s_spm_dist_genrhs_check( const spmatrix_t *spm, spm_int_t nrhs,
const float *bloc, const float *bdst, int root );
int s_spm_dist_matvec_check( spm_int_t baseval, spm_trans_t trans, const spmatrix_t *spm );
+int s_spm_sort_check_values( const spmatrix_t *spm1, const spmatrix_t *spm2 );
void p_spm_print_check( char *filename, const spmatrix_t *spm );
diff --git a/tests/z_spm_sort_tests.c b/tests/z_spm_sort_tests.c
new file mode 100644
index 0000000000000000000000000000000000000000..8f8e399763b8577c722e553b372d0776850ea734
--- /dev/null
+++ b/tests/z_spm_sort_tests.c
@@ -0,0 +1,292 @@
+/**
+ *
+ * @file z_spm_sort_tests.c
+ *
+ * Tests and validate the spm_sort routines.
+ *
+ * @copyright 2015-2020 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+ * Univ. Bordeaux. All rights reserved.
+ *
+ * @version 1.0.0
+ * @author Mathieu Faverge
+ * @author Tony Delarue
+ * @date 2020-09-14
+ *
+ * @precisions normal z -> c d s
+ *
+ **/
+#include <stdint.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <spm_tests.h>
+#include <z_spm.h>
+
+/**
+ * @brief Indicates if the col/row couple is already in the memory array
+ *
+ * @param[inout] memory
+ * The array that contains the already inspect col/row.
+ * If a value has already been seen, it's value is set to -1.
+ *
+ * @param[in] col
+ * The column index.
+ *
+ * @param[in] row
+ * The row index.
+ *
+ * @param[in] size
+ * The size of the memory array.
+ *
+ * @param[in] ijv
+ * Indicates spm format type, ijv or csx.
+ */
+static inline spm_int_t
+z_spm_sort_is_in( spm_int_t *memory,
+ spm_int_t col,
+ spm_int_t row,
+ spm_int_t size,
+ int ijv )
+{
+ spm_int_t i;
+
+ if( !ijv ) {
+ for ( i = 0; i < size; i++)
+ {
+ if( memory[i] == row ) {
+ memory[i] = -1;
+ return 1;
+ }
+ }
+ }
+ else
+ {
+ for ( i = 0; i < size; i++)
+ {
+ if( (memory[2*i] == col) && (memory[2*i+1] == row) ) {
+ memory[2*i] = -1;
+ memory[2*i+1] = -1;
+ return 1;
+ }
+ }
+ }
+
+ return 0;
+}
+
+/**
+ * @brief Check that the value array of this two spm match
+ * for csx formatted matrices.
+ *
+ * We have to be aware that doublons may exists in the spm, so a memory will be
+ * set up in this routine to store the already seen col/row index. It's
+ * necessary to be sure we've checked the same information.
+ *
+ * @param[in] spm1
+ * The unsorted spm.
+ *
+ * @param[in] spm2
+ * The sorted spm.
+ */
+static inline int
+z_spm_sort_check_values_csx( const spmatrix_t *spm1,
+ const spmatrix_t *spm2 )
+{
+ /* The spm2 is the one which is sorted */
+ spm_int_t i, j, krow, kval, n;
+ spm_int_t count, size;
+ spm_int_t dofi, dofj, dof2, baseval;
+ spm_int_t *col1, *col2;
+ spm_int_t *row1, *row2;
+ spm_int_t *dofs;
+ spm_complex64_t *val1, *val2, *valtmp;
+
+ col1 = (spm1->fmttype == SpmCSC) ? spm1->colptr : spm1->rowptr;
+ col2 = (spm2->fmttype == SpmCSC) ? spm2->colptr : spm2->rowptr;
+ row1 = (spm1->fmttype == SpmCSC) ? spm1->rowptr : spm1->colptr;
+ row2 = (spm2->fmttype == SpmCSC) ? spm2->rowptr : spm2->colptr;
+ val1 = spm1->values;
+ val2 = spm2->values;
+
+ n = spm1->n;
+ dofs = spm1->dofs;
+ baseval = spmFindBase(spm1);
+ for ( i = 0; i < n; i++, col1++, col2++)
+ {
+ /* Make sure we're on the same col */
+ if( *col1 != *col2 ) {
+ return 1;
+ }
+
+ size = col1[1] - col1[0];
+ count = 0;
+ dofj = (spm1->dof > 0) ? spm1->dof : dofs[i+1] - dofs[i];
+ valtmp = val1;
+ while( count < size )
+ {
+ spm_int_t row_stored[size];
+ memset( row_stored, 0xff, size * sizeof(spm_int_t) );
+ krow = 0;
+ kval = 0;
+ /*
+ * Get the same row to do the comparison
+ * row2 is sorted.
+ */
+ while( *row1 > row2[krow] )
+ {
+ j = row2[krow] - baseval;
+ dofi = (spm2->dof > 0) ? spm2->dof : dofs[j+1] - dofs[j];
+ kval += dofj * dofi;
+ krow++;
+ }
+
+ j = *row1 - baseval;
+ dofi = (spm1->dof > 0) ? spm1->dof : dofs[j+1] - dofs[j];
+ dof2 = dofj * dofi;
+
+ while( z_spm_sort_is_in( row_stored, -1, *row1, size, 0 ) )
+ {
+ krow++;
+ kval += dof2;
+ }
+ /* Make sure we're on the same row */
+ assert( *row1 == row2[krow] );
+
+ row_stored[count] = *row1;
+
+ /* Check the values array with multidof */
+ for ( j = 0; j < dof2; j++)
+ {
+ if ( val1[j] != val2[kval + j] ) {
+ return 1;
+ }
+ }
+ count++;
+ row1++;
+ val1 += dof2;
+ }
+ row2 += size;
+ val2 += (val1 - valtmp);
+ }
+
+ return 0;
+}
+
+/**
+ * @brief Check that the value array of this two spm match
+ * for ijv formatted matrices.
+ *
+ * We have to be aware that doublons may exists in the spm, so a memory will be
+ * set up in this routine to store the already seen col/row index. It's
+ * necessary to be sure we've checked the same information.
+ *
+ * @param[in] spm1
+ * The unsorted spm.
+ *
+ * @param[in] spm2
+ * The sorted spm.
+ */
+static inline int
+z_spm_sort_check_values_ijv( const spmatrix_t *spm1,
+ const spmatrix_t *spm2 )
+{
+ spm_int_t i, j, k, kptr, kval, n;
+ spm_int_t dofi, dofj, dof2, baseval;
+ spm_int_t *col1, *col2;
+ spm_int_t *row1, *row2;
+ spm_int_t *dofs;
+ spm_int_t *memory;
+ spm_complex64_t *val1, *val2;
+
+ col1 = spm1->colptr;
+ col2 = spm2->colptr;
+ row1 = spm1->rowptr;
+ row2 = spm2->rowptr;
+ val1 = spm1->values;
+ val2 = spm2->values;
+
+ n = spm1->nnz;
+ dofs = spm1->dofs;
+ baseval = spmFindBase(spm1);
+ memory = malloc( 2* n * sizeof(spm_int_t) );
+ memset( memory, 0xff, 2 * n * sizeof(spm_int_t) );
+ for ( k = 0; k < n; k++, col1++, row1++ )
+ {
+ kptr = 0;
+ kval = 0;
+ while( (*col1 != col2[kptr]) || (*row1 != row2[kptr]) )
+ {
+ i = row2[kptr] - baseval;
+ j = col2[kptr] - baseval;
+ dofi = (spm2->dof > 0) ? spm2->dof : dofs[i+1] - dofs[i];
+ dofj = (spm2->dof > 0) ? spm2->dof : dofs[j+1] - dofs[j];
+ kval += dofj * dofi;
+
+ kptr++;
+ assert(kptr < n);
+ }
+
+ i = row2[kptr] - baseval;
+ j = col2[kptr] - baseval;
+ dofi = (spm2->dof > 0) ? spm2->dof : dofs[i+1] - dofs[i];
+ dofj = (spm2->dof > 0) ? spm2->dof : dofs[j+1] - dofs[j];
+ dof2 = dofj * dofi;
+
+ /* Check for doublons */
+ while( z_spm_sort_is_in( memory, *col1, *row1, k, 1 ) )
+ {
+ kptr++;
+ kval += dof2;
+ }
+ assert( (*col1 == col2[kptr]) && (*row1 == row2[kptr]) );
+
+ memory[2*k] = *col1;
+ memory[2*k+1] = *row1;
+
+ for ( i = 0; i < dof2; i++)
+ {
+ if ( val1[i] != val2[kval + i] ) {
+ free(memory);
+ return 1;
+ }
+ }
+
+ val1 += dof2;
+ }
+ free(memory);
+
+ return 0;
+}
+
+/**
+ *******************************************************************************
+ *
+ * @brief Check that the value array corresponds after the sort routine.
+ * These routines should not be called with multiple MPI processes.
+ *
+ *******************************************************************************
+ *
+ * @param[in] spm1
+ * A pointer to the unsorted spm.
+ *
+ * @param[in] spm2
+ * A pointer to the sorted spm.
+ *
+ *******************************************************************************/
+int
+z_spm_sort_check_values( const spmatrix_t *spm1,
+ const spmatrix_t *spm2 )
+{
+ int rc;
+
+ assert(spm1->fmttype == spm2->fmttype);
+
+ if( spm1->fmttype != SpmIJV ) {
+ rc = z_spm_sort_check_values_csx(spm1, spm2);
+ }
+ else
+ {
+ rc = z_spm_sort_check_values_ijv(spm1, spm2);
+ }
+ return rc;
+}