diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1a0d7de78590ef8ca16602b60bac1e0b30cf94f9..3f897fcdd6718df2df7fb1caef63200e90c20644 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -196,7 +196,7 @@ set(SOURCES
src/z_spm_dof_extend.c
src/z_spm_norm.c
src/z_spm_scal.c
-
+ src/z_spm_reduce_rhs.c
src/z_spm_convert_to_csc.c
src/z_spm_convert_to_csr.c
src/z_spm_convert_to_ijv.c
diff --git a/src/z_spm.h b/src/z_spm.h
index 526985f3f32611252154749ba0a25579fbd15d40..4fbe1b06122721419c0bff646b3df1ccce6f075b 100644
--- a/src/z_spm.h
+++ b/src/z_spm.h
@@ -73,8 +73,9 @@ void z_spmSort( spmatrix_t *spm );
spm_int_t z_spmMergeDuplicate( spmatrix_t *spm );
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 );
+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_spmReduceRHS( const spmatrix_t *spm, int nrhs, const spm_complex64_t *x, spm_int_t ldx, int root );
/**
* Output routines
diff --git a/src/z_spm_reduce_rhs.c b/src/z_spm_reduce_rhs.c
new file mode 100644
index 0000000000000000000000000000000000000000..69e96fb3dfd0df5e1cec86b81c2a51500659c993
--- /dev/null
+++ b/src/z_spm_reduce_rhs.c
@@ -0,0 +1,155 @@
+/**
+ *
+ * @file z_spm_reduce_rhs.c
+ *
+ * SParse Matrix package right hand side reduction routine.
+ *
+ * @copyright 2020-2020 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+ * Univ. Bordeaux. All rights reserved.
+ *
+ * @version 1.0.0
+ * @author Delarue Tony
+ * @author Faverge Mathieu
+ * @date 2020-06-11
+ *
+ * @precisions normal z -> c d s
+ *
+ **/
+#include "common.h"
+#include "z_spm.h"
+
+/**
+ *******************************************************************************
+ *
+ * @brief Reduce all the global C coefficients and store the good ones in local
+ *
+ *******************************************************************************
+ *
+ * @param[in] spm
+ * The sparse matrix spm
+ *
+ * @param[inout] glob
+ * The global vector to reduce
+ *
+ * @param[inout] loc
+ * Local vector
+ *
+ *******************************************************************************/
+spm_complex64_t *
+z_spmReduceRHS( const spmatrix_t *spm,
+ int nrhs,
+ const spm_complex64_t *x,
+ spm_int_t ldx,
+ int root )
+{
+ spm_complex64_t *outptr, *out = NULL;
+ int c;
+ int n, nmax = 0;
+ int nexp, nexpmax = 0;
+ spm_complex64_t *tmp_out, *current_out;
+ spm_int_t *tmp_l2g, *current_l2g;
+ spm_int_t current_n, current_nexp;
+ spm_int_t i, j, k, ig, dofi, row;
+ int local = ( ( root == -1 ) || (root == spm->clustnum) );
+ spm_int_t baseval = spmFindBase( spm );
+
+ /* We do not handle cases where ldx is different from spm->n */
+ assert( ldx == spm->nexp );
+
+ if ( spm->loc2glob == NULL ) {
+ if ( local ) {
+ out = malloc( spm->gNexp * nrhs * sizeof( spm_complex64_t ) );
+ memcpy( out, x, spm->gNexp * nrhs * sizeof( spm_complex64_t ) );
+ }
+ return out;
+ }
+
+#if !defined(SPM_WITH_MPI)
+ assert( 0 );
+#else
+ n = spm->n;
+ nexp = spm->nexp;
+ if ( root == -1 ) {
+ MPI_Allreduce( &n, &nmax, 1, MPI_INT, MPI_MAX, spm->comm );
+ MPI_Allreduce( &nexp, &nexpmax, 1, MPI_INT, MPI_MAX, spm->comm );
+ }
+ else {
+ MPI_Reduce( &n, &nmax, 1, MPI_INT, MPI_MAX, root, spm->comm );
+ MPI_Reduce( &nexp, &nexpmax, 1, MPI_INT, MPI_MAX, root, spm->comm );
+ }
+
+ if ( local ) {
+ out = malloc( spm->gNexp * nrhs * sizeof( spm_complex64_t ) );
+ tmp_out = malloc( nexpmax * nrhs * sizeof( spm_complex64_t ) );
+ tmp_l2g = malloc( nmax * sizeof( spm_int_t ) );
+ }
+
+ for( c=0; c<spm->clustnbr; c++ ) {
+ MPI_Status status;
+
+ if ( c == spm->clustnum ) {
+ current_n = spm->n;
+ current_nexp = spm->nexp;
+ current_l2g = spm->loc2glob;
+ current_out = (spm_complex64_t*)x;
+ }
+ else {
+ current_n = 0;
+ current_nexp = 0;
+ current_l2g = tmp_l2g;
+ current_out = tmp_out;
+ }
+
+ if ( root == -1 ) {
+ MPI_Bcast( ¤t_n, 1, SPM_MPI_INT, c, spm->comm );
+ MPI_Bcast( ¤t_nexp, 1, SPM_MPI_INT, c, spm->comm );
+ if ( current_n > 0 ) {
+ MPI_Bcast( current_l2g, current_n, SPM_MPI_INT, c, spm->comm );
+ MPI_Bcast( current_out, current_nexp * nrhs, SPM_MPI_COMPLEX64, c, spm->comm );
+ }
+ }
+ else if ( root != c ) { /* No communication if c == root */
+ if ( c == spm->clustnum ) {
+ MPI_Send( ¤t_n, 1, SPM_MPI_INT, root, 0, spm->comm );
+ MPI_Send( ¤t_nexp, 1, SPM_MPI_INT, root, 1, spm->comm );
+ if ( current_n > 0 ) {
+ MPI_Send( current_l2g, current_n, SPM_MPI_INT, root, 2, spm->comm );
+ MPI_Send( current_out, current_nexp * nrhs, SPM_MPI_COMPLEX64, root, 3, spm->comm );
+ }
+ }
+ if ( root == spm->clustnum ) {
+ MPI_Recv( ¤t_n, 1, SPM_MPI_INT, c, 0, spm->comm, &status );
+ MPI_Recv( ¤t_nexp, 1, SPM_MPI_INT, c, 1, spm->comm, &status );
+ if ( current_n > 0 ) {
+ MPI_Recv( current_l2g, current_n, SPM_MPI_INT, c, 2, spm->comm, &status );
+ MPI_Recv( current_out, current_nexp * nrhs, SPM_MPI_COMPLEX64, c, 3, spm->comm, &status );
+ }
+ }
+ }
+
+ if ( !local ) {
+ continue;
+ }
+
+ outptr = out;
+ for( i=0; i<current_n; i++, current_l2g++ ) {
+ ig = *current_l2g - baseval;
+ dofi = ( spm->dof > 0 ) ? spm->dof : spm->dofs[ig+1] - spm->dofs[ig];
+ row = ( spm->dof > 0 ) ? spm->dof * ig : spm->dofs[ig] - baseval;
+ for( j=0; j<nrhs; j++ ) {
+ for( k=0; k<dofi; k++ ) {
+ outptr[ row + j * spm->gNexp + k ] = current_out[ j * current_nexp + k ];
+ }
+ }
+ current_out += dofi;
+ }
+ }
+
+ if ( local ) {
+ free( tmp_out );
+ free( tmp_l2g );
+ }
+
+#endif
+ return out;
+}