[MaL] PETSc matrix stores pointer to matrix

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)

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);
 };
 
 /*!