diff --git a/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp b/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
index a563e4e47c9f9cdb117f7aed1ba991216c3f9853..566d775d5c565a6f5978421f5a1012c75468d05c 100644
--- a/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
+++ b/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
@@ -14,7 +14,6 @@
http://www.opengeosys.org/project/license
*/
-
#include "PETScLinearSolver.h"
#include "BaseLib/RunTime.h"
#include "MathLib/LinAlg/LinearSolverOptions.h"
@@ -24,13 +23,15 @@ namespace MathLib
PETScLinearSolver::PETScLinearSolver(const std::string /*prefix*/,
BaseLib::ConfigTree const* const option)
{
- // Insert options into petsc database. Default options are given in the string below.
- std::string petsc_options
- = "-ksp_type cg -pc_type bjacobi -ksp_rtol 1e-16 -ksp_max_it 10000";
+ // Insert options into petsc database. Default options are given in the
+ // string below.
+ std::string petsc_options =
+ "-ksp_type cg -pc_type bjacobi -ksp_rtol 1e-16 -ksp_max_it 10000";
std::string prefix;
- if (option) {
+ if (option)
+ {
ignoreOtherLinearSolvers(*option, "petsc");
//! \ogs_file_param{prj__linear_solvers__linear_solver__petsc}
@@ -38,18 +39,21 @@ PETScLinearSolver::PETScLinearSolver(const std::string /*prefix*/,
{
if (auto const parameters =
//! \ogs_file_param{prj__linear_solvers__linear_solver__petsc__parameters}
- subtree->getConfigParameterOptional<std::string>("parameters")) {
+ subtree->getConfigParameterOptional<std::string>("parameters"))
+ {
petsc_options = *parameters;
}
- //! \ogs_file_param{prj__linear_solvers__linear_solver__petsc__prefix}
- if (auto const pre = subtree->getConfigParameterOptional<std::string>("prefix")) {
+ if (auto const pre =
+ //! \ogs_file_param{prj__linear_solvers__linear_solver__petsc__prefix}
+ subtree->getConfigParameterOptional<std::string>("prefix"))
+ {
if (!pre->empty())
prefix = *pre + "_";
}
}
}
-#if PETSC_VERSION_LT(3,7,0)
+#if PETSC_VERSION_LT(3, 7, 0)
PetscOptionsInsertString(petsc_options.c_str());
#else
PetscOptionsInsertString(nullptr, petsc_options.c_str());
@@ -59,14 +63,15 @@ PETScLinearSolver::PETScLinearSolver(const std::string /*prefix*/,
KSPGetPC(_solver, &_pc);
- if (!prefix.empty()) {
+ if (!prefix.empty())
+ {
KSPSetOptionsPrefix(_solver, prefix.c_str());
}
- KSPSetFromOptions(_solver); // set running time option
+ KSPSetFromOptions(_solver); // set running time option
}
-bool PETScLinearSolver::solve(PETScMatrix& A, PETScVector &b, PETScVector &x)
+bool PETScLinearSolver::solve(PETScMatrix& A, PETScVector& b, PETScVector& x)
{
BaseLib::RunTime wtimer;
wtimer.start();
@@ -80,7 +85,8 @@ bool PETScLinearSolver::solve(PETScMatrix& A, PETScVector &b, PETScVector &x)
#if (PETSC_VERSION_NUMBER > 3040)
KSPSetOperators(_solver, A.getRawMatrix(), A.getRawMatrix());
#else
- KSPSetOperators(_solver, A.getRawMatrix(), A.getRawMatrix(), DIFFERENT_NONZERO_PATTERN);
+ KSPSetOperators(_solver, A.getRawMatrix(), A.getRawMatrix(),
+ DIFFERENT_NONZERO_PATTERN);
#endif
KSPSolve(_solver, b.getRawVector(), x.getRawVector());
@@ -89,58 +95,76 @@ bool PETScLinearSolver::solve(PETScMatrix& A, PETScVector &b, PETScVector &x)
KSPGetConvergedReason(_solver, &reason);
bool converged = true;
- if(reason > 0)
+ if (reason > 0)
{
- const char *ksp_type;
- const char *pc_type;
+ const char* ksp_type;
+ const char* pc_type;
KSPGetType(_solver, &ksp_type);
PCGetType(_pc, &pc_type);
- PetscPrintf(PETSC_COMM_WORLD,"\n================================================");
- PetscPrintf(PETSC_COMM_WORLD, "\nLinear solver %s with %s preconditioner",
- ksp_type, pc_type);
+ PetscPrintf(PETSC_COMM_WORLD,
+ "\n================================================");
+ PetscPrintf(PETSC_COMM_WORLD,
+ "\nLinear solver %s with %s preconditioner", ksp_type,
+ pc_type);
PetscInt its;
KSPGetIterationNumber(_solver, &its);
- PetscPrintf(PETSC_COMM_WORLD,"\nconverged in %d iterations.", its);
- PetscPrintf(PETSC_COMM_WORLD,"\n================================================\n");
+ PetscPrintf(PETSC_COMM_WORLD, "\nconverged in %d iterations.", its);
+ PetscPrintf(PETSC_COMM_WORLD,
+ "\n================================================\n");
}
- else if(reason == KSP_DIVERGED_ITS)
+ else if (reason == KSP_DIVERGED_ITS)
{
- const char *ksp_type;
- const char *pc_type;
+ const char* ksp_type;
+ const char* pc_type;
KSPGetType(_solver, &ksp_type);
PCGetType(_pc, &pc_type);
- PetscPrintf(PETSC_COMM_WORLD, "\nLinear solver %s with %s preconditioner",
- ksp_type, pc_type);
- PetscPrintf(PETSC_COMM_WORLD, "\nWarning: maximum number of iterations reached.\n");
+ PetscPrintf(PETSC_COMM_WORLD,
+ "\nLinear solver %s with %s preconditioner", ksp_type,
+ pc_type);
+ PetscPrintf(PETSC_COMM_WORLD,
+ "\nWarning: maximum number of iterations reached.\n");
}
else
{
converged = false;
- if(reason == KSP_DIVERGED_INDEFINITE_PC)
+ if (reason == KSP_DIVERGED_INDEFINITE_PC)
{
- PetscPrintf(PETSC_COMM_WORLD, "\nDivergence because of indefinite preconditioner,");
- PetscPrintf(PETSC_COMM_WORLD, "\nTry to run again with -pc_factor_shift_positive_definite option.\n");
+ PetscPrintf(PETSC_COMM_WORLD,
+ "\nDivergence because of indefinite preconditioner,");
+ PetscPrintf(PETSC_COMM_WORLD,
+ "\nTry to run again with "
+ "-pc_factor_shift_positive_definite option.\n");
}
- else if(reason == KSP_DIVERGED_BREAKDOWN_BICG)
+ else if (reason == KSP_DIVERGED_BREAKDOWN_BICG)
{
- PetscPrintf(PETSC_COMM_WORLD, "\nKSPBICG method was detected so the method could not continue to enlarge the Krylov space.");
- PetscPrintf(PETSC_COMM_WORLD, "\nTry to run again with another solver.\n");
+ PetscPrintf(PETSC_COMM_WORLD,
+ "\nKSPBICG method was detected so the method could not "
+ "continue to enlarge the Krylov space.");
+ PetscPrintf(PETSC_COMM_WORLD,
+ "\nTry to run again with another solver.\n");
}
- else if(reason == KSP_DIVERGED_NONSYMMETRIC)
+ else if (reason == KSP_DIVERGED_NONSYMMETRIC)
{
- PetscPrintf(PETSC_COMM_WORLD, "\nMatrix or preconditioner is unsymmetric but KSP requires symmetric.\n");
+ PetscPrintf(PETSC_COMM_WORLD,
+ "\nMatrix or preconditioner is unsymmetric but KSP "
+ "requires symmetric.\n");
}
else
{
- PetscPrintf(PETSC_COMM_WORLD, "\nDivergence detected, use command option -ksp_monitor or -log_summary to check the details.\n");
+ PetscPrintf(PETSC_COMM_WORLD,
+ "\nDivergence detected, use command option "
+ "-ksp_monitor or -log_summary to check the details.\n");
}
}
#ifdef TEST_MEM_PETSC
PetscMemoryGetCurrentUsage(&mem2);
- PetscPrintf(PETSC_COMM_WORLD, "###Memory usage by solver. Before: %f After: %f Increase: %d\n", mem1, mem2, (int)(mem2 - mem1));
+ PetscPrintf(
+ PETSC_COMM_WORLD,
+ "###Memory usage by solver. Before: %f After: %f Increase: %d\n", mem1,
+ mem2, (int)(mem2 - mem1));
#endif
_elapsed_ctime += wtimer.elapsed();
@@ -148,5 +172,4 @@ bool PETScLinearSolver::solve(PETScMatrix& A, PETScVector &b, PETScVector &x)
return converged;
}
-} //end of namespace
-
+} // end of namespace
diff --git a/MathLib/LinAlg/PETSc/PETScLinearSolver.h b/MathLib/LinAlg/PETSc/PETScLinearSolver.h
index d3c69583775aef946377df8abe3a5f4ac04c6d01..c22f00759b9d8a71563ac05bfa05a599d2e9ce73 100644
--- a/MathLib/LinAlg/PETSc/PETScLinearSolver.h
+++ b/MathLib/LinAlg/PETSc/PETScLinearSolver.h
@@ -16,7 +16,7 @@
#pragma once
-#include<string>
+#include <string>
#include <petscksp.h>
@@ -29,7 +29,6 @@
namespace MathLib
{
-
/*! A class of linear solver based on PETSc rountines.
All command-line options that are not recognized by OGS are passed on to
@@ -39,46 +38,37 @@ namespace MathLib
*/
class PETScLinearSolver
{
- public:
-
- /*!
- Constructor
- \param prefix Name used to distinguish the options in the command
- line for this solver. It can be the name of the PDE
- that owns an instance of this class.
- \param option Petsc options, which will be inserted into the global
- petsc options database.
- */
- PETScLinearSolver(const std::string prefix,
- BaseLib::ConfigTree const* const option);
-
- ~PETScLinearSolver()
- {
- KSPDestroy(&_solver);
- }
-
- // TODO check if some args in LinearSolver interface can be made const&.
- bool solve(PETScMatrix& A, PETScVector &b, PETScVector &x);
-
- /// Get number of iterations.
- PetscInt getNumberOfIterations() const
- {
- PetscInt its = 0;
- KSPGetIterationNumber(_solver, &its);
- return its;
- }
-
- /// Get elapsed wall clock time.
- double getElapsedTime() const
- {
- return _elapsed_ctime;
- }
-
- private:
- KSP _solver; ///< Solver type.
- PC _pc; ///< Preconditioner type.
-
- double _elapsed_ctime = 0.0; ///< Clock time
+public:
+ /*!
+ Constructor
+ \param prefix Name used to distinguish the options in the command
+ line for this solver. It can be the name of the PDE
+ that owns an instance of this class.
+ \param option Petsc options, which will be inserted into the global
+ petsc options database.
+ */
+ PETScLinearSolver(const std::string prefix,
+ BaseLib::ConfigTree const* const option);
+
+ ~PETScLinearSolver() { KSPDestroy(&_solver); }
+ // TODO check if some args in LinearSolver interface can be made const&.
+ bool solve(PETScMatrix& A, PETScVector& b, PETScVector& x);
+
+ /// Get number of iterations.
+ PetscInt getNumberOfIterations() const
+ {
+ PetscInt its = 0;
+ KSPGetIterationNumber(_solver, &its);
+ return its;
+ }
+
+ /// Get elapsed wall clock time.
+ double getElapsedTime() const { return _elapsed_ctime; }
+private:
+ KSP _solver; ///< Solver type.
+ PC _pc; ///< Preconditioner type.
+
+ double _elapsed_ctime = 0.0; ///< Clock time
};
-} // end namespace
+} // end namespace
diff --git a/MathLib/LinAlg/PETSc/PETScMatrix.cpp b/MathLib/LinAlg/PETSc/PETScMatrix.cpp
index 4a35c649e54206f52f46dc578c24f14cf71ec15e..4da7f29ab02113abe0e9abadcdd097c6b13590e4 100644
--- a/MathLib/LinAlg/PETSc/PETScMatrix.cpp
+++ b/MathLib/LinAlg/PETSc/PETScMatrix.cpp
@@ -1,6 +1,7 @@
/*!
\file PETScMatrix.cpp
- \brief Definition of member functions of class PETScMatrix, which provides an interface to
+ \brief Definition of member functions of class PETScMatrix, which provides an
+ interface to
PETSc matrix routines.
\author Wenqing Wang
@@ -17,12 +18,13 @@
namespace MathLib
{
-
-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)
+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)
+ if (!mat_opt.is_global_size)
{
_n_loc_rows = nrows;
_n_loc_cols = nrows;
@@ -33,11 +35,14 @@ PETScMatrix::PETScMatrix (const PetscInt nrows, const PETScMatrixOption &mat_opt
create(mat_opt.d_nz, mat_opt.o_nz);
}
-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)
+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)
{
- if(!mat_opt.is_global_size)
+ if (!mat_opt.is_global_size)
{
_nrows = PETSC_DECIDE;
_ncols = PETSC_DECIDE;
@@ -48,19 +53,18 @@ PETScMatrix::PETScMatrix (const PetscInt nrows, const PetscInt ncols, const PETS
create(mat_opt.d_nz, mat_opt.o_nz);
}
-PETScMatrix::PETScMatrix(const PETScMatrix &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)
+PETScMatrix::PETScMatrix(const PETScMatrix& 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)
{
MatConvert(A._A, MATSAME, MAT_INITIAL_MATRIX, &_A);
}
-PETScMatrix&
-PETScMatrix::operator=(PETScMatrix const& A)
+PETScMatrix& PETScMatrix::operator=(PETScMatrix const& A)
{
_nrows = A._nrows;
_ncols = A._ncols;
@@ -69,10 +73,13 @@ PETScMatrix::operator=(PETScMatrix const& A)
_start_rank = A._start_rank;
_end_rank = A._end_rank;
- if (_A) {
+ if (_A)
+ {
// TODO this is the slowest option for copying
MatCopy(A._A, _A, DIFFERENT_NONZERO_PATTERN);
- } else {
+ }
+ else
+ {
destroy();
MatConvert(A._A, MATSAME, MAT_INITIAL_MATRIX, &_A);
}
@@ -82,9 +89,10 @@ PETScMatrix::operator=(PETScMatrix const& A)
void PETScMatrix::setRowsColumnsZero(std::vector<PetscInt> const& row_pos)
{
- // Each rank (compute core) processes only the rows that belong to the rank itself.
+ // Each rank (compute core) processes only the rows that belong to the rank
+ // itself.
const PetscScalar one = 1.0;
- const PetscInt nrows = static_cast<PetscInt> (row_pos.size());
+ const PetscInt nrows = static_cast<PetscInt>(row_pos.size());
// Each process will only zero its own rows.
// This avoids all reductions in the zero row routines
@@ -93,15 +101,16 @@ void PETScMatrix::setRowsColumnsZero(std::vector<PetscInt> const& row_pos)
MatSetOption(_A, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
// Keep the non-zero pattern for the assignment operator.
- MatSetOption(_A, MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);
+ MatSetOption(_A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
- if(nrows>0)
+ if (nrows > 0)
MatZeroRows(_A, nrows, &row_pos[0], one, PETSC_NULL, PETSC_NULL);
else
MatZeroRows(_A, 0, PETSC_NULL, one, PETSC_NULL, PETSC_NULL);
}
-void PETScMatrix::viewer(const std::string &file_name, const PetscViewerFormat vw_format)
+void PETScMatrix::viewer(const std::string& file_name,
+ const PetscViewerFormat vw_format)
{
PetscViewer viewer;
PetscViewerASCIIOpen(PETSC_COMM_WORLD, file_name.c_str(), &viewer);
@@ -109,17 +118,17 @@ 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.
+// This preprocessor is only for debugging, e.g. dump the matrix and exit the
+// program.
//#define EXIT_TEST
#ifdef EXIT_TEST
MatDestroy(_A);
PetscFinalize();
exit(0);
#endif
-
}
void PETScMatrix::create(const PetscInt d_nz, const PetscInt o_nz)
@@ -136,15 +145,14 @@ void PETScMatrix::create(const PetscInt d_nz, const PetscInt o_nz)
// slower.
MatGetOwnershipRange(_A, &_start_rank, &_end_rank);
- MatGetSize(_A, &_nrows, &_ncols);
+ MatGetSize(_A, &_nrows, &_ncols);
MatGetLocalSize(_A, &_n_loc_rows, &_n_loc_cols);
}
-bool finalizeMatrixAssembly(PETScMatrix &mat, const MatAssemblyType asm_type)
+bool finalizeMatrixAssembly(PETScMatrix& mat, const MatAssemblyType asm_type)
{
mat.finalizeAssembly(asm_type);
return true;
}
-} //end of namespace
-
+} // end of namespace
diff --git a/MathLib/LinAlg/PETSc/PETScMatrix.h b/MathLib/LinAlg/PETSc/PETScMatrix.h
index 90fc96ac00c967d276a0a061aba22e6e976695dd..d388febec0375723b7f33784e51738d1dda7d015 100644
--- a/MathLib/LinAlg/PETSc/PETScMatrix.h
+++ b/MathLib/LinAlg/PETSc/PETScMatrix.h
@@ -26,232 +26,210 @@ typedef Mat PETSc_Mat;
namespace MathLib
{
-
/*!
\brief Wrapper class for PETSc matrix routines for matrix.
*/
class PETScMatrix
{
- public:
- using IndexType = PetscInt;
-
- public:
- PETScMatrix() {}
-
- /*!
- \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 nrows, const PETScMatrixOption &mat_op = PETScMatrixOption() );
-
- /*!
- \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 nrows, const PetscInt ncols,
- const PETScMatrixOption &mat_op = PETScMatrixOption() );
-
- ~PETScMatrix()
- {
- 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
- or MAT_FINAL_ASSEMBLY
- */
- void finalizeAssembly(const MatAssemblyType asm_type = MAT_FINAL_ASSEMBLY)
- {
- MatAssemblyBegin(_A, asm_type);
- MatAssemblyEnd(_A, asm_type);
- }
-
- /// Get the number of rows.
- PetscInt getNumberOfRows() const
- {
- return _nrows;
- }
-
- /// Get the number of columns.
- PetscInt getNumberOfColumns() const
- {
- return _ncols;
- }
-
-
- /// Get the number of local rows.
- PetscInt getNumberOfLocalRows() const
- {
- return _n_loc_rows;
- }
-
- /// Get the number of local columns.
- PetscInt getNumberOfLocalColumns() const
- {
- 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.
- Mat& getRawMatrix() { return _A; }
-
- /*! Get a matrix reference.
- *
- * \warning
- * This method is dangerous insofar as you can do arbitrary things also
- * with a const PETSc matrix.
- */
- Mat const& getRawMatrix() const { return _A; }
-
- /// Set all entries to zero.
- void setZero()
- {
- MatZeroEntries(_A);
- }
-
- /*
- \brief Set the specified rows to zero except diagonal entries, i.e.
- \f$A(k, j) = \begin{cases}
- 0.0, &j\not=k, j=1,2,\dots,k-1, k+1, \dots, n
- 1.0, &j = k
- \end{cases}f$, where \f$k \in \mbox{row\_pos}\f$
- This function must be called by all ranks.
- \param row_pos The row indicies of the specified rows.
- */
- void setRowsColumnsZero(std::vector<PetscInt> const& row_pos);
-
- /*!
- \brief Set a single entry with a value.
- \param i The row index.
- \param j The column index.
- \param value The entry value.
- */
- void set(const PetscInt i, const PetscInt j, const PetscScalar value)
- {
- MatSetValue(_A, i, j, value, INSERT_VALUES);
- }
-
- /*!
- \brief Add value to a single entry.
- \param i The row index.
- \param j The column index.
- \param value The entry value.
- */
- void add(const PetscInt i, const PetscInt j, const PetscScalar value)
- {
- MatSetValue(_A, i, j, value, ADD_VALUES);
- }
-
- /// Add sub-matrix at positions given by \c indices.
- template<class T_DENSE_MATRIX>
- void add(RowColumnIndices<PetscInt> const& indices,
- const T_DENSE_MATRIX &sub_matrix)
+public:
+ using IndexType = PetscInt;
+
+public:
+ PETScMatrix() {}
+ /*!
+ \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 nrows,
+ const PETScMatrixOption& mat_op = PETScMatrixOption());
+
+ /*!
+ \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 nrows, const PetscInt ncols,
+ const PETScMatrixOption& mat_op = PETScMatrixOption());
+
+ ~PETScMatrix() { 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
+ or MAT_FINAL_ASSEMBLY
+ */
+ void finalizeAssembly(const MatAssemblyType asm_type = MAT_FINAL_ASSEMBLY)
+ {
+ MatAssemblyBegin(_A, asm_type);
+ MatAssemblyEnd(_A, asm_type);
+ }
+
+ /// Get the number of rows.
+ PetscInt getNumberOfRows() const { return _nrows; }
+ /// Get the number of columns.
+ PetscInt getNumberOfColumns() const { return _ncols; }
+ /// Get the number of local rows.
+ PetscInt getNumberOfLocalRows() const { return _n_loc_rows; }
+ /// Get the number of local columns.
+ PetscInt getNumberOfLocalColumns() const { 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.
+ Mat& getRawMatrix() { return _A; }
+ /*! Get a matrix reference.
+ *
+ * \warning
+ * This method is dangerous insofar as you can do arbitrary things also
+ * with a const PETSc matrix.
+ */
+ Mat const& getRawMatrix() const { return _A; }
+ /// Set all entries to zero.
+ void setZero() { MatZeroEntries(_A); }
+ /*
+ \brief Set the specified rows to zero except diagonal entries, i.e.
+ \f$A(k, j) = \begin{cases}
+ 0.0, &j\not=k, j=1,2,\dots,k-1, k+1, \dots, n
+ 1.0, &j = k
+ \end{cases}f$, where \f$k \in \mbox{row\_pos}\f$
+ This function must be called by all ranks.
+ \param row_pos The row indicies of the specified rows.
+ */
+ void setRowsColumnsZero(std::vector<PetscInt> const& row_pos);
+
+ /*!
+ \brief Set a single entry with a value.
+ \param i The row index.
+ \param j The column index.
+ \param value The entry value.
+ */
+ void set(const PetscInt i, const PetscInt j, const PetscScalar value)
+ {
+ MatSetValue(_A, i, j, value, INSERT_VALUES);
+ }
+
+ /*!
+ \brief Add value to a single entry.
+ \param i The row index.
+ \param j The column index.
+ \param value The entry value.
+ */
+ void add(const PetscInt i, const PetscInt j, const PetscScalar value)
+ {
+ MatSetValue(_A, i, j, value, ADD_VALUES);
+ }
+
+ /// Add sub-matrix at positions given by \c indices.
+ template <class T_DENSE_MATRIX>
+ void add(RowColumnIndices<PetscInt> const& indices,
+ const T_DENSE_MATRIX& sub_matrix)
+ {
+ std::vector<PetscInt> cols;
+ cols.reserve(indices.columns.size());
+ for (auto col : indices.columns)
{
- 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);
+ // Ghost entries, and its original index is 0.
+ if (col == -_ncols)
+ cols.push_back(0);
+ else
+ cols.push_back(std::abs(col));
}
- /*!
- \brief Add a submatrix to this.
- \param row_pos The row indices of the entries of the submatrix.
- \param col_pos The column indices of the entries of the submatrix.
- \param sub_mat A dense matrix to be added on.
- */
- 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 );
-
- /*! View the global vector for test purpose. Do not use it for output a big vector.
- \param file_name File name for output
- \param vw_format File format listed as:
- PETSC_VIEWER_DEFAULT Default format
- PETSC_VIEWER_ASCII_MATLAB MATLAB format
- PETSC_VIEWER_ASCII_DENSE Print matrix as dense
- PETSC_VIEWER_ASCII_IMPL Implementation-specific format
- (which is in many cases the same as the default)
- PETSC_VIEWER_ASCII_INFO Basic information about object
- PETSC_VIEWER_ASCII_INFO_DETAIL More detailed info about object
- PETSC_VIEWER_ASCII_COMMON Identical output format for all objects of a particular type
- PETSC_VIEWER_ASCII_INDEX (for vectors) Prints the vector element number next to
- each vector entry
- PETSC_VIEWER_ASCII_SYMMODU Print parallel vectors without indicating the processor ranges
- PETSC_VIEWER_ASCII_VTK Outputs the object to a VTK file
- PETSC_VIEWER_NATIVE Store the object to the binary file in its native format
- (for example, dense matrices are stored as dense),
- DMDA vectors are dumped directly to the file instead of
- being first put in the natural ordering
- PETSC_VIEWER_DRAW_BASIC Views the vector with a simple 1d plot
- PETSC_VIEWER_DRAW_LG Views the vector with a line graph
- PETSC_VIEWER_DRAW_CONTOUR Views the vector with a contour plot
- */
- void viewer(const std::string &file_name,
- const PetscViewerFormat vw_format = PETSC_VIEWER_ASCII_MATLAB );
-
- private:
- void destroy() { if (_A) MatDestroy(&_A); _A = nullptr; }
-
- /// PETSc matrix
- Mat _A = nullptr;
-
- /// 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;
-
- /*!
- \brief Create the matrix, configure memory allocation and set the
- related member data.
- \param d_nz Number of nonzeros per row in the diagonal portion of
- local submatrix (same value is used for all local rows),
- \param o_nz Number of nonzeros per row in the off-diagonal portion of
- local submatrix (same value is used for all local rows)
- */
- void create(const PetscInt d_nz, const PetscInt o_nz);
-
- friend bool finalizeMatrixAssembly(PETScMatrix &mat, const MatAssemblyType asm_type);
+ add(indices.rows, cols, sub_matrix);
+ }
+
+ /*!
+ \brief Add a submatrix to this.
+ \param row_pos The row indices of the entries of the submatrix.
+ \param col_pos The column indices of the entries of the submatrix.
+ \param sub_mat A dense matrix to be added on.
+ */
+ 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);
+
+ /*! View the global vector for test purpose. Do not use it for output a big
+ vector.
+ \param file_name File name for output
+ \param vw_format File format listed as:
+ PETSC_VIEWER_DEFAULT Default format
+ PETSC_VIEWER_ASCII_MATLAB MATLAB format
+ PETSC_VIEWER_ASCII_DENSE Print matrix as dense
+ PETSC_VIEWER_ASCII_IMPL Implementation-specific format
+ (which is in many cases the same as
+ the default)
+ PETSC_VIEWER_ASCII_INFO Basic information about object
+ PETSC_VIEWER_ASCII_INFO_DETAIL More detailed info about object
+ PETSC_VIEWER_ASCII_COMMON Identical output format for all objects
+ of a particular type
+ PETSC_VIEWER_ASCII_INDEX (for vectors) Prints the vector element
+ number next to each vector entry
+ PETSC_VIEWER_ASCII_SYMMODU Print parallel vectors without
+ indicating the processor ranges
+ PETSC_VIEWER_ASCII_VTK Outputs the object to a VTK file
+ PETSC_VIEWER_NATIVE Store the object to the binary file in
+ its native format (for example,
+ dense matrices are stored as dense),
+ DMDA vectors are dumped directly to
+ the file instead of being first put in
+ the natural ordering
+ PETSC_VIEWER_DRAW_BASIC Views the vector with a simple 1d plot
+ PETSC_VIEWER_DRAW_LG Views the vector with a line graph
+ PETSC_VIEWER_DRAW_CONTOUR Views the vector with a contour plot
+ */
+ void viewer(const std::string& file_name,
+ const PetscViewerFormat vw_format = PETSC_VIEWER_ASCII_MATLAB);
+
+private:
+ void destroy()
+ {
+ if (_A)
+ MatDestroy(&_A);
+ _A = nullptr;
+ }
+
+ /// PETSc matrix
+ Mat _A = nullptr;
+
+ /// 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;
+
+ /*!
+ \brief Create the matrix, configure memory allocation and set the
+ related member data.
+ \param d_nz Number of nonzeros per row in the diagonal portion of
+ local submatrix (same value is used for all local rows),
+ \param o_nz Number of nonzeros per row in the off-diagonal portion of
+ local submatrix (same value is used for all local rows)
+ */
+ void create(const PetscInt d_nz, const PetscInt o_nz);
+
+ friend bool finalizeMatrixAssembly(PETScMatrix& mat,
+ const MatAssemblyType asm_type);
};
/*!
@@ -260,15 +238,16 @@ class PETScMatrix
\param col_pos The global indices of the colums of the dense sub-matrix.
\param sub_mat A dense sub-matrix to be added.
*/
-template<class T_DENSE_MATRIX>
+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 T_DENSE_MATRIX& sub_mat)
{
- const PetscInt nrows = static_cast<PetscInt> (row_pos.size());
- const PetscInt ncols = static_cast<PetscInt> (col_pos.size());
+ 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);
};
/*!
@@ -277,6 +256,7 @@ void PETScMatrix::add(std::vector<PetscInt> const& row_pos,
\param asm_type Assmebly type, either MAT_FLUSH_ASSEMBLY
or MAT_FINAL_ASSEMBLY.
*/
-bool finalizeMatrixAssembly(PETScMatrix &mat, const MatAssemblyType asm_type = MAT_FINAL_ASSEMBLY);
+bool finalizeMatrixAssembly(
+ PETScMatrix& mat, const MatAssemblyType asm_type = MAT_FINAL_ASSEMBLY);
-} // end namespace
+} // end namespace
diff --git a/MathLib/LinAlg/PETSc/PETScMatrixOption.h b/MathLib/LinAlg/PETSc/PETScMatrixOption.h
index 66a6750cc4c25b82f0bbdece319e003324810f76..5c08a7020d89e48eeb69560df33f0faba549cf7e 100644
--- a/MathLib/LinAlg/PETSc/PETScMatrixOption.h
+++ b/MathLib/LinAlg/PETSc/PETScMatrixOption.h
@@ -18,13 +18,18 @@
namespace MathLib
{
/*!
- \brief This a struct data containing the configuration information to create a PETSc type matrix
+ \brief This a struct data containing the configuration information to create
+ a PETSc type matrix
*/
struct PETScMatrixOption
{
- PETScMatrixOption() : is_global_size(true), n_local_cols(PETSC_DECIDE),
- d_nz(PETSC_DECIDE), o_nz(PETSC_DECIDE)
- { }
+ PETScMatrixOption()
+ : is_global_size(true),
+ n_local_cols(PETSC_DECIDE),
+ d_nz(PETSC_DECIDE),
+ o_nz(PETSC_DECIDE)
+ {
+ }
/*!
\brief Flag for the type of size, which is one of arguments of
@@ -39,16 +44,18 @@ struct PETScMatrixOption
PetscInt n_local_cols;
/*!
- \brief Number of nonzeros per row in the diagonal portion of local submatrix
+ \brief Number of nonzeros per row in the diagonal portion of local
+ submatrix
(same value is used for all local rows), the default is PETSC_DECIDE
*/
PetscInt d_nz;
/*!
- \brief Number of nonzeros per row in the off-diagonal portion of local submatrix
+ \brief Number of nonzeros per row in the off-diagonal portion of local
+ submatrix
(same value is used for all local rows), the default is PETSC_DECIDE
*/
PetscInt o_nz;
};
-} // end namespace
+} // end namespace
diff --git a/MathLib/LinAlg/PETSc/PETScTools.cpp b/MathLib/LinAlg/PETSc/PETScTools.cpp
index 68a0a2059a8957dd69d4c2d8adf7c7d59cc25ca5..4b177a35543b6a64a97b9f7b977d2ffc6ef1f6a0 100644
--- a/MathLib/LinAlg/PETSc/PETScTools.cpp
+++ b/MathLib/LinAlg/PETSc/PETScTools.cpp
@@ -18,10 +18,9 @@
namespace MathLib
{
-
-void applyKnownSolution(PETScMatrix &A, PETScVector &b, PETScVector &x,
- const std::vector<PetscInt> &vec_knownX_id,
- const std::vector<PetscScalar> &vec_knownX_x)
+void applyKnownSolution(PETScMatrix& A, PETScVector& b, PETScVector& x,
+ const std::vector<PetscInt>& vec_knownX_id,
+ const std::vector<PetscScalar>& vec_knownX_x)
{
A.finalizeAssembly();
@@ -30,7 +29,7 @@ void applyKnownSolution(PETScMatrix &A, PETScVector &b, PETScVector &x,
x.finalizeAssembly();
b.finalizeAssembly();
- if(vec_knownX_id.size() > 0)
+ if (vec_knownX_id.size() > 0)
{
x.set(vec_knownX_id, vec_knownX_x);
b.set(vec_knownX_id, vec_knownX_x);
@@ -40,6 +39,4 @@ void applyKnownSolution(PETScMatrix &A, PETScVector &b, PETScVector &x,
b.finalizeAssembly();
}
-} // end of namespace MathLib
-
-
+} // end of namespace MathLib
diff --git a/MathLib/LinAlg/PETSc/PETScTools.h b/MathLib/LinAlg/PETSc/PETScTools.h
index 9c0a20db04a95b090ac107053bee626bb401045e..0499064cf582dbd494069391fc623af720b10033 100644
--- a/MathLib/LinAlg/PETSc/PETScTools.h
+++ b/MathLib/LinAlg/PETSc/PETScTools.h
@@ -32,7 +32,7 @@ namespace MathLib
\param vec_knownX_id A vector of known solution entry IDs
\param vec_knownX_x A vector of known solutions
*/
-void applyKnownSolution(PETScMatrix &A, PETScVector &b, PETScVector &x,
- const std::vector<PetscInt> &vec_knownX_id,
- const std::vector<PetscScalar> &vec_knownX_x);
-} // end of namespace MathLib
+void applyKnownSolution(PETScMatrix& A, PETScVector& b, PETScVector& x,
+ const std::vector<PetscInt>& vec_knownX_id,
+ const std::vector<PetscScalar>& vec_knownX_x);
+} // end of namespace MathLib
diff --git a/MathLib/LinAlg/PETSc/PETScVector.cpp b/MathLib/LinAlg/PETSc/PETScVector.cpp
index 12e1a2b10e109965edfed5965cb456dddb26dcc3..071b90dcef17346acda8831ceec253d68da06bca 100644
--- a/MathLib/LinAlg/PETSc/PETScVector.cpp
+++ b/MathLib/LinAlg/PETSc/PETScVector.cpp
@@ -18,7 +18,6 @@
*
*/
-
#include "PETScVector.h"
#include <algorithm>
@@ -28,7 +27,7 @@ namespace MathLib
{
PETScVector::PETScVector(const PetscInt vec_size, const bool is_global_size)
{
- if( is_global_size )
+ if (is_global_size)
{
VecCreate(PETSC_COMM_WORLD, &_v);
VecSetSizes(_v, PETSC_DECIDE, vec_size);
@@ -46,11 +45,10 @@ PETScVector::PETScVector(const PetscInt vec_size, const bool is_global_size)
PETScVector::PETScVector(const PetscInt vec_size,
const std::vector<PetscInt>& ghost_ids,
const bool is_global_size)
- : _size_ghosts {static_cast<PetscInt>(ghost_ids.size())},
- _has_ghost_id {true}
+ : _size_ghosts{static_cast<PetscInt>(ghost_ids.size())}, _has_ghost_id{true}
{
- PetscInt nghosts = static_cast<PetscInt>( ghost_ids.size() );
- if ( is_global_size )
+ PetscInt nghosts = static_cast<PetscInt>(ghost_ids.size());
+ if (is_global_size)
{
VecCreateGhost(PETSC_COMM_WORLD, PETSC_DECIDE, vec_size, nghosts,
ghost_ids.data(), &_v);
@@ -65,33 +63,33 @@ PETScVector::PETScVector(const PetscInt vec_size,
config();
- for (PetscInt i=0; i<nghosts; i++)
- _global_ids2local_ids_ghost.emplace(ghost_ids[i],
- _size_loc + i);
+ for (PetscInt i = 0; i < nghosts; i++)
+ _global_ids2local_ids_ghost.emplace(ghost_ids[i], _size_loc + i);
}
-PETScVector::PETScVector(const PETScVector &existing_vec, const bool deep_copy)
+PETScVector::PETScVector(const PETScVector& existing_vec, const bool deep_copy)
{
shallowCopy(existing_vec);
// Copy values
- if(deep_copy)
+ if (deep_copy)
{
VecCopy(existing_vec._v, _v);
}
}
-PETScVector::PETScVector(PETScVector &&other)
- : _v {std::move(other._v)},
- _v_loc {std::move(other._v_loc)},
- _start_rank {other._start_rank},
- _end_rank {other._end_rank},
- _size {other._size},
- _size_loc {other._size_loc},
- _size_ghosts {other._size_ghosts},
- _has_ghost_id {other._has_ghost_id},
- _global_ids2local_ids_ghost {other._global_ids2local_ids_ghost}
-{}
+PETScVector::PETScVector(PETScVector&& other)
+ : _v{std::move(other._v)},
+ _v_loc{std::move(other._v_loc)},
+ _start_rank{other._start_rank},
+ _end_rank{other._end_rank},
+ _size{other._size},
+ _size_loc{other._size_loc},
+ _size_ghosts{other._size_ghosts},
+ _has_ghost_id{other._has_ghost_id},
+ _global_ids2local_ids_ghost{other._global_ids2local_ids_ghost}
+{
+}
void PETScVector::config()
{
@@ -111,36 +109,35 @@ void PETScVector::finalizeAssembly()
VecAssemblyEnd(_v);
}
-void PETScVector::gatherLocalVectors( PetscScalar local_array[],
- PetscScalar global_array[]) const
+void PETScVector::gatherLocalVectors(PetscScalar local_array[],
+ PetscScalar global_array[]) const
{
// Collect vectors from processors.
int size_rank;
MPI_Comm_size(PETSC_COMM_WORLD, &size_rank);
// number of elements to be sent for each rank
- std::vector<PetscInt> i_cnt(size_rank);
+ std::vector<PetscInt> i_cnt(size_rank);
- MPI_Allgather(&_size_loc, 1, MPI_INT, &i_cnt[0], 1, MPI_INT, PETSC_COMM_WORLD);
+ MPI_Allgather(&_size_loc, 1, MPI_INT, &i_cnt[0], 1, MPI_INT,
+ PETSC_COMM_WORLD);
// collect local array
PetscInt offset = 0;
// offset in the receive vector of the data from each rank
- std::vector<PetscInt> i_disp(size_rank);
- for(PetscInt i=0; i<size_rank; i++)
+ std::vector<PetscInt> i_disp(size_rank);
+ for (PetscInt i = 0; i < size_rank; i++)
{
i_disp[i] = offset;
offset += i_cnt[i];
}
- MPI_Allgatherv(local_array, _size_loc, MPI_DOUBLE,
- global_array, &i_cnt[0], &i_disp[0], MPI_DOUBLE, PETSC_COMM_WORLD);
-
+ MPI_Allgatherv(local_array, _size_loc, MPI_DOUBLE, global_array, &i_cnt[0],
+ &i_disp[0], MPI_DOUBLE, PETSC_COMM_WORLD);
}
void PETScVector::getGlobalVector(std::vector<PetscScalar>& u) const
{
-
#ifdef TEST_MEM_PETSC
PetscLogDouble mem1, mem2;
PetscMemoryGetCurrentUsage(&mem1);
@@ -148,21 +145,24 @@ void PETScVector::getGlobalVector(std::vector<PetscScalar>& u) const
assert(u.size() == _size);
- PetscScalar *xp = nullptr;
+ PetscScalar* xp = nullptr;
VecGetArray(_v, &xp);
gatherLocalVectors(xp, u.data());
- //This following line may be needed late on
+ // This following line may be needed late on
// for a communication load balance:
- //MPI_Barrier(PETSC_COMM_WORLD);
+ // MPI_Barrier(PETSC_COMM_WORLD);
VecRestoreArray(_v, &xp);
- //TEST
+// TEST
#ifdef TEST_MEM_PETSC
PetscMemoryGetCurrentUsage(&mem2);
- PetscPrintf(PETSC_COMM_WORLD, "### Memory usage by Updating. Before :%f After:%f Increase:%d\n", mem1, mem2, (int)(mem2 - mem1));
+ PetscPrintf(
+ PETSC_COMM_WORLD,
+ "### Memory usage by Updating. Before :%f After:%f Increase:%d\n", mem1,
+ mem2, (int)(mem2 - mem1));
#endif
}
@@ -170,17 +170,16 @@ void PETScVector::setLocalAccessibleVector() const
{
if (_entry_array.empty())
{
- const PetscInt array_size
- = _global_ids2local_ids_ghost.size() > 0 ?
- _size_loc + _size_ghosts: _size;
+ const PetscInt array_size = _global_ids2local_ids_ghost.size() > 0
+ ? _size_loc + _size_ghosts
+ : _size;
_entry_array.resize(array_size);
}
- if ( !_global_ids2local_ids_ghost.empty() )
+ if (!_global_ids2local_ids_ghost.empty())
{
PetscScalar* loc_x = getLocalVector();
- std::copy_n(loc_x, _size_loc + _size_ghosts,
- _entry_array.begin());
+ std::copy_n(loc_x, _size_loc + _size_ghosts, _entry_array.begin());
restoreArray(loc_x);
}
else
@@ -189,7 +188,7 @@ void PETScVector::setLocalAccessibleVector() const
void PETScVector::copyValues(std::vector<PetscScalar>& u) const
{
- assert(u.size() == (std::size_t) (getLocalSize() + getGhostSize()));
+ assert(u.size() == (std::size_t)(getLocalSize() + getGhostSize()));
PetscScalar* loc_x = getLocalVector();
std::copy_n(loc_x, getLocalSize() + getGhostSize(), u.begin());
@@ -204,12 +203,12 @@ PetscScalar PETScVector::get(const PetscInt idx) const
}
// Ghost entries, and its original index is 0.
- const PetscInt id_p = (idx == -_size) ? 0 : std::abs(idx);
+ const PetscInt id_p = (idx == -_size) ? 0 : std::abs(idx);
return _entry_array[id_p];
}
-
-std::vector<PetscScalar> PETScVector::get(std::vector<IndexType> const& indices) const
+std::vector<PetscScalar> PETScVector::get(
+ std::vector<IndexType> const& indices) const
{
std::vector<PetscScalar> local_x(indices.size());
// If VecGetValues can get values from different processors,
@@ -218,18 +217,18 @@ std::vector<PetscScalar> PETScVector::get(std::vector<IndexType> const& indices)
if (_global_ids2local_ids_ghost.size() > 0)
{
- for (std::size_t i=0; i<indices.size(); i++)
+ for (std::size_t i = 0; i < indices.size(); i++)
{
local_x[i] = _entry_array[getLocalIndex(indices[i])];
}
}
else
{
- for (std::size_t i=0; i<indices.size(); i++)
+ for (std::size_t i = 0; i < indices.size(); i++)
{
// Ghost entries, and its original index is 0.
- const IndexType id_p = (indices[i] == -_size)
- ? 0 : std::abs(indices[i]);
+ const IndexType id_p =
+ (indices[i] == -_size) ? 0 : std::abs(indices[i]);
local_x[i] = _entry_array[id_p];
}
}
@@ -239,7 +238,7 @@ std::vector<PetscScalar> PETScVector::get(std::vector<IndexType> const& indices)
PetscScalar* PETScVector::getLocalVector() const
{
- PetscScalar *loc_array;
+ PetscScalar* loc_array;
if (_has_ghost_id)
{
VecGhostUpdateBegin(_v, INSERT_VALUES, SCATTER_FORWARD);
@@ -254,13 +253,12 @@ PetscScalar* PETScVector::getLocalVector() const
PetscInt PETScVector::getLocalIndex(const PetscInt global_index) const
{
- if (global_index >= 0) // non-ghost entry.
+ if (global_index >= 0) // non-ghost entry.
return global_index - _start_rank;
// A special case for a ghost location with global index equal to
// the size of the local vector:
- PetscInt real_global_index =
- (-global_index == _size) ? 0 : -global_index;
+ PetscInt real_global_index = (-global_index == _size) ? 0 : -global_index;
return _global_ids2local_ids_ghost.at(real_global_index);
}
@@ -276,7 +274,8 @@ void PETScVector::restoreArray(PetscScalar* array) const
VecRestoreArray(_v, &array);
}
-void PETScVector::viewer(const std::string &file_name, const PetscViewerFormat vw_format) const
+void PETScVector::viewer(const std::string& file_name,
+ const PetscViewerFormat vw_format) const
{
PetscViewer viewer;
PetscViewerASCIIOpen(PETSC_COMM_WORLD, file_name.c_str(), &viewer);
@@ -285,35 +284,34 @@ void PETScVector::viewer(const std::string &file_name, const PetscViewerFormat v
PetscObjectSetName((PetscObject)_v, file_name.c_str());
VecView(_v, viewer);
-#define nEXIT_TEST
+#define nEXIT_TEST
#ifdef EXIT_TEST
VecDestroy(_v);
PetscFinalize();
exit(0);
#endif
-
}
-void PETScVector::shallowCopy(const PETScVector &v)
+void PETScVector::shallowCopy(const PETScVector& v)
{
destroy();
VecDuplicate(v.getRawVector(), &_v);
- _start_rank = v._start_rank;
- _end_rank = v._end_rank;
- _size = v._size;
- _size_loc = v._size_loc;
- _size_ghosts = v._size_ghosts;
+ _start_rank = v._start_rank;
+ _end_rank = v._end_rank;
+ _size = v._size;
+ _size_loc = v._size_loc;
+ _size_ghosts = v._size_ghosts;
_has_ghost_id = v._has_ghost_id;
_global_ids2local_ids_ghost = v._global_ids2local_ids_ghost;
VecSetOption(_v, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
}
-void finalizeVectorAssembly(PETScVector &vec)
+void finalizeVectorAssembly(PETScVector& vec)
{
vec.finalizeAssembly();
}
-} //end of namespace
+} // end of namespace
diff --git a/MathLib/LinAlg/PETSc/PETScVector.h b/MathLib/LinAlg/PETSc/PETScVector.h
index bb5193109d521e3a668c5c0339ca03f7e39ab9c0..5c91637fa1c6a8ecaecade9e97069b17acc70b2a 100644
--- a/MathLib/LinAlg/PETSc/PETScVector.h
+++ b/MathLib/LinAlg/PETSc/PETScVector.h
@@ -31,276 +31,261 @@ namespace MathLib
*/
class PETScVector
{
- public:
- using IndexType = PetscInt;
- // TODO make this class opaque, s.t. the PETSc symbols are not scattered all
- // over the global namespace
- using PETSc_Vec = Vec;
-
- public:
- PETScVector() {}
-
- /*!
- \brief Constructor
- \param vec_size The size of the vector, either global or local
- \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.
- If is_global_size is true, the vector
- is created by the global size, the local size
- of the vector is determined by PETSc,
- and vice versa is the same.
- */
- PETScVector(const PetscInt vec_size, const bool is_global_size = true);
-
- /*!
- \brief Constructor
- \param vec_size The size of the vector, either global or local
- \param ghost_ids The global indices of ghost entries
- \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.
- */
- PETScVector(const PetscInt vec_size, const std::vector<PetscInt>& ghost_ids,
- const bool is_global_size = true);
-
- /*!
- \brief Copy constructor
- \param existing_vec The vector to be copied
- \param deep_copy The flag for a deep copy, which means to copy the values as well,
- the default is true
- */
- explicit PETScVector(const PETScVector &existing_vec, const bool deep_copy = true);
-
- PETScVector(PETScVector&& other);
-
- ~PETScVector()
- {
- destroy();
- }
-
- /// Perform MPI collection of assembled entries in buffer
- void finalizeAssembly();
-
- /// Get the global size of the vector
- PetscInt size() const
- {
- return _size;
- }
-
- /// Get the number of entries in the same rank
- PetscInt getLocalSize() const
- {
- return _size_loc;
- }
-
- /// Get the number of ghost entries in the same rank
- PetscInt getGhostSize() const
- {
- return _size_ghosts;
- }
-
- /// Get the start index of the local vector
- PetscInt getRangeBegin() const
- {
- return _start_rank;
- }
-
- /// Get the end index of the local vector
- PetscInt getRangeEnd() const
- {
- return _end_rank;
- }
-
- /*!
- Insert a single entry with value.
- \param i Entry index
- \param value Entry value
-
- */
- void set(const PetscInt i, const PetscScalar value)
- {
- VecSetValue(_v, i, value, INSERT_VALUES);
- }
-
- /*!
- Add a value to an entry.
- \param i Number of the entry
- \param value Value.
- */
- void add(const PetscInt i, const PetscScalar value)
- {
- VecSetValue(_v, i, value, ADD_VALUES);
- }
-
- /*!
- Add values to several entries
- \param e_idxs Indicies of entries to be added
- Note: std::size_t cannot be the type of e_idxs template argument
- \param sub_vec Entries to be added
- */
- template<class T_SUBVEC> void add(const std::vector<PetscInt> &e_idxs,
- const T_SUBVEC &sub_vec)
- {
- VecSetValues(_v, e_idxs.size(), &e_idxs[0], &sub_vec[0], ADD_VALUES);
- }
-
- /*!
- Add values to several entries
- \param e_idxs Indicies of entries to be added.
- Note: std::size_t cannot be the type of e_idxs template argument
- \param sub_vec Entries to be added
- */
- template<class T_SUBVEC> void set(const std::vector<PetscInt> &e_idxs,
- const T_SUBVEC &sub_vec)
- {
- VecSetValues(_v, e_idxs.size(), &e_idxs[0], &sub_vec[0], INSERT_VALUES);
- }
-
- // TODO preliminary
- void setZero() { VecSet(_v, 0.0); }
-
- /// Disallow moving.
- /// \todo This operator should be implemented properly when doing a
- /// general cleanup of all matrix and vector classes.
- PETScVector& operator = (PETScVector &&) = delete;
-
- /// Set local accessible vector in order to get entries.
- /// Call this before call operator[] or get(...).
- void setLocalAccessibleVector() const;
-
- /// Get several entries. setLocalAccessibleVector() must be
- /// called beforehand.
- std::vector<PetscScalar> get(std::vector<IndexType> const& indices) const;
-
- /// Get the value of an entry by [] operator.
- /// setLocalAccessibleVector() must be called beforehand.
- PetscScalar operator[] (PetscInt idx) const
- {
- return get(idx);
- }
-
- /*!
- Get global vector
- \param u Array to store the global vector. Memory allocation is needed in advance
- */
- void getGlobalVector(std::vector<PetscScalar>& u) const;
-
- /* Get an entry value. This is an expensive operation,
- and it only get local value. Use it for only test purpose
- Get the value of an entry by [] operator.
- setLocalAccessibleVector() must be called beforehand.
- */
- PetscScalar get(const PetscInt idx) const;
-
- //! Exposes the underlying PETSc vector.
- PETSc_Vec getRawVector() { return _v; }
-
- /*! Exposes the underlying PETSc vector.
- *
- * \warning
- * This method is dangerous insofar as you can do arbitrary things also
- * with a const PETSc vector!
- */
- PETSc_Vec getRawVector() const { return _v; }
-
- /*!
- Copy local entries including ghost ones to an array
- \param u Preallocated vector for the values of local entries.
- */
- void copyValues(std::vector<PetscScalar>& u) const;
-
- /*! View the global vector for test purpose. Do not use it for output a big vector.
- \param file_name File name for output
- \param vw_format File format listed as:
- PETSC_VIEWER_DEFAULT Default format
- PETSC_VIEWER_ASCII_MATLAB MATLAB format
- PETSC_VIEWER_ASCII_DENSE Print matrix as dense
- PETSC_VIEWER_ASCII_IMPL Implementation-specific format
- (which is in many cases the same as the default)
- PETSC_VIEWER_ASCII_INFO Basic information about object
- PETSC_VIEWER_ASCII_INFO_DETAIL More detailed info about object
- PETSC_VIEWER_ASCII_COMMON Identical output format for all objects of a particular type
- PETSC_VIEWER_ASCII_INDEX (for vectors) Prints the vector element number next to
- each vector entry
- PETSC_VIEWER_ASCII_SYMMODU Print parallel vectors without indicating the processor ranges
- PETSC_VIEWER_ASCII_VTK Outputs the object to a VTK file
- PETSC_VIEWER_NATIVE Store the object to the binary file in its native format
- (for example, dense matrices are stored as dense),
- DMDA vectors are dumped directly to the file instead of
- being first put in the natural ordering
- PETSC_VIEWER_DRAW_BASIC Views the vector with a simple 1d plot
- PETSC_VIEWER_DRAW_LG Views the vector with a line graph
- PETSC_VIEWER_DRAW_CONTOUR Views the vector with a contour plot
- */
- void viewer(const std::string &file_name,
- const PetscViewerFormat vw_format = PETSC_VIEWER_ASCII_MATLAB ) const;
-
- void shallowCopy(const PETScVector &v);
-
- private:
- void destroy() { if (_v) VecDestroy(&_v); _v = nullptr; }
-
- PETSc_Vec _v = nullptr;
- /// Local vector, which is only for the case that _v is created
- /// with ghost entries.
- mutable PETSc_Vec _v_loc = nullptr;
-
- /// Starting index in a rank
- PetscInt _start_rank;
- /// Ending index in a rank
- PetscInt _end_rank;
-
- /// Size of the vector
- PetscInt _size;
- /// Size of local entries
- PetscInt _size_loc;
- /// Size of local ghost entries
- PetscInt _size_ghosts = 0;
-
- /// Flag to indicate whether the vector is created with ghost entry indices
- bool _has_ghost_id = false;
-
- /*!
- \brief Array containing the entries of the vector.
- If the vector is created without given ghost IDs, the array
- contains all entries of the global vector, _v. Otherwise it
- only contains the entries of the global vector owned by the
- current rank.
- */
- mutable std::vector<PetscScalar> _entry_array;
-
- /// Map global indices of ghost enrties to local indices
- mutable std::map<PetscInt, PetscInt> _global_ids2local_ids_ghost;
-
- /*!
- \brief Collect local vectors
- \param local_array Local array
- \param global_array Global array
- */
- void gatherLocalVectors(PetscScalar local_array[],
- PetscScalar global_array[]) const;
-
- /*!
- Get local vector, i.e. entries in the same rank
- */
- PetscScalar* getLocalVector() const;
-
- /// Get local index by a global index
- PetscInt getLocalIndex(const PetscInt global_index) const;
-
- /*!
- Restore array after finish access local array
- \param array Pointer to the local array fetched by VecGetArray
- */
- inline void restoreArray(PetscScalar* array) const;
-
- /// A funtion called by constructors to configure members
- void config();
-
- friend void finalizeVectorAssembly(PETScVector &vec);
+public:
+ using IndexType = PetscInt;
+ // TODO make this class opaque, s.t. the PETSc symbols are not scattered all
+ // over the global namespace
+ using PETSc_Vec = Vec;
+
+public:
+ PETScVector() {}
+ /*!
+ \brief Constructor
+ \param vec_size The size of the vector, either global or local
+ \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.
+ If is_global_size is true, the vector
+ is created by the global size, the local size
+ of the vector is determined by PETSc,
+ and vice versa is the same.
+ */
+ PETScVector(const PetscInt vec_size, const bool is_global_size = true);
+
+ /*!
+ \brief Constructor
+ \param vec_size The size of the vector, either global or local
+ \param ghost_ids The global indices of ghost entries
+ \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.
+ */
+ PETScVector(const PetscInt vec_size, const std::vector<PetscInt>& ghost_ids,
+ const bool is_global_size = true);
+
+ /*!
+ \brief Copy constructor
+ \param existing_vec The vector to be copied
+ \param deep_copy The flag for a deep copy, which means to copy the
+ values as well, the default is true
+ */
+ explicit PETScVector(const PETScVector& existing_vec,
+ const bool deep_copy = true);
+
+ PETScVector(PETScVector&& other);
+
+ ~PETScVector() { destroy(); }
+ /// Perform MPI collection of assembled entries in buffer
+ void finalizeAssembly();
+
+ /// Get the global size of the vector
+ PetscInt size() const { return _size; }
+ /// Get the number of entries in the same rank
+ PetscInt getLocalSize() const { return _size_loc; }
+ /// Get the number of ghost entries in the same rank
+ PetscInt getGhostSize() const { return _size_ghosts; }
+ /// Get the start index of the local vector
+ PetscInt getRangeBegin() const { return _start_rank; }
+ /// Get the end index of the local vector
+ PetscInt getRangeEnd() const { return _end_rank; }
+ /*!
+ Insert a single entry with value.
+ \param i Entry index
+ \param value Entry value
+
+ */
+ void set(const PetscInt i, const PetscScalar value)
+ {
+ VecSetValue(_v, i, value, INSERT_VALUES);
+ }
+
+ /*!
+ Add a value to an entry.
+ \param i Number of the entry
+ \param value Value.
+ */
+ void add(const PetscInt i, const PetscScalar value)
+ {
+ VecSetValue(_v, i, value, ADD_VALUES);
+ }
+
+ /*!
+ Add values to several entries.
+ \param e_idxs Indicies of entries to be added
+ Note: std::size_t cannot be the type of e_idxs template
+ argument.
+ \param sub_vec Entries to be added.
+ */
+ template <class T_SUBVEC>
+ void add(const std::vector<PetscInt>& e_idxs, const T_SUBVEC& sub_vec)
+ {
+ VecSetValues(_v, e_idxs.size(), &e_idxs[0], &sub_vec[0], ADD_VALUES);
+ }
+
+ /*!
+ Add values to several entries
+ \param e_idxs Indicies of entries to be added.
+ Note: std::size_t cannot be the type of e_idxs template
+ argument
+ \param sub_vec Entries to be added
+ */
+ template <class T_SUBVEC>
+ void set(const std::vector<PetscInt>& e_idxs, const T_SUBVEC& sub_vec)
+ {
+ VecSetValues(_v, e_idxs.size(), &e_idxs[0], &sub_vec[0], INSERT_VALUES);
+ }
+
+ // TODO preliminary
+ void setZero() { VecSet(_v, 0.0); }
+ /// Disallow moving.
+ /// \todo This operator should be implemented properly when doing a
+ /// general cleanup of all matrix and vector classes.
+ PETScVector& operator=(PETScVector&&) = delete;
+
+ /// Set local accessible vector in order to get entries.
+ /// Call this before call operator[] or get(...).
+ void setLocalAccessibleVector() const;
+
+ /// Get several entries. setLocalAccessibleVector() must be
+ /// called beforehand.
+ std::vector<PetscScalar> get(std::vector<IndexType> const& indices) const;
+
+ /// Get the value of an entry by [] operator.
+ /// setLocalAccessibleVector() must be called beforehand.
+ PetscScalar operator[](PetscInt idx) const { return get(idx); }
+ /*!
+ Get global vector
+ \param u Array to store the global vector. Memory allocation is needed in
+ advance
+ */
+ void getGlobalVector(std::vector<PetscScalar>& u) const;
+
+ /* Get an entry value. This is an expensive operation,
+ and it only get local value. Use it for only test purpose
+ Get the value of an entry by [] operator.
+ setLocalAccessibleVector() must be called beforehand.
+ */
+ PetscScalar get(const PetscInt idx) const;
+
+ //! Exposes the underlying PETSc vector.
+ PETSc_Vec getRawVector() { return _v; }
+ /*! Exposes the underlying PETSc vector.
+ *
+ * \warning
+ * This method is dangerous insofar as you can do arbitrary things also
+ * with a const PETSc vector!
+ */
+ PETSc_Vec getRawVector() const { return _v; }
+ /*!
+ Copy local entries including ghost ones to an array
+ \param u Preallocated vector for the values of local entries.
+ */
+ void copyValues(std::vector<PetscScalar>& u) const;
+
+ /*! View the global vector for test purpose. Do not use it for output a big
+ vector.
+ \param file_name File name for output
+ \param vw_format File format listed as:
+ PETSC_VIEWER_DEFAULT Default format
+ PETSC_VIEWER_ASCII_MATLAB MATLAB format
+ PETSC_VIEWER_ASCII_DENSE Print matrix as dense
+ PETSC_VIEWER_ASCII_IMPL Implementation-specific format
+ (which is in many cases the same as
+ the default)
+ PETSC_VIEWER_ASCII_INFO Basic information about object
+ PETSC_VIEWER_ASCII_INFO_DETAIL More detailed info about object
+ PETSC_VIEWER_ASCII_COMMON Identical output format for all objects
+ of a particular type
+ PETSC_VIEWER_ASCII_INDEX (for vectors) Prints the vector element
+ number next to each vector entry
+ PETSC_VIEWER_ASCII_SYMMODU Print parallel vectors without
+ indicating the processor ranges
+ PETSC_VIEWER_ASCII_VTK Outputs the object to a VTK file
+ PETSC_VIEWER_NATIVE Store the object to the binary file in
+ its native format (for example,
+ dense matrices are stored as dense),
+ DMDA vectors are dumped directly to
+ the file instead of being first put in
+ the natural ordering
+ PETSC_VIEWER_DRAW_BASIC Views the vector with a simple 1d plot
+ PETSC_VIEWER_DRAW_LG Views the vector with a line graph
+ PETSC_VIEWER_DRAW_CONTOUR Views the vector with a contour plot
+ */
+ void viewer(
+ const std::string& file_name,
+ const PetscViewerFormat vw_format = PETSC_VIEWER_ASCII_MATLAB) const;
+
+ void shallowCopy(const PETScVector& v);
+
+private:
+ void destroy()
+ {
+ if (_v)
+ VecDestroy(&_v);
+ _v = nullptr;
+ }
+
+ PETSc_Vec _v = nullptr;
+ /// Local vector, which is only for the case that _v is created
+ /// with ghost entries.
+ mutable PETSc_Vec _v_loc = nullptr;
+
+ /// Starting index in a rank
+ PetscInt _start_rank;
+ /// Ending index in a rank
+ PetscInt _end_rank;
+
+ /// Size of the vector
+ PetscInt _size;
+ /// Size of local entries
+ PetscInt _size_loc;
+ /// Size of local ghost entries
+ PetscInt _size_ghosts = 0;
+
+ /// Flag to indicate whether the vector is created with ghost entry indices
+ bool _has_ghost_id = false;
+
+ /*!
+ \brief Array containing the entries of the vector.
+ If the vector is created without given ghost IDs, the array
+ contains all entries of the global vector, _v. Otherwise it
+ only contains the entries of the global vector owned by the
+ current rank.
+ */
+ mutable std::vector<PetscScalar> _entry_array;
+
+ /// Map global indices of ghost enrties to local indices
+ mutable std::map<PetscInt, PetscInt> _global_ids2local_ids_ghost;
+
+ /*!
+ \brief Collect local vectors
+ \param local_array Local array
+ \param global_array Global array
+ */
+ void gatherLocalVectors(PetscScalar local_array[],
+ PetscScalar global_array[]) const;
+
+ /*!
+ Get local vector, i.e. entries in the same rank
+ */
+ PetscScalar* getLocalVector() const;
+
+ /// Get local index by a global index
+ PetscInt getLocalIndex(const PetscInt global_index) const;
+
+ /*!
+ Restore array after finish access local array
+ \param array Pointer to the local array fetched by VecGetArray
+ */
+ inline void restoreArray(PetscScalar* array) const;
+
+ /// A funtion called by constructors to configure members
+ void config();
+
+ friend void finalizeVectorAssembly(PETScVector& vec);
};
/// Function to finalize the vector assembly
-void finalizeVectorAssembly(PETScVector &vec);
+void finalizeVectorAssembly(PETScVector& vec);
-} // end namespace
+} // end namespace