From b37166fbd7b4c9f2e14cf81cff0db81091f7ffba Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Wed, 26 Mar 2014 17:41:59 +0100
Subject: [PATCH] Rectangular PETSc matrix
---
MathLib/LinAlg/PETSc/PETScMatrix.cpp | 80 ++++++++++-------
MathLib/LinAlg/PETSc/PETScMatrix.h | 97 ++++++++++++---------
Tests/MathLib/TestGlobalMatrixInterface.cpp | 89 +++++++++++++++++--
3 files changed, 183 insertions(+), 83 deletions(-)
diff --git a/MathLib/LinAlg/PETSc/PETScMatrix.cpp b/MathLib/LinAlg/PETSc/PETScMatrix.cpp
index aa5d482f8da..14ce090f443 100644
--- a/MathLib/LinAlg/PETSc/PETScMatrix.cpp
+++ b/MathLib/LinAlg/PETSc/PETScMatrix.cpp
@@ -18,55 +18,55 @@
namespace MathLib
{
-PETScMatrix::PETScMatrix (const PetscInt size)
- :_size(size), _n_loc_rows(PETSC_DECIDE), _n_loc_cols(PETSC_DECIDE)
+PETScMatrix::PETScMatrix(const PetscInt dim)
+ :_nrows(dim), _ncols(dim), _n_loc_rows(PETSC_DECIDE), _n_loc_cols(PETSC_DECIDE)
{
create();
+ config(PETSC_DECIDE, PETSC_DECIDE);
}
-PETScMatrix::PETScMatrix (const PetscInt size, const PetscInt n_loc_cols,
- const bool is_global_size)
- :_size(size), _n_loc_rows(PETSC_DECIDE), _n_loc_cols(n_loc_cols)
-{
- if(!is_global_size)
- {
- _size = PETSC_DECIDE;
- _n_loc_rows = size;
- }
-
- create();
-}
-
-PETScMatrix::PETScMatrix (const PetscInt size, const PETScMatrixOption &mat_opt)
- :_size(size), _n_loc_rows(PETSC_DECIDE), _n_loc_cols(mat_opt._n_local_cols)
+PETScMatrix::PETScMatrix (const PetscInt nrows, const PETScMatrixOption &mat_opt)
+ :_nrows(nrows), _ncols(nrows), _n_loc_rows(PETSC_DECIDE),
+ _n_loc_cols(mat_opt._n_local_cols)
{
if(!mat_opt._is_global_size)
{
- _size = PETSC_DECIDE;
- _n_loc_rows = size;
+ _nrows = PETSC_DECIDE;
+ _ncols = PETSC_DECIDE;
+ _n_loc_rows = nrows;
+
+ // Make the matrix be square.
+ MPI_Allreduce(&_n_loc_rows, &_nrows, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
+ _ncols = _nrows;
}
create();
- config(mat_opt);
+ config(mat_opt._d_nz, mat_opt._o_nz);
}
-void PETScMatrix::create()
+PETScMatrix::PETScMatrix(const PetscInt nrows_global, const PetscInt ncols_global)
+ :_nrows(nrows_global), _ncols(ncols_global), _n_loc_rows(PETSC_DECIDE), _n_loc_cols(PETSC_DECIDE)
{
- MatCreate(PETSC_COMM_WORLD, &_A);
- MatSetSizes(_A, _n_loc_rows, _n_loc_cols, _size, _size);
-
- MatSetFromOptions(_A);
+ create();
+ config(PETSC_DECIDE, PETSC_DECIDE);
}
-void PETScMatrix::config(const PETScMatrixOption &mat_opt)
+PETScMatrix::PETScMatrix (const PetscInt nrows, const PetscInt ncols, const PETScMatrixOption &mat_opt)
+ :_nrows(nrows), _ncols(ncols), _n_loc_rows(PETSC_DECIDE),
+ _n_loc_cols(mat_opt._n_local_cols)
{
- // for a dense matrix: MatSeqAIJSetPreallocation(_A, d_nz, PETSC_NULL);
- MatMPIAIJSetPreallocation(_A, mat_opt._d_nz, PETSC_NULL, mat_opt._o_nz, PETSC_NULL);
+ if(!mat_opt._is_global_size)
+ {
+ _nrows = PETSC_DECIDE;
+ _ncols = PETSC_DECIDE;
+ _n_loc_rows = nrows;
+ _n_loc_cols = ncols;
+ }
- MatGetOwnershipRange(_A, &_start_rank, &_end_rank);
- MatGetSize(_A, &_size, PETSC_NULL);
- MatGetLocalSize(_A, &_n_loc_rows, &_n_loc_cols);
+ create();
+
+ config(mat_opt._d_nz, mat_opt._o_nz);
}
void PETScMatrix::setRowsColumnsZero(std::vector<PetscInt> const& row_pos)
@@ -101,6 +101,24 @@ void PETScMatrix::viewer(const std::string &file_name, const PetscViewerFormat v
}
+void PETScMatrix::create()
+{
+ MatCreate(PETSC_COMM_WORLD, &_A);
+ MatSetSizes(_A, _n_loc_rows, _n_loc_cols, _nrows, _ncols);
+
+ MatSetFromOptions(_A);
+}
+
+void PETScMatrix::config(const PetscInt d_nz, const PetscInt o_nz)
+{
+ // for a dense matrix: MatSeqAIJSetPreallocation(_A, d_nz, PETSC_NULL);
+ MatMPIAIJSetPreallocation(_A, d_nz, PETSC_NULL, o_nz, PETSC_NULL);
+
+ MatGetOwnershipRange(_A, &_start_rank, &_end_rank);
+ MatGetSize(_A, &_nrows, &_ncols);
+ MatGetLocalSize(_A, &_n_loc_rows, &_n_loc_cols);
+}
+
bool finalizeMatrixAssembly(PETScMatrix &mat, const MatAssemblyType asm_type)
{
mat.finalizeAssembly(asm_type);
diff --git a/MathLib/LinAlg/PETSc/PETScMatrix.h b/MathLib/LinAlg/PETSc/PETScMatrix.h
index 40eb50c94e5..83bf3a91336 100644
--- a/MathLib/LinAlg/PETSc/PETScMatrix.h
+++ b/MathLib/LinAlg/PETSc/PETScMatrix.h
@@ -27,53 +27,48 @@ namespace MathLib
{
/*!
- \brief Wrapper class for PETSc matrix routines for a square matrix for
- the system of linear equations.
+ \brief Wrapper class for PETSc matrix routines for matrix.
*/
class PETScMatrix
{
public:
/*!
- \brief Constructor for the matrix partitioning with default options:
- The size of the glabal size, the numbers of local rows
- and columns have the value of PETSC_DECIDE.
- \param size The dimension of the matrix.
+ \brief Constructor for a square matrix partitioning with default options:
+ The numbers of local rows and columns have the value of PETSC_DECIDE.
+ \param dim The dimension of the matrix.
*/
- explicit PETScMatrix(const PetscInt size);
+ explicit PETScMatrix(const PetscInt dim);
/*!
- \brief Constructor for the user determined matrix partitioning
- with two options.
- \param size The number of rows of the matrix or the local matrix.
- \param n_loc_cols Number of the columns of the rank.
- If it is -1, it is PETSC_DECIDE.
- \param is_global_size The flag of the type of vec_size,
- i.e. whether it is a global size or local size.
- The default is true.
+ \brief Constructor for a square matrix partitioning with more options
+ \param nrows The number of rows of the matrix or the local matrix.
+ \param mat_op The configuration information for creating a matrix.
*/
- PETScMatrix(const PetscInt size, const PetscInt n_loc_cols,
- const bool is_global_size);
+ PETScMatrix(const PetscInt nrows, const PETScMatrixOption &mat_op = PETScMatrixOption() );
/*!
- \brief Constructor for the user determined partitioning with more options
- \param size The number of rows of the matrix or the local matrix.
+ \brief Constructor for a rectangular matrix partitioning with default options:
+ The numbers of local rows and columns have the value of PETSC_DECIDE.
+ \param nrows_global The number of global rows.
+ \param ncols_global The number of global columns.
+ */
+ PETScMatrix(const PetscInt nrows_global, const PetscInt ncols_global);
+
+ /*!
+ \brief Constructor for a rectangular matrix partitioning with more options
+ \param nrows The number of global or local rows.
+ \param ncols The number of global or local columns.
\param mat_op The configuration information for creating a matrix.
*/
- PETScMatrix(const PetscInt size, const PETScMatrixOption &mat_op = PETScMatrixOption() );
+ PETScMatrix(const PetscInt nrows, const PetscInt ncols,
+ const PETScMatrixOption &mat_op = PETScMatrixOption() );
+
~PETScMatrix()
{
MatDestroy(&_A);
}
- /*!
- \brief Config memory allocation and set the related member data.
- Only called after an intance is created with the constructor
- wthout PETScMatrixOption argument.
- \param mat_op The configuration information for creating a matrix.
- */
- void config(const PETScMatrixOption &mat_op);
-
/*!
\brief Perform MPI collection of assembled entries in buffer
\param asm_type Assmebly type, either MAT_FLUSH_ASSEMBLY
@@ -88,26 +83,15 @@ class PETScMatrix
/// Get the number of rows.
PetscInt getNRows() const
{
- return _size;
+ return _nrows;
}
/// Get the number of columns.
PetscInt getNCols() const
{
- return _size;
- }
-
- /// Get the start global index of the rows of the same rank.
- PetscInt getRangeBegin() const
- {
- return _start_rank;
+ return _ncols;
}
- /// Get the end global index of the rows in the same rank.
- PetscInt getRangeEnd() const
- {
- return _end_rank;
- }
/// Get the number of local rows.
PetscInt getNLocalRows() const
@@ -121,6 +105,18 @@ class PETScMatrix
return _n_loc_cols;
}
+ /// Get the start global index of the rows of the same rank.
+ PetscInt getRangeBegin() const
+ {
+ return _start_rank;
+ }
+
+ /// Get the end global index of the rows in the same rank.
+ PetscInt getRangeEnd() const
+ {
+ return _end_rank;
+ }
+
/// Get matrix reference.
PETSc_Mat &getRawMatrix()
{
@@ -215,21 +211,36 @@ class PETScMatrix
/// PETSc matrix
PETSc_Mat _A;
- /// Dimension of matrix
- PetscInt _size;
+ /// Number of the global rows
+ PetscInt _nrows;
+
+ /// Number of the global columns
+ PetscInt _ncols;
+
/// Number of the local rows
PetscInt _n_loc_rows;
+
/// Number of the local columns
PetscInt _n_loc_cols;
/// Starting index in a rank
PetscInt _start_rank;
+
/// Ending index in a rank
PetscInt _end_rank;
/// Create the matrix
void create();
+ /*!
+ \brief Config memory allocation and set the related member data.
+ \param Number of nonzeros per row in the diagonal portion of local submatrix
+ (same value is used for all local rows),
+ \param Number of nonzeros per row in the off-diagonal portion of local submatrix
+ (same value is used for all local rows)
+ */
+ void config(const PetscInt d_nz, const PetscInt o_nz);
+
friend bool finalizeMatrixAssembly(PETScMatrix &mat, const MatAssemblyType asm_type);
};
diff --git a/Tests/MathLib/TestGlobalMatrixInterface.cpp b/Tests/MathLib/TestGlobalMatrixInterface.cpp
index bae80573d6d..48d789731cf 100644
--- a/Tests/MathLib/TestGlobalMatrixInterface.cpp
+++ b/Tests/MathLib/TestGlobalMatrixInterface.cpp
@@ -61,13 +61,17 @@ void checkGlobalMatrixInterfaceMPI(T_MATRIX &m, T_VECTOR &v)
MPI_Comm_rank(PETSC_COMM_WORLD, &mrank);
ASSERT_EQ(3u, msize);
- ASSERT_EQ(6u, m.getNRows());
ASSERT_EQ(m.getRangeEnd()-m.getRangeBegin(), m.getNLocalRows());
+ int gathered_rows;
+ int local_rows = m.getNLocalRows();
+ MPI_Allreduce(&local_rows, &gathered_rows, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
+ ASSERT_EQ(m.getNRows(), gathered_rows);
+
int gathered_cols;
int local_cols = m.getNLocalColumns();
MPI_Allreduce(&local_cols, &gathered_cols, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
- ASSERT_EQ(6u, gathered_cols);
+ ASSERT_EQ(m.getNCols(), gathered_cols);
// Add entries
MathLib::DenseMatrix<double> loc_m(2, 2);
@@ -106,6 +110,61 @@ void checkGlobalMatrixInterfaceMPI(T_MATRIX &m, T_VECTOR &v)
ASSERT_EQ(sqrt((2*3*3 + 8*8 + 3*7*7)), y.getNorm());
}
+
+// Rectanglular matrix
+template <class T_MATRIX, class T_VECTOR>
+void checkGlobalRectangularMatrixInterfaceMPI(T_MATRIX &m, T_VECTOR &v)
+{
+ int mrank;
+ MPI_Comm_rank(PETSC_COMM_WORLD, &mrank);
+
+ ASSERT_EQ(m.getRangeEnd()-m.getRangeBegin(), m.getNLocalRows());
+
+ int gathered_rows;
+ int local_rows = m.getNLocalRows();
+ MPI_Allreduce(&local_rows, &gathered_rows, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
+ ASSERT_EQ(m.getNRows(), gathered_rows);
+
+ int gathered_cols;
+ int local_cols = m.getNLocalColumns();
+ MPI_Allreduce(&local_cols, &gathered_cols, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
+ ASSERT_EQ(m.getNCols(), gathered_cols);
+
+ printf("m.getNRows()%d m.getNRows() %d %d", m.getNRows(), m.getNCols(), v.size() );
+
+ // Add entries
+ MathLib::DenseMatrix<double> loc_m(2, 3);
+ loc_m(0, 0) = 1.;
+ loc_m(0, 1) = 2.;
+ loc_m(0, 1) = 3.;
+ loc_m(1, 0) = 1.;
+ loc_m(1, 1) = 2.;
+ loc_m(1, 1) = 3.;
+
+ std::vector<int> row_pos(2);
+ std::vector<int> col_pos(3);
+ row_pos[0] = 2 * mrank;
+ row_pos[1] = 2 * mrank + 1;
+ col_pos[0] = 3 * mrank;
+ col_pos[1] = 3 * mrank + 1;
+ col_pos[2] = 3 * mrank + 2;
+
+ m.add(row_pos, col_pos, loc_m);
+
+ MathLib::finalizeMatrixAssembly(m);
+
+ /*
+ Cannot perform the following multiplication due to PETSc bug.
+ // Multiply by a vector
+ v = 1.;
+ const bool deep_copy = false;
+ T_VECTOR y(v, deep_copy);
+ m.multi(v, y);
+
+ ASSERT_EQ(sqrt(24.), y.getNorm());
+ */
+}
+
#endif // end of: ifdef USE_PETSC // or MPI
} // end namespace
@@ -152,20 +211,32 @@ TEST(Math, CheckInterface_PETScMatrix_Global_Size)
checkGlobalMatrixInterfaceMPI(A, x);
}
-TEST(Math, CheckInterface_PETScMatrix_Global_local_late_config)
+// Test rectangular matrix
+TEST(Math, CheckInterface_PETSc_Rectangular_Matrix_Local_Size)
{
MathLib::PETScMatrixOption opt;
- opt._d_nz = 2;
+ opt._d_nz = 3;
opt._o_nz = 0;
opt._is_global_size = false;
- opt._n_local_cols = 2;
- MathLib::PETScMatrix A(2, 2, false);
+ opt._n_local_cols = -1;
+ MathLib::PETScMatrix A(2, 3, opt);
const bool is_gloabal_size = false;
- MathLib::PETScVector x(2, is_gloabal_size);
+ MathLib::PETScVector x(3, is_gloabal_size);
- A.config(opt);
- checkGlobalMatrixInterfaceMPI(A, x);
+ checkGlobalRectangularMatrixInterfaceMPI(A, x);
+}
+
+TEST(Math, CheckInterface_PETSc_Rectangular_Matrix_Global_Size)
+{
+ MathLib::PETScMatrixOption opt;
+ opt._d_nz = 3;
+ opt._o_nz = 0;
+ MathLib::PETScMatrix A(6, 9, opt);
+
+ MathLib::PETScVector x(9);
+
+ checkGlobalRectangularMatrixInterfaceMPI(A, x);
}
#endif // end of: ifdef USE_PETSC // or MPI
--
GitLab