diff --git a/src/z_spm.c b/src/z_spm.c
index 0b0f96c112a609fbe354af4146715a24b58b70df..8b6c4617561dc99bd4a2d1e4e7e1526f83d99127 100644
--- a/src/z_spm.c
+++ b/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
+ }
}
/**
diff --git a/src/z_spm_convert_to_csc.c b/src/z_spm_convert_to_csc.c
index 971ddf0f5cbd1f9c24cdc4140963f560ee287be3..0d7946890a1a5a84506490355356fbcd2571b693 100644
--- a/src/z_spm_convert_to_csc.c
+++ b/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;
diff --git a/src/z_spm_convert_to_csr.c b/src/z_spm_convert_to_csr.c
index 9e412ff68701e8662583bc10208a0d6f44cc2aee..1cae234dd77859572e3b21b4bb923b8535e46a6e 100644
--- a/src/z_spm_convert_to_csr.c
+++ b/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;