Matrix.hpp 24 KB
Newer Older
1
2
//! \file 
//
3
4
//
//  Petsc.h
5
//  HappyHeart
6
7
8
9
10
//
//  Created by Sébastien Gilles on 20/09/13.
//  Copyright (c) 2013 Inria. All rights reserved.
//

11
12
#ifndef HAPPY_HEART_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_MATRIX_x_MATRIX_HPP_
# define HAPPY_HEART_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_MATRIX_x_MATRIX_HPP_
13
14
15
16
17


# include <iostream>
# include <vector>

18
19
# include "ThirdParty/Wrappers/Mpi/Mpi.hpp"

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

23
24
# include "ThirdParty/Wrappers/Petsc/Matrix/Private/BaseMatrix.hpp"

25
# include "Utilities/SparseMatrix/CSRPattern.hpp"
26

27
28
29
30


namespace HappyHeart
{
31
32
    
    
33
    // ============================
34
35
    //! \cond IGNORE_BLOCK_IN_DOXYGEN
    // Forward declarations.    
36
37
38
    // ============================
    
    
39
40
    namespace Wrappers
    {
41
42
        
        
43
44
45
46
        // Forward declaration.
        class Mpi;
        
        
47
48
        namespace Petsc
        {
49
50
51
            
            
            class Vector;
52
            class MatrixPattern;
53
54
            
            
55
56
57
58
59
60
61
62
63
        } // namespace Petsc
        
        
    } // namespace Wrappers
    
    
    
    // ============================
    // End of forward declarations.
64
    //! \endcond IGNORE_BLOCK_IN_DOXYGEN
65
66
67
68
69
70
71
72
73
74
75
76
    // ============================

    
    
    namespace Wrappers
    {
        
        
        namespace Petsc
        {
            
     
77
            
78
            /*!
79
             * \brief A wrapper class over Petsc Mat objects.
80
             *
81
82
             * Most of the Petsc functions used over Petsc matrices have been included as methods in this class, which
             * also acts as RAII over Petsc Mat object.
83
             *
84
             * \internal <b><tt>[internal]</tt></b> This class is way more trickier to implement that it might seems because of Petsc's
85
86
             * internal structure: Mat objects are in fact pointers over an internal Petsc type. So copy and
             * destruction operations must be made with great care! That's why several choices have been made:
87
             *
88
89
90
91
92
             * - No implicit conversion to the internal Mat. It seems alluring at first sight to allow it but can
             * lead very fastly to runtime problems (covered by asserts in debug mode).
             * - No copy semantic for this class.
             * - The most usual way to proceed is to construct with the default constructor and then init it either'
             * by a 'Init***()' method or by using 'DuplicateLayout', 'Copy' or 'CompleteCopy'methods.
93
94
             *
             */
95
            class Matrix : private Private::BaseMatrix
96
            {
97
98
99
100
101
            private:
                
                //! Alias to parent.
                using parent = Private::BaseMatrix;
                
102
            public:
103
104
105
106
107
108
109
                
                //! Constructor.
                explicit Matrix();
                
                //! Destructor.
                ~Matrix();
                
110
                //! Disable move.
111
                Matrix(Matrix&&) = delete;
112
                
113
114
                //! Copy constructor: copy both layout and data.
                Matrix(const Matrix&);
115
116
117
118
                
                //! Disable affectation operator.
                Matrix& operator=(const Matrix&) = delete;
                
119
120
121
                //! Disable affectation operator.
                Matrix& operator=(Matrix&&) = delete;
                
122
123
124
125
126
127
128
                
                /*!
                 * \brief Create a sequential sparse matrix.
                 *
                 *
                 * \param[in] Nrow Number of rows.
                 * \param[in] Ncolumn Number of columns.
129
                 * \param[in] matrix_pattern Pattern of the matrix (number of elements expected on each row_).
130
                 * \copydetails doxygen_hide_mpi_param
131
132
                 * \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__.
133
134
135
136
                 *
                 * The matrix hence created will keep the non-zero structure: even if a value becomes 0 it will
                 * go on being held as a non-zero emplacement.
                 */
137
                void InitSequentialMatrix(unsigned int Nrow, unsigned int Ncolumn,
138
                                          const MatrixPattern& matrix_pattern,
139
                                          const Mpi& mpi, const char* invoking_file, int invoking_line);
140
141
142
143
144
145
146
147
148
                
                
                /*!
                 * \brief Create a parallel sparse matrix.
                 *
                 * \param[in] Nlocal_row Number of rows on the local processor.
                 * \param[in] Nlocal_column Number of columns on the local processor.
                 * \param[in] Nglobal_row Number of rows on the global processor.
                 * \param[in] Nglobal_column Number of columns on the global processor.
149
                 * \param[in] matrix_pattern Pattern of the matrix (number of elements expected on each row_).
150
                 * \copydetails doxygen_hide_mpi_param
151
152
                 * \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__.
153
154
                 *
                 */
155
156
                void InitParallelMatrix(unsigned int Nlocal_row, unsigned int Nlocal_column,
                                        unsigned int Nglobal_row, unsigned int Nglobal_column,
157
                                        const MatrixPattern& matrix_pattern,
158
                                        const Mpi& mpi, const char* invoking_file, int invoking_line);
159
                
160
                                
161
162
163
                /*!
                 * \brief Handle over the internal Mat object.
                 *
164
165
                 * \return Internal Mat object, which is indeed a pointer in Petsc.
                 *
166
167
168
                 * Ideally it shouldn't be used at all except in the implementation of the Petsc Wrapper: a wrapper
                 * method should be implemented over the function that might need access to the Vec internal object.
                 */
169
                Mat Internal() const noexcept;
170
                
171
172
173
174
175
176
177
178
179
180
                
                /*!
                 * \brief Init a Matrix from a Petsc Mat object.
                 *
                 * Ideally it shouldn't be used at all except in the implementation of few of the associated 
                 * function such as MatMatMult.
                 */
                void SetFromPetscMat(Mat petsc_matrix);
                
                
181
182
183
184
185
186
187
188
                /*!
                 * \brief Duplicate the layout of the matrix.
                 *
                 * \param[in] original Original matrix, which layout is to be duplicated.
                 * \param[in] option Fine tune the behaviour of the duplication. Three values are possible:
                 * - MAT_DO_NOT_COPY_VALUES
                 * - MAT_COPY_VALUES
                 * - MAT_SHARE_NONZERO_PATTERN (share the nonzero patterns with the previous matrix but do not copy them.)
189
190
                 * \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__.
191
                 *
192
193
194
195
196
197
                 * To make the third option work, the \a original matrix must have been assembled (i.e. where are the
                 * non-zero values must already be determined accurately, not just how many of them there are on
                 * each line - which is the information we get when a matrix is created). So technically
                 * it means you can't just duplicate the layout of a matrix you've just initialized but in which
                 * you haven't put any value yet.
                 *
198
                 */
199
                void DuplicateLayout(const Matrix& original, MatDuplicateOption option, const char* invoking_file, int invoking_line);
200
201
202
203
204
205
                
                
                
                /*!
                 * \brief A wrapper over MatCopy, which assumed target already gets the right layout.
                 *
206
207
                 * \param[in] source Original matrix which content is copied. It is assumed the object
                 * that called the method already gets the same layout as \a source.
208
209
210
211
212
                 * \param[in] structure SAME_NONZERO_PATTERN or DIFFERENT_NONZERO_PATTERN. SAME_NONZERO_PATTERN 
                 * should apply for almost all cases in HappyHeart; however if you're using raw Petsc library
                 * you'll see there are more prerequisite to this choice (the pattern must be fully built). 
                 * Within the frame of this wrapper there shouldn't be any problem provided one of the InitXXX() 
                 * method has been called.
213
214
                 * \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__.
215
                 */
216
                void Copy(const Matrix& source, const char* invoking_file, int invoking_line,
217
218
219
220
221
                          const MatStructure& structure = SAME_NONZERO_PATTERN);
                
                
                /*!
                 * \brief A complete copy: layout is copied first and then the values.
222
                 *
223
224
225
                 * \param[in] source Original matrix which layout AND content is copied.
                 * Beware: this method requires that \a source has already been deeply initialized (SetValues()
                 * and Assembly() should already have been called). See DuplicateLayout() for a lenghtier explanation.
226
227
                 * \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__.
228
                 *
229
                 */
230
                void CompleteCopy(const Matrix& source, const char* invoking_file, int invoking_line);
231
                
232
                /*!
233
                 * \brief Get the number of elements in the local matrix.
234
235
236
                 *
                 * \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__.
237
238
                 *
                 * \return First index is number of rows, second one the number of columns.
239
                 */
240
                std::pair<PetscInt, PetscInt> GetProcessorWiseSize(const char* invoking_file, int invoking_line) const;
241
                
242
243
244
245
246
                /*!
                 * \brief Get the number of elements in the global matrix (first is number of rows, second number of columns).
                 *
                 * \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__.
247
248
                 *
                 * \return First index is number of rows, second one the number of columns.
249
                 */
250
                std::pair<PetscInt, PetscInt> GetProgramWiseSize(const char* invoking_file, int invoking_line) const;
251
252
                
                
253
254
255
256
257
258
                /*!
                 * \brief Set all the entries to zero.
                 *
                 * \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__.
                 */
259
                void ZeroEntries(const char* invoking_file, int invoking_line);
260
261
262
263
264
265
266
                
                
                /*!
                 * \brief Zeros all entries (except possibly the main diagonal) of a set of rows of a matrix.
                 *
                 * \param[in] row_indexes Indexes of the row considered (C - numbering).
                 * \param[in] diagonal_value Value fior the diagonal term. Put 0. if you want it zeroed too.
267
268
269
                 * \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__.
                 
270
271
272
273
                 */
                
                void ZeroRows(const std::vector<PetscInt>& row_indexes,
                              PetscScalar diagonal_value,
274
                              const char* invoking_file, int invoking_line);
275
276
277
278
                
                
                
                /*!
279
                 * \brief Petsc Assembling of the matrix.
280
                 *
281
                 * \param[in] invoking_file File that invoked the function or class; usually __FILE__.
282
                 * \param[in] invoking_line File that invoked the function or class; usually __LINE__.
283
                 */
284
                void Assembly(const char* invoking_file, int invoking_line);
285
286
                
                
287
288
289
290
291
292
293
294
295
                /*!
                 * \brief Indicates if a matrix has been assembled and is ready for use; for example, in matrix-vector product.
                 *
                 * \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__.
                 */
                bool IsAssembled(const char* invoking_file, int invoking_line) const;
                
                
296
                /*!
297
298
                 * \brief Wrapper over MatScale, that performs Y = a * Y.
                 *
299
300
                 * \param[in] a Factor by which the vector is scaled.
                 * \param[in] invoking_file File that invoked the function or class; usually __FILE__.
301
                 * \param[in] invoking_line File that invoked the function or class; usually __LINE__.
302
                 */
303
                void Scale(PetscScalar a, const char* invoking_file, int invoking_line);
304
305
306
307
308
                
                
                /*!
                 * \brief Add or modify values inside a Petsc matrix.
                 *
309
310
                 * \warning Assembly() must be called afterwards!
                 *
GILLES Sebastien's avatar
GILLES Sebastien committed
311
312
                 * \param[in] row_indexing Program-wise index of the rows which values will be set.
                 * \param[in] column_indexing Program-wise index of the columns which values will be set.
313
314
315
                 * Should be the same size as \a row_indexing.
                 * \param[in] values Values to put in the matrix. This array should be the same size as \a row_indexing
                 * (unfortunately we can't check that here as it is a C array).
316
                 * \param [in] insertOrAppend Petsc ADD_VALUES or INSERT_VALUES (see Petsc documentation).
317
318
                 * \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__.
319
320
321
322
                 */
                void SetValues(const std::vector<PetscInt>& row_indexing,
                               const std::vector<PetscInt>& column_indexing,
                               const PetscScalar* values, InsertMode insertOrAppend,
323
                               const char* invoking_file, int invoking_line);
324
                
325
326
327
328
329
330
331
332
                
                /*!
                 * \brief Set the values of a row.
                 *
                 * \warning Pattern must have been set beforehand!
                 * \warning Assembly() must be called afterwards!
                 *
                 * \param[in] row_index Program-wise index of the row which values will be set.
333
                 * \param[in] values Non-zero values to report in the matrix.
334
335
336
337
338
339
                 * \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__.
                 */
                void SetValuesRow(PetscInt row_index,
                                  const PetscScalar* values,
                                  const char* invoking_file, int invoking_line);
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
                
                
                /*!
                 * \brief Add or modify a block of values into a Petsc matrix.
                 *
                 * According to Petsc documentation, it is likely more efficient than \a SetValues().
                 *
                 * \warning Assembly() must be called afterwards!
                 *
                 * \param[in] row_indexing Program-wise index of the rows which values will be set.
                 * \param[in] column_indexing Program-wise index of the columns which values will be set.
                 * Should be the same size as \a row_indexing.
                 * \param[in] values Values to put in the matrix. This array should be the same size as \a row_indexing
                 * (unfortunately we can't check that here as it is a C array).
                 * \param [in] insertOrAppend Petsc ADD_VALUES or INSERT_VALUES (see Petsc documentation).
                 * \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__.
                 */
                void SetValuesBlocked(const std::vector<PetscInt>& row_indexing,
                                      const std::vector<PetscInt>& column_indexing,
                                      const PetscScalar* values,
                                      InsertMode insertOrAppend,
                                      const char* invoking_file, int invoking_line);
363
364

                
365
         
366
367
368
                /*!
                 * \brief Set a single entry into a Petsc matrix.
                 *
369
370
                 * \warning Assembly() must be called afterwards!
                 *
371
372
373
374
375
376
                 * As noted in Petsc manual pages, if several values the SetValues() method should be preferred.
                 *
                 * \param[in] row_index Index of the row which value will be set.
                 * \param[in] column_index Index of the column which value will be set.
                 * \param[in] value Value to put in the matrix.
                 * \param [in] insertOrAppend Petsc ADD_VALUES or INSERT_VALUES (see Petsc documentation).
377
378
                 * \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__.
379
380
381
382
                 */
                void SetValue(PetscInt row_index,
                              PetscInt column_index,
                              PetscScalar value, InsertMode insertOrAppend,
383
                              const char* invoking_file, int invoking_line);
384
                
385
386
387
388
389
390
391
392
                
                /*!
                 * \brief Get the value of an element of the matrix.
                 *
                 * The underlying Petsc function requires that the row belong to the current processor; however for 
                 * the column it is much more liberal: if there is an attempt to get a value that doesn't belong
                 * to the CSR pattern, the function simply returns 0.
                 *
393
                 * \param[in] row_index Program-wise index of the row. The row MUST be
394
                 * one of the row handled by the current processor.
395
                 * \param[in] column_index Program-wise index of the column.
396
397
398
399
400
401
402
                 * \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__.
                 */
                PetscScalar GetValue(PetscInt row_index, PetscInt column_index,
                                     const char* invoking_file, int invoking_line) const;
                
                
403
404
405
406
407
408
409
410
411
412
                /*!
                 * \brief Get a row for the matrix.
                 *
                 * This is a wrapper over MatGetRow; it should be noticed that Petsc people hope this function is 
                 * actually not used...
                 *
                 * \param[in] row_index Program-wise index of the row. The row MUST be
                 * one of the row handled by the current processor.
                 * \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__.
413
414
415
                 * \param[out] row_content Key is position of the non-zero values, value the actual value.
                 * \todo At the moment row_content is allocated each time; it might be useful not to do this... But this
                 * means structure has to be kept for each line of the matrix considered.
416
                 */
417
418
                void GetRow(PetscInt row_index,
                            std::vector<std::pair<PetscInt, PetscScalar>>& row_content,
419
420
                            const char* invoking_file, int invoking_line) const;
                
421
                
422
423
424
                /*!
                 * \brief Wrapper over MatView.
                 *
425
                 * \copydetails doxygen_hide_mpi_param
426
427
428
                 * \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__.
                 */
429
                void View(const Mpi& mpi, const char* invoking_file, int invoking_line) const;
430
                
431
                
432
433
434
435
436
                /*!
                 * \brief Wrapper over MatView in the case the viewer is a file.
                 *
                 * \param[in] format Format in which the matrix is written. See Petsc manual pages to get all the
                 * formats available; relevant ones so far are PETSC_VIEWER_DEFAULT	and PETSC_VIEWER_ASCII_MATLAB.
437
                 * \copydetails doxygen_hide_mpi_param
438
                 * \param[in] output_file File into which the vector content will be written.
439
440
                 * \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__.
441
                 */
442
                void View(const Mpi& mpi, const std::string& output_file, const char* invoking_file, int invoking_line,
443
                          PetscViewerFormat format = PETSC_VIEWER_DEFAULT) const;
444
                
445
                /*!
446
                 * \brief Wrapper over MatNorm.
447
                 *
448
                 * Available norms are NORM_1, NORM_FROBENIUS and NORM_INFINITY.
449
450
451
                 */
                double Norm(NormType type, const char* invoking_file, int invoking_line) const;
                
452
                
453
454
            private:
                
455
                // ===========================================================================
456
                // \attention Do not forget to update Swap() if a new data member is added!
457
458
                // =============================================================================
                
459
                //! Underlying Petsc matrix.
460
461
462
                Mat petsc_matrix_;
                
                /*!
463
                 * \brief Whether the underlying Petsc matrix will be destroyed upon destruction of the object.
464
                 *
465
                 * Default behaviour is to do so, but in some cases it is wiser not to.
466
                 */
467
                bool do_petsc_destroy_;
468
469
470
                
                //! Friendship.
                friend void Swap(Matrix& A, Matrix& B);
471
472
473
            };
            
            
474
475
476
477
478
            /*!
             * \brief Swap two matrices.
             */
            void Swap(Matrix& A, Matrix& B);
            
479
                  
480
481
482
483
484
485
486
487
488
        } // namespace Petsc
        
        
    } // namespace Wrappers
    
    
} // namespace HappyHeart


489
# include "ThirdParty/Wrappers/Petsc/Matrix/Matrix.hxx"
490
491


492
#endif // HAPPY_HEART_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_MATRIX_x_MATRIX_HPP_