-
GILLES Sebastien authoredGILLES Sebastien authored
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_