Matrix.hpp 25.9 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/LinearAlgebra/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
        class Mpi;
        
        
46
47
        namespace Petsc
        {
48
49
50
            
            
            class Vector;
51
            class MatrixPattern;
52
53
            
            
54
55
56
57
58
59
60
61
62
        } // namespace Petsc
        
        
    } // namespace Wrappers
    
    
    
    // ============================
    // End of forward declarations.
63
    //! \endcond IGNORE_BLOCK_IN_DOXYGEN
64
65
66
67
68
69
70
71
72
73
74
75
    // ============================

    
    
    namespace Wrappers
    {
        
        
        namespace Petsc
        {
            
     
76
            
77
            /*!
78
             * \brief A wrapper class over Petsc Mat objects.
79
             *
80
81
             * 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.
82
             *
83
             * \internal <b><tt>[internal]</tt></b> This class is way more trickier to implement that it might seems because of Petsc's
84
85
             * 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:
86
             *
87
88
89
90
91
             * - 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.
92
93
             *
             */
94
            class Matrix : private Private::BaseMatrix
95
            {
96
97
98
99
100
            private:
                
                //! Alias to parent.
                using parent = Private::BaseMatrix;
                
101
            public:
102
103
104
105
                
                //! Constructor.
                explicit Matrix();
                
106
107
108
109
110
111
112
113
114
115
                /*!
                 * \brief Constructor from a file: load a matrix dumped with View() method.
                 *
                 * \param[in] binary_file File from which the matrix must be loaded. This file must be in binary format.
                 * \copydoc doxygen_hide_invoking_file_and_line
                 */
                explicit Matrix(const Mpi& mpi,
                                const std::string& binary_file,
                                const char* invoking_file, int invoking_line);
                
116
117
118
                //! Destructor.
                ~Matrix();
                
119
120
                //! Move constructor.
                Matrix(Matrix&&) = default;
121
                
122
123
                //! Copy constructor: copy both layout and data.
                Matrix(const Matrix&);
124
125
126
127
                
                //! Disable affectation operator.
                Matrix& operator=(const Matrix&) = delete;
                
128
129
130
                //! Disable affectation operator.
                Matrix& operator=(Matrix&&) = delete;
                
131
132
133
134
135
136
137
                
                /*!
                 * \brief Create a sequential sparse matrix.
                 *
                 *
                 * \param[in] Nrow Number of rows.
                 * \param[in] Ncolumn Number of columns.
138
                 * \param[in] matrix_pattern Pattern of the matrix (number of elements expected on each row_).
139
                 * \copydetails doxygen_hide_mpi_param
140
141
                 * \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__.
142
143
144
145
                 *
                 * 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.
                 */
146
                void InitSequentialMatrix(unsigned int Nrow, unsigned int Ncolumn,
147
                                          const MatrixPattern& matrix_pattern,
148
                                          const Mpi& mpi, const char* invoking_file, int invoking_line);
149
150
151
152
153
154
155
156
157
                
                
                /*!
                 * \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.
158
                 * \param[in] matrix_pattern Pattern of the matrix (number of elements expected on each row_).
159
                 * \copydetails doxygen_hide_mpi_param
160
161
                 * \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__.
162
163
                 *
                 */
164
165
                void InitParallelMatrix(unsigned int Nlocal_row, unsigned int Nlocal_column,
                                        unsigned int Nglobal_row, unsigned int Nglobal_column,
166
                                        const MatrixPattern& matrix_pattern,
167
                                        const Mpi& mpi, const char* invoking_file, int invoking_line);
168
                
169
                                
170
171
172
                /*!
                 * \brief Handle over the internal Mat object.
                 *
173
174
                 * \return Internal Mat object, which is indeed a pointer in Petsc.
                 *
175
176
177
                 * 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.
                 */
178
                Mat Internal() const noexcept;
179
                
180
181
182
183
184
185
186
187
188
189
                
                /*!
                 * \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);
                
                
190
191
192
193
194
195
196
197
                /*!
                 * \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.)
198
199
                 * \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__.
200
                 *
201
202
203
204
205
206
                 * 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.
                 *
207
                 */
208
                void DuplicateLayout(const Matrix& original, MatDuplicateOption option, const char* invoking_file, int invoking_line);
209
210
211
212
213
214
                
                
                
                /*!
                 * \brief A wrapper over MatCopy, which assumed target already gets the right layout.
                 *
215
216
                 * \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.
217
218
219
220
221
                 * \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.
222
223
                 * \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__.
224
                 */
225
                void Copy(const Matrix& source, const char* invoking_file, int invoking_line,
226
227
228
229
230
                          const MatStructure& structure = SAME_NONZERO_PATTERN);
                
                
                /*!
                 * \brief A complete copy: layout is copied first and then the values.
231
                 *
232
233
234
                 * \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.
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
                 */
239
                void CompleteCopy(const Matrix& source, const char* invoking_file, int invoking_line);
240
                
241
                /*!
242
                 * \brief Get the number of elements in the local matrix.
243
244
245
                 *
                 * \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__.
246
247
                 *
                 * \return First index is number of rows, second one the number of columns.
248
                 */
249
                std::pair<PetscInt, PetscInt> GetProcessorWiseSize(const char* invoking_file, int invoking_line) const;
250
                
251
252
253
254
255
                /*!
                 * \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__.
256
257
                 *
                 * \return First index is number of rows, second one the number of columns.
258
                 */
259
                std::pair<PetscInt, PetscInt> GetProgramWiseSize(const char* invoking_file, int invoking_line) const;
260
261
                
                
262
263
264
265
266
267
                /*!
                 * \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__.
                 */
268
                void ZeroEntries(const char* invoking_file, int invoking_line);
269
270
271
272
273
274
275
                
                
                /*!
                 * \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.
276
277
278
                 * \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__.
                 
279
280
281
282
                 */
                
                void ZeroRows(const std::vector<PetscInt>& row_indexes,
                              PetscScalar diagonal_value,
283
                              const char* invoking_file, int invoking_line);
284
285
286
287
                
                
                
                /*!
288
                 * \brief Petsc Assembling of the matrix.
289
                 *
290
                 * \param[in] invoking_file File that invoked the function or class; usually __FILE__.
291
                 * \param[in] invoking_line File that invoked the function or class; usually __LINE__.
292
                 */
293
                void Assembly(const char* invoking_file, int invoking_line);
294
295
                
                
296
297
298
299
300
301
302
303
304
                /*!
                 * \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;
                
                
305
                /*!
306
307
                 * \brief Wrapper over MatScale, that performs Y = a * Y.
                 *
308
309
                 * \param[in] a Factor by which the vector is scaled.
                 * \param[in] invoking_file File that invoked the function or class; usually __FILE__.
310
                 * \param[in] invoking_line File that invoked the function or class; usually __LINE__.
311
                 */
312
                void Scale(PetscScalar a, const char* invoking_file, int invoking_line);
313
314
315
316
317
                
                
                /*!
                 * \brief Add or modify values inside a Petsc matrix.
                 *
318
319
320
321
322
323
324
                 * \note Petsc documentation says you should use this over MatSetValue, but usage is actually not that 
                 * wide: you actually define here a block of values. So for instance you choose \a row_indexing = { 1, 2 }
                 * and \a col_indexing = { 3, 5 }, you modify the four values (1, 3), (1, 5), (2, 3), (2, 5)... (which
                 * might fail if some of it are actually not non zero values!)
                 * To my knowledge there are no way to modify just (1, 3) and (2, 5) in one call! I skimmed a bit through
                 * petsc-user help list, and they seem to suggest MatSetValues should be called for one row at a time.
                 *
325
326
                 * \warning Assembly() must be called afterwards!
                 *
GILLES Sebastien's avatar
GILLES Sebastien committed
327
328
                 * \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.
329
                 * \param[in] values Values to put in the matrix. This array should be the same size as \a row_indexing
330
                 * x \a column_indexing: it is a block of values that is actually introduced!
331
                 * \param [in] insertOrAppend Petsc ADD_VALUES or INSERT_VALUES (see Petsc documentation).
332
333
                 * \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__.
334
335
336
337
                 */
                void SetValues(const std::vector<PetscInt>& row_indexing,
                               const std::vector<PetscInt>& column_indexing,
                               const PetscScalar* values, InsertMode insertOrAppend,
338
                               const char* invoking_file, int invoking_line);
339
                
340
341
342
343
344
345
346
347
                
                /*!
                 * \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.
348
                 * \param[in] values Non-zero values to report in the matrix.
349
350
351
352
353
354
                 * \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);
355
                
356
357
358
                /*!
                 * \brief Set a single entry into a Petsc matrix.
                 *
359
360
                 * \warning Assembly() must be called afterwards!
                 *
361
362
                 * As noted in Petsc manual pages, if several values the SetValues() method should be preferred;
                 * however see \a SetValues() documentation that specify the cases in which it can be applied.
363
364
365
366
367
                 *
                 * \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).
368
369
                 * \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__.
370
371
372
373
                 */
                void SetValue(PetscInt row_index,
                              PetscInt column_index,
                              PetscScalar value, InsertMode insertOrAppend,
374
                              const char* invoking_file, int invoking_line);
375
                
376
377
378
379
380
381
382
383
                
                /*!
                 * \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.
                 *
384
                 * \param[in] row_index Program-wise index of the row. The row MUST be
385
                 * one of the row handled by the current processor.
386
                 * \param[in] column_index Program-wise index of the column.
387
388
389
390
391
392
393
                 * \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;
                
                
394
395
396
397
398
399
400
401
402
403
                /*!
                 * \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__.
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
                 * \param[out] row_content_position_list Position of the non-zero values in the row.
                 * \param[out] row_content_value_list Values of non-zero values in the row. This container is the same 
                 * size as \a row_content_position
                 * \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.
                 */
                void GetRow(PetscInt row_index,
                            std::vector<PetscInt>& row_content_position_list,
                            std::vector<PetscScalar>& row_content_value_list,
                            const char* invoking_file, int invoking_line) const;
                
                
                /*!
                 * \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__.
426
427
428
                 * \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.
429
                 */
430
431
                void GetRow(PetscInt row_index,
                            std::vector<std::pair<PetscInt, PetscScalar>>& row_content,
432
433
                            const char* invoking_file, int invoking_line) const;
                
434
                
435
                
436
437
438
                /*!
                 * \brief Wrapper over MatView.
                 *
439
                 * \copydetails doxygen_hide_mpi_param
440
441
442
                 * \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__.
                 */
443
                void View(const Mpi& mpi, const char* invoking_file, int invoking_line) const;
444
                
445
                
446
447
448
449
450
                /*!
                 * \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.
451
                 * \copydetails doxygen_hide_mpi_param
452
                 * \param[in] output_file File into which the vector content will be written.
453
454
                 * \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__.
455
                 */
456
                void View(const Mpi& mpi, const std::string& output_file, const char* invoking_file, int invoking_line,
457
                          PetscViewerFormat format = PETSC_VIEWER_DEFAULT) const;
458
                
459
460
461
462
463
464
465
466
467
468
469
470
471
                
                /*!
                 * \brief Wrapper over MatView in the case the viewer is a binary file.
                 *
                 * \param[in] output_file File into which the vector content will be written.
                 * \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 ViewBinary(const Mpi& mpi,
                                const std::string& output_file,
                                const char* invoking_file, int invoking_line) const;

                
472
                /*!
473
                 * \brief Wrapper over MatNorm.
474
                 *
475
                 * Available norms are NORM_1, NORM_FROBENIUS and NORM_INFINITY.
476
477
478
                 */
                double Norm(NormType type, const char* invoking_file, int invoking_line) const;
                
479
                
480
481
            private:
                
482
                // ===========================================================================
483
                // \attention Do not forget to update Swap() if a new data member is added!
484
485
                // =============================================================================
                
486
                //! Underlying Petsc matrix.
487
488
489
                Mat petsc_matrix_;
                
                /*!
490
                 * \brief Whether the underlying Petsc matrix will be destroyed upon destruction of the object.
491
                 *
492
                 * Default behaviour is to do so, but in some cases it is wiser not to.
493
                 */
494
                bool do_petsc_destroy_;
495
496
497
                
                //! Friendship.
                friend void Swap(Matrix& A, Matrix& B);
498
499
500
            };
            
            
501
502
503
504
505
            /*!
             * \brief Swap two matrices.
             */
            void Swap(Matrix& A, Matrix& B);
            
506
                  
507
508
509
510
511
512
513
514
515
        } // namespace Petsc
        
        
    } // namespace Wrappers
    
    
} // namespace HappyHeart


516
# include "ThirdParty/Wrappers/Petsc/Matrix/Matrix.hxx"
517
518


519
#endif // HAPPY_HEART_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_MATRIX_x_MATRIX_HPP_