Mentions légales du service

Skip to content
Snippets Groups Projects
Vector.hpp 32.51 KiB
//! \file 
//
//
//  Petsc.h
//  HappyHeart
//
//  Created by Sébastien Gilles on 20/09/13.
//  Copyright (c) 2013 Inria. All rights reserved.
//

#ifndef HAPPY_HEART_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_VECTOR_x_VECTOR_HPP_
# define HAPPY_HEART_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_VECTOR_x_VECTOR_HPP_


# include <vector>
# include <memory>

# include "Utilities/Miscellaneous.hpp"

# include "ThirdParty/Wrappers/Mpi/Mpi.hpp"
# include "ThirdParty/Wrappers/Mpi/MpiScale.hpp"

# include "ThirdParty/IncludeWithoutWarning/Petsc/PetscVec.hpp"
# include "ThirdParty/Wrappers/Petsc/Vector/Private/VectorHelper.hpp"
# include "ThirdParty/Wrappers/Mpi/Mpi.hpp"
# include "ThirdParty/Wrappers/Mpi/MpiScale.hpp"


namespace HappyHeart
{
    
    
    // ============================
    //! \cond IGNORE_BLOCK_IN_DOXYGEN
    // Forward declarations.    
    // ============================
    
    
    namespace Wrappers
    {
        
        
        class Mpi;
        
        
        namespace Petsc
        {
            
            
            template<Utilities::Access AccessT>
            class AccessVectorContent;
            
            
        } // namespace Petsc
        
        
    } // namespace Wrappers

    
    // ============================
    // End of forward declarations.
    //! \endcond IGNORE_BLOCK_IN_DOXYGEN
    // ============================

    
    
    namespace Wrappers
    {
        
        namespace Petsc
        {
            
            
            /*!
             * \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.
             *
             * \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:
             *
             * - 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.
             * - 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.
             * - There is a constructor that takes as argument a Petsc Vec object. It should be avoided as much
             * as possible (current class should oversee most of vector operations) but is nonetheless required 
             * 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.
             *
             */
            class Vector
            {
            public:
                
                //! Alias to unique_ptr.
                using unique_ptr = std::unique_ptr<Vector>;
                
            public:
                
                
                /// \name Special members.
                
                ///@{
                
                
                //! Constructor.
                explicit Vector();
                
                /*!
                 * \brief Constructor from an existing Petsc Vec.
                 *
                 * \internal <b><tt>[internal]</tt></b> Avoid this as much as possible (current class should hold all Petsc Vector information)
                 * but in some very specific case (for instance defining a SNES function) it is easier at the moment 
                 * to resort to this one.
                 *
                 * \param[in] petsc_vector The petsc vector encapsulated within the object.
                 * \param[in] do_destroy_petsc Whether the underlying Petsc Vec must be destroyed or not in the 
                 * destructor.
                 *
                 *
                 */
                explicit Vector(const Vec& petsc_vector, bool do_destroy_petsc);
                
                /*!
                 * \brief Constructor from a 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.
                 * \copydoc doxygen_hide_invoking_file_and_line
                 */
                explicit Vector(const Mpi& mpi,
                                const std::string& binary_file,
                                const char* invoking_file, int invoking_line);

       
                //! Destructor.
                ~Vector();
                
                //! Copy constructor: copy both layout and data.
                Vector(const Vector&);
                
                //! Move constructor.
                Vector(Vector&&) = default;
                
                //! Disable affectation operator.
                Vector& operator=(const Vector&) = delete;
                
                //! Disable affectation operator.
                Vector& operator=(Vector&&) = delete;

                ///@}
                
                
            public:
                
                /*!
                 * \brief Set the vector as sequential.
                 *
                 * \param[in] size Number of elements in the vector.
                 * \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__.
                 * \copydetails doxygen_hide_mpi_param
                 */
                void InitSequentialVector(unsigned int size, const Mpi& mpi,
                                          const char* invoking_file, int invoking_line);
                
                
                /*!
                 * \brief Set the vector as parallel.
                 *
                 * \param[in] size_processor_wise Size of the vector on the current processor (the ghost are not counted
                 * here).
                 * \param[in] size_program_wise Size of the vector on the program (i.e. the one before partitioning).
                 * \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__.
                 * \copydetails doxygen_hide_mpi_param
                 */
                void InitMpiVector(unsigned int size_processor_wise , unsigned int size_program_wise, const Mpi& mpi,
                                   const char* invoking_file, int invoking_line);
                
                
                /*!
                 * \brief Create a parallel vector with ghost padding.
                 *
                 * \param[in] size_processor_wise Size of the vector on the current processor (the ghost are not counted
                 * here).
                 * \param[in] size_program_wise Size of the vector on the program (i.e. the one before partitioning).
                 * \param[in] ghost_padding List of program-wise index of values that are ghosted (i.e. required
                 * processor-wise but owned by another 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__.
                 * \copydetails doxygen_hide_mpi_param
                 */
                void InitMpiVectorWithGhost(unsigned int size_processor_wise , unsigned int size_program_wise,
                                            const std::vector<PetscInt>& ghost_padding, const Mpi& mpi,
                                            const char* invoking_file, int invoking_line);
                
                
                /*!
                 * \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.
                 *
                 * Parallel case is not handled at all at the moment!
                 *
                 * \param[in] file File from which vector content is read.
                 * \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__.
                 * \copydetails doxygen_hide_mpi_param
                 */
                void InitSequentialFromFile(const std::string& file,
                                            const Mpi& mpi,
                                            const char* invoking_file, int invoking_line);
                
                
                /*!
                 * \brief Handle over the internal Vec object.
                 *
                 * \return Internal Vec object, which is indeed a pointer in Petsc.
                 *
                 * 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.
                 */
                Vec Internal() const;
                
                /*!
                 * \brief Same as internal, except the pointer might be NULL (no check on it).
                 *
                 * \return Internal Vec object, which is indeed a pointer in Petsc.
                 *
                 * \internal <b><tt>[internal]</tt></b> I need this for instance for AccessGhostContent: I need to pass the address 
                 * of the internal pointer to the Petsc function that will initialize it properly.
                 */
                Vec InternalWithoutCheck() const;
                

                /*!
                 * \brief Duplicate layout (in memory, processor repartition, etc...)
                 *
                 * \param[in] original Vector which layout is copied.
                 * \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 DuplicateLayout(const Vector& original, const char* invoking_file, int invoking_line);
                
                
                /*!
                 * \brief Get the number of elements in the local vector.
                 *
                 * \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__.
                 *
                 * \return Number of processor-wise elements of the vector (ghost excluded).
                 */
                PetscInt GetProcessorWiseSize(const char* invoking_file, int invoking_line) const;
                
                /*!
                 * \brief Get the number of elements in the global vector.
                 *
                 * \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__.
                 *
                 * \return Number of program-wise elements of the vector.
                 */
                PetscInt GetProgramWiseSize(const char* invoking_file, int invoking_line) const;
                
                /*!
                 * \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__.
                 */
                void ZeroEntries(const char* invoking_file, int invoking_line);
                
                /*!
                 * \brief Set all the entries to alpha.
                 *
                 * \param[in] alpha value which will fill wll the vector entries.
                 * \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 Set(const PetscScalar alpha, const char* invoking_file, int invoking_line);
                
                
                //! Convenient enum class for Assembly().
                enum class DoUpdateGhosts { yes, no };
                
                /*!
                 * \brief Petsc Assembling of the vector.
                 *
                 * \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] do_update_ghosts If 'yes', update the values of ghosts just after the assembling is done.
                 * 'yes' is therefore more secure but lead to slower code if you don't actually need the updated values
                 * of the ghosts. 'yes' is set by default.
                 *
                 */
                void Assembly(const char* invoking_file, int invoking_line,
                              DoUpdateGhosts do_update_ghosts = DoUpdateGhosts::yes);
                                
                /*!
                 * \brief Add or modify values inside a Petsc vector.
                 *
                 * \param[in] indexing All indexes (program-wise) that have to be modified in the vector are stored here.
                 * \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).
                 * \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 SetValues(const std::vector<PetscInt>& indexing,
                               const PetscScalar* values,
                               InsertMode insertOrAppend,
                               const char* invoking_file, int invoking_line);
                
                
                           
                /*!
                 * \brief Add or modify values inside a Petsc vector.
                 *
                 * \param[in] indexing All indexes (program-wise) (program-wise) that have to be modified in the vector are stored here.
                 * \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).
                 * \param[in] invoking_file File that invoked the function or class; usually __FILE__.
                 * \param[in] invoking_line File that invoked the function or class; usually __LINE__.
                 */
                template<Utilities::Access AccessT>
                void SetValues(const std::vector<PetscInt>& indexing,
                               const AccessVectorContent<AccessT>& local_vec,
                               InsertMode insertOrAppend,
                               const char* invoking_file, int invoking_line);
                
                /*!
                 * \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.
                 *
                 */
                void SetFromPetscVec(const Vec& petsc_vector, const char* invoking_file, int invoking_line);

                
                /*!
                 * \brief Get the values of ta vctor on a current process.
                 *
                 * Used to get the values of a sequential vector (see AccessVectorContent for mpi vectors).
                 *
                 * \param[in] indexing All indexes which are required from the vector.
                 * \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__.
                 *
                 * \return The values for all the indexes given in input.
                 */
                std::vector<PetscScalar> GetValues(const std::vector<PetscInt>& indexing,
                                                   const char* invoking_file, int invoking_line) const;
                
                
                /*!
                 * \brief Same as GetValues() for a unique index.
                 *
                 * \param[in] index Index which associated value is seeked.
                 * \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__.
                 *
                 * \return Associated value.
                 */

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

                
                
                /*!
                 * \brief Add or modify one value inside a Petsc vector.
                 *
                 * \param[in] index Index (program-wise) to be modified.
                 * \param[in] value Value to set.
                 * \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 SetValue(PetscInt index, PetscScalar value, InsertMode insertOrAppend, const char* invoking_file, int invoking_line);
                
                           
                /*!
                 * \brief Set the same value to all entries of the vector.
                 *
                 * \param[in] value Value to set.
                 * \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 SetUniformValue(PetscScalar value, const char* invoking_file, int invoking_line);

                
                
                
                /*!
                 * \brief A wrapper over VecCopy, which assumed target already gets the right layout.
                 *
                 * Doesn't do much except check the return value.
                 * \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.
                 * \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 Copy(const Vector& source, const char* invoking_file, int invoking_line);
                
                
                /*!
                 * \brief A complete copy: layout is copied first and then the values.
                 *
                 * \param[in] source Original vector which layout AND content is copied.
                 * \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 CompleteCopy(const Vector& source, const char* invoking_file, int invoking_line);
                
                
                /*!
                 * \brief Wrapper over VecScale, that performs Y = a * Y.
                 *
                 * \param[in] a Factor by which the vector is scaled.
                 * \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 Scale(PetscScalar a, const char* invoking_file, int invoking_line);
                
                
                /*!
                 * \brief Wrapper over VecView.
                 *
                 * \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__.
                 * \copydetails doxygen_hide_mpi_param
                 */
                void View(const Mpi& mpi, const char* invoking_file, int invoking_line) const;
                
                /*!
                 * \brief Wrapper over VecView 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.
                 * \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__.
                 * \copydetails doxygen_hide_mpi_param
                 * \param[in] output_file File into which the vector content will be written.
                 */
                void View(const Mpi& mpi, const std::string& output_file, const char* invoking_file, int invoking_line,
                          PetscViewerFormat format = PETSC_VIEWER_DEFAULT) const;
                
                
                /*!
                 * \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;
                

                
                
                /*!
                 * \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.
                 * If processor-wise is chosen, this function does not rely on VecView: what I want to achieve is write 
                 * to a different file for each processor and VecView doesn't seem to be able to do so.
                 *
                 * \copydetails doxygen_hide_mpi_param
                 * \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] output_file File into which the vector content will be written.
                 */
                template<MpiScale MpiScaleT>
                void Print(const Mpi& mpi,
                           const std::string& output_file,
                           const char* invoking_file, int invoking_line) const;
                
                /*!
                 * \brief Get the minimum.
                 *
                 * \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__.
                 *
                 * \return First element is the position of the minimum found, second is its value.
                 */
                std::pair<PetscInt, PetscReal> Min(const char* invoking_file, int invoking_line) const;
                
                /*!
                 * \brief Get the maximum.
                 *
                 * \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__.
                 * \return First element is the position of the maximum found, second is its value.
                 */
                std::pair<PetscInt, PetscReal> Max(const char* invoking_file, int invoking_line) const;
                
                /*!
                 * \brief Wrapper over VecNorm. 
                 *
                 * Available norms are NORM_1, NORM_2 and NORM_INFINITY.
                 *
                 */
                double Norm(NormType type, const char* invoking_file, int invoking_line) const;

                /*!
                 * \brief Update the ghost values.
                 *
                 * \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__.
                 * 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...
                 */
                void UpdateGhosts(const char* invoking_file, int invoking_line);
                
                //! Number of ghosts.
                unsigned int Nghost() const;
                
            
                /*!
                 * \brief Tells the class not to destroy the underlying vector through a call to VecDestroy().
                 *
                 * \attention This method should be avoided most of the time; the only case in which it's relevant
                 * is the one exposed in the documentation of member function \a Swap.
                 */
                void SetDoNotDestroyPetscVector();
                
                
            private:
                
                /// \name Methods reserved for friends.
                ///@{
                
                //! Returns the position of the given global index in the ghost vector.
                unsigned int GetPositionInGhostVector(unsigned int global_index) const;                
                 
                
                ///@}
                
                //! Friendship!
                template<Utilities::Access AccessT>
                friend class AccessVectorContent;
                
                friend class AccessGhostContent;
                
                
            private:
                
                
                // ===========================================================================
                // \attention Do not forget to update Swap() if a new data member is added!
                // =============================================================================
                
                //! Underlying Petsc vector.
                Vec petsc_vector_;
                
                /*!
                 * \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_;
                
                //! Friendship.
                friend void Swap(Vector& A, Vector& B);
                
            };
            
            
            
            
            
            //! Swap two Vectors.
            void Swap(Wrappers::Petsc::Vector& A, Wrappers::Petsc::Vector& B);
            
            

         

            /*!
             * \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.
             */
            void DisplaySomeValues(std::ostream& stream, const Vector& vector, PetscInt firstIndex, PetscInt lastIndex,
                                   int rankProc, const char* invoking_file, int invoking_line);


            /*
             * \brief Checks whether Petsc vectors are (almost) equal
             *
             * Note: Petsc proposes VecEqual, but it is not designed to compare two vectors obtained independantly (their
             * documentation has been updated following mail exchanges I had with them).
             *
             * \param[in] epsilon Precision actually required to get an equality (a == b if fabs(a - b) < epsilon)
             * \param[out] inequalityDescription Description of when the discrepancy happened. Empty if true is returned.
             */
            bool AreEqual(const Vector& vec1, const Vector& vec2, double epsilon, std::string& inequalityDescription,
                          const char* invoking_file, int invoking_line);
            

            /*!
             * \brief Wrapper over VecAXPY, that performs Y = alpha * X + Y.
             *
             */
            void AXPY(PetscScalar alpha, const Vector& x, Vector& y, const char* invoking_file, int invoking_line);
            
            
            /*!
             * \brief Wrapper over VecDot, that performs Y = alpha * X + Y.
             *
             */
            void DotProduct(const Vector& x, const Vector& y, PetscScalar& val, const char* invoking_file, int invoking_line);
            

            /*!
             * \brief Gather several mpi vectors into a sequential one.
             */
            void GatherVector(const Mpi& mpi,
                              const Wrappers::Petsc::Vector& local_parallel_vector,
                              Wrappers::Petsc::Vector& sequential_vector,
                              const char* invoking_file, int invoking_line);
            
            
            /*!
             * \brief Perform v3 = tV1 x V2.
             */
            double ScalarProduct(const Mpi& mpi,
                                 const Vector& v1,
                                 const Vector& v2,
                                 const char* invoking_file, int invoking_line);
            
            /*!
             * \brief Compute Norm L2. For development purposes only. Another method exists to compute the norm.
             */
            double NormL2(const Mpi& mpi,
                          const Vector& vector,
                          const char* invoking_file, int invoking_line);

   
            template<Utilities::Access AccessT>
            struct VectorForAccess;
            
            
            template<>
            struct VectorForAccess<Utilities::Access::read_and_write>
            {
                using Type = Vector;
                
                using scalar_array_type = PetscScalar*;
                
                using scalar_type = PetscScalar;
            };

            
            template<>
            struct VectorForAccess<Utilities::Access::read_only>
            {
                using Type = const Vector;
                
                using scalar_array_type = const PetscScalar*;
                
                using scalar_type = const PetscScalar;
            };
            
            
            # ifndef NDEBUG
            
            /*!
             * \brief Assert none of the values is inf or nan.
             */
            void AssertNumericValues(const Vector& vector,
                                     const char* invoking_file, int invoking_line);
            
            
            # endif // NDEBUG


            
        } // namespace Petsc
        
        
    } // namespace Wrappers
    
    
} // namespace HappyHeart


# include "ThirdParty/Wrappers/Petsc/Vector/Vector.hxx"


#endif // HAPPY_HEART_x_THIRD_PARTY_x_WRAPPERS_x_PETSC_x_VECTOR_x_VECTOR_HPP_