Attention une mise à jour du serveur va être effectuée le lundi 17 mai entre 13h et 13h30. Cette mise à jour va générer une interruption du service de quelques minutes.

MatrixOperations.hpp 20.1 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
            /*!
             * \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.
68 69 70
             * If is set to in_place, the input matrix will be replaced by the result of the function.
             * This is not applicable to all functions but can be useful in some, for instance with
             * MatTranspose which replaces the input matrix by its transposed matrix.
71
             */
72

73 74 75 76 77
            /*!
             * \class doxygen_hide_petsc_do_reuse_matrix_arg
             *
             * \param[in] do_reuse_matrix \copydetails doxygen_hide_petsc_do_reuse_matrix
             */
78 79


80 81 82
            /*!
             * \brief \copydoc doxygen_hide_petsc_do_reuse_matrix
             */
83
            enum class DoReuseMatrix { yes, no, in_place };
84

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


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


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

149 150 151
            /*!
             * \brief Wrapper over MatTranspose, that performs in-place or out-of-place transpose of a matrix.
             *
152 153
             * \param[in] matrix1 Matrix to transpose.
             * \param[in] matrix2 Matrix to contain the result of the transpose.
154 155
             * \copydoc doxygen_hide_invoking_file_and_line
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
156 157 158
             */
            template
            <
159 160
                class MatrixT,
                class MatrixU
161 162 163
            >
            std::enable_if_t
            <
164 165 166
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value,
                void
167 168 169
            >
            MatTranspose(MatrixT& matrix1,
                         MatrixU& matrix2,
170 171
                         const char* invoking_file, int invoking_line,
                         DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
172 173


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



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


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

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

             * \endinternal
252
             */
253 254 255



256
            /*!
257
             * \brief Wrapper over MatMatMult, that performs m3 = m1 * m2.
258
             *
259 260
             * \copydetails doxygen_hide_matmatmult_warning
             *
261 262
             * 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.
263

264
             * One can also reuse a pattern for several matrices, but so far I do not need that in MoReFEM so
265
             * I equally use the default setting.
266 267 268
             * \param[in] m1 See formula above.
             * \param[in] m2 See formula above.
             * \param[in] m3 See formula above.
269
             *
270
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
271
             * \copydoc doxygen_hide_invoking_file_and_line
272
             */
273 274
            template
            <
275
            class MatrixT,
276 277
            class MatrixU,
            class MatrixV
278 279 280
            >
            std::enable_if_t
            <
281
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
282 283
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixV>::value,
284 285
                void
            >
286
            MatMatMult(const MatrixT& m1,
287 288
                       const MatrixU& m2,
                       MatrixV& 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
            class MatrixT,
306 307 308
            class MatrixU,
            class MatrixV,
            class MatrixW
309 310 311
            >
            std::enable_if_t
            <
312
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
313 314 315
                && 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,
316 317
                void
            >
318
            MatMatMatMult(const MatrixT& m1,
319 320 321
                          const MatrixU& m2,
                          const MatrixV& m3,
                          MatrixW& m4,
322 323
                          const char* invoking_file, int invoking_line,
                          DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
324

325 326 327 328 329
            /*!
             * \brief Creates a new matrix object that behaves like A'.
             *
             *
             * \param[in] A matrix to transpose.
330
             * \param[out] transpose The matrix that figuratively represents A'. This matrix must not have been
331
             * allocated!
332
             * \copydoc doxygen_hide_invoking_file_and_line
333 334 335 336
             *
             * \attention The transpose A' is NOT actually formed! Rather the new matrix object performs the matrix-vector product
             * by using the MatMultTranspose() on the original matrix.
             *
337
             */
338 339
            template
            <
340 341
                class MatrixT,
                class MatrixU
342
            >
343 344
            std::enable_if_t
            <
345 346
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value,
347 348 349
                void
            >
            MatCreateTranspose(const MatrixT& A,
350
                               MatrixU& transpose,
351
                               const char* invoking_file, int invoking_line);
352 353


354
            /*!
355
             * \class doxygen_hide_mat_transpose_mat_mult
356
             *
357 358
             * Formula is:
             * \verbatim
359
             m3 = m1^T * m2.
360 361 362 363 364 365 366
             * \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
367
             */
368 369


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



396 397 398
            /*!
             * \brief Performs Matrix-Matrix Multiplication C = A * B^T.
             *
399
             * \copydetails doxygen_hide_matmatmult_warning
400
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
401
             *
402 403 404
             * \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.
405
             * \copydoc doxygen_hide_invoking_file_and_line
406 407 408
             *
             * \attention It is not advised to use this wrapper (at least currently: as of PETSc v1.13.2, there is no support for this operation
             * for MPIAIJ and MPIAIJ matrices - so the code you would write with it will not run in parallel.
409 410 411
             */
            template
            <
412
            class MatrixT,
413 414
            class MatrixU,
            class MatrixV
415 416 417
            >
            std::enable_if_t
            <
418
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
419 420
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixV>::value,
421 422 423
                void
            >
            MatMatTransposeMult(const MatrixT& matrix1,
424 425
                                const MatrixU& matrix2,
                                MatrixV& matrix3,
426 427
                                const char* invoking_file, int invoking_line,
                                DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
428 429


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

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


478
        } // namespace Petsc
479 480


481
    } // namespace Wrappers
482 483 484 485



    /*!
486 487
     * \brief Convenient alias to avoid repeating the namespaces in each call.
     *
488
     * It should be used anyway only as template arguments of functions within this namespaces, such as
489 490 491
     * Wrappers::Petsc::AXPY.
     */
    using NonZeroPattern = Wrappers::Petsc::NonZeroPattern;
492 493


494
} // namespace MoReFEM
495 496


497 498 499
/// @} // addtogroup ThirdPartyGroup


500
# include "ThirdParty/Wrappers/Petsc/Matrix/MatrixOperations.hxx"
501 502


503
#endif // MOREFEM_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_MATRIX_x_MATRIX_OPERATIONS_HPP_