Skip to content

GitLab

Commit 787f66b0 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#859 Petsc: I wasn't at ease with the Numeric/Symbolic hack for... 

#859 Petsc: I wasn't at ease with the Numeric/Symbolic hack for MatTransposeMatMult, so I rewrote it with explicit reference to the function argument used.
parent 9fd2c1bd
Hide whitespace changes
Inline Side-by-side
Showing with 46 additions and 131 deletions
Sources/ThirdParty/Wrappers/Petsc/Matrix/MatrixOperations.hpp
View file @ 787f66b0
......@@ -34,6 +34,28 @@ namespace HappyHeart
{
/*!
* \class doxygen_hide_petsc_do_reuse_matrix
*
* Whether matrix must be reused or initialized from scratch in functions such
* as MatTransposeMatMult.
* If yes, it must have been called once with no beforehand, and one has to be sure the pattern
* of the matrices involved is still the same.
*/
/*!
* \class doxygen_hide_petsc_do_reuse_matrix_arg
*
* \param[in] do_reuse_matrix \copydetails doxygen_hide_petsc_do_reuse_matrix
*/
/*!
* \brief \copydoc doxygen_hide_petsc_do_reuse_matrix
*/
enum class DoReuseMatrix { yes, no };
/*!
* \brief Wrapper over MatAXPY, that performs Y = a * X + Y.
*
......@@ -357,7 +379,7 @@ namespace HappyHeart
* \brief Performs Matrix-Matrix Multiplication C = A^T * B.
*
* \copydetails doxygen_hide_mat_transpose_mat_mult
*
* \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
*/
template
<
......@@ -373,69 +395,16 @@ namespace HappyHeart
MatTransposeMatMult(const MatrixT& matrix1,
const MatrixT& matrix2,
MatrixU& matrix3,
const char* invoking_file, int invoking_line);
/*!
* \class doxygen_hide_matrix_transpose_matrix_mult_hack
*
* \warning Petsc doesn't propose this handy function, so I hacked it with MatReuse = MAT_INITIAL_MATRIX
* in first (Symbolic) call and MatReuse = MAT_REUSE_MATRIX in subsequent calls. This works in my case...
* but is still clear ly hack.
*/
/*!
* \brief Performs Matrix-Matrix Multiplication C = A^T * B in order to compute structure of C
*
* \copydetails doxygen_hide_matrix_transpose_matrix_mult_hack
* \copydetails doxygen_hide_mat_transpose_mat_mult
*/
template
<
class MatrixT,
class MatrixU
>
std::enable_if_t
<
std::is_base_of<Private::BaseMatrix, MatrixT>::value
&& std::is_base_of<Private::BaseMatrix, MatrixU>::value,
void
>
MatTransposeMatMultSymbolic(const MatrixT& matrix1,
const MatrixT& matrix2,
MatrixU& matrix3,
const char* invoking_file, int invoking_line);
/*!
* \brief Performs Matrix-Matrix Multiplication C = A^T * B assuming C pattern is already in place.
*
* \copydetails doxygen_hide_matrix_transpose_matrix_mult_hack
* \copydetails doxygen_hide_mat_transpose_mat_mult
*/
template
<
class MatrixT,
class MatrixU
>
std::enable_if_t
<
std::is_base_of<Private::BaseMatrix, MatrixT>::value
&& std::is_base_of<Private::BaseMatrix, MatrixU>::value,
void
>
MatTransposeMatMultNumeric(const MatrixT& matrix1,
const MatrixT& matrix2,
MatrixU& matrix3,
const char* invoking_file, int invoking_line);
const char* invoking_file, int invoking_line,
DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
/*!
* \brief Performs Matrix-Matrix Multiplication C = A * B^T.
*
* \copydetails doxygen_hide_matmatmult_warning
* \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
*
* \param[in] matrix1 A in C = A * B^T.
* \param[in] matrix2 B in C = A * B^T.
......@@ -457,7 +426,8 @@ namespace HappyHeart
MatMatTransposeMult(const MatrixT& matrix1,
const MatrixT& matrix2,
MatrixU& matrix3,
const char* invoking_file, int invoking_line);
const char* invoking_file, int invoking_line,
DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
/*!
......
Sources/ThirdParty/Wrappers/Petsc/Matrix/MatrixOperations.hxx
View file @ 787f66b0
......@@ -325,92 +325,32 @@ namespace HappyHeart
MatTransposeMatMult(const MatrixT& matrix1,
const MatrixT& matrix2,
MatrixU& matrix3,
const char* invoking_file, int invoking_line)
const char* invoking_file, int invoking_line,
DoReuseMatrix do_reuse_matrix)
{
Mat result;
int error_code = ::MatTransposeMatMult(matrix1.Internal(),
matrix2.Internal(),
MAT_INITIAL_MATRIX,
PETSC_DEFAULT,
&result);
matrix3.SetFromPetscMat(result);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatTransposeMatMult", invoking_file, invoking_line);
}
template
<
class MatrixT,
class MatrixU
>
inline std::enable_if_t
<
std::is_base_of<Private::BaseMatrix, MatrixT>::value
&& std::is_base_of<Private::BaseMatrix, MatrixU>::value,
void
>
MatTransposeMatMultSymbolic(const MatrixT& matrix1,
const MatrixT& matrix2,
MatrixU& matrix3,
const char* invoking_file, int invoking_line)
{
MatTransposeMatMult(matrix1,
matrix2,
matrix3,
invoking_file, invoking_line);
}
template
<
class MatrixT,
class MatrixU
>
std::enable_if_t
<
std::is_base_of<Private::BaseMatrix, MatrixT>::value
&& std::is_base_of<Private::BaseMatrix, MatrixU>::value,
void
>
MatTransposeMatMultNumeric(const MatrixT& matrix1,
const MatrixT& matrix2,
MatrixU& matrix3,
const char* invoking_file, int invoking_line)
{
Mat result = matrix3.Internal();
if (do_reuse_matrix == DoReuseMatrix::yes)
result = matrix3.Internal();
int error_code = ::MatTransposeMatMult(matrix1.Internal(),
matrix2.Internal(),
MAT_REUSE_MATRIX,
MAT_INITIAL_MATRIX,
PETSC_DEFAULT,
&result);
if (do_reuse_matrix == DoReuseMatrix::no)
matrix3.SetFromPetscMat(result);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatTransposeMatMult", invoking_file, invoking_line);
}
/*!
* \brief Performs Matrix-Matrix Multiplication C = A * B^T.
*
* \todo #684 Investigate to use the argument fill, which provides an estimation of the non zero of the
* resulting matrix. Currently PETSC_DEFAULT is used.
*
* \param[in] matrix1 A in C = A * B^T.
* \param[in] matrix2 B in C = A * B^T.
* \param[out] matrix3 C in B in C = A * B^T. The matrix must be not allocated when this function is called.
* \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,
class MatrixU
class MatrixU
>
std::enable_if_t
<
......@@ -421,17 +361,22 @@ namespace HappyHeart
MatMatTransposeMult(const MatrixT& matrix1,
const MatrixT& matrix2,
MatrixU& matrix3,
const char* invoking_file, int invoking_line)
const char* invoking_file, int invoking_line,
DoReuseMatrix do_reuse_matrix)
{
Mat result;
if (do_reuse_matrix == DoReuseMatrix::yes)
result = matrix3.Internal();
int error_code = ::MatMatTransposeMult(matrix1.Internal(),
matrix2.Internal(),
MAT_INITIAL_MATRIX,
PETSC_DEFAULT,
&result);
matrix3.SetFromPetscMat(result);
if (do_reuse_matrix == DoReuseMatrix::no)
matrix3.SetFromPetscMat(result);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatMatTransposeMult", invoking_file, invoking_line);
......
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