Une MAJ de sécurité est nécessaire sur notre version actuelle. Elle sera effectuée lundi 02/08 entre 12h30 et 13h. L'interruption de service devrait durer quelques minutes (probablement moins de 5 minutes).

Commit e7480106 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#859 Petsc: add wrappers to use MatPtAP. However it doesn't seem to work at...

#859 Petsc: add wrappers to use MatPtAP. However it doesn't seem to work at the moment: PETSC_DEFAULT doesn't work for fill parameter, and the one I compute doesn't respect their constraint fill >= 1. I have contacted Petsc developers to help me on that topic; in the time being MatMatMult followed by MatTransposeMatMult do the job, albeit probably much less efficiently.
parent 6db9d6ed
......@@ -426,6 +426,8 @@ namespace HappyHeart
row_content_position_list.resize(size);
row_content_value_list.resize(size);
assert(values != nullptr);
for (auto i = 0ul; i < size; ++i)
{
row_content_position_list[i] = columns[i];
......
......@@ -115,6 +115,16 @@ namespace HappyHeart
* in the richer HappyHeart interface you should call in debug mode \a AssertMatrixRespectPattern()
* to make sure the resulting matrix respects the pattern defined for the \a GlobalMatrix pair
* of \a NumberingSubset.
*
* \internal <b><tt>[internal]</tt></b> If you know Petsc, you might see I didn't give access to
* argument MatReuse, setting it each time to MAT_INITIAL_MATRIX and skipping entirely MAT_REUSE_MATRIX.
* This is because at the time being MatMatMult operations are seldom in the code (only Poromechanics so
* far) and using Symbolic/Numeric seems more elegant. Of course, reintroducing the argument is really easy;
* feel free to do so if you need it (for instance for MatMatMatMult support: Symbolic/Numeric doesn't
* work for them and Petsc guys seemed unlikely to fix that in our exchanges).
*
* \internal <b><tt>[internal]</tt></b> \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.
*/
......@@ -126,14 +136,10 @@ namespace HappyHeart
*
* 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.
*
* \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.
* 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.
*
* \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
<
......@@ -215,13 +221,7 @@ namespace HappyHeart
* If the operation is performed many times with each time the same non zero pattern for the matrices,
* rather use MatMatMatMultSymbolic/MatMatMatMultNumeric to improve efficiency.
*
* \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.
* 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.
*
* \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__.
* \copydetails
*/
template
<
......@@ -244,7 +244,7 @@ namespace HappyHeart
/*!
* \brief Wrapper over MatMatMultSymbolic, that performs D = A * B * C.
*
* Actual computation is done with a call to NatMatMatMultNumeric, which might be called many times.
* Actual computation is done with a call to MatMatMatMultNumeric, which might be called many times.
*
* \param[in] matrix1 A in D = A * B * C.
* \param[in] matrix2 B in D = A * B * C.
......@@ -346,9 +346,6 @@ namespace HappyHeart
*
* \copydetails doxygen_hide_matmatmult_warning
*
* \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^T * B.
* \param[in] matrix2 B in C = A^T * B.
* \param[out] matrix3 C in B in C = A^T * B. The matrix must be not allocated when this function is called.
......@@ -377,9 +374,6 @@ namespace HappyHeart
*
* \copydetails doxygen_hide_matmatmult_warning
*
* \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.
......@@ -403,6 +397,107 @@ namespace HappyHeart
const char* invoking_file, int invoking_line);
/*!
* \class doxygen_hide_mat_pt_a_p
*
* \warning Unfortunately a simple test with P = I leads to error
* \verbatim
* [0]PETSC ERROR: Nonconforming object sizes
* [0]PETSC ERROR: Expected fill=-2 must be >= 1.0
* \endverbatim
* The reason is that PETSC_DEFAULT may not be supported (I've asked Petsc developers); but even with
* hand-tailored fill it doesn't seem to work...
* So unfortunately I have to advise instead MatMatMult followed by MatTransposeMatMult instead...
*
* \copydetails doxygen_hide_matmatmult_warning
*
* \param[in] A A in C = P^T * A * P.
* \param[in] P P in C = P^T * A * P
* \param[out] out C in C = P^T * A * P. The matrix must be not allocated when this function is called.
* \copydetails doxygen_hide_invoking_file_and_line
*/
/*!
* \brief Performs the matrix product C = P^T * A * P
*
* \copydoc doxygen_hide_mat_pt_a_p
*
* \attention If the operation is performed many times with each time the same non zero pattern for the matrices,
* rather use PtAPSymbolic/PtAPNumeric to improve efficiency.
*
*/
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
>
PtAP(const MatrixT& A,
const MatrixT& P,
MatrixU& out,
const char* invoking_file, int invoking_line);
/*!
* \brief Wrapper over PtAPSymbolic, that computes the structure of C = P^T * A * P.
*
* Actual computations must be performed through calls to PtAPNumeric.
* \see http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatPtAPSymbolic.html#MatPtAPSymbolic
* for more details.
*
* \copydoc doxygen_hide_mat_pt_a_p
*
*/
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
>
PtAPSymbolic(const MatrixT& A,
const MatrixT& P,
MatrixU& out,
const char* invoking_file, int invoking_line);
/*!
* \brief Wrapper over PtAPNumeric, that computes C = P^T * A * P. in the case C has already been allocated.
*
* \see http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatPtAPNumeric.html#MatPtAPNumeric
* for more details.
*
* \copydoc doxygen_hide_mat_pt_a_p
*
*/
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
>
PtAPNumeric(const MatrixT& A,
const MatrixT& P,
MatrixU& out,
const char* invoking_file, int invoking_line);
} //namespace Petsc
......
......@@ -11,6 +11,7 @@
#ifndef HAPPY_HEART_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_MATRIX_x_MATRIX_OPERATIONS_HXX_
# define HAPPY_HEART_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_MATRIX_x_MATRIX_OPERATIONS_HXX_
#include <iostream>
namespace HappyHeart
{
......@@ -385,6 +386,94 @@ namespace HappyHeart
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
>
PtAP(const MatrixT& A,
const MatrixT& P,
MatrixU& out,
const char* invoking_file, int invoking_line)
{
Mat result;
int error_code = ::MatPtAP(A.Internal(),
P.Internal(),
MAT_INITIAL_MATRIX,
PETSC_DEFAULT,
&result);
out.SetFromPetscMat(result);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatPtAP", 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
>
PtAPSymbolic(const MatrixT& A,
const MatrixT& P,
MatrixU& out,
const char* invoking_file, int invoking_line)
{
Mat result;
int error_code = ::MatPtAPSymbolic(A.Internal(),
P.Internal(),
PETSC_DEFAULT,
&result);
out.SetFromPetscMat(result);
if (error_code)
throw ExceptionNS::Exception(error_code, "MatPtAPSymbolic", 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
>
PtAPNumeric(const MatrixT& A,
const MatrixT& P,
MatrixU& out,
const char* invoking_file, int invoking_line)
{
int error_code = ::MatPtAPNumeric(A.Internal(),
P.Internal(),
out.Internal());
if (error_code)
throw ExceptionNS::Exception(error_code, "MatPtAPNumeric", invoking_file, invoking_line);
}
} //namespace Petsc
......
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