#888 Petsc/Wrappers: add wrapper for MatMatMultSymbolic/NUmeric and MatMatMatMultSymbolic/Numeric.

Sources/ThirdParty/Wrappers/Petsc/Matrix/Matrix.hpp

@@ -139,7 +139,6 @@ namespace HappyHeart
                                           const Mpi& mpi, const char* invoking_file, int invoking_line);
                 
                 
-                
                 /*!
                  * \brief Create a parallel sparse matrix.
                  *

Sources/ThirdParty/Wrappers/Petsc/Matrix/Matrix.hxx

@@ -24,7 +24,6 @@ namespace HappyHeart
         namespace Petsc
         {
 
-
             
             inline Mat Matrix::Internal() const noexcept
             {
@@ -33,8 +32,6 @@ namespace HappyHeart
             }
             
             
-              
-            
         } // namespace Petsc
         
         

Sources/ThirdParty/Wrappers/Petsc/Matrix/MatrixOperations.hpp

@@ -111,6 +111,9 @@ namespace HappyHeart
             /*!
              * \brief Wrapper over MatMatMult, that performs C = A * B.
              *
+             * 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
@@ -136,9 +139,67 @@ namespace HappyHeart
                        const char* invoking_file, int invoking_line);
             
             
+            /*!
+             * \brief Wrapper over MatMatMultSymbolic, that computes the structure of C = A * B. 
+             *
+             * Actual computation is done with a call to MatMatMultNumeric, which might be called many times.
+             *
+             * \param[in] matrix1 A in C = A * B.
+             * \param[in] matrix2 B in C = A * B.
+             * \param[out] matrix3 C in C = A * B. Its pattern is correctly allocated after the call.
+             * \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
+            >
+            std::enable_if_t
+            <
+                std::is_base_of<Private::BaseMatrix, MatrixT>::value
+                && std::is_base_of<Private::BaseMatrix, MatrixU>::value,
+                void
+            >
+            MatMatMultSymbolic(const MatrixT& matrix1,
+                               const MatrixT& matrix2,
+                               MatrixU& matrix3,
+                               const char* invoking_file, int invoking_line);
+            
+            
+            /*!
+             * \brief Wrapper over MatMatMultNumeric, that computes C = A * B in the case MatMatMultSymbolic
+             * has been used to allocate structure of C.
+             *
+             * \param[in] matrix1 A in C = A * B.
+             * \param[in] matrix2 B in C = A * B.
+             * \param[out] matrix3 C in C = A * B.
+             * \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
+            >
+            std::enable_if_t
+            <
+                std::is_base_of<Private::BaseMatrix, MatrixT>::value
+                && std::is_base_of<Private::BaseMatrix, MatrixU>::value,
+                void
+            >
+            MatMatMultNumeric(const MatrixT& matrix1,
+                              const MatrixT& matrix2,
+                              MatrixU& matrix3,
+                              const char* invoking_file, int invoking_line);
+            
+            
             /*!
              * \brief Wrapper over MatMatMatMult, that performs D = A * B * C.
              *
+             * 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
@@ -165,6 +226,77 @@ namespace HappyHeart
                           const char* invoking_file, int invoking_line);
             
             
+            /*!
+             * \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.
+             *
+             * \param[in] matrix1 A in D = A * B * C.
+             * \param[in] matrix2 B in D = A * B * C.
+             * \param[in] matrix3 C in D = A * B * C.
+             * \param[out] matrix4 D in D = A * B * C. Its pattern is correctly allocated after the call.
+
+             * \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
+            <
+                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
+            >
+            MatMatMatMultSymbolic(const MatrixT& matrix1,
+                                  const MatrixT& matrix2,
+                                  const MatrixT& matrix3,
+                                  MatrixU& matrix4,
+                                  const char* invoking_file, int invoking_line);
+            
+            
+            /*!
+             * \brief Wrapper over MatMatMatMultNumeric, that performs D = A * B * C in the case MatMatMatMultSymbolic
+             * has been used to allocate structure of D.
+             *
+             * \param[in] matrix1 A in D = A * B * C.
+             * \param[in] matrix2 B in D = A * B * C.
+             * \param[in] matrix3 C in D = A * B * C.
+             * \param[out] matrix4 D in D = A * B * C.
+             *
+             * \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
+            <
+                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
+            >
+            MatMatMatMultNumeric(const MatrixT& matrix1,
+                                 const MatrixT& matrix2,
+                                 const MatrixT& matrix3,
+                                 MatrixU& matrix4,
+                                 const char* invoking_file, int invoking_line);
+            
+            
             /*!
              * \brief Creates a new matrix object that behaves like A'.
              *

Sources/ThirdParty/Wrappers/Petsc/Matrix/MatrixOperations.hxx

@@ -82,6 +82,64 @@ 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
+            >
+            MatMatMultSymbolic(const MatrixT& matrix1,
+                               const MatrixT& matrix2,
+                               MatrixU& matrix3,
+                               const char* invoking_file, int invoking_line)
+            {
+                Mat result;
+                
+                int error_code = ::MatMatMultSymbolic(matrix1.Internal(),
+                                                      matrix2.Internal(),
+                                                      PETSC_DEFAULT,
+                                                      &result);
+                
+                
+                matrix3.SetFromPetscMat(result);
+                
+                if (error_code)
+                    throw ExceptionNS::Exception(error_code, "MatMatMultSymbolic", 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
+            >
+            MatMatMultNumeric(const MatrixT& matrix1,
+                              const MatrixT& matrix2,
+                              MatrixU& matrix3,
+                              const char* invoking_file, int invoking_line)
+            {
+                int error_code = ::MatMatMultNumeric(matrix1.Internal(),
+                                                     matrix2.Internal(),
+                                                     matrix3.Internal());
+                
+                if (error_code)
+                    throw ExceptionNS::Exception(error_code, "MatMatMultNumeric", invoking_file, invoking_line);
+            }
+
+
+            
+            
             template
             <
                 class MatrixT,
@@ -185,7 +243,70 @@ namespace HappyHeart
                 
                 transpose.SetFromPetscMat(result);
             }
-
+            
+            
+            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
+            >
+            MatMatMatMultSymbolic(const MatrixT& matrix1,
+                                  const MatrixT& matrix2,
+                                  const MatrixT& matrix3,
+                                  MatrixU& matrix4,
+                                  const char* invoking_file, int invoking_line)
+            
+            
+            
+            {
+                Mat result;
+                
+                int error_code = ::MatMatMatMultSymbolic(matrix1.Internal(),
+                                                         matrix2.Internal(),
+                                                         matrix3.Internal(),
+                                                         matrix4.Internal(),
+                                                         PETSC_DEFAULT,
+                                                         &result);
+                
+                
+                matrix4.SetFromPetscMat(result);
+                
+                if (error_code)
+                    throw ExceptionNS::Exception(error_code, "MatMatMatMultSymbolic", 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
+            >
+            MatMatMatMultNumeric(const MatrixT& matrix1,
+                                 const MatrixT& matrix2,
+                                 const MatrixT& matrix3,
+                                 MatrixU& matrix4,
+                                 const char* invoking_file, int invoking_line)
+            {
+                int error_code = ::MatMatMatMultNumeric(matrix1.Internal(),
+                                                        matrix2.Internal(),
+                                                        matrix3.Internal(),
+                                                        matrix4.Internal());
+                
+                if (error_code)
+                    throw ExceptionNS::Exception(error_code, "MatMatMatMultNumeric", invoking_file, invoking_line);
+            }
             
             
         } //namespace Petsc