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 8c92c26a66acda08a6e694bd8ca7f009bb7e6a11..82a7fff45721e5b5dddfdf1e778f209874bd7337 100644
--- a/src/spm.c
+++ b/src/spm.c
@@ -540,11 +540,7 @@ spmNorm( spm_normtype_t   ntype,
 /**
  *******************************************************************************
  *
- * @brief Sort the subarray of edges of each vertex in a CSC or CSR format.
- *
- * Nothing is performed if IJV format is used.
- *
- * @warning This function should NOT be called if dof is greater than 1.
+ * @brief Sort the subarray of edges of each vertex.
  *
  *******************************************************************************
  *
@@ -562,10 +558,6 @@ spmNorm( spm_normtype_t   ntype,
 int
 spmSort( spmatrix_t *spm )
 {
-    if ( (spm->dof != 1) && (spm->flttype != SpmPattern) ) {
-        assert( 0 );
-        fprintf(stderr, "ERROR: spmSort should not be called with non expanded matrices including values\n");
-    }
     switch (spm->flttype) {
     case SpmPattern:
         p_spmSort( spm );
@@ -1687,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/spm_gather.c b/src/spm_gather.c
index 82cdce679c6808cf692b1f74b81f627f2e19a9ed..93bd05b7fbd7ffe0f3f7ca5bceea17b6c3dfbbd9 100644
--- a/src/spm_gather.c
+++ b/src/spm_gather.c
@@ -115,7 +115,7 @@ spm_gather_csx_update( const spmatrix_t *spm,
     for ( c=1; c<spm->clustnbr; c++ ) {
         /* Let's start shifting the value after the first array */
         spm_int_t shift = recvdispls[c];
-        spm_int_t end   = ( c < (spm->clustnbr-1) ) ? recvdispls[c+1] : spm->gN+1;
+        spm_int_t end   = ( c < (spm->clustnbr-1) ) ? recvdispls[c+1] : spm->gN;
         spm_int_t i;
 
         to_add += recvcounts[c-1];
@@ -123,6 +123,9 @@ spm_gather_csx_update( const spmatrix_t *spm,
             colptr[i] += to_add;
         }
     }
+    assert( to_add + recvcounts[spm->clustnbr-1] == spm->gnnz );
+    /* Set the last value */
+    colptr[spm->gN] = colptr[0] + spm->gnnz;
 }
 
 /**
@@ -177,7 +180,6 @@ spm_gather_csx_continuous( const spmatrix_t *oldspm,
                 recvdispls[c] = recvdispls[c-1] + recvcounts[c-1];
                 recvcounts[c] = allcounts[ 3 * c ];
             }
-            recvcounts[ oldspm->clustnbr - 1 ] += 1; /* Add the extra elements */
 
             /* Initialize the new pointers */
             assert( newspm );
@@ -186,11 +188,6 @@ spm_gather_csx_continuous( const spmatrix_t *oldspm,
             newval =  newspm->values;
         }
 
-        /* Add the first value from the first node */
-        if ( oldspm->clustnum == oldspm->clustnbr-1 ) {
-            n++;
-        }
-
         /* Gather the colptr */
         if ( root == -1 ) {
             MPI_Allgatherv( oldcol, n, SPM_MPI_INT,
diff --git a/src/z_spm.c b/src/z_spm.c
index cc64cfefdc713f78e7be73e8f148d1b83cf6b6f5..273593dbd204562b12dd4fd34496626f2b2b2fe6 100644
--- a/src/z_spm.c
+++ b/src/z_spm.c
@@ -18,91 +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( rowptr, 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
-    }
-}
-
 /**
  *******************************************************************************
  *
diff --git a/src/z_spm_expand.c b/src/z_spm_expand.c
index f57aef7b5dfe022bb3bf7c754bcef8b16c23432a..07a38e8ce4def2c48d8a7e2ee8f69f365f4cf27e 100644
--- a/src/z_spm_expand.c
+++ b/src/z_spm_expand.c
@@ -117,9 +117,9 @@ z_spmCSCExpand( const spmatrix_t *spm_in, spmatrix_t *spm_out )
 
         /* Sum the heights of the elements in the column */
         newcol[1] = newcol[0];
-        for(k=oldcol[0]; k<oldcol[1]; k++)
+        for(k=oldcol[0]; k<oldcol[1]; k++, oldrow++ )
         {
-            ig   = oldrow[k-baseval] - baseval;
+            ig   = *oldrow - baseval;
             dofi = (spm_in->dof > 0 ) ? spm_in->dof : dofs[ig+1] - dofs[ig];
             newcol[1] += dofi;
 
@@ -280,9 +280,9 @@ z_spmCSRExpand( const spmatrix_t *spm_in, spmatrix_t *spm_out )
 
         /* Sum the width of the elements in the row */
         newrow[1] = newrow[0];
-        for(k=oldrow[0]; k<oldrow[1]; k++)
+        for(k=oldrow[0]; k<oldrow[1]; k++, oldcol++)
         {
-            jg = oldcol[k-baseval] - baseval;
+            jg   = *oldcol - baseval;
             dofj = (spm_in->dof > 0 ) ? spm_in->dof : dofs[jg+1] - dofs[jg];
             newrow[1] += dofj;
 
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 0040b2ac9603e75bf62a68b5124cce7e234e60d9..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
@@ -46,6 +47,7 @@ set (TESTS
   spm_dof_expand_tests.c
   spm_dof_norm_tests.c
   spm_dof_matvec_tests.c
+  spm_dof_sort_tests.c
   )
 
 if ( SPM_WITH_MPI )
@@ -54,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()
 
@@ -65,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_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
new file mode 100644
index 0000000000000000000000000000000000000000..146be200d0be9cb612421e0cfe95d156e5520427
--- /dev/null
+++ b/tests/spm_dof_sort_tests.c
@@ -0,0 +1,310 @@
+/**
+ *
+ * @file spm_dof_sort_tests.c
+ *
+ * Tests and validate the spm_sort routines when the spm contains constant and/or variadic dofs.
+ *
+ * @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>
+
+#define PRINT_RES(_ret_)                        \
+    if(_ret_) {                                 \
+        printf("FAILED(%d)\n", _ret_);          \
+    }                                           \
+    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_csx( const spmatrix_t *spm )
+{
+    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;
+
+    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 );
+
+    /* 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;
+
+    case SpmDouble:
+        rc2 = d_spm_sort_check_values(spm, spmcpy);
+        break;
+
+    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(spmcpy);
+    free(spmcpy);
+
+    return rc1 + rc2;
+}
+
+int main (int argc, char **argv)
+{
+    spmatrix_t    original, *spm;
+    spm_driver_t  driver;
+    char         *filename;
+    spm_mtxtype_t spmtype, mtxtype;
+    spm_fmttype_t fmttype;
+    int baseval;
+    int rc = SPM_SUCCESS;
+    int err = 0;
+    int i, dofmax = 4;
+
+#if defined(SPM_WITH_MPI)
+    MPI_Init( &argc, &argv );
+#endif
+
+    /**
+     * 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;
+    }
+
+    spmtype = original.mtxtype;
+    printf(" -- SPM Sort Dof Test --\n");
+
+    for( i=0; i<2; i++ )
+    {
+        for( mtxtype=SpmGeneral; mtxtype<=SpmHermitian; mtxtype++ )
+        {
+            if ( (mtxtype == SpmHermitian) &&
+                 ( ((original.flttype != SpmComplex64) && (original.flttype != SpmComplex32)) ||
+                   (spmtype != SpmHermitian) ) )
+            {
+                continue;
+            }
+            if ( (mtxtype != SpmGeneral) &&
+                 (spmtype == SpmGeneral) )
+            {
+                continue;
+            }
+            original.mtxtype = mtxtype;
+
+            for( baseval=0; baseval<2; baseval++ )
+            {
+                spmBase( &original, baseval );
+
+                for( fmttype=SpmCSC; fmttype<=SpmIJV; fmttype++ )
+                {
+                    spmConvert( fmttype, &original );
+                    spm = spmDofExtend( &original, i, dofmax );
+                    if ( spm == NULL ) {
+                        fprintf( stderr, "FAILED to extend matrix\n" );
+                        PRINT_RES(1);
+                        continue;
+                    }
+
+                    printf( " Case: %s / %s / %s / %d / %s\n",
+                            fltnames[spm->flttype],
+                            dofname[i+1],
+                            mtxnames[mtxtype - SpmGeneral],
+                            baseval,
+                            fmtnames[spm->fmttype] );
+
+                    rc = spm_sort_check( spm );
+                    err = (rc == 0) ? err : err + 1;
+                    PRINT_RES(rc);
+
+                    spmExit( spm );
+                    free(spm);
+                    spm = NULL;
+                }
+            }
+        }
+    }
+    spmExit( &original );
+
+#if defined(SPM_WITH_MPI)
+    MPI_Finalize();
+#endif
+
+    if( err == 0 ) {
+        printf(" -- All tests PASSED --\n");
+        return EXIT_SUCCESS;
+    }
+    else
+    {
+        printf(" -- %d tests FAILED --\n", err);
+        return EXIT_FAILURE;
+    }
+}
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;
+}