MatrixOperations.hpp 23.5 KB
Newer Older
GILLES Sebastien's avatar
GILLES Sebastien committed
1 2 3 4 5 6 7 8 9 10 11 12 13
/*!
//
// \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
// \{
*/

14

15 16
#ifndef MOREFEM_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_MATRIX_x_MATRIX_OPERATIONS_HPP_
# define MOREFEM_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_MATRIX_x_MATRIX_OPERATIONS_HPP_
17 18 19 20

# include <memory>
# include <vector>

21
# include "Utilities/Containers/EnumClass.hpp"
22

23 24
# include "ThirdParty/IncludeWithoutWarning/Petsc/PetscSys.hpp"
# include "ThirdParty/IncludeWithoutWarning/Petsc/PetscMat.hpp"
25

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


32
namespace MoReFEM
33
{
34 35


36 37 38 39
    // ============================
    //! \cond IGNORE_BLOCK_IN_DOXYGEN
    // Forward declarations.
    // ============================
40 41 42



43
    // Matrix actually used in MoReFEM, that wraps a Wrappers::Petsc::Matrix.
44
    class GlobalMatrix;
45

46 47 48 49 50 51 52

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


53 54
    namespace Wrappers
    {
55 56


57 58
        namespace Petsc
        {
59 60


61 62 63 64 65 66 67 68
            /*!
             * \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.
             */
69

70 71 72 73 74
            /*!
             * \class doxygen_hide_petsc_do_reuse_matrix_arg
             *
             * \param[in] do_reuse_matrix \copydetails doxygen_hide_petsc_do_reuse_matrix
             */
75 76


77 78 79 80
            /*!
             * \brief \copydoc doxygen_hide_petsc_do_reuse_matrix
             */
            enum class DoReuseMatrix { yes, no };
81

82 83 84
            /*!
             * \brief Wrapper over MatAXPY, that performs Y = a * X + Y.
             *
85
             *
86
             * \tparam NonZeroPatternT This value indicates the level of similarity between X and Y non zero patterns.
87
             * \warning Beware: if X and Y aren't following the same pattern, 'Same' won't yield what
88
             * you expect! In the case you're considering adding a part computed by a transfert matrix, you should
89
             * use 'Subset'. In MoReFEM, I advise you to call in debug mode AssertSameNumberingSubset()
90
             * (defined in Core) just before or after the function call; if the assert is raised it means you can't
91
             * use 'Same'. \todo This could be enforced by a proper overload (but it mix Utilities and Core
92
             * libraries...), but for the time being the separate assert will do.
93 94 95 96 97
             * \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.
98
             */
99 100
            template
            <
101
                NonZeroPattern NonZeroPatternT,
102 103 104 105 106
                class MatrixT,
                class MatrixU
            >
            std::enable_if_t
            <
107 108
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
             && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value,
109 110
                void
            >
111
            AXPY(PetscScalar a,
112
                 const MatrixT& X,
113
                 MatrixU& Y,
114
                 const char* invoking_file, int invoking_line);
115 116


117 118 119
            /*!
             * \brief Wrapper over MatShift, that performs M = M + a I, where a is a PetscScalar and I is the identity matrix.
             *
120 121 122
             * \param[in] matrix See formula above.
             * \param[in] a See formula above.
             * \copydoc doxygen_hide_invoking_file_and_line
123 124 125 126
             */
            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);
127 128


129 130 131
            /*!
             * \brief Wrapper over MatMult, that performs v2 = matrix * v1.
             *
132 133 134 135
             * \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.
136
             * \copydoc doxygen_hide_do_update_ghost_arg
137 138
             */
            template<class MatrixT>
139
            std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
140 141 142
            MatMult(const MatrixT& matrix,
                    const Vector& v1,
                    Vector& v2,
143 144
                    const char* invoking_file, int invoking_line,
                    update_ghost do_update_ghost = update_ghost::yes);
145

146 147 148
            /*!
             * \brief Wrapper over MatTranspose, that performs in-place or out-of-place transpose of a matrix.
             *
149 150 151
             * \copydoc doxygen_hide_invoking_file_and_line
             * \param[in] matrix1 Matrix to transpose.
             * \param[in] matrix2 Matrix to contain the result of the transpose.
152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167
             *
             */
            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);
168 169


170 171 172
            /*!
             * \brief Wrapper over MatMultAdd, that performs v3 = v2 + matrix * v1.
             *
173 174 175 176 177
             * \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.
178
             * \copydoc doxygen_hide_do_update_ghost_arg
179 180
             */
            template<class MatrixT>
181
            std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
182 183 184 185
            MatMultAdd(const MatrixT& matrix,
                       const Vector& v1,
                       const Vector& v2,
                       Vector& v3,
186 187
                       const char* invoking_file, int invoking_line,
                       update_ghost do_update_ghost = update_ghost::yes);
188 189 190



191 192 193
            /*!
             * \brief Wrapper over MatMultTranspose, that performs v2 = transpose(matrix) * v1.
             *
194
             * \copydoc doxygen_hide_invoking_file_and_line
195
             * \copydoc doxygen_hide_do_update_ghost_arg
196 197 198
             * \param[in] matrix See formula above.
             * \param[in] v1 See formula above.
             * \param[in] v2 See formula above.
199 200
             */
            template<class MatrixT>
201
            std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
202 203 204
            MatMultTranspose(const MatrixT& matrix,
                             const Vector& v1,
                             Vector& v2,
205 206
                             const char* invoking_file, int invoking_line,
                             update_ghost do_update_ghost = update_ghost::yes);
207 208


209 210 211
            /*!
             * \brief Wrapper over MatMultTransposeAdd, that performs v3 = v2 + transpose(matrix) * v1.
             *
212
             * \copydoc doxygen_hide_invoking_file_and_line
213
             * \copydoc doxygen_hide_do_update_ghost_arg
214 215 216 217
             * \param[in] matrix See formula above.
             * \param[in] v1 See formula above.
             * \param[in] v2 See formula above.
             * \param[in] v3 See formula above.
218 219
             */
            template<class MatrixT>
220
            std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
221 222 223 224
            MatMultTransposeAdd(const MatrixT& matrix,
                                const Vector& v1,
                                const Vector& v2,
                                Vector& v3,
225 226
                                const char* invoking_file, int invoking_line,
                                update_ghost do_update_ghost = update_ghost::yes);
227

228 229 230 231 232
            /*!
             * \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;
233
             * in the richer MoReFEM interface you should call in debug mode \a AssertMatrixRespectPattern()
234
             * to make sure the resulting matrix respects the pattern defined for the \a GlobalMatrix pair
235
             * of \a NumberingSubset.
236 237 238 239 240
             *
             * \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;
241
             * feel free to do so if you need it (for instance for MatMatMatMult support: Symbolic/Numeric doesn't
242 243
             * work for them and Petsc guys seemed unlikely to fix that in our exchanges).
             *
244
             * <b><tt>[internal]</tt></b> \todo #684 Investigate to use the argument fill, which provides an
245
             * estimation of the non zero of the resulting matrix. Currently PETSC_DEFAULT is used.
246 247

             * \endinternal
248
             */
249 250 251



252
            /*!
253
             * \brief Wrapper over MatMatMult, that performs m3 = m1 * m2.
254
             *
255 256
             * \copydetails doxygen_hide_matmatmult_warning
             *
257 258
             * 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.
259

260
             * One can also reuse a pattern for several matrices, but so far I do not need that in MoReFEM so
261
             * I equally use the default setting.
262 263 264
             * \param[in] m1 See formula above.
             * \param[in] m2 See formula above.
             * \param[in] m3 See formula above.
265
             *
266
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
267
             * \copydoc doxygen_hide_invoking_file_and_line
268
             */
269 270
            template
            <
271
            class MatrixT,
272 273
            class MatrixU,
            class MatrixV
274 275 276
            >
            std::enable_if_t
            <
277
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
278 279
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixV>::value,
280 281
                void
            >
282
            MatMatMult(const MatrixT& m1,
283 284
                       const MatrixU& m2,
                       MatrixV& m3,
285 286
                       const char* invoking_file, int invoking_line,
                       DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
287 288


289
            /*!
290
             * \brief Wrapper over MatMatMatMult, that performs m4 = m1 * m2 * m3.
291
             *
292 293 294 295
             * \param[in] m1 See formula above.
             * \param[in] m2 See formula above.
             * \param[in] m3 See formula above.
             * \param[in] m4 See formula above.
296
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
297
             * \copydoc doxygen_hide_invoking_file_and_line
298
             */
299 300
            template
            <
301
            class MatrixT,
302 303 304
            class MatrixU,
            class MatrixV,
            class MatrixW
305 306 307
            >
            std::enable_if_t
            <
308
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
309 310 311
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixV>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixW>::value,
312 313
                void
            >
314
            MatMatMatMult(const MatrixT& m1,
315 316 317
                          const MatrixU& m2,
                          const MatrixV& m3,
                          MatrixW& 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 335 336
            template
            <
            class MatrixT,
            class MatrixU
            >
337 338
            std::enable_if_t
            <
339 340
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value,
341 342 343
                void
            >
            MatCreateTranspose(const MatrixT& A,
344
                               MatrixU& transpose,
345
                               const char* invoking_file, int invoking_line);
346 347


348
            /*!
349
             * \class doxygen_hide_mat_transpose_mat_mult
350
             *
351 352
             * Formula is:
             * \verbatim
353
             m3 = m1^T * m2.
354 355 356 357 358 359 360
             * \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
361
             */
362 363


364
            /*!
365
             * \brief Performs Matrix-Matrix Multiplication m3 = m1^T * m2.
366 367
             *
             * \copydetails doxygen_hide_mat_transpose_mat_mult
368
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
369
             */
370 371
            template
            <
372 373 374
                class MatrixT,
                class MatrixU,
                class MatrixV
375 376 377
            >
            std::enable_if_t
            <
378
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
379
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixV>::value,
380 381
                void
            >
382
            MatTransposeMatMult(const MatrixT& m1,
383 384
                                const MatrixU& m2,
                                MatrixV& m3,
385 386
                                const char* invoking_file, int invoking_line,
                                DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
387 388 389



390 391 392
            /*!
             * \brief Performs Matrix-Matrix Multiplication C = A * B^T.
             *
393
             * \copydetails doxygen_hide_matmatmult_warning
394
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
395
             *
396 397 398
             * \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.
399
             * \copydoc doxygen_hide_invoking_file_and_line
400 401 402
             */
            template
            <
403
            class MatrixT,
404 405
            class MatrixU,
            class MatrixV
406 407 408
            >
            std::enable_if_t
            <
409
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
410 411
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixV>::value,
412 413 414
                void
            >
            MatMatTransposeMult(const MatrixT& matrix1,
415 416
                                const MatrixU& matrix2,
                                MatrixV& matrix3,
417 418
                                const char* invoking_file, int invoking_line,
                                DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
419 420


421 422 423
            /*!
             * \class doxygen_hide_mat_pt_a_p
             *
424
             * \warning Unfortunately a simple test with P = I leads to error
425
             * \verbatim
426 427 428
             [0;39m[0;49m[0]PETSC ERROR: Nonconforming object sizes
             [0]PETSC ERROR: Expected fill=-2 must be >= 1.0
             \endverbatim
429
             * The reason is that PETSC_DEFAULT may not be supported (I've asked Petsc developers); but even with
430 431 432 433 434 435 436 437 438
             * 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
439
             * \copydoc doxygen_hide_petsc_do_reuse_matrix_arg
440
             */
441

442 443 444 445 446 447 448 449 450
            /*!
             * \brief Performs the matrix product C = P^T * A * P
             *
             * \copydoc doxygen_hide_mat_pt_a_p
             *
             *
             */
            template
            <
451
            class MatrixT,
452 453
            class MatrixU,
            class MatrixV
454 455 456
            >
            std::enable_if_t
            <
457
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
458 459
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixV>::value,
460 461 462
                void
            >
            PtAP(const MatrixT& A,
463 464
                 const MatrixU& P,
                 MatrixV& out,
465 466
                 const char* invoking_file, int invoking_line,
                 DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
467 468 469



470 471 472 473 474 475
            /*!
             * \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.
             *
476 477 478 479
             * \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.
480
             * \copydetails doxygen_hide_invoking_file_and_line
481 482 483 484 485 486 487 488 489
             *
             */
            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);
490 491


492
            /*!
493
             * \brief Wrapper over MatLUFactor, that performs in-place LU factorization of matrix.
494 495 496 497
             *
             * \see http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatLUFactor.html#MatLUFactor
             * for more details.
             *
498 499 500 501 502
             * \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
503 504 505 506 507 508 509 510
             *
             */
            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);
511 512


513 514 515 516 517 518
            /*!
             * \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.
             *
519 520 521 522
             * \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
523 524 525 526
             *
             */
            template
            <
GILLES Sebastien's avatar
GILLES Sebastien committed
527 528 529
                class MatrixT,
                class MatrixU,
                class MatrixV
530 531 532
            >
            std::enable_if_t
            <
GILLES Sebastien's avatar
GILLES Sebastien committed
533 534 535 536
                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
537 538 539 540 541
            >
            MatMatSolve(const MatrixT& A,
                        const MatrixU& B,
                        MatrixV& X,
                        const char* invoking_file, int invoking_line);
542 543


544 545 546 547 548 549
            /*!
             * \brief Wrapper over MatCholeskyFactor, that performs in-place Cholesky factorization of a symmetric matrix.
             *
             * \see http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCholeskyFactor.html
             * for more details.
             *
550
             * \attention Not used yet, but might be useful in midterm developments. If finally too unwieldy for our
551 552
             * purposes it will be removed.
             *
553 554 555 556 557 558 559 560 561 562 563 564
             * \param[in] mat Matrix to factor.
             * \param[in] perm Result of the GetOrdering, row and column permutations.
             * \param[in] info Info of CholeskyFactor.
             * \copydetails doxygen_hide_invoking_file_and_line
             *
             */
            template<class MatrixT>
            std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
            CholeskyFactor(MatrixT& mat,
                           IS perm,
                           const MatFactorInfo *info,
                           const char* invoking_file, int invoking_line);
565 566


567 568


569
        } //namespace Petsc
570 571


572
    } // namespace Wrappers
573 574 575 576



    /*!
577 578
     * \brief Convenient alias to avoid repeating the namespaces in each call.
     *
579
     * It should be used anyway only as template arguments of functions within this namespaces, such as
580 581 582
     * Wrappers::Petsc::AXPY.
     */
    using NonZeroPattern = Wrappers::Petsc::NonZeroPattern;
583 584


585
} // namespace MoReFEM
586 587


588 589 590
/// @} // addtogroup ThirdPartyGroup


591
# include "ThirdParty/Wrappers/Petsc/Matrix/MatrixOperations.hxx"
592 593


594
#endif // MOREFEM_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_MATRIX_x_MATRIX_OPERATIONS_HPP_