From 43d8d8cb7c944f1eb5bde059b5704608fdf1d30a Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Tue, 8 Mar 2016 11:18:45 +0100
Subject: [PATCH] [MaL] PETSc matrix stores pointer to matrix
---
MathLib/LinAlg/PETSc/PETScMatrix.cpp | 69 +++++++++++++++++++++-------
MathLib/LinAlg/PETSc/PETScMatrix.h | 44 ++++++++++--------
2 files changed, 79 insertions(+), 34 deletions(-)
diff --git a/MathLib/LinAlg/PETSc/PETScMatrix.cpp b/MathLib/LinAlg/PETSc/PETScMatrix.cpp
index d04845c6d4a..a0c56f6b21b 100644
--- a/MathLib/LinAlg/PETSc/PETScMatrix.cpp
+++ b/MathLib/LinAlg/PETSc/PETScMatrix.cpp
@@ -48,6 +48,41 @@ PETScMatrix::PETScMatrix (const PetscInt nrows, const PetscInt ncols, const PETS
create(mat_opt.d_nz, mat_opt.o_nz);
}
+PETScMatrix::PETScMatrix(const PETScMatrix &A)
+ : _A(new PETSc_Mat)
+ , _nrows(A._nrows)
+ , _ncols(A._ncols)
+ , _n_loc_rows(A._n_loc_rows)
+ , _n_loc_cols(A._n_loc_cols)
+ , _start_rank(A._start_rank)
+ , _end_rank(A._end_rank)
+{
+ _A.reset(new PETSc_Mat);
+ MatConvert(*A._A, MATSAME, MAT_INITIAL_MATRIX, _A.get());
+}
+
+PETScMatrix&
+PETScMatrix::operator=(PETScMatrix const& A)
+{
+ _nrows = A._nrows;
+ _ncols = A._ncols;
+ _n_loc_rows = A._n_loc_rows;
+ _n_loc_cols = A._n_loc_cols;
+ _start_rank = A._start_rank;
+ _end_rank = A._end_rank;
+
+ if (_A) {
+ // TODO this is the slowest option for copying
+ MatCopy(*A._A, *_A, DIFFERENT_NONZERO_PATTERN);
+ } else {
+ destroy();
+ _A.reset(new PETSc_Mat);
+ MatConvert(*A._A, MATSAME, MAT_INITIAL_MATRIX, _A.get());
+ }
+
+ return *this;
+}
+
void PETScMatrix::setRowsColumnsZero(std::vector<PetscInt> const& row_pos)
{
// Each rank (compute core) processes only the rows that belong to the rank itself.
@@ -58,11 +93,11 @@ void PETScMatrix::setRowsColumnsZero(std::vector<PetscInt> const& row_pos)
// This avoids all reductions in the zero row routines
// and thus improves performance for very large process counts.
// See PETSc doc about MAT_NO_OFF_PROC_ZERO_ROWS.
- MatSetOption(_A, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
+ MatSetOption(*_A, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
if(nrows>0)
- MatZeroRows(_A, nrows, &row_pos[0], one, PETSC_NULL, PETSC_NULL);
+ MatZeroRows(*_A, nrows, &row_pos[0], one, PETSC_NULL, PETSC_NULL);
else
- MatZeroRows(_A, 0, PETSC_NULL, one, PETSC_NULL, PETSC_NULL);
+ MatZeroRows(*_A, 0, PETSC_NULL, one, PETSC_NULL, PETSC_NULL);
}
void PETScMatrix::viewer(const std::string &file_name, const PetscViewerFormat vw_format)
@@ -73,13 +108,13 @@ void PETScMatrix::viewer(const std::string &file_name, const PetscViewerFormat v
finalizeAssembly();
- PetscObjectSetName((PetscObject)_A,"Stiffness_matrix");
- MatView(_A,viewer);
+ PetscObjectSetName((PetscObject)*_A,"Stiffness_matrix");
+ MatView(*_A,viewer);
// This preprocessor is only for debugging, e.g. dump the matrix and exit the program.
//#define EXIT_TEST
#ifdef EXIT_TEST
- MatDestroy(&_A);
+ MatDestroy(_A);
PetscFinalize();
exit(0);
#endif
@@ -88,20 +123,22 @@ void PETScMatrix::viewer(const std::string &file_name, const PetscViewerFormat v
void PETScMatrix::create(const PetscInt d_nz, const PetscInt o_nz)
{
- MatCreate(PETSC_COMM_WORLD, &_A);
- MatSetSizes(_A, _n_loc_rows, _n_loc_cols, _nrows, _ncols);
+ _A.reset(new PETSc_Mat);
+
+ MatCreate(PETSC_COMM_WORLD, _A.get());
+ MatSetSizes(*_A, _n_loc_rows, _n_loc_cols, _nrows, _ncols);
- MatSetFromOptions(_A);
+ MatSetFromOptions(*_A);
- MatSetType(_A, MATMPIAIJ);
- MatSeqAIJSetPreallocation(_A, d_nz, PETSC_NULL);
- MatMPIAIJSetPreallocation(_A, d_nz, PETSC_NULL, o_nz, PETSC_NULL);
- // If pre-allocation does not work one can use MatSetUp(_A), which is much
+ MatSetType(*_A, MATMPIAIJ);
+ MatSeqAIJSetPreallocation(*_A, d_nz, PETSC_NULL);
+ MatMPIAIJSetPreallocation(*_A, d_nz, PETSC_NULL, o_nz, PETSC_NULL);
+ // If pre-allocation does not work one can use MatSetUp(*_A), which is much
// slower.
- MatGetOwnershipRange(_A, &_start_rank, &_end_rank);
- MatGetSize(_A, &_nrows, &_ncols);
- MatGetLocalSize(_A, &_n_loc_rows, &_n_loc_cols);
+ 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)
diff --git a/MathLib/LinAlg/PETSc/PETScMatrix.h b/MathLib/LinAlg/PETSc/PETScMatrix.h
index 091facf9eb8..d063d585037 100644
--- a/MathLib/LinAlg/PETSc/PETScMatrix.h
+++ b/MathLib/LinAlg/PETSc/PETScMatrix.h
@@ -15,6 +15,7 @@
#ifndef PETSCMATRIX_H_
#define PETSCMATRIX_H_
+#include <memory>
#include <string>
#include <vector>
@@ -37,10 +38,7 @@ class PETScMatrix
using IndexType = PetscInt;
public:
- // TODO preliminary
- PETScMatrix() {
- // TODO implement
- }
+ PETScMatrix() {}
/*!
\brief Constructor for a square matrix partitioning with more options
@@ -60,9 +58,13 @@ class PETScMatrix
~PETScMatrix()
{
- MatDestroy(&_A);
+ destroy();
}
+ PETScMatrix(PETScMatrix const& A);
+
+ PETScMatrix& operator=(PETScMatrix const& A);
+
/*!
\brief Perform MPI collection of assembled entries in buffer
\param asm_type Assmebly type, either MAT_FLUSH_ASSEMBLY
@@ -70,8 +72,8 @@ class PETScMatrix
*/
void finalizeAssembly(const MatAssemblyType asm_type = MAT_FINAL_ASSEMBLY)
{
- MatAssemblyBegin(_A, asm_type);
- MatAssemblyEnd(_A, asm_type);
+ MatAssemblyBegin(*_A, asm_type);
+ MatAssemblyEnd(*_A, asm_type);
}
/// Get the number of rows.
@@ -114,20 +116,24 @@ class PETScMatrix
/// Get matrix reference.
PETSc_Mat &getRawMatrix()
{
- return _A;
+ return *_A;
}
- // TODO preliminary
- // this method is dangerous insofar you can do arbitrary things also
- // with a const PETSc matrix.
+
+ /*! Get a matrix reference.
+ *
+ * \warning
+ * This method is dangerous insofar as you can do arbitrary things also
+ * with a const PETSc matrix.
+ */
PETSc_Mat const& getRawMatrix() const
{
- return _A;
+ return *_A;
}
/// Set all entries to zero.
void setZero()
{
- MatZeroEntries(_A);
+ MatZeroEntries(*_A);
}
/*
@@ -149,7 +155,7 @@ class PETScMatrix
*/
void multiply(const PETScVector &vec, PETScVector &vec_r)
{
- MatMult(_A, vec.getData(), vec_r.getData() );
+ MatMult(*_A, vec.getData(), vec_r.getData() );
}
/*!
@@ -160,7 +166,7 @@ class PETScMatrix
*/
void set(const PetscInt i, const PetscInt j, const PetscScalar value)
{
- MatSetValue(_A, i, j, value, INSERT_VALUES);
+ MatSetValue(*_A, i, j, value, INSERT_VALUES);
}
/*!
@@ -171,7 +177,7 @@ class PETScMatrix
*/
void add(const PetscInt i, const PetscInt j, const PetscScalar value)
{
- MatSetValue(_A, i, j, value, ADD_VALUES);
+ MatSetValue(*_A, i, j, value, ADD_VALUES);
}
/// Add sub-matrix at positions given by \c indices.
@@ -231,8 +237,10 @@ class PETScMatrix
const PetscViewerFormat vw_format = PETSC_VIEWER_ASCII_MATLAB );
private:
+ void destroy() { if (_A) MatDestroy(_A.get()); _A.reset(nullptr); }
+
/// PETSc matrix
- PETSc_Mat _A;
+ std::unique_ptr<PETSc_Mat> _A;
/// Number of the global rows
PetscInt _nrows;
@@ -278,7 +286,7 @@ void PETScMatrix::add(std::vector<PetscInt> const& row_pos,
const PetscInt nrows = static_cast<PetscInt> (row_pos.size());
const PetscInt ncols = static_cast<PetscInt> (col_pos.size());
- MatSetValues(_A, nrows, &row_pos[0], ncols, &col_pos[0], &sub_mat(0,0), ADD_VALUES);
+ MatSetValues(*_A, nrows, &row_pos[0], ncols, &col_pos[0], &sub_mat(0,0), ADD_VALUES);
};
/*!
--
GitLab