MatrixOperations.hpp 23 KB
Newer Older
1 2 3 4 5 6 7 8 9 10
///
////// \file
///
///
/// Created by Sebastien Gilles <sebastien.gilles@inria.fr> on the Fri, 30 Oct 2015 12:41:42 +0100
/// Copyright (c) Inria. All rights reserved.
///
/// \ingroup ThirdPartyGroup
/// \addtogroup ThirdPartyGroup
/// \{
11 12 13 14 15 16 17

#ifndef HAPPY_HEART_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_MATRIX_x_MATRIX_OPERATIONS_HPP_
# define HAPPY_HEART_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_MATRIX_x_MATRIX_OPERATIONS_HPP_

# include <memory>
# include <vector>

18
# include "Utilities/Containers/EnumClass.hpp"
19

20 21
# include "ThirdParty/IncludeWithoutWarning/Petsc/PetscSys.hpp"
# include "ThirdParty/IncludeWithoutWarning/Petsc/PetscMat.hpp"
22

23 24
# include "ThirdParty/Wrappers/Petsc/Matrix/NonZeroPattern.hpp"
# include "ThirdParty/Wrappers/Petsc/Vector/Vector.hpp"
25
# include "ThirdParty/Wrappers/Petsc/Matrix/Internal/BaseMatrix.hpp"
26
# include "ThirdParty/Wrappers/Petsc/Exceptions/Petsc.hpp"
27 28 29 30 31 32


namespace HappyHeart
{
    
    
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
    // ============================
    //! \cond IGNORE_BLOCK_IN_DOXYGEN
    // Forward declarations.
    // ============================
    
    
    
    // Matrix actually used in HappyHeart, that wraps a Wrappers::Petsc::Matrix.
    class GlobalMatrix;
    

    // ============================
    // End of forward declarations.
    //! \endcond IGNORE_BLOCK_IN_DOXYGEN
    // ============================


50 51 52 53 54 55
    namespace Wrappers
    {
        
        
        namespace Petsc
        {
56 57
            
            
58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
            /*!
             * \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 };
78
                   
79 80 81
            /*!
             * \brief Wrapper over MatAXPY, that performs Y = a * X + Y.
             *
82
             *
83
             * \tparam NonZeroPatternT This value indicates the level of similarity between X and Y non zero patterns.
84
             * \warning Beware: if X and Y aren't following the same pattern, 'Same' won't yield what
85
             * you expect! In the case you're considering adding a part computed by a transfert matrix, you should
86
             * use 'Subset'. In HappyHeart, I advise you to call in debug mode AssertSameNumberingSubset()
87
             * (defined in Core) just before or after the function call; if the assert is raised it means you can't
88
             * use 'Same'. \todo This could be enforced by a proper overload (but it mix Utilities and Core
89
             * libraries...), but for the time being the separate assert will do.
90 91 92 93 94
             * \copydoc doxygen_hide_invoking_file_and_line
             *
             * \param[in] a See formula above.
             * \param[in] X See formula above.
             * \param[in] Y See formula above.
95
             */
96 97
            template
            <
98
                NonZeroPattern NonZeroPatternT,
99 100 101 102 103
                class MatrixT,
                class MatrixU
            >
            std::enable_if_t
            <
104 105
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
             && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value,
106 107
                void
            >
108
            AXPY(PetscScalar a,
109
                 const MatrixT& X,
110
                 MatrixU& Y,
111
                 const char* invoking_file, int invoking_line);
112
            
113 114 115
            /*!
             * \brief Wrapper over MatShift, that performs M = M + a I, where a is a PetscScalar and I is the identity matrix.
             *
116 117 118
             * \param[in] matrix See formula above.
             * \param[in] a See formula above.
             * \copydoc doxygen_hide_invoking_file_and_line
119 120 121 122 123
             */
            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);
            
124
            
125 126 127
            /*!
             * \brief Wrapper over MatShift, that performs M = M + a I, where a is a PetscScalar and I is the identity matrix.
             *
128 129 130
             * \param[in] matrix See formula above.
             * \param[in] a See formula above.
             * \copydoc doxygen_hide_invoking_file_and_line
131 132 133 134 135 136
             */
            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);
            
            
137 138 139
            /*!
             * \brief Wrapper over MatMult, that performs v2 = matrix * v1.
             *
140 141 142 143
             * \copydoc doxygen_hide_invoking_file_and_line
             * \param[in] matrix See formula above.
             * \param[in] v1 See formula above.
             * \param[in] v2 See formula above.
144
             * \copydoc doxygen_hide_do_update_ghost_arg
145 146
             */
            template<class MatrixT>
147
            std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
148 149 150
            MatMult(const MatrixT& matrix,
                    const Vector& v1,
                    Vector& v2,
151 152
                    const char* invoking_file, int invoking_line,
                    update_ghost do_update_ghost = update_ghost::yes);
153
            
154 155 156
            /*!
             * \brief Wrapper over MatTranspose, that performs in-place or out-of-place transpose of a matrix.
             *
157 158 159
             * \copydoc doxygen_hide_invoking_file_and_line
             * \param[in] matrix1 Matrix to transpose.
             * \param[in] matrix2 Matrix to contain the result of the transpose.
160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176
             *
             */
            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);
            
177 178 179 180
            
            /*!
             * \brief Wrapper over MatMultAdd, that performs v3 = v2 + matrix * v1.
             *
181 182 183 184 185
             * \copydoc doxygen_hide_invoking_file_and_line
             * \param[in] matrix See formula above.
             * \param[in] v1 See formula above.
             * \param[in] v2 See formula above.
             * \param[in] v3 See formula above.
186
             * \copydoc doxygen_hide_do_update_ghost_arg
187 188
             */
            template<class MatrixT>
189
            std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
190 191 192 193
            MatMultAdd(const MatrixT& matrix,
                       const Vector& v1,
                       const Vector& v2,
                       Vector& v3,
194 195
                       const char* invoking_file, int invoking_line,
                       update_ghost do_update_ghost = update_ghost::yes);
196 197 198 199 200 201
            
            
            
            /*!
             * \brief Wrapper over MatMultTranspose, that performs v2 = transpose(matrix) * v1.
             *
202
             * \copydoc doxygen_hide_invoking_file_and_line
203
             * \copydoc doxygen_hide_do_update_ghost_arg
204 205 206
             * \param[in] matrix See formula above.
             * \param[in] v1 See formula above.
             * \param[in] v2 See formula above.
207 208
             */
            template<class MatrixT>
209
            std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
210 211 212
            MatMultTranspose(const MatrixT& matrix,
                             const Vector& v1,
                             Vector& v2,
213 214
                             const char* invoking_file, int invoking_line,
                             update_ghost do_update_ghost = update_ghost::yes);
215 216 217 218 219
            
            
            /*!
             * \brief Wrapper over MatMultTransposeAdd, that performs v3 = v2 + transpose(matrix) * v1.
             *
220
             * \copydoc doxygen_hide_invoking_file_and_line
221
             * \copydoc doxygen_hide_do_update_ghost_arg
222 223 224 225
             * \param[in] matrix See formula above.
             * \param[in] v1 See formula above.
             * \param[in] v2 See formula above.
             * \param[in] v3 See formula above.
226 227
             */
            template<class MatrixT>
228
            std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
229 230 231 232
            MatMultTransposeAdd(const MatrixT& matrix,
                                const Vector& v1,
                                const Vector& v2,
                                Vector& v3,
233 234
                                const char* invoking_file, int invoking_line,
                                update_ghost do_update_ghost = update_ghost::yes);
235
            
236 237 238 239 240 241
            /*!
             * \class doxygen_hide_matmatmult_warning
             *
             * \attention In Petsc, matrix-matrix multiplication functions compute on the fly the
             * pattern for the resulting matrix. However, it is possible this pattern isn't the one expected;
             * in the richer HappyHeart interface you should call in debug mode \a AssertMatrixRespectPattern()
242
             * to make sure the resulting matrix respects the pattern defined for the \a GlobalMatrix pair
243
             * of \a NumberingSubset.
244 245 246 247 248
             *
             * \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;
249
             * feel free to do so if you need it (for instance for MatMatMatMult support: Symbolic/Numeric doesn't
250 251
             * work for them and Petsc guys seemed unlikely to fix that in our exchanges).
             *
252
             * \internal <b><tt>[internal]</tt></b> \todo #684 Investigate to use the argument fill, which provides an
253
             * estimation of the non zero of the resulting matrix. Currently PETSC_DEFAULT is used.
254 255 256
             */
            
            
257 258
            
            /*!
259
             * \brief Wrapper over MatMatMult, that performs m3 = m1 * m2.
260
             *
261 262
             * \copydetails doxygen_hide_matmatmult_warning
             *
263 264
             * 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.
265
             
266 267
             * 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.
268 269 270
             * \param[in] m1 See formula above.
             * \param[in] m2 See formula above.
             * \param[in] m3 See formula above.
271
             *
272
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
273
             * \copydoc doxygen_hide_invoking_file_and_line
274
             */
275 276
            template
            <
277 278
            class MatrixT,
            class MatrixU
279 280 281
            >
            std::enable_if_t
            <
282 283
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value,
284 285
                void
            >
286 287 288
            MatMatMult(const MatrixT& m1,
                       const MatrixT& m2,
                       MatrixU& m3,
289 290
                       const char* invoking_file, int invoking_line,
                       DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
291
            
292
                      
293
            /*!
294
             * \brief Wrapper over MatMatMatMult, that performs m4 = m1 * m2 * m3.
295
             *
296 297 298 299
             * \param[in] m1 See formula above.
             * \param[in] m2 See formula above.
             * \param[in] m3 See formula above.
             * \param[in] m4 See formula above.
300
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
301
             * \copydoc doxygen_hide_invoking_file_and_line
302
             */
303 304
            template
            <
305 306
            class MatrixT,
            class MatrixU
307 308 309
            >
            std::enable_if_t
            <
310 311
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value,
312 313
                void
            >
314 315 316 317
            MatMatMatMult(const MatrixT& m1,
                          const MatrixT& m2,
                          const MatrixT& m3,
                          MatrixU& m4,
318 319
                          const char* invoking_file, int invoking_line,
                          DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
320
            
321 322 323
            /*!
             * \brief Creates a new matrix object that behaves like A'.
             *
324
             * The transpose A' is NOT actually formed! Rather the new matrix object performs the matrix-vector product
325 326 327
             * by using the MatMultTranspose() on the original matrix.
             *
             * \param[in] A matrix to transpose.
328
             * \param[out] transpose The matrix that figuratively represents A'. This matrix must not have been
329
             * allocated!
330
             * \copydoc doxygen_hide_invoking_file_and_line
331 332 333 334
             */
            template<class MatrixT>
            std::enable_if_t
            <
335
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value,
336 337 338 339 340 341
                void
            >
            MatCreateTranspose(const MatrixT& A,
                               MatrixT& transpose,
                               const char* invoking_file, int invoking_line);
            
342
            
343
            /*!
344
             * \class doxygen_hide_mat_transpose_mat_mult
345
             *
346 347 348 349 350 351 352 353 354 355
             * Formula is:
             * \verbatim
             * m3 = m1^T * m2.
             * \endverbatim
             *
             * \copydetails doxygen_hide_matmatmult_warning
             * \param[in] m1 See formula above.
             * \param[in] m2 See formula above.
             * \param[out] m3 See formula above. The matrix must be not allocated when this function is called.
             * \copydoc doxygen_hide_invoking_file_and_line
356
             */
357 358 359
            
            
            /*!
360
             * \brief Performs Matrix-Matrix Multiplication m3 = m1^T * m2.
361 362
             *
             * \copydetails doxygen_hide_mat_transpose_mat_mult
363
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
364 365 366
             * \param[in] m1 See formula above.
             * \param[in] m2 See formula above.
             * \param[in] m3 See formula above.
367
             */
368 369
            template
            <
370 371
            class MatrixT,
            class MatrixU
372 373 374
            >
            std::enable_if_t
            <
375 376
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value,
377 378
                void
            >
379 380 381
            MatTransposeMatMult(const MatrixT& m1,
                                const MatrixT& m2,
                                MatrixU& m3,
382 383
                                const char* invoking_file, int invoking_line,
                                DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
384 385
            
            
386
            
387 388 389
            /*!
             * \brief Performs Matrix-Matrix Multiplication C = A * B^T.
             *
390
             * \copydetails doxygen_hide_matmatmult_warning
391
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
392
             *
393 394 395
             * \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.
396
             * \copydoc doxygen_hide_invoking_file_and_line
397 398 399
             */
            template
            <
400 401
            class MatrixT,
            class MatrixU
402 403 404
            >
            std::enable_if_t
            <
405 406
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value,
407 408 409 410 411
                void
            >
            MatMatTransposeMult(const MatrixT& matrix1,
                                const MatrixT& matrix2,
                                MatrixU& matrix3,
412 413
                                const char* invoking_file, int invoking_line,
                                DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
414 415
            
            
416 417 418
            /*!
             * \class doxygen_hide_mat_pt_a_p
             *
419
             * \warning Unfortunately a simple test with P = I leads to error
420 421 422 423
             * \verbatim
             * [0]PETSC ERROR: Nonconforming object sizes
             * [0]PETSC ERROR: Expected fill=-2 must be >= 1.0
             * \endverbatim
424
             * The reason is that PETSC_DEFAULT may not be supported (I've asked Petsc developers); but even with
425 426 427 428 429 430 431 432 433
             * 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
434
             * \copydoc doxygen_hide_petsc_do_reuse_matrix_arg
435 436 437 438 439 440 441 442 443 444 445
             */
            
            /*!
             * \brief Performs the matrix product C = P^T * A * P
             *
             * \copydoc doxygen_hide_mat_pt_a_p
             *
             *
             */
            template
            <
446 447
            class MatrixT,
            class MatrixU
448 449 450
            >
            std::enable_if_t
            <
451 452
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value,
453 454 455 456 457
                void
            >
            PtAP(const MatrixT& A,
                 const MatrixT& P,
                 MatrixU& out,
458 459
                 const char* invoking_file, int invoking_line,
                 DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
460 461
            
            
462 463 464 465 466 467 468
            
            /*!
             * \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.
             *
469 470 471 472 473
             * \param[in] A Matrix to get the ordering.
             * \param[in] type Type of the ordering.
             * \param[in] rperm Row permutation for the ordering.
             * \param[in] cperm Column permutation for the ordering.
             * \copydetails doxygen_hide_invoking_file_and_line             
474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490
             *
             */
            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.
             *
491 492 493 494 495
             * \param[in] A Matrix to factor.
             * \param[in] row Result of the GetOrdering.
             * \param[in] col Result of the GetOrdering.
             * \param[in] info Info of LUFactor.
             * \copydetails doxygen_hide_invoking_file_and_line
496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511
             *
             */
            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.
             *
512 513 514 515
             * \param[in] A See formula above.
             * \param[in] B See formula above.
             * \param[in] X See formula above.
             * \copydetails doxygen_hide_invoking_file_and_line
516 517 518 519
             *
             */
            template
            <
GILLES Sebastien's avatar
GILLES Sebastien committed
520 521 522
                class MatrixT,
                class MatrixU,
                class MatrixV
523 524 525
            >
            std::enable_if_t
            <
GILLES Sebastien's avatar
GILLES Sebastien committed
526 527 528 529
                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
530 531 532 533 534 535 536 537 538 539
            >
            MatMatSolve(const MatrixT& A,
                        const MatrixU& B,
                        MatrixV& X,
                        const char* invoking_file, int invoking_line);


            
            
        } //namespace Petsc
540 541
        
        
542
    } // namespace Wrappers
543 544
    
    
545 546 547 548 549 550 551 552 553 554
    
    /*! 
     * \brief Convenient alias to avoid repeating the namespaces in each call.
     *
     * It should be used anyway only as template arguments of functions within this namespaces, such as 
     * Wrappers::Petsc::AXPY.
     */
    using NonZeroPattern = Wrappers::Petsc::NonZeroPattern;
    
    
555 556 557
} // namespace HappyHeart


558 559 560
/// @} // addtogroup ThirdPartyGroup


561
# include "ThirdParty/Wrappers/Petsc/Matrix/MatrixOperations.hxx"
562 563 564


#endif // HAPPY_HEART_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_MATRIX_x_MATRIX_OPERATIONS_HPP_