PETScMatrix.h

Go to the documentation of this file.

 
#pragma once
 
#include <string>
#include <vector>
 
#include "MathLib/LinAlg/RowColumnIndices.h"
#include "PETScMatrixOption.h"
#include "PETScVector.h"
 
typedef Mat PETSc_Mat;
 
namespace MathLib
{
class PETScMatrix
{
public:
    using IndexType = PetscInt;
 
public:
    PETScMatrix() {}
    PETScMatrix(const PetscInt nrows,
                const PETScMatrixOption& mat_op = PETScMatrixOption());
 
    PETScMatrix(const PetscInt nrows, const PetscInt ncols,
                const PETScMatrixOption& mat_op = PETScMatrixOption());
 
    ~PETScMatrix() { destroy(); }
    PETScMatrix(PETScMatrix const& A);
 
    PETScMatrix& operator=(PETScMatrix const& A);
 
    void finalizeAssembly(const MatAssemblyType asm_type = MAT_FINAL_ASSEMBLY)
    {
        MatAssemblyBegin(A_, asm_type);
        MatAssemblyEnd(A_, asm_type);
    }
 
    PetscInt getNumberOfRows() const { return nrows_; }
    PetscInt getNumberOfColumns() const { return ncols_; }
    PetscInt getNumberOfLocalRows() const { return n_loc_rows_; }
    PetscInt getNumberOfLocalColumns() const { return n_loc_cols_; }
    PetscInt getRangeBegin() const { return start_rank_; }
    PetscInt getRangeEnd() const { return end_rank_; }
    Mat& getRawMatrix() { return A_; }
    Mat const& getRawMatrix() const { return A_; }
    void setZero() { MatZeroEntries(A_); }
    void setRowsColumnsZero(std::vector<PetscInt> const& row_pos);
 
    void set(const PetscInt i, const PetscInt j, const PetscScalar value)
    {
        MatSetValue(A_, i, j, value, INSERT_VALUES);
    }
 
    void add(const PetscInt i, const PetscInt j, const PetscScalar value)
    {
        MatSetValue(A_, i, j, value, ADD_VALUES);
    }
 
    template <class T_DENSE_MATRIX>
    void add(RowColumnIndices<PetscInt> const& indices,
             const T_DENSE_MATRIX& sub_matrix)
    {
        // Set global column indices to positive to allow all entries of columns
        // to be added to the global matrix. For the ghost columns, only the
        // off diagonal entries are added due to the negative indices of the
        // corresponding rows.
        std::vector<PetscInt> cols;
        cols.reserve(indices.columns.size());
        for (auto col : indices.columns)
        {
            // Ghost entries, and its original index is 0.
            if (col == -ncols_)
            {
                cols.push_back(0);
            }
            else
            {
                cols.push_back(std::abs(col));
            }
        }
 
        add(indices.rows, cols, sub_matrix);
    }
 
    template <class T_DENSE_MATRIX>
    void add(std::vector<PetscInt> const& row_pos,
             std::vector<PetscInt> const& col_pos,
             const T_DENSE_MATRIX& sub_mat);
 
    void viewer(const std::string& file_name,
                const PetscViewerFormat vw_format = PETSC_VIEWER_ASCII_MATLAB);
 
private:
    void destroy()
    {
        if (A_ != nullptr)
        {
            MatDestroy(&A_);
        }
        A_ = nullptr;
    }
 
    Mat A_ = nullptr;
 
    PetscInt nrows_;
 
    PetscInt ncols_;
 
    PetscInt n_loc_rows_;
 
    PetscInt n_loc_cols_;
 
    PetscInt start_rank_;
 
    PetscInt end_rank_;
 
    void create(const PetscInt d_nz, const PetscInt o_nz);
 
    friend bool finalizeMatrixAssembly(PETScMatrix& mat,
                                       const MatAssemblyType asm_type);
};
 
template <class T_DENSE_MATRIX>
void PETScMatrix::add(std::vector<PetscInt> const& row_pos,
                      std::vector<PetscInt> const& col_pos,
                      const T_DENSE_MATRIX& sub_mat)
{
    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);
};
 
bool finalizeMatrixAssembly(
    PETScMatrix& mat, const MatAssemblyType asm_type = MAT_FINAL_ASSEMBLY);
 
}  // namespace MathLib