MatrixOperations.hpp 17.4 KB
Newer Older
1 2
//! \file 
//
3 4 5 6 7 8 9 10 11 12 13 14 15 16
//
//  MatrixTOperations.hpp
//  HappyHeart
//
//  Created by Sebastien Gilles on 30/10/15.
//  Copyright © 2015 Inria. All rights reserved.
//

#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>

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

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

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


namespace HappyHeart
{
    
    
32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
    // ============================
    //! \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
    // ============================


49 50 51 52 53 54
    namespace Wrappers
    {
        
        
        namespace Petsc
        {
55 56
            
            
57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76
            /*!
             * \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 };
77
                   
78 79 80
            /*!
             * \brief Wrapper over MatAXPY, that performs Y = a * X + Y.
             *
81
             *
82
             * \tparam NonZeroPatternT This value indicates the level of similarity between X and Y non zero patterns.
83
             * \warning Beware: if X and Y aren't following the same pattern, 'Same' won't yield what
84
             * you expect! In the case you're considering adding a part computed by a transfert matrix, you should
85
             * use 'Subset'. In HappyHeart, I advise you to call in debug mode AssertSameNumberingSubset()
86
             * (defined in Core) just before or after the function call; if the assert is raised it means you can't
87
             * use 'Same'. \todo This could be enforced by a proper overload (but it mix Utilities and Core
88
             * libraries...), but for the time being the separate assert will do.
89 90 91 92 93
             * \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.
94
             */
95 96
            template
            <
97
                NonZeroPattern NonZeroPatternT,
98 99 100 101 102 103 104 105 106
                class MatrixT,
                class MatrixU
            >
            std::enable_if_t
            <
                std::is_base_of<Private::BaseMatrix, MatrixT>::value
             && std::is_base_of<Private::BaseMatrix, MatrixU>::value,
                void
            >
107
            AXPY(PetscScalar a,
108
                 const MatrixT& X,
109
                 MatrixU& Y,
110
                 const char* invoking_file, int invoking_line);
111 112 113 114 115
            
            
            /*!
             * \brief Wrapper over MatMult, that performs v2 = matrix * v1.
             *
116 117 118 119
             * \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.
120
             * \copydoc doxygen_hide_do_update_ghost_arg
121 122 123
             */
            template<class MatrixT>
            std::enable_if_t<std::is_base_of<Private::BaseMatrix, MatrixT>::value, void>
124 125 126
            MatMult(const MatrixT& matrix,
                    const Vector& v1,
                    Vector& v2,
127 128
                    const char* invoking_file, int invoking_line,
                    update_ghost do_update_ghost = update_ghost::yes);
129 130 131 132 133
            
            
            /*!
             * \brief Wrapper over MatMultAdd, that performs v3 = v2 + matrix * v1.
             *
134 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.
             * \param[in] v3 See formula above.
139
             * \copydoc doxygen_hide_do_update_ghost_arg
140 141 142
             */
            template<class MatrixT>
            std::enable_if_t<std::is_base_of<Private::BaseMatrix, MatrixT>::value, void>
143 144 145 146
            MatMultAdd(const MatrixT& matrix,
                       const Vector& v1,
                       const Vector& v2,
                       Vector& v3,
147 148
                       const char* invoking_file, int invoking_line,
                       update_ghost do_update_ghost = update_ghost::yes);
149 150 151 152 153 154
            
            
            
            /*!
             * \brief Wrapper over MatMultTranspose, that performs v2 = transpose(matrix) * v1.
             *
155
             * \copydoc doxygen_hide_invoking_file_and_line
156
             * \copydoc doxygen_hide_do_update_ghost_arg
157 158 159
             * \param[in] matrix See formula above.
             * \param[in] v1 See formula above.
             * \param[in] v2 See formula above.
160 161 162
             */
            template<class MatrixT>
            std::enable_if_t<std::is_base_of<Private::BaseMatrix, MatrixT>::value, void>
163 164 165
            MatMultTranspose(const MatrixT& matrix,
                             const Vector& v1,
                             Vector& v2,
166 167
                             const char* invoking_file, int invoking_line,
                             update_ghost do_update_ghost = update_ghost::yes);
168 169 170 171 172
            
            
            /*!
             * \brief Wrapper over MatMultTransposeAdd, that performs v3 = v2 + transpose(matrix) * v1.
             *
173
             * \copydoc doxygen_hide_invoking_file_and_line
174
             * \copydoc doxygen_hide_do_update_ghost_arg
175 176 177 178
             * \param[in] matrix See formula above.
             * \param[in] v1 See formula above.
             * \param[in] v2 See formula above.
             * \param[in] v3 See formula above.
179 180 181
             */
            template<class MatrixT>
            std::enable_if_t<std::is_base_of<Private::BaseMatrix, MatrixT>::value, void>
182 183 184 185
            MatMultTransposeAdd(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 194 195 196
            /*!
             * \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()
             * to make sure the resulting matrix respects the pattern defined for the \a GlobalMatrix pair 
             * of \a NumberingSubset.
197 198 199 200 201 202 203 204 205 206
             *
             * \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;
             * feel free to do so if you need it (for instance for MatMatMatMult support: Symbolic/Numeric doesn't 
             * work for them and Petsc guys seemed unlikely to fix that in our exchanges).
             *
             * \internal <b><tt>[internal]</tt></b> \todo #684 Investigate to use the argument fill, which provides an 
             * estimation of the non zero of the resulting matrix. Currently PETSC_DEFAULT is used.
207 208 209
             */
            
            
210 211
            
            /*!
212
             * \brief Wrapper over MatMatMult, that performs m3 = m1 * m2.
213
             *
214 215
             * \copydetails doxygen_hide_matmatmult_warning
             *
216 217
             * 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.
218
            
219 220
             * 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.
221 222 223
             * \param[in] m1 See formula above.
             * \param[in] m2 See formula above.
             * \param[in] m3 See formula above.
224
             *
225
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
226
             * \copydoc doxygen_hide_invoking_file_and_line
227
             */
228 229 230 231 232 233 234 235 236 237 238
            template
            <
                class MatrixT,
                class MatrixU
            >
            std::enable_if_t
            <
                std::is_base_of<Private::BaseMatrix, MatrixT>::value
                && std::is_base_of<Private::BaseMatrix, MatrixU>::value,
                void
            >
239 240 241
            MatMatMult(const MatrixT& m1,
                       const MatrixT& m2,
                       MatrixU& m3,
242 243
                       const char* invoking_file, int invoking_line,
                       DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
244
            
245
                      
246
            /*!
247
             * \brief Wrapper over MatMatMatMult, that performs m4 = m1 * m2 * m3.
248
             *
249 250 251 252
             * \param[in] m1 See formula above.
             * \param[in] m2 See formula above.
             * \param[in] m3 See formula above.
             * \param[in] m4 See formula above.
253
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
254
             * \copydoc doxygen_hide_invoking_file_and_line
255
             */
256 257 258 259 260 261 262 263 264 265 266
            template
            <
                class MatrixT,
                class MatrixU
            >
            std::enable_if_t
            <
                std::is_base_of<Private::BaseMatrix, MatrixT>::value
                && std::is_base_of<Private::BaseMatrix, MatrixU>::value,
                void
            >
267 268 269 270
            MatMatMatMult(const MatrixT& m1,
                          const MatrixT& m2,
                          const MatrixT& m3,
                          MatrixU& m4,
271 272
                          const char* invoking_file, int invoking_line,
                          DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
273
            
274
                  
275 276 277 278 279 280 281 282 283
            /*!
             * \brief Creates a new matrix object that behaves like A'.
             *
             * 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.
             *
             * \param[in] A matrix to transpose.
             * \param[out] transpose The matrix that figuratively represents A'. This matrix must not have been 
             * allocated!
284
             * \copydoc doxygen_hide_invoking_file_and_line
285 286 287 288 289 290 291 292 293 294 295
             */
            template<class MatrixT>
            std::enable_if_t
            <
                std::is_base_of<Private::BaseMatrix, MatrixT>::value,
                void
            >
            MatCreateTranspose(const MatrixT& A,
                               MatrixT& transpose,
                               const char* invoking_file, int invoking_line);
            
296
            
297
            /*!
298
             * \class doxygen_hide_mat_transpose_mat_mult
299
             *
300 301 302 303 304 305 306 307 308 309
             * 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
310
             */
311 312 313
            
            
            /*!
314
             * \brief Performs Matrix-Matrix Multiplication m3 = m1^T * m2.
315 316
             *
             * \copydetails doxygen_hide_mat_transpose_mat_mult
317
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
318 319 320
             * \param[in] m1 See formula above.
             * \param[in] m2 See formula above.
             * \param[in] m3 See formula above.
321
             */
322 323 324 325 326 327 328 329 330 331 332
            template
            <
                class MatrixT,
                class MatrixU
            >
            std::enable_if_t
            <
                std::is_base_of<Private::BaseMatrix, MatrixT>::value
                && std::is_base_of<Private::BaseMatrix, MatrixU>::value,
                void
            >
333 334 335
            MatTransposeMatMult(const MatrixT& m1,
                                const MatrixT& m2,
                                MatrixU& m3,
336 337
                                const char* invoking_file, int invoking_line,
                                DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
338 339
            
            
340
                     
341 342 343
            /*!
             * \brief Performs Matrix-Matrix Multiplication C = A * B^T.
             *
344
             * \copydetails doxygen_hide_matmatmult_warning
345
             * \copydetails doxygen_hide_petsc_do_reuse_matrix_arg
346
             *
347 348 349
             * \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.
350
             * \copydoc doxygen_hide_invoking_file_and_line
351 352 353 354 355 356 357 358 359 360 361 362 363 364 365
             */
            template
            <
                class MatrixT,
                class MatrixU
            >
            std::enable_if_t
            <
                std::is_base_of<Private::BaseMatrix, MatrixT>::value
                && std::is_base_of<Private::BaseMatrix, MatrixU>::value,
                void
            >
            MatMatTransposeMult(const MatrixT& matrix1,
                                const MatrixT& matrix2,
                                MatrixU& matrix3,
366 367
                                const char* invoking_file, int invoking_line,
                                DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
368 369
            
            
370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387
            /*!
             * \class doxygen_hide_mat_pt_a_p
             *
             * \warning Unfortunately a simple test with P = I leads to error  
             * \verbatim
             * [0]PETSC ERROR: Nonconforming object sizes
             * [0]PETSC ERROR: Expected fill=-2 must be >= 1.0
             * \endverbatim
             * The reason is that PETSC_DEFAULT may not be supported (I've asked Petsc developers); but even with 
             * 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
388
             * \copydoc doxygen_hide_petsc_do_reuse_matrix_arg
389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411
             */
            
            /*!
             * \brief Performs the matrix product C = P^T * A * P
             *
             * \copydoc doxygen_hide_mat_pt_a_p
             *
             *
             */
            template
            <
                class MatrixT,
                class MatrixU
            >
            std::enable_if_t
            <
                std::is_base_of<Private::BaseMatrix, MatrixT>::value
                && std::is_base_of<Private::BaseMatrix, MatrixU>::value,
                void
            >
            PtAP(const MatrixT& A,
                 const MatrixT& P,
                 MatrixU& out,
412 413
                 const char* invoking_file, int invoking_line,
                 DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
414 415
            
            
416 417 418 419 420 421
        } //namespace Petsc
        
        
    } //namespace Wrappers
    
    
422 423 424 425 426 427 428 429 430 431
    
    /*! 
     * \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;
    
    
432 433 434
} // namespace HappyHeart


435
# include "ThirdParty/Wrappers/Petsc/Matrix/MatrixOperations.hxx"
436 437 438


#endif // HAPPY_HEART_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_MATRIX_x_MATRIX_OPERATIONS_HPP_