spmm and spmv now support multidof

include/spm.h

@@ -116,7 +116,6 @@ spmatrix_t *spmScatter( const spmatrix_t *spm,
                         SPM_Comm          comm );
 spmatrix_t *spmGather ( const spmatrix_t *spm,
                               int         root );
-int spmGetDistribution( const spmatrix_t *spm );
 
 /**
  * @}

include/spm/const.h

@@ -32,6 +32,12 @@ BEGIN_C_DECLS
 #define CBLAS_SADDR( a_ ) (&(a_))
 #endif
 
+/**
+ * @brief Distribution of the matrix storage
+ */
+#define SpmDistByColumn (0x1 << 0) /**< Storage in column distributed */
+#define SpmDistByRow    (0x1 << 1) /**< Storage in row distributed    */
+
 /**
  * @brief Verbose modes
  */

src/common.h

@@ -54,6 +54,7 @@ spm_get_datatype( const spmatrix_t *spm )
 #endif
 
 spm_int_t *spm_get_glob2loc( spmatrix_t *spm, spm_int_t baseval );
+int        spm_get_distribution( const spmatrix_t *spm );
 
 /********************************************************************
  * Conjuguate/Id functions

src/spm.c

@@ -128,70 +128,6 @@ spmInitDist( spmatrix_t *spm, SPM_Comm comm )
 #endif /* defined(SPM_WITH_MPI) */
 }
 
-/**
- *******************************************************************************
- *
- * @brief Search the distribution pattern used in the spm structure.
- *
- *******************************************************************************
- *
- * @param[in] spm
- *          The sparse matrix structure.
- *
- ********************************************************************************
- *
- * @return  1 if the distribution is column based.
- *          0 otherwise.
- *
- *******************************************************************************/
-int
-spmGetDistribution( const spmatrix_t *spm )
-{
-    int distribution = 1;
-
-    if( spm->fmttype == SpmCSC ){
-        distribution = 1;
-    }
-    else if ( spm->fmttype == SpmCSR ) {
-        distribution = 0;
-    }
-    else {
-        spm_int_t  i, baseval;
-        spm_int_t *colptr   = spm->colptr;
-        spm_int_t *glob2loc = spm->glob2loc;
-
-        baseval = spmFindBase( spm );
-        assert( glob2loc != NULL );
-        for ( i = 0; i < spm->nnz; i++, colptr++ )
-        {
-            /*
-             * If the global index is not in the local colptr
-             * -> row distribution
-             */
-            if( glob2loc[ *colptr - baseval  ] < 0 ) {
-                distribution = 0;
-                break;
-            }
-        }
-
-#if defined(SPM_WITH_MPI)
-        {
-            int check = 0;
-            MPI_Allreduce( &distribution, &check, 1, MPI_INT,
-                           MPI_SUM, spm->comm );
-            if( distribution == 0) {
-                assert( check == 0 );
-            }
-            else {
-                assert( check == spm->clustnbr );
-            }
-        }
-#endif
-    }
-
-    return distribution;
-}
-
 /**
  *******************************************************************************
  *
@@ -1308,10 +1244,6 @@ spmMatMat(       spm_trans_t trans,
     spmatrix_t *espm = (spmatrix_t*)A;
     int rc = SPM_SUCCESS;
 
-    if ( A->dof != 1 ) {
-        espm = malloc( sizeof(spmatrix_t) );
-        spmExpand( A, espm );
-    }
     switch (A->flttype) {
     case SpmFloat:
         rc = spm_sspmm( SpmLeft, trans, SpmNoTrans, n, alpha, espm, B, ldb, beta, C, ldc );
@@ -1682,3 +1614,76 @@ spm_get_glob2loc( spmatrix_t *spm,
     (void) baseval;
     return spm->glob2loc;
 }
+
+/**
+ *******************************************************************************
+ *
+ * @ingroup spm_mpi_dev
+ *
+ * @brief Search the distribution pattern used in the spm structure.
+ *
+ *******************************************************************************
+ *
+ * @param[in] spm
+ *          The sparse matrix structure.
+ *
+ ********************************************************************************
+ *
+ * @return  SpmDistByColumn if the distribution is column based.
+ *          SpmDistByRow if the distribution is row based.
+ *          (SpmDistByColumn|SpmDistByRow) if the matrix is not distributed.
+ *
+ *******************************************************************************/
+int
+spm_get_distribution( const spmatrix_t *spm )
+{
+    int distribution = 0;
+
+    if( (spm->loc2glob == NULL) || (spm->n == spm->gN) ) {
+        distribution = ( SpmDistByColumn | SpmDistByRow );
+    }
+    else {
+        if( spm->fmttype == SpmCSC ){
+            distribution = SpmDistByColumn;
+        }
+        else if ( spm->fmttype == SpmCSR ) {
+            distribution = SpmDistByRow;
+        }
+        else {
+            spm_int_t  i, baseval;
+            spm_int_t *colptr   = spm->colptr;
+            spm_int_t *glob2loc = spm->glob2loc;
+
+            baseval = spmFindBase( spm );
+            distribution = 1;
+            assert( glob2loc != NULL );
+            for ( i = 0; i < spm->nnz; i++, colptr++ )
+            {
+                /*
+                * If the global index is not in the local colptr
+                * -> row distribution
+                */
+                if( glob2loc[ *colptr - baseval  ] < 0 ) {
+                    distribution = SpmDistByRow;
+                    break;
+                }
+            }
+
+    #if defined(SPM_WITH_MPI)
+            {
+                int check = 0;
+                MPI_Allreduce( &distribution, &check, 1, MPI_INT,
+                               MPI_BOR, spm->comm );
+                /*
+                 * If a matrix is distributed
+                 * it cannot be distributed by row AND column
+                 */
+                assert( check != ( SpmDistByColumn | SpmDistByRow ) );
+                assert( distribution == check );
+            }
+    #endif
+        }
+    }
+    assert(distribution > 0);
+    return distribution;
+}

src/z_spm.h

@@ -76,7 +76,7 @@ spm_int_t z_spmSymmetrize( spmatrix_t *spm );
 int              z_spmGenRHS(spm_rhstype_t type, int nrhs, const spmatrix_t *spm, void *x, int ldx, void *b, int ldb );
 int              z_spmCheckAxb( spm_fixdbl_t eps, int nrhs, const spmatrix_t *spm, void *x0, int ldx0, void *b, int ldb, const void *x, int ldx );
 spm_complex64_t *z_spmGatherRHS( const spmatrix_t *spm, int nrhs, const spm_complex64_t *x, spm_int_t ldx, int root );
-void             z_spmReduceRhs( const spmatrix_t *spm, int nrhs, spm_complex64_t *bglob, spm_complex64_t *b, spm_int_t ldb );
+void             z_spmReduceRHS( const spmatrix_t *spm, int nrhs, spm_complex64_t *bglob, spm_int_t ldbglob, spm_complex64_t *b, spm_int_t ldb );
 
 /**
  * Output routines

src/z_spm_convert_to_csc.c

@@ -65,30 +65,13 @@ z_spmConvertIJV2CSC( spmatrix_t *spm )
 
 #if defined(SPM_WITH_MPI)
     if ( spm->loc2glob != NULL ) {
-        /*
-         * Check if the distribution is by column or row by exploiting the fact
-         * that the array is sorted.
-         * This is not completely safe, but that avoids going through the full
-         * matrix.
-         */
         const spm_int_t *glob2loc;
-        spm_int_t m = spm->rowptr[spm->nnz-1] - spm->rowptr[0] + 1; /* This may be not correct */
-        spm_int_t n = spm->colptr[spm->nnz-1] - spm->colptr[0] + 1;
         spm_int_t jg;
-        int distribution = 0;
+        int distribution = spm_get_distribution( spm );
 
-        if ( m <= spm->n ) { /* By row */
-            distribution |= 1;
-        }
-        if ( n <= spm->n ) { /* By column */
-            distribution |= 2;
-        }
-        MPI_Allreduce( MPI_IN_PLACE, &distribution, 1, MPI_INT,
-                       MPI_BAND, spm->comm );
-
-        if ( !(distribution & 2) ) {
-            //fprintf( stderr, "spmConvert: Conversion of column distributed matrices to CSC is not yet implemented\n");
-            return SPM_ERR_NOTIMPLEMENTED;
+        if ( !(distribution & SpmDistByColumn) ) {
+            fprintf( stderr, "spmConvert: Conversion of non column distributed matrices to CSC is not yet implemented\n");
+            return SPM_ERR_BADPARAMETER;
         }
 
         /* Allocate and compute the glob2loc array */

src/z_spm_convert_to_csr.c

@@ -70,29 +70,12 @@ z_spmConvertIJV2CSR( spmatrix_t *spm )
 
 #if defined(SPM_WITH_MPI)
     if ( spm->loc2glob != NULL ) {
-        /*
-         * Check if the distribution is by column or row by exploiting the fact
-         * that the array is sorted.
-         * This is not completely safe, but that avoids going through the full
-         * matrix.
-         */
         const spm_int_t *glob2loc;
-        spm_int_t m = spm->rowptr[spm->nnz-1] - spm->rowptr[0] + 1; /* This may be not correct */
-        spm_int_t n = spm->colptr[spm->nnz-1] - spm->colptr[0] + 1;
         spm_int_t ig;
-        int distribution = 0;
+        int distribution = spm_get_distribution( spm );
 
-        if ( m <= spm->n ) { /* By row */
-            distribution |= 1;
-        }
-        if ( n <= spm->n ) { /* By column */
-            distribution |= 2;
-        }
-        MPI_Allreduce( MPI_IN_PLACE, &distribution, 1, MPI_INT,
-                       MPI_BAND, spm->comm );
-
-        if ( !(distribution & 1) ) {
-            //fprintf( stderr, "spmConvert: Conversion of column distributed matrices to CSC is not yet implemented\n");
+        if ( !(distribution & SpmDistByRow) ) {
+            fprintf( stderr, "spmConvert: Conversion of non row distributed matrices to CSR is not yet implemented\n");
             return SPM_ERR_NOTIMPLEMENTED;
         }
 

src/z_spm_reduce_rhs.c

@@ -41,9 +41,10 @@
  *
  *******************************************************************************/
 void
-z_spmReduceRhs( const spmatrix_t      *spm,
+z_spmReduceRHS( const spmatrix_t      *spm,
                       int              nrhs,
                       spm_complex64_t *bglob,
+                      spm_int_t        ldbglob,
                       spm_complex64_t *b,
                       spm_int_t        ldb )
 {
@@ -57,7 +58,7 @@ z_spmReduceRhs( const spmatrix_t      *spm,
         return;
     }
 
-    MPI_Allreduce( MPI_IN_PLACE, bglob, ldb * nrhs, SPM_MPI_COMPLEX64, MPI_SUM, spm->comm );
+    MPI_Allreduce( MPI_IN_PLACE, bglob, ldbglob * nrhs, SPM_MPI_COMPLEX64, MPI_SUM, spm->comm );
 
     baseval  = spmFindBase( spm );
     loc2glob = spm->loc2glob;
@@ -67,7 +68,7 @@ z_spmReduceRhs( const spmatrix_t      *spm,
         row  = ( spm->dof > 0 ) ? spm->dof * ig : spm->dofs[ig] - baseval;
         for( j=0; j<nrhs; j++ ) {
             for( k=0; k<dofi; k++ ) {
-                rhs[ j * spm->nexp + k ] = bglob[ row + j * ldb + k ];
+                rhs[ j * ldb + k ] = bglob[ row + j * ldbglob + k ];
             }
         }
         rhs += dofi;
@@ -76,6 +77,7 @@ z_spmReduceRhs( const spmatrix_t      *spm,
     (void)spm;
     (void)nrhs;
     (void)bglob;
+    (void)ldbglob;
     (void)b;
     (void)ldb;
 #endif