MatrixOperations.hpp 22.7 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 128 129 130 131 132 133 134 135
            /*!
             * \brief Wrapper over MatShift, that performs M = M + a I, where a is a PetscScalar and I is the identity matrix.
             *
             * \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>
            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);
            
            
136 137 138
            /*!
             * \brief Wrapper over MatMult, that performs v2 = matrix * v1.
             *
139 140 141 142
             * \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.
143
             * \copydoc doxygen_hide_do_update_ghost_arg
144 145
             */
            template<class MatrixT>
146
            std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
147 148 149
            MatMult(const MatrixT& matrix,
                    const Vector& v1,
                    Vector& v2,
150 151
                    const char* invoking_file, int invoking_line,
                    update_ghost do_update_ghost = update_ghost::yes);
152
            
153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181
            /*!
             * \brief Wrapper over MatTranspose, that performs in-place or out-of-place transpose of a matrix.
             *
             * \copydetails doxygen_hide_matranpose_warning
             *
             * 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.
             
             * 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.
             *
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
             *
             */
            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);
            
182 183 184 185
            
            /*!
             * \brief Wrapper over MatMultAdd, that performs v3 = v2 + matrix * v1.
             *
186 187 188 189 190
             * \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.
191
             * \copydoc doxygen_hide_do_update_ghost_arg
192 193
             */
            template<class MatrixT>
194
            std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
195 196 197 198
            MatMultAdd(const MatrixT& matrix,
                       const Vector& v1,
                       const Vector& v2,
                       Vector& v3,
199 200
                       const char* invoking_file, int invoking_line,
                       update_ghost do_update_ghost = update_ghost::yes);
201 202 203 204 205 206
            
            
            
            /*!
             * \brief Wrapper over MatMultTranspose, that performs v2 = transpose(matrix) * v1.
             *
207
             * \copydoc doxygen_hide_invoking_file_and_line
208
             * \copydoc doxygen_hide_do_update_ghost_arg
209 210 211
             * \param[in] matrix See formula above.
             * \param[in] v1 See formula above.
             * \param[in] v2 See formula above.
212 213
             */
            template<class MatrixT>
214
            std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
215 216 217
            MatMultTranspose(const MatrixT& matrix,
                             const Vector& v1,
                             Vector& v2,
218 219
                             const char* invoking_file, int invoking_line,
                             update_ghost do_update_ghost = update_ghost::yes);
220 221 222 223 224
            
            
            /*!
             * \brief Wrapper over MatMultTransposeAdd, that performs v3 = v2 + transpose(matrix) * v1.
             *
225
             * \copydoc doxygen_hide_invoking_file_and_line
226
             * \copydoc doxygen_hide_do_update_ghost_arg
227 228 229 230
             * \param[in] matrix See formula above.
             * \param[in] v1 See formula above.
             * \param[in] v2 See formula above.
             * \param[in] v3 See formula above.
231 232
             */
            template<class MatrixT>
233
            std::enable_if_t<std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value, void>
234 235 236 237
            MatMultTransposeAdd(const MatrixT& matrix,
                                const Vector& v1,
                                const Vector& v2,
                                Vector& v3,
238 239
                                const char* invoking_file, int invoking_line,
                                update_ghost do_update_ghost = update_ghost::yes);
240
            
241 242 243 244 245 246
            /*!
             * \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()
247
             * to make sure the resulting matrix respects the pattern defined for the \a GlobalMatrix pair
248
             * of \a NumberingSubset.
249 250 251 252 253
             *
             * \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;
254
             * feel free to do so if you need it (for instance for MatMatMatMult support: Symbolic/Numeric doesn't
255 256
             * work for them and Petsc guys seemed unlikely to fix that in our exchanges).
             *
257
             * \internal <b><tt>[internal]</tt></b> \todo #684 Investigate to use the argument fill, which provides an
258
             * estimation of the non zero of the resulting matrix. Currently PETSC_DEFAULT is used.
259 260 261
             */
            
            
262 263
            
            /*!
264
             * \brief Wrapper over MatMatMult, that performs m3 = m1 * m2.
265
             *
266 267
             * \copydetails doxygen_hide_matmatmult_warning
             *
268 269
             * 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.
270
             
271 272
             * 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.
273 274 275
             * \param[in] m1 See formula above.
             * \param[in] m2 See formula above.
             * \param[in] m3 See formula above.
276
             *
277
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
278
             * \copydoc doxygen_hide_invoking_file_and_line
279
             */
280 281
            template
            <
282 283
            class MatrixT,
            class MatrixU
284 285 286
            >
            std::enable_if_t
            <
287 288
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value,
289 290
                void
            >
291 292 293
            MatMatMult(const MatrixT& m1,
                       const MatrixT& m2,
                       MatrixU& m3,
294 295
                       const char* invoking_file, int invoking_line,
                       DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
296
            
297
                      
298
            /*!
299
             * \brief Wrapper over MatMatMatMult, that performs m4 = m1 * m2 * m3.
300
             *
301 302 303 304
             * \param[in] m1 See formula above.
             * \param[in] m2 See formula above.
             * \param[in] m3 See formula above.
             * \param[in] m4 See formula above.
305
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
306
             * \copydoc doxygen_hide_invoking_file_and_line
307
             */
308 309
            template
            <
310 311
            class MatrixT,
            class MatrixU
312 313 314
            >
            std::enable_if_t
            <
315 316
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value,
317 318
                void
            >
319 320 321 322
            MatMatMatMult(const MatrixT& m1,
                          const MatrixT& m2,
                          const MatrixT& m3,
                          MatrixU& m4,
323 324
                          const char* invoking_file, int invoking_line,
                          DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
325
            
326 327 328
            /*!
             * \brief Creates a new matrix object that behaves like A'.
             *
329
             * The transpose A' is NOT actually formed! Rather the new matrix object performs the matrix-vector product
330 331 332
             * by using the MatMultTranspose() on the original matrix.
             *
             * \param[in] A matrix to transpose.
333
             * \param[out] transpose The matrix that figuratively represents A'. This matrix must not have been
334
             * allocated!
335
             * \copydoc doxygen_hide_invoking_file_and_line
336 337 338 339
             */
            template<class MatrixT>
            std::enable_if_t
            <
340
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value,
341 342 343 344 345 346
                void
            >
            MatCreateTranspose(const MatrixT& A,
                               MatrixT& transpose,
                               const char* invoking_file, int invoking_line);
            
347
            
348
            /*!
349
             * \class doxygen_hide_mat_transpose_mat_mult
350
             *
351 352 353 354 355 356 357 358 359 360
             * 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
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
             * \param[in] m1 See formula above.
             * \param[in] m2 See formula above.
             * \param[in] m3 See formula above.
372
             */
373 374
            template
            <
375 376
            class MatrixT,
            class MatrixU
377 378 379
            >
            std::enable_if_t
            <
380 381
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value,
382 383
                void
            >
384 385 386
            MatTransposeMatMult(const MatrixT& m1,
                                const MatrixT& m2,
                                MatrixU& m3,
387 388
                                const char* invoking_file, int invoking_line,
                                DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
389 390
            
            
391
            
392 393 394
            /*!
             * \brief Performs Matrix-Matrix Multiplication C = A * B^T.
             *
395
             * \copydetails doxygen_hide_matmatmult_warning
396
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
397
             *
398 399 400
             * \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.
401
             * \copydoc doxygen_hide_invoking_file_and_line
402 403 404
             */
            template
            <
405 406
            class MatrixT,
            class MatrixU
407 408 409
            >
            std::enable_if_t
            <
410 411
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value,
412 413 414 415 416
                void
            >
            MatMatTransposeMult(const MatrixT& matrix1,
                                const MatrixT& matrix2,
                                MatrixU& 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 426 427 428
             * \verbatim
             * [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 452
            class MatrixT,
            class MatrixU
453 454 455
            >
            std::enable_if_t
            <
456 457
                std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixT>::value
                && std::is_base_of<Internal::Wrappers::Petsc::BaseMatrix, MatrixU>::value,
458 459 460 461 462
                void
            >
            PtAP(const MatrixT& A,
                 const MatrixT& P,
                 MatrixU& out,
463 464
                 const char* invoking_file, int invoking_line,
                 DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
465 466
            
            
467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513
            
            /*!
             * \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.
             *
             * \copydoc doxygen_hide_mat_get_ordering
             *
             */
            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.
             *
             * \copydoc doxygen_hide_mat_lu_factor
             *
             */
            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.
             *
             * \copydoc doxygen_hide_mat_mat_solve
             *
             */
            template
            <
GILLES Sebastien's avatar
GILLES Sebastien committed
514 515 516
                class MatrixT,
                class MatrixU,
                class MatrixV
517 518 519
            >
            std::enable_if_t
            <
GILLES Sebastien's avatar
GILLES Sebastien committed
520 521 522 523
                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
524 525 526 527 528 529 530 531 532 533
            >
            MatMatSolve(const MatrixT& A,
                        const MatrixU& B,
                        MatrixV& X,
                        const char* invoking_file, int invoking_line);


            
            
        } //namespace Petsc
534 535
        
        
536
    } // namespace Wrappers
537 538
    
    
539 540 541 542 543 544 545 546 547 548
    
    /*! 
     * \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;
    
    
549 550 551
} // namespace HappyHeart


552 553 554
/// @} // addtogroup ThirdPartyGroup


555
# include "ThirdParty/Wrappers/Petsc/Matrix/MatrixOperations.hxx"
556 557 558


#endif // HAPPY_HEART_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_MATRIX_x_MATRIX_OPERATIONS_HPP_