Commit 582a81cd authored by Gautier Bureau's avatar Gautier Bureau Committed by GILLES Sebastien
Browse files

#410 Added some new operations in Matrix and the possibility to create a DenseMatrix.

parent 0ddff7ea
......@@ -208,6 +208,39 @@ namespace HappyHeart
do_petsc_destroy_ = true;
}
void Matrix::InitSequentialDenseMatrix(unsigned int Nrow, unsigned int Ncolumn,
const Mpi& mpi, const char* invoking_file, int invoking_line)
{
const auto& communicator = mpi.GetCommunicator();
assert(petsc_matrix_ == PETSC_NULL && "Should not be initialized when this method is called!");
// Follow Petsc's advice and build the matrix step by step (see MatCreateSeqAIJ for this advice).
int error_code = MatCreate(communicator, &petsc_matrix_);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatCreate", invoking_file, invoking_line);
const Mat& internal = Internal(); // const is indeed ignored below as Mat is truly a pointer...
error_code = MatSetType(internal, MATSEQDENSE);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatSetType", invoking_file, invoking_line);
PetscInt Nrow_int = static_cast<PetscInt>(Nrow);
PetscInt Ncol_int = static_cast<PetscInt>(Ncolumn);
error_code = MatSetSizes(internal, Nrow_int, Ncol_int, Nrow_int, Ncol_int);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatSetSizes", invoking_file, invoking_line);
error_code = MatSetUp(internal);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatSetUp", invoking_file, invoking_line);
do_petsc_destroy_ = true;
}
void Matrix::ZeroEntries(const char* invoking_file, int invoking_line)
......@@ -513,7 +546,16 @@ namespace HappyHeart
for (auto i = 0ul; i < size; ++i)
row_content.emplace_back(std::make_pair(row_content_position_list[i], row_content_value_list[i]));
}
void Matrix::GetColumnVector(PetscInt column_index,
Vector& column,
const char* invoking_file, int invoking_line) const
{
int error_code = MatGetColumnVector(Internal(), column.Internal(), column_index);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatGetColumnVector", invoking_file, invoking_line);
}
void Matrix::SetFromPetscMat(Mat petsc_matrix)
......
......@@ -107,6 +107,11 @@ namespace HappyHeart
//! Alias to parent.
using parent = Internal::Wrappers::Petsc::BaseMatrix;
public:
//! Alias to unique pointer.
using unique_ptr = std::unique_ptr<Matrix>;
public:
//! Constructor.
......@@ -171,6 +176,23 @@ namespace HappyHeart
const MatrixPattern& matrix_pattern,
const Mpi& mpi, const char* invoking_file, int invoking_line);
/*!
* \brief Create a sequential sparse matrix.
*
*
* \param[in] Nrow Number of rows.
* \param[in] Ncolumn Number of columns.
* \copydetails doxygen_hide_mpi_param
* \param[in] invoking_file File that invoked the function or class; usually __FILE__.
* \param[in] invoking_line File that invoked the function or class; usually __LINE__.
*
* The matrix hence created will keep the non-zero structure: even if a value becomes 0 it will
* go on being held as a non-zero emplacement.
*/
void InitSequentialDenseMatrix(unsigned int Nrow, unsigned int Ncolumn,
const Mpi& mpi, const char* invoking_file, int invoking_line);
/*!
* \brief Handle over the internal Mat object.
......@@ -436,6 +458,21 @@ namespace HappyHeart
const char* invoking_file, int invoking_line) const;
/*!
* \brief Get the column of the matrix.
*
* This is a wrapper over MatGetColumnVector; it should be noticed that Petsc people hope this function is
* actually not used...
*
* \param[in] column_index Program-wise index of the column.
* \param[in] invoking_file File that invoked the function or class; usually __FILE__.
* \param[in] invoking_line File that invoked the function or class; usually __LINE__.
* \param[out] column Vector containing the column.
*/
void GetColumnVector(PetscInt column_index,
Vector& column,
const char* invoking_file, int invoking_line) const;
/*!
* \brief Wrapper over MatView.
......
......@@ -122,6 +122,17 @@ namespace HappyHeart
MatShift(const PetscScalar a, MatrixT& matrix, const char* invoking_file, int invoking_line);
/*!
* \brief Wrapper over MatShift, that performs M = M + a I, where a is a PetscScalar and I is the identity matrix.
*
* \param[in] invoking_file File that invoked the function or class; usually __FILE__.
* \param[in] invoking_line File that invoked the function or class; usually __LINE__.
*/
template<class MatrixT>
std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
MatShift(const PetscScalar a, MatrixT& matrix, const char* invoking_file, int invoking_line);
/*!
* \brief Wrapper over MatMult, that performs v2 = matrix * v1.
*
......@@ -139,6 +150,35 @@ namespace HappyHeart
const char* invoking_file, int invoking_line,
update_ghost do_update_ghost = update_ghost::yes);
/*!
* \brief Wrapper over MatTranspose, that performs in-place or out-of-place transpose of a matrix.
*
* \copydetails doxygen_hide_matranpose_warning
*
* If the operation is performed many times with each time the same non zero pattern for the matrices,
* rather use MatMatMultSymbolic/MatMatMultNumeric to improve efficiency.
* One can also reuse a pattern for several matrices, but so far I do not need that in HappyHeart so
* I equally use the default setting.
*
* \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
*
*/
template
<
class MatrixT,
class MatrixU
>
std::enable_if_t
<
std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
&& std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value,
void
>
MatTranspose(MatrixT& matrix1,
MatrixU& matrix2,
const char* invoking_file, int invoking_line);
/*!
* \brief Wrapper over MatMultAdd, that performs v3 = v2 + matrix * v1.
......@@ -424,7 +464,73 @@ namespace HappyHeart
DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
} // namespace Petsc
/*!
* \brief Wrapper over MatGetOrdering, gets a reordering for a matrix to reduce fill or to improve numerical stability of LU factorization.
*
* \see http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/MatOrderings/MatGetOrdering.html#MatGetOrdering
* for more details.
*
* \copydoc doxygen_hide_mat_get_ordering
*
*/
template<class MatrixT>
std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
GetOrdering(MatrixT& A,
MatOrderingType type,
IS *rperm,
IS *cperm,
const char* invoking_file, int invoking_line);
/*!
* \brief Wrapper over MatLUFactor, that performs in-place LU factorization of matrix. .
*
* \see http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatLUFactor.html#MatLUFactor
* for more details.
*
* \copydoc doxygen_hide_mat_lu_factor
*
*/
template<class MatrixT>
std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
LUFactor(MatrixT& A,
IS row, IS col,
const MatFactorInfo *info,
const char* invoking_file, int invoking_line);
/*!
* \brief Wrapper over MatMatSolve, solves A X = B, given a factored matrix.
*
* \see http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatMatSolve.html
* for more details.
*
* \copydoc doxygen_hide_mat_mat_solve
*
*/
template
<
class MatrixT,
class MatrixU,
class MatrixV
>
std::enable_if_t
<
std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
&& std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value
&& std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixV>::value,
void
>
MatMatSolve(const MatrixT& A,
const MatrixU& B,
MatrixV& X,
const char* invoking_file, int invoking_line);
} //namespace Petsc
} // namespace Wrappers
......
......@@ -191,11 +191,10 @@ namespace HappyHeart
template<class MatrixT>
std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
MatShift(const PetscScalar a,
MatrixT& matrix,
const char* invoking_file, int invoking_line)
const char* invoking_file, int invoking_line)
{
int error_code = ::MatShift(matrix.Internal(), a);
if (error_code)
......@@ -219,6 +218,34 @@ namespace HappyHeart
}
template
<
class MatrixT,
class MatrixU
>
std::enable_if_t
<
std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
&& std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value,
void
>
MatTranspose(MatrixT& matrix1,
MatrixU& matrix2,
const char* invoking_file, int invoking_line)
{
Mat result;
int error_code = ::MatTranspose(matrix1.Internal(), MAT_INITIAL_MATRIX, &result);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatTranspose", invoking_file, invoking_line);
error_code = ::MatCopy(result, matrix2.Internal(), SAME_NONZERO_PATTERN);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatCopy", invoking_file, invoking_line);
}
template<class MatrixT>
std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
MatMultAdd(const MatrixT& matrix,
......@@ -410,8 +437,70 @@ namespace HappyHeart
throw ExceptionNS::Exception(error_code, "MatPtAP", invoking_file, invoking_line);
}
} // namespace Petsc
template<class MatrixT>
std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
GetOrdering(MatrixT& A,
MatOrderingType type,
IS *rperm,
IS *cperm,
const char* invoking_file, int invoking_line)
{
int error_code = ::MatGetOrdering(A.Internal(),
type,
rperm,
cperm);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatGetOrdering", invoking_file, invoking_line);
}
template<class MatrixT>
std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
LUFactor(MatrixT& A,
IS row, IS col,
const MatFactorInfo *info,
const char* invoking_file, int invoking_line)
{
int error_code = ::MatLUFactor(A.Internal(),
row,
col,
info);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatLUFactor", invoking_file, invoking_line);
}
template
<
class MatrixT,
class MatrixU,
class MatrixV
>
std::enable_if_t
<
std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
&& std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value
&& std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixV>::value,
void
>
MatMatSolve(const MatrixT& A,
const MatrixU& B,
MatrixV& X,
const char* invoking_file, int invoking_line)
{
int error_code = ::MatMatSolve(A.Internal(),
B.Internal(),
X.Internal());
if (error_code)
throw ExceptionNS::Exception(error_code, "MatMatSolve", invoking_file, invoking_line);
}
} //namespace Petsc
} // namespace Wrappers
......
......@@ -566,7 +566,17 @@ namespace HappyHeart
throw ExceptionNS::Exception(error_code, "VecView", invoking_file, invoking_line);
}
void Vector::Load(const Mpi& mpi, const std::string& input_file, const char* invoking_file, int invoking_line,
PetscViewerFormat format) const
{
Viewer viewer(mpi, input_file, format, invoking_file, invoking_line);
int error_code = VecLoad(Internal(), viewer.GetUnderlyingPetscObject());
if (error_code)
throw ExceptionNS::Exception(error_code, "VecLoad", invoking_file, invoking_line);
}
void GatherVector(const Mpi& mpi,
......
......@@ -480,6 +480,20 @@ namespace HappyHeart
const std::string& output_file,
const char* invoking_file, int invoking_line) const;
/*!
* \brief Wrapper over VecLoad in the case the viewer is a file.
*
* \param[in] format Format in which the matrix is loaded. See Petsc manual pages to get all the
* formats available; relevant ones so far are PETSC_VIEWER_DEFAULT and PETSC_VIEWER_ASCII_MATLAB.
* \param[in] invoking_file File that invoked the function or class; usually __FILE__.
* \param[in] invoking_line File that invoked the function or class; usually __LINE__.
* \param[in] mpi Mpi object which knows the rank of the processor, the total number of processors, etc...
* \param[in] input_file File into which the vector content will be loaded.
*/
void Load(const Mpi& mpi, const std::string& input_file, const char* invoking_file, int invoking_line,
PetscViewerFormat format = PETSC_VIEWER_DEFAULT) const;
/*!
* \brief Get the minimum.
*
......
......@@ -18,9 +18,6 @@
# include "ThirdParty/Wrappers/Petsc/Vector/Vector.hpp"
# include "ThirdParty/Wrappers/Petsc/Matrix/Matrix.hpp"
# include "ThirdParty/Wrappers/Seldon/SubVector.hpp"
# include "ThirdParty/Wrappers/Seldon/HappyHeartPetscVector.hpp"
namespace HappyHeart
......@@ -36,12 +33,6 @@ namespace HappyHeart
//! Convenient alias for the local symmetrix matrices, used in some specific operators.
using LocalSymmetricMatrix = ::Seldon::Matrix<double, Seldon::General, Seldon::RowSym, Seldon::MallocAlloc<double>>;
//! Alias used when coupling HappyHeart with Verdandi.
using SeldonWrapperVector = ::Seldon::Vector<PetscScalar, Seldon::HappyHeartPetsc, Seldon::PetscPseudoAllocator>;
/*!
* \brief Enum to distinguish matrix and vectors.
*/
enum class IsMatrixOrVector
{
vector,
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment