Vector.hpp 41.9 KB
Newer Older
GILLES Sebastien's avatar
GILLES Sebastien committed
1 2 3 4 5 6 7 8 9 10 11 12 13
/*!
//
// \file
//
//
// Created by Sebastien Gilles <sebastien.gilles@inria.fr> on the Fri, 4 Oct 2013 11:00:51 +0200
// Copyright (c) Inria. All rights reserved.
//
// \ingroup ThirdPartyGroup
// \addtogroup ThirdPartyGroup
// \{
*/

14

15 16
#ifndef MOREFEM_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_VECTOR_x_VECTOR_HPP_
# define MOREFEM_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_VECTOR_x_VECTOR_HPP_
17 18 19


# include <vector>
20
# include <memory>
21

22
# include "Utilities/Miscellaneous.hpp"
23
# include "Utilities/OutputFormat/OutputFormat.hpp"
24

25 26
# include "ThirdParty/Wrappers/Mpi/Mpi.hpp"
# include "ThirdParty/Wrappers/Mpi/MpiScale.hpp"
27

28
# include "ThirdParty/IncludeWithoutWarning/Petsc/PetscVec.hpp"
29
# include "ThirdParty/Wrappers/Petsc/Vector/BinaryOrAscii.hpp"
30
# include "ThirdParty/Wrappers/Petsc/Vector/Internal/VectorHelper.hpp"
31

32

33

34
namespace MoReFEM
35
{
36 37


38
    // ============================
39
    //! \cond IGNORE_BLOCK_IN_DOXYGEN
40
    // Forward declarations.
41
    // ============================
42 43


44 45
    namespace Wrappers
    {
46 47


48
        class Mpi;
49 50


51 52
        namespace Petsc
        {
53 54


55
            template<Utilities::Access AccessT>
56
            class AccessVectorContent;
57 58


59
        } // namespace Petsc
60 61


62 63
    } // namespace Wrappers

64

65 66
    // ============================
    // End of forward declarations.
67
    //! \endcond IGNORE_BLOCK_IN_DOXYGEN
68 69
    // ============================

70

71 72
    //! Enum class to specify if in a copy ghosts are updated or not.
    enum class update_ghost { yes, no };
73

74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
    /*!
     * \class doxygen_hide_parallel_size_arg
     *
     * \param[in] processor_wise_size Number of elements processor-wise (ghosts excluded).
     * \param[in] program_wise_size Number of elements program-wise.
     */


    /*!
     * \class doxygen_hide_parallel_with_ghosts_arg
     *
     * \copydetails doxygen_hide_parallel_size_arg
     * \param[in] ghost_padding List of program-wise index of values that are ghosted (i.e. required
     * processor-wise but owned by another processor).
     */

90

91 92
    namespace Wrappers
    {
93

94 95 96

        namespace Petsc
        {
97 98


99
            /*!
100 101 102 103 104
             * \brief A wrapper class over Petsc Vec objects.
             *
             * Most of the Petsc functions used over Petsc vectors have been included as methods in this class, which
             * also acts as RAII over Petsc Vec object.
             *
105 106 107
             * \internal <b><tt>[internal]</tt></b> This class is way more trickier to implement that it might seems
             * because of Petsc's internal structure: Vec 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:
108 109 110 111
             *
             * - No implicit conversion to the internal Vec. 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.
112
             * - The most usual way to proceed is to construct with the default constructor and then init it either
GILLES Sebastien's avatar
GILLES Sebastien committed
113
             * by a 'Init***()' method or by using 'DuplicateLayout', 'Copy' or 'CompleteCopy' methods.
114
             * - There is a constructor that takes as argument a Petsc Vec object. It should be avoided as much
115
             * as possible (current class should oversee most of vector operations) but is nonetheless required
116 117
             * by user-defined Snes functions. When this constructor is used VecDestroy() is NOT called upon
             * destruction, so that the vector taken as argument still lives.
118
             * \endinternal
119 120 121 122
             *
             */
            class Vector
            {
123
            public:
124

125 126
                //! Alias to unique_ptr.
                using unique_ptr = std::unique_ptr<Vector>;
127

128

129
            public:
130 131


132
                /// \name Special members.
133

134
                ///@{
135 136


137
                //! Constructor.
138
                explicit Vector();
139

140 141 142
                /*!
                 * \brief Constructor from an existing Petsc Vec.
                 *
143
                 * \internal <b><tt>[internal]</tt></b> Avoid this as much as possible (current class should hold all Petsc Vector information)
144
                 * but in some very specific case (for instance defining a SNES function) it is easier at the moment
145
                 * to resort to this one.
146
                 *
147
                 * \endinternal
148
                 * \param[in] petsc_vector The petsc vector encapsulated within the object.
149
                 * \param[in] do_destroy_petsc Whether the underlying Petsc Vec must be destroyed or not in the
150
                 * destructor.
151
                 *
152
                 */
153
                explicit Vector(const Vec& petsc_vector, bool do_destroy_petsc);
154

155
                //! Destructor.
156
                virtual ~Vector();
157

158 159 160
                //! \copydoc doxygen_hide_copy_constructor
                //! Both layout and data are copied.
                Vector(const Vector& rhs);
161

162 163
                //! \copydoc doxygen_hide_move_constructor
                Vector(Vector&& rhs) = default;
164

165 166
                //! \copydoc doxygen_hide_copy_affectation
                Vector& operator=(const Vector& rhs) = delete;
167

168 169
                //! \copydoc doxygen_hide_move_affectation
                Vector& operator=(Vector&& rhs) = delete;
170

171
                ///@}
172 173


174
            public:
175

176 177
                /*!
                 * \brief Set the vector as sequential.
178
                 *
179 180 181
                 * \attention This method is to be called just after creation of the object, whcih should still
                 * be a blank state.
                 *
182
                 * \param[in] size Number of elements in the vector.
183
                 * \copydoc doxygen_hide_invoking_file_and_line
184
                 * \copydetails doxygen_hide_mpi_param
185
                 */
186 187
                void InitSequentialVector(const Mpi& mpi,
                                          unsigned int size,
188
                                          const char* invoking_file, int invoking_line);
189

190
                /*!
191
                 * \brief Set the vector as sequential.
192
                 *
193 194 195 196
                 * \attention This method is to be called just after creation of the object, whcih should still
                 * be a blank state.
                 *
                 * \copydetails doxygen_hide_parallel_size_arg
197
                 * \copydoc doxygen_hide_invoking_file_and_line
198
                 * \copydetails doxygen_hide_mpi_param
199
                 */
200 201 202
                void InitMpiVector(const Mpi& mpi,
                                   unsigned int processor_wise_size,
                                   unsigned int program_wise_size,
203
                                   const char* invoking_file, int invoking_line);
204

205 206
                /*!
                 * \brief Create a parallel vector with ghost padding.
207
                 *
208
                 * \copydetails doxygen_hide_parallel_with_ghosts_arg
209
                 * \copydoc doxygen_hide_invoking_file_and_line
210
                 * \copydetails doxygen_hide_mpi_param
211
                 */
212 213 214 215
                void InitMpiVectorWithGhost(const Mpi& mpi,
                                            unsigned int processor_wise_size ,
                                            unsigned int program_wise_size,
                                            const std::vector<PetscInt>& ghost_padding,
216
                                            const char* invoking_file, int invoking_line);
217 218


219 220 221 222 223 224
                /*!
                 * \brief Init a sequential vector with the data read in the file.
                 *
                 * This file is assumed to have been created with Print() method for a sequential vector.
                 *
                 * \param[in] file File from which vector content is read.
225
                 * \copydoc doxygen_hide_invoking_file_and_line
226
                 * \copydetails doxygen_hide_mpi_param
227
                 */
228 229
                void InitSequentialFromFile(const Mpi& mpi,
                                            const std::string& file,
230
                                            const char* invoking_file, int invoking_line);
231

232 233 234 235 236 237 238 239 240

                /*!
                 * \brief Init a vector from the data read in the file.
                 *
                 * This file is assumed to have been created with Print() method.
                 *
                 * Current method is in fact able to create a sequential vector as well, but use rather
                 * InitSequentialFromFile() with its more friendly API if you need only the sequential case.
                 *
241
                 * \param[in] ascii_file File from which vector content is read.
242 243 244 245 246 247 248 249
                 * \copydetails doxygen_hide_parallel_with_ghosts_arg
                 * \copydoc doxygen_hide_invoking_file_and_line
                 * \copydetails doxygen_hide_mpi_param
                 */
                void InitParallelFromProcessorWiseAsciiFile(const Mpi& mpi,
                                                            unsigned int processor_wise_size,
                                                            unsigned int program_wise_size,
                                                            const std::vector<PetscInt>& ghost_padding,
250
                                                            const std::string& ascii_file,
251 252
                                                            const char* invoking_file, int invoking_line);

253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269
                /*!
                 * \brief Init a vector from the data read in the file.
                 *
                 * This file is assumed to have been created with Print() method.
                 *
                 * \param[in] binary_file File from which vector content is read.
                 * \copydetails doxygen_hide_parallel_with_ghosts_arg
                 * \copydoc doxygen_hide_invoking_file_and_line
                 * \copydetails doxygen_hide_mpi_param
                 */
                void InitParallelFromProcessorWiseBinaryFile(const Mpi& mpi,
                                                             unsigned int processor_wise_size,
                                                             unsigned int program_wise_size,
                                                             const std::vector<PetscInt>& ghost_padding,
                                                             const std::string& binary_file,
                                                             const char* invoking_file, int invoking_line);

270 271 272 273 274 275 276 277 278 279 280 281 282 283
                /*!
                 * \brief Init from a program-wise binary file: load a vector dumped with View() method.
                 *
                 * \param[in] binary_file File from which the vector must be loaded. This file must be in binary format.
                 * \copydetails doxygen_hide_parallel_with_ghosts_arg
                 * \copydoc doxygen_hide_invoking_file_and_line
                 * \copydetails doxygen_hide_mpi_param
                 */
                void InitFromProgramWiseBinaryFile(const Mpi& mpi,
                                                   unsigned int processor_wise_size,
                                                   unsigned int program_wise_size,
                                                   const std::vector<PetscInt>& ghost_padding,
                                                   const std::string& binary_file,
                                                   const char* invoking_file, int invoking_line);
284

285 286 287
                /*!
                 * \brief Handle over the internal Vec object.
                 *
288 289
                 * \return Internal Vec object, which is indeed a pointer in Petsc.
                 *
290 291 292
                 * 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.
                 */
293
                Vec Internal() const;
294

295 296 297
                /*!
                 * \brief Same as internal, except the pointer might be NULL (no check on it).
                 *
298
                 * \return Internal Vec object, which is indeed a pointer in Petsc.
299
                 *
300
                 * \internal <b><tt>[internal]</tt></b> I need this for instance for AccessGhostContent: I need to pass the address
301
                 * of the internal pointer to the Petsc function that will initialize it properly.
302
                 * \endinternal
303 304
                 */
                Vec InternalWithoutCheck() const;
305

306

307 308 309
                /*!
                 * \brief Duplicate layout (in memory, processor repartition, etc...)
                 *
310
                 * \param[in] original Vector which layout is copied.
311
                 * \copydoc doxygen_hide_invoking_file_and_line
312
                 */
313
                void DuplicateLayout(const Vector& original, const char* invoking_file, int invoking_line);
314 315


316 317 318
                /*!
                 * \brief Get the number of elements in the local vector.
                 *
319
                 * \copydoc doxygen_hide_invoking_file_and_line
320 321 322 323
                 *
                 * \return Number of processor-wise elements of the vector (ghost excluded).
                 */
                PetscInt GetProcessorWiseSize(const char* invoking_file, int invoking_line) const;
324

325 326 327
                /*!
                 * \brief Get the number of elements in the global vector.
                 *
328
                 * \copydoc doxygen_hide_invoking_file_and_line
329 330 331 332
                 *
                 * \return Number of program-wise elements of the vector.
                 */
                PetscInt GetProgramWiseSize(const char* invoking_file, int invoking_line) const;
333

334 335 336
                /*!
                 * \brief Set all the entries to zero.
                 *
337
                 * \copydoc doxygen_hide_invoking_file_and_line
338
                 */
339
                void ZeroEntries(const char* invoking_file, int invoking_line);
340 341


342
                /*!
343 344 345 346 347 348 349 350 351
                 * \class doxygen_hide_petsc_vec_assembly
                 *
                 * You have to call Assembly() method after you're done with all your SetXXX() calls; otherwise you
                 * will in all likelihood get error message in parallel about wrong state of the vector (but nothing in
                 * sequential run).
                 * It has not been put directly inside these SetXXX() methods as usually several of them are called in
                 * a row.
                 * For more details, see original explanation in Petsc page about vecSetValues:
                 * http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecSetValues.html#VecSetValues
352 353
                 *
                 */
354 355


356 357
                /*!
                 * \brief Petsc Assembling of the vector.
358
                 *
359 360
                 * \copydoc doxygen_hide_petsc_vec_assembly
                 *
361
                 * \copydoc doxygen_hide_invoking_file_and_line
362
                 * \copydoc doxygen_hide_do_update_ghost_arg
363
                 *
364 365 366
                 * \internal There is a hidden call to it in boundary condition appliance; so maybe you
                 * \endinternal
                 *
367
                 */
368
                void Assembly(const char* invoking_file, int invoking_line,
369
                              update_ghost do_update_ghost = update_ghost::yes);
370

371 372 373
                /*!
                 * \brief Add or modify values inside a Petsc vector.
                 *
374 375
                 * \copydoc doxygen_hide_petsc_vec_assembly
                 *
376
                 * \param[in] indexing All indexes (program-wise) that have to be modified in the vector are stored here.
377 378 379
                 * \param[in] values Values to put in the vector. This array should be the same size as \a 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).
380
                 * \copydoc doxygen_hide_invoking_file_and_line
381
                 *
382
                 */
GILLES Sebastien's avatar
GILLES Sebastien committed
383 384
                void SetValues(const std::vector<PetscInt>& indexing,
                               const PetscScalar* values,
385
                               InsertMode insertOrAppend,
386
                               const char* invoking_file, int invoking_line);
387 388 389



390 391 392
                /*!
                 * \brief Add or modify values inside a Petsc vector.
                 *
393 394
                 * \copydoc doxygen_hide_petsc_vec_assembly
                 *
395
                 * \param[in] indexing All indexes (program-wise) (program-wise) that have to be modified in the vector are stored here.
396 397
                 * \param[in] local_vec Local vector which values will be put inside vector.
                 * \param [in] insertOrAppend Petsc ADD_VALUES or INSERT_VALUES (see Petsc documentation).
398
                 * \copydoc doxygen_hide_invoking_file_and_line
399
                 */
400
                template<Utilities::Access AccessT>
GILLES Sebastien's avatar
GILLES Sebastien committed
401 402
                void SetValues(const std::vector<PetscInt>& indexing,
                               const AccessVectorContent<AccessT>& local_vec,
403
                               InsertMode insertOrAppend,
404
                               const char* invoking_file, int invoking_line);
405

406 407 408 409 410 411 412
                /*!
                 * \brief Set the values of a vector from a Petsc Vec object.
                 *
                 * \attention This method should be used as little as possible: the purpose of current class
                 * is to avoid interacting at all with native Petsc objects. However, in some cases we do not have
                 * the choice: for instance in the definition of a Snes function (for Newton method) we get
                 * a Vec argument.
413 414
                 * \param[in] petsc_vector Petsc \a Vec object.
                 * \copydoc doxygen_hide_invoking_file_and_line
415 416 417 418
                 *
                 */
                void SetFromPetscVec(const Vec& petsc_vector, const char* invoking_file, int invoking_line);

419

420
                /*!
421
                 * \brief Get the values of a vector on a current process.
422
                 *
423 424 425 426
                 * This method allocates a vector at each call, another method with the same name does the same thing
                 * without the allocation at each time.
                 *
                 * Used to get the values of a sequential vector (see \a AccessVectorContent for mpi vectors).
427
                 *
428
                 * \param[in] indexing All indexes which are required from the vector.
429
                 * \copydoc doxygen_hide_invoking_file_and_line
430 431 432 433
                 *
                 * \return The values for all the indexes given in input.
                 */
                std::vector<PetscScalar> GetValues(const std::vector<PetscInt>& indexing,
434
                                                   const char* invoking_file, int invoking_line) const;
435

Gautier Bureau's avatar
Gautier Bureau committed
436
                /*!
437
                 * \brief Get the values of a vector on a current process.
438

439 440 441
                 * Same as the other method with the same name but the vector to contain the values is not allocated
                 * at each call. This method should be called if you need to get the values of a vector with the same
                 * size a lot of time.
Gautier Bureau's avatar
Gautier Bureau committed
442
                 *
443
                 * Used to get the values of a sequential vector (see \a AccessVectorContent for mpi vectors).
Gautier Bureau's avatar
Gautier Bureau committed
444 445 446
                 *
                 * \param[in] indexing All indexes which are required from the vector.
                 * \copydoc doxygen_hide_invoking_file_and_line
447
                 * \param[in] values Vector used to contain the values. Should have been allocated before the call.                 
Gautier Bureau's avatar
Gautier Bureau committed
448 449 450 451 452
                 *
                 */
                void GetValues(const std::vector<PetscInt>& indexing,
                               std::vector<PetscScalar>& values,
                               const char* invoking_file, int invoking_line) const;
453 454


455 456 457 458
                /*!
                 * \brief Same as GetValues() for a unique index.
                 *
                 * \param[in] index Index which associated value is seeked.
459
                 * \copydoc doxygen_hide_invoking_file_and_line
460 461 462 463 464
                 *
                 * \return Associated value.
                 */

                PetscScalar GetValue(PetscInt index,
465
                                     const char* invoking_file, int invoking_line) const;
466

467 468


469 470 471
                /*!
                 * \brief Add or modify one value inside a Petsc vector.
                 *
472 473
                 * \copydoc doxygen_hide_petsc_vec_assembly
                 *
474
                 * \param[in] index Index (program-wise) to be modified.
475 476
                 * \param[in] value Value to set.
                 * \param[in] insertOrAppend Petsc ADD_VALUES or INSERT_VALUES (see Petsc documentation).
477
                 * \copydoc doxygen_hide_invoking_file_and_line
478
                 */
479 480
                void SetValue(PetscInt index, PetscScalar value, InsertMode insertOrAppend,
                              const char* invoking_file, int invoking_line);
481 482


483 484 485
                /*!
                 * \brief Set the same value to all entries of the vector.
                 *
486 487
                 * \copydoc doxygen_hide_petsc_vec_assembly
                 *
488
                 * \param[in] value Value to set.
489
                 * \copydoc doxygen_hide_invoking_file_and_line
490 491 492
                 */
                void SetUniformValue(PetscScalar value, const char* invoking_file, int invoking_line);

493

494 495 496 497 498 499
                /*!
                 * \class doxygen_hide_do_update_ghost_arg
                 *
                 * \param[in] do_update_ghost Whether the target gets its ghost automatically updated or not.
                 * Default is yes.
                 */
500

501 502 503 504
                /*!
                 * \brief A wrapper over VecCopy, which assumed target already gets the right layout.
                 *
                 * Doesn't do much except check the return value.
505 506
                 * \param[in] source Original vector which content is copied. Layout is assumed to be already
                 * the same between object for which the method is called and \a source.
507
                 * \copydoc doxygen_hide_invoking_file_and_line
508
                 * \copydoc doxygen_hide_do_update_ghost_arg
509
                 */
510 511
                void Copy(const Vector& source, const char* invoking_file, int invoking_line,
                          update_ghost do_update_ghost = update_ghost::yes);
512 513


514 515
                /*!
                 * \brief A complete copy: layout is copied first and then the values.
516
                 *
517
                 * \param[in] source Original vector which layout AND content is copied.
518
                 * \copydoc doxygen_hide_invoking_file_and_line
519
                 * \copydoc doxygen_hide_do_update_ghost_arg
520
                 */
521 522
                void CompleteCopy(const Vector& source, const char* invoking_file, int invoking_line,
                                  update_ghost do_update_ghost = update_ghost::yes);
523 524


525
                /*!
526 527
                 * \brief Wrapper over VecScale, that performs Y = a * Y.
                 *
528
                 * \param[in] a Factor by which the vector is scaled.
529
                 * \copydoc doxygen_hide_invoking_file_and_line
530
                 * \copydoc doxygen_hide_do_update_ghost_arg
531
                 */
532 533
                void Scale(PetscScalar a, const char* invoking_file, int invoking_line,
                           update_ghost do_update_ghost = update_ghost::yes);
534 535


Gautier Bureau's avatar
Gautier Bureau committed
536 537 538 539 540 541 542 543 544 545
                /*!
                 * \brief Wrapper over VecShift, that performs Y = Y + a * 1.
                 *
                 * \param[in] a Factor by which the vector is shifted.
                 * \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__.
                 * \copydoc doxygen_hide_do_update_ghost_arg
                 */
                void Shift(PetscScalar a, const char* invoking_file, int invoking_line,
                           update_ghost do_update_ghost = update_ghost::yes);
546 547


548 549 550
                /*!
                 * \brief Wrapper over VecView.
                 *
551
                 * \copydoc doxygen_hide_invoking_file_and_line
552
                 * \copydetails doxygen_hide_mpi_param
553
                 */
554
                void View(const Mpi& mpi, const char* invoking_file, int invoking_line) const;
555

556
                /*!
557
                 * \brief Wrapper over VecView in the case the viewer is a file.
558 559
                 *
                 * \param[in] format Format in which the matrix is written. See Petsc manual pages to get all the
GILLES Sebastien's avatar
GILLES Sebastien committed
560
                 * formats available; relevant ones so far are PETSC_VIEWER_DEFAULT and PETSC_VIEWER_ASCII_MATLAB.
561
                 * \copydoc doxygen_hide_invoking_file_and_line
562
                 * \copydetails doxygen_hide_mpi_param
563
                 * \param[in] output_file File into which the vector content will be written.
564 565
                 *
                 * \attention To my knowledge there are no way to reload a PETSc vector from the file; if you need
566
                 * to do so rather use binary format. 
567
                 *
568
                 */
569
                void View(const Mpi& mpi, const std::string& output_file, const char* invoking_file, int invoking_line,
570
                          PetscViewerFormat format = PETSC_VIEWER_DEFAULT) const;
571 572


573 574 575 576
                /*!
                 * \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.
577
                 * \copydoc doxygen_hide_invoking_file_and_line
578
                 * \copydetails doxygen_hide_mpi_param
579 580 581 582
                 *
                 * \attention If you want to be able to reload the vector later with the exact same structure,
                 * you have to also store its local and global sizes (called processor- and program-wise sizes
                 * within MoReFEM) and the ghost padding.
583 584 585 586 587
                 */
                void ViewBinary(const Mpi& mpi,
                                const std::string& output_file,
                                const char* invoking_file, int invoking_line) const;

588

589 590 591 592 593
                /*!
                 * \brief Print the content of a vector in a file.
                 *
                 * \tparam MpiScaleT Whether we want to print program-wise (in which case View function above is called)
                 * or processor-wise data.
594
                 * If processor-wise is chosen, this function does not rely on VecView: what I want to achieve is write
GILLES Sebastien's avatar
GILLES Sebastien committed
595
                 * to a different file for each processor and VecView doesn't seem to be able to do so.
596
                 *
597
                 * \copydetails doxygen_hide_mpi_param
598
                 * \copydoc doxygen_hide_invoking_file_and_line
599
                 * \param[in] output_file File into which the vector content will be written.
600 601
                 * \param[in] binary_or_ascii_choice Whether the vector should be printed as binary or ascii. Default
                 * value takes its cue from the choice written in the input data file.
602 603 604 605 606 607 608 609
                 *
                 * \attention If the purpose is to be able to reload the file from disk later on, format MUST be
                 * binary (PETSc [VecLoad](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecLoad.html)
                 * works only with binary or HDF5). The file is not enough to rebuild the vector: you will also need
                 * to know:
                 * - The number of elements processor-wise.
                 * - The number of elements program-wise.
                 * - The program-wise indexes on the ghosted elements.
610 611 612 613
                 */
                template<MpiScale MpiScaleT>
                void Print(const Mpi& mpi,
                           const std::string& output_file,
614 615
                           const char* invoking_file, int invoking_line,
                           binary_or_ascii binary_or_ascii_choice = binary_or_ascii::from_input_data) const;
616 617


618 619 620 621
                /*!
                 * \brief Wrapper over VecLoad in the case the viewer is a file.
                 *
                 * \param[in] format Format in which the matrix is loaded. See Petsc manual pages to get all the
GILLES Sebastien's avatar
GILLES Sebastien committed
622
                 * formats available; relevant ones so far are PETSC_VIEWER_DEFAULT and PETSC_VIEWER_ASCII_MATLAB.
623 624 625 626 627 628 629
                 * \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__.
                 * \param[in] mpi Mpi object which knows the rank of the processor, the total number of processors, etc...
                 * \param[in] input_file File into which the vector content will be loaded.
                 */
                void Load(const Mpi& mpi, const std::string& input_file, const char* invoking_file, int invoking_line,
                          PetscViewerFormat format = PETSC_VIEWER_DEFAULT) const;
630

631 632 633
                /*!
                 * \brief Get the minimum.
                 *
634
                 * \copydoc doxygen_hide_invoking_file_and_line
635
                 *
636 637
                 * \return First element is the position of the minimum found, second is its value.
                 */
638
                std::pair<PetscInt, PetscReal> Min(const char* invoking_file, int invoking_line) const;
639

640 641 642
                /*!
                 * \brief Get the maximum.
                 *
643
                 * \copydoc doxygen_hide_invoking_file_and_line
644 645
                 * \return First element is the position of the maximum found, second is its value.
                 */
646
                std::pair<PetscInt, PetscReal> Max(const char* invoking_file, int invoking_line) const;
647

648
                /*!
649
                 * \brief Wrapper over VecNorm.
650 651
                 *
                 * Available norms are NORM_1, NORM_2 and NORM_INFINITY.
652
                 *
653 654 655 656
                 * \param[in] type NORM_1, NORM_2 or NORM_INFINITY.
                 * \copydoc doxygen_hide_invoking_file_and_line
                 *
                 * \return Value of the norm.
657 658
                 */
                double Norm(NormType type, const char* invoking_file, int invoking_line) const;
659

660 661 662
                /*!
                 * \brief Update the ghost values.
                 *
663
                 * \copydoc doxygen_hide_invoking_file_and_line
664 665 666 667 668
                 * Behaviour when the vector is without ghost is up to Petsc (if it results in an error
                 * an exception will be thrown).
                 *
                 * Beware: all processor must make the call to this method! If for instance you're in a loop
                 * and one processor leaves it before the other, the code will be deadlocked...
669
                 */
670
                void UpdateGhosts(const char* invoking_file, int invoking_line);
671 672


673 674 675
                /*!
                 * \brief Update the ghost values if do_update_ghost is set to yes.
                 *
676
                 * This convenient method should be used only in Vector or Matrix related functions that
677 678 679
                 * provides the possibility to automatically update ghosts through an ad hoc argument (see
                 * for instance Copy() method).
                 *
680
                 * \copydoc doxygen_hide_invoking_file_and_line
681 682 683 684 685 686 687 688 689 690
                 * Behaviour when the vector is without ghost is up to Petsc (if it results in an error
                 * an exception will be thrown).
                 *
                 * Beware: all processor must make the call to this method! If for instance you're in a loop
                 * and one processor leaves it before the other, the code will be deadlocked...
                 *
                 * \copydoc doxygen_hide_do_update_ghost_arg
                 */
                void UpdateGhosts(const char* invoking_file, int invoking_line,
                                  update_ghost do_update_ghost);
691

692
                //! Number of ghosts.
693
                unsigned int Nghost() const;
694 695


696 697 698
                /*!
                 * \brief Tells the class not to destroy the underlying vector through a call to VecDestroy().
                 *
699 700 701
                 * \attention This method should be avoided most of the time; the only cases in which it's relevant are:
                 * - the one exposed in the documentation of member function \a Swap.
                
702 703
                 */
                void SetDoNotDestroyPetscVector();
704 705


706 707 708 709 710 711 712 713 714
                /*!
                 * \brief Accessor to the list of program-wise index of values that are ghosted
                 *
                 * (i.e. required processor-wise but owned by another processor).
                 *
                 * Shouldn't be called if not a vector with ghost (an assert checks that).
                 */
                const std::vector<PetscInt>& GetGhostPadding() const noexcept;

715 716 717 718 719 720 721 722 723 724 725 726 727 728 729
            protected:

                /*!
                 * \brief Change the underlying PETSc pointer .
                 *
                 * \param[in] new_petsc_vector The new underlying PETSc \a Vec (which is a typedef to a pointer).
                 *
                 * This functionality was introduced after version 3.12 of PETSc: I realized after the non linear code failed that I didn't set it up properly.
                 * The functions given to \a SNESSetFunction and \a SNESSetJacobian must work on \a Vec and \a Mat objects handled
                 * completely by PETSc; it just happened in prior versions of the library it used up the system matrices and vectors given to
                 * \a SolveNonLinear. So these functions now start by creating on the stack respectively a \a GlobalVector and a \a GlobalMatrix;
                 * \a ChangeInternal is then called so that these temporary objects handle internally the objects provided by PETSc API.
                 *
                 */
                void ChangeInternal(Vec new_petsc_vector);
730

731
            private:
732

733 734 735 736 737
                // ============================
                //! \cond IGNORE_BLOCK_IN_DOXYGEN
                // Friendship.
                // ============================

738

739
                template<Utilities::Access AccessT>
740
                friend class AccessVectorContent;
741

742
                friend class AccessGhostContent;
743

744
                friend void Swap(Vector& lhs, Vector& rhs);
745

746 747 748 749
                // ============================
                //! \endcond IGNORE_BLOCK_IN_DOXYGEN
                // ============================

750

751
            private:
752 753


754
                // ===========================================================================
755
                // \attention Do not forget to update Swap() if a new data member is added!
756
                // =============================================================================
757

758
                //! Underlying Petsc vector.
759
                Vec petsc_vector_;
760

761 762 763 764 765 766 767 768
                /*!
                 * \brief List of program-wise index of values that are ghosted (i.e. required
                 * processor-wise but owned by another processor).
                 *
                 * Left empty if sequential or mpi without ghost.
                 */
                std::vector<PetscInt> ghost_padding_;

769 770 771 772 773 774 775
                /*!
                 * \brief Whether the underlying Petsc vector will be destroyed upon destruction of the object.
                 *
                 * Default behaviour is to do so, but in some cases (for instance when vector has been built from
                 * an existing Petsc Vec) it is wiser not to.
                 */
                bool do_petsc_destroy_;
776

777
            };
778 779


780 781 782 783 784 785
            /*!
             * \brief Swap two Vectors.
             *
             * \copydoc doxygen_hide_lhs_rhs_arg
             */
            void Swap(Vector& lhs, Vector& rhs);
786

787 788 789 790 791

            /*!
             * \brief A quick and dirty way to display some values of a petsc vector (for debug purposes)
             *
             * I'm not even sure it works as intended in parallelism context, but it is quite useful in sequential.
792 793 794 795 796 797 798 799
             *
             *
             * \param[in,out] stream Stream onto which the values are written.
             * \param[in] vector Vector being investigated.
             * \param[in] first_index First (processor-wise) index to be printed.
             * \param[in] last_index Last (processor-wise) index to be printed.
             * \param[in] rank Mpi rank of the current processor.
             * \copydoc doxygen_hide_invoking_file_and_line
800
             */
801 802 803 804 805 806
            void DisplaySomeValues(std::ostream& stream,
                                   const Vector& vector,
                                   PetscInt first_index,
                                   PetscInt last_index,
                                   int rank,
                                   const char* invoking_file, int invoking_line);
807 808


809
            /*!