Merge branch 'feature/inplace_ijv_to_csx_recup' into 'master'

Feature/inplace ijv to csx recup

See merge request !3

src/z_spm.c

@@ -23,10 +23,12 @@
  *
  * @ingroup spm_dev_check
  *
- * @brief This routine sorts the subarray of edges of each vertex in a
- * centralized spm stored in CSC or CSR format.
+ * @brief This routine sorts the spm matrix.
  *
- * Nothing is performed if IJV  format is used.
+ * 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.
  *
@@ -44,13 +46,13 @@ 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[2];
+    void *sortptr[3];
     spm_int_t n = spm->n;
     spm_int_t i, size;
     (void)sortptr;
 
     if (spm->dof > 1){
-        fprintf(stderr, "z_spmSort: Number of dof (%d) greater than one not supported\n", (int)spm->dof);
+        fprintf(stderr, "z_spmSort: Number of dof (%d) different from one not supported\n", (int)spm->dof);
         exit(1);
     }
 
@@ -87,6 +89,19 @@ z_spmSort( spmatrix_t *spm )
             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
+    }
 }
 
 /**

src/z_spm_convert_to_csc.c

@@ -40,12 +40,8 @@
 int
 z_spmConvertIJV2CSC( spmatrix_t *spm )
 {
-#if !defined(PRECISION_p)
-    spm_complex64_t *navals = NULL;
-    spm_complex64_t *oavals = NULL;
-#endif
-    spm_int_t       *spmptx, *otmp;
-    spm_int_t i, j, tmp, baseval, total;
+    spm_int_t *newcol, *oldcol;
+    spm_int_t  i, tmp, baseval, total;
     spmatrix_t oldspm;
 
     /* Backup the input */
@@ -56,66 +52,35 @@ z_spmConvertIJV2CSC( spmatrix_t *spm )
      */
     baseval = spmFindBase( spm );
 
-    /* Compute the new colptr */
+    /*
+     * Sort the IJV structure by column/row indexes
+     */
+    z_spmSort( spm );
+
+    /* Allocate and compute the new colptr */
     spm->colptr = (spm_int_t *) calloc(spm->n+1,sizeof(spm_int_t));
 
     /* Compute the number of edges per row */
-    spmptx = spm->colptr - baseval;
-    otmp   = oldspm.colptr;
-    for (i=0; i<spm->nnz; i++, otmp++)
+    newcol = spm->colptr - baseval;
+    oldcol = oldspm.colptr;
+    for (i=0; i<spm->nnz; i++, oldcol++)
     {
-        spmptx[ *otmp ] ++;
+        newcol[ *oldcol ] ++;
     }
 
-    /* Compute the indexes in C numbering for the following sort */
-    total = 0;
-    spmptx = spm->colptr;
-    for (i=0; i<(spm->n+1); i++, spmptx++)
+    /* Update the colptr */
+    total  = baseval;
+    newcol = spm->colptr;
+    for (i=0; i<(spm->n+1); i++, newcol++)
     {
-        tmp = *spmptx;
-        *spmptx = total;
+        tmp = *newcol;
+        *newcol = total;
         total += tmp;
     }
-    assert( total == spm->nnz );
-
-    /* Sort the rows and avals arrays by column */
-    spm->rowptr = malloc(spm->nnz * sizeof(spm_int_t));
-
-#if defined(PRECISION_p)
-    spm->values = NULL;
-#else
-    spm->values = malloc(spm->nnz * sizeof(spm_complex64_t));
-    navals = (spm_complex64_t*)(spm->values);
-    oavals = (spm_complex64_t*)(oldspm.values);
-#endif
-
-    for (j=0; j<spm->nnz; j++)
-    {
-        i = oldspm.colptr[j] - baseval;
-
-        spm->rowptr[ spm->colptr[i] ] = oldspm.rowptr[j];
-
-#if !defined(PRECISION_p)
-        navals[ spm->colptr[i] ] = oavals[j];
-#endif
-        (spm->colptr[i])++;
-
-        assert( spm->colptr[i] <= spm->colptr[i+1] );
-    }
-
-    /* Rebuild the colptr with the correct baseval */
-    tmp = spm->colptr[0];
-    spm->colptr[0] = baseval;
-
-    spmptx = spm->colptr + 1;
-    for (i=1; i<(spm->n+1); i++, spmptx++)
-    {
-        total = *spmptx;
-        *spmptx = tmp + baseval;
-        tmp = total;
-    }
-    assert( spm->colptr[ spm->n ] == (spm->nnz+baseval) );
+    assert( (total-baseval) == spm->nnz );
 
+    oldspm.rowptr = NULL;
+    oldspm.values = NULL;
     spmExit( &oldspm );
 
     spm->fmttype = SpmCSC;

src/z_spm_convert_to_csr.c

@@ -13,6 +13,7 @@
  * @date 2015-01-01
  *
  * @precisions normal z -> c d s p
+ *
  **/
 #include "common.h"
 #include "z_spm.h"
@@ -123,12 +124,8 @@ z_spmConvertCSC2CSR( spmatrix_t *spm )
 int
 z_spmConvertIJV2CSR( spmatrix_t *spm )
 {
-#if !defined(PRECISION_p)
-    spm_complex64_t *navals = NULL;
-    spm_complex64_t *oavals = NULL;
-#endif
-    spm_int_t       *spmptx, *otmp;
-    spm_int_t i, j, tmp, baseval, total;
+    spm_int_t *newrow, *oldrow;
+    spm_int_t  i, tmp, baseval, total;
     spmatrix_t oldspm;
 
     /* Backup the input */
@@ -139,66 +136,40 @@ z_spmConvertIJV2CSR( spmatrix_t *spm )
      */
     baseval = spmFindBase( spm );
 
+    /*
+     * Sort the IJV structure by column/row indexes
+     */
+    newrow = spm->colptr;
+    spm->colptr = spm->rowptr;
+    spm->rowptr = newrow;
+    z_spmSort( spm );
+    spm->rowptr = spm->colptr;
+    spm->colptr = newrow;
+
     /* Compute the new rowptr */
     spm->rowptr = (spm_int_t *) calloc(spm->n+1,sizeof(spm_int_t));
 
     /* Compute the number of edges per row */
-    spmptx = spm->rowptr - baseval;
-    otmp   = oldspm.rowptr;
-    for (i=0; i<spm->nnz; i++, otmp++)
+    newrow = spm->rowptr - baseval;
+    oldrow = oldspm.rowptr;
+    for (i=0; i<spm->nnz; i++, oldrow++)
     {
-        spmptx[ *otmp ] ++;
+        newrow[ *oldrow ] ++;
     }
 
-    /* Compute the indexes in C numbering for the following sort */
-    total = 0;
-    spmptx = spm->rowptr;
-    for (i=0; i<(spm->n+1); i++, spmptx++)
+    /* Update the rowptr */
+    total  = baseval;
+    newrow = spm->rowptr;
+    for (i=0; i<(spm->n+1); i++, newrow++)
     {
-        tmp = *spmptx;
-        *spmptx = total;
+        tmp = *newrow;
+        *newrow = total;
         total += tmp;
     }
-    assert( total == spm->nnz );
-
-    /* Sort the colptr and avals arrays by rows */
-    spm->colptr  = malloc(spm->nnz * sizeof(spm_int_t));
-
-#if defined(PRECISION_p)
-    spm->values = NULL;
-#else
-    spm->values = malloc(spm->nnz * sizeof(spm_complex64_t));
-    navals = (spm_complex64_t*)(spm->values);
-    oavals = (spm_complex64_t*)(oldspm.values);
-#endif
-
-    for (j=0; j<spm->nnz; j++)
-    {
-        i = oldspm.rowptr[j] - baseval;
-
-        spm->colptr[ spm->rowptr[i] ] = oldspm.colptr[j];
-
-#if !defined(PRECISION_p)
-        navals[ spm->rowptr[i] ] = oavals[j];
-#endif
-        (spm->rowptr[i])++;
-
-        assert( spm->rowptr[i] <= spm->rowptr[i+1] );
-    }
-
-    /* Rebuild the rows (rowptr) with the correct baseval */
-    tmp = spm->rowptr[0];
-    spm->rowptr[0] = baseval;
-
-    spmptx = spm->rowptr + 1;
-    for (i=1; i<(spm->n+1); i++, spmptx++)
-    {
-        total = *spmptx;
-        *spmptx = tmp + baseval;
-        tmp = total;
-    }
-    assert( spm->rowptr[ spm->n ] == (spm->nnz+baseval) );
+    assert( (total-baseval) == spm->nnz );
 
+    oldspm.colptr = NULL;
+    oldspm.values = NULL;
     spmExit( &oldspm );
 
     spm->fmttype = SpmCSR;