read_write_vector.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2015 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_read_write_vector_h
#define dealii_read_write_vector_h

#include <deal.II/base/config.h>

#include <deal.II/base/index_set.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/base/types.h>
#include <deal.II/base/utilities.h>

#include <deal.II/lac/vector_operation.h>

#include <cstdlib>
#include <cstring>
#include <iomanip>

#ifdef DEAL_II_WITH_TRILINOS
#  include <deal.II/lac/trilinos_epetra_communication_pattern.h>
#  include <deal.II/lac/trilinos_epetra_vector.h>

#  include <Epetra_MultiVector.h>

#endif

DEAL_II_NAMESPACE_OPEN

namespace LinearAlgebra
{
  class CommunicationPatternBase;
  namespace distributed
  {
    template <typename>
    class Vector;
  }
} // namespace LinearAlgebra

#ifdef DEAL_II_WITH_PETSC
namespace PETScWrappers
{
  namespace MPI
  {
    class Vector;
  }
} // namespace PETScWrappers
#endif

#ifdef DEAL_II_WITH_TRILINOS
namespace TrilinosWrappers
{
  namespace MPI
  {
    class Vector;
  }
} // namespace TrilinosWrappers
#endif

#ifdef DEAL_II_WITH_CUDA
namespace LinearAlgebra
{
  namespace CUDAWrappers
  {
    template <typename>
    class Vector;
  }
} // namespace LinearAlgebra
#endif

namespace LinearAlgebra
{
  template <typename Number>
  class ReadWriteVector : public Subscriptor
  {
  public:
    using value_type      = Number;
    using pointer         = value_type *;
    using const_pointer   = const value_type *;
    using iterator        = value_type *;
    using const_iterator  = const value_type *;
    using reference       = value_type &;
    using const_reference = const value_type &;
    using size_type       = types::global_dof_index;
    using real_type       = typename numbers::NumberTraits<Number>::real_type;

    ReadWriteVector();

    ReadWriteVector(const ReadWriteVector<Number> &in_vector);

    explicit ReadWriteVector(const size_type size);

    explicit ReadWriteVector(const IndexSet &locally_stored_indices);

    ~ReadWriteVector() override = default;

    virtual void
    reinit(const size_type size, const bool omit_zeroing_entries = false);

    template <typename Number2>
    void
    reinit(const ReadWriteVector<Number2> &in_vector,
           const bool                      omit_zeroing_entries = false);

    virtual void
    reinit(const IndexSet &locally_stored_indices,
           const bool      omit_zeroing_entries = false);


#ifdef DEAL_II_WITH_TRILINOS
#  ifdef DEAL_II_WITH_MPI

    void
    reinit(const TrilinosWrappers::MPI::Vector &trilinos_vec);
#  endif
#endif

    template <typename Functor>
    void
    apply(const Functor &func);

    void
    swap(ReadWriteVector<Number> &v);

    ReadWriteVector<Number> &
    operator=(const ReadWriteVector<Number> &in_vector);

    template <typename Number2>
    ReadWriteVector<Number> &
    operator=(const ReadWriteVector<Number2> &in_vector);

    ReadWriteVector<Number> &
    operator=(const Number s);

    void
    import(const distributed::Vector<Number> &vec,
           VectorOperation::values            operation,
           const std::shared_ptr<const CommunicationPatternBase>
             &communication_pattern =
               std::shared_ptr<const CommunicationPatternBase>());

#ifdef DEAL_II_WITH_PETSC

    void
    import(const PETScWrappers::MPI::Vector &petsc_vec,
           VectorOperation::values           operation,
           const std::shared_ptr<const CommunicationPatternBase>
             &communication_pattern =
               std::shared_ptr<const CommunicationPatternBase>());
#endif

#ifdef DEAL_II_WITH_TRILINOS

    void
    import(
      const TrilinosWrappers::MPI::Vector &           trilinos_vec,
      VectorOperation::values                         operation,
      std::shared_ptr<const CommunicationPatternBase> communication_pattern =
        std::shared_ptr<const CommunicationPatternBase>());

#  ifdef DEAL_II_WITH_MPI

    void
    import(
      const EpetraWrappers::Vector &                  epetra_vec,
      VectorOperation::values                         operation,
      std::shared_ptr<const CommunicationPatternBase> communication_pattern =
        std::shared_ptr<const CommunicationPatternBase>());
#  endif
#endif

#ifdef DEAL_II_WITH_CUDA

    void
    import(
      const CUDAWrappers::Vector<Number> &            cuda_vec,
      VectorOperation::values                         operation,
      std::shared_ptr<const CommunicationPatternBase> communication_pattern =
        std::shared_ptr<const CommunicationPatternBase>());
#endif

    size_type
    size() const;

    size_type
    n_elements() const;

    const IndexSet &
    get_stored_elements() const;

    iterator
    begin();

    const_iterator
    begin() const;

    iterator
    end();

    const_iterator
    end() const;



    Number
    operator()(const size_type global_index) const;

    Number &
    operator()(const size_type global_index);

    Number operator[](const size_type global_index) const;

    Number &operator[](const size_type global_index);

    template <typename Number2>
    void
    extract_subvector_to(const std::vector<size_type> &indices,
                         std::vector<Number2> &        values) const;

    template <typename ForwardIterator, typename OutputIterator>
    void
    extract_subvector_to(ForwardIterator       indices_begin,
                         const ForwardIterator indices_end,
                         OutputIterator        values_begin) const;

    Number
    local_element(const size_type local_index) const;

    Number &
    local_element(const size_type local_index);



    template <typename Number2>
    void
    add(const std::vector<size_type> &indices,
        const std::vector<Number2> &  values);

    template <typename Number2>
    void
    add(const std::vector<size_type> &  indices,
        const ReadWriteVector<Number2> &values);

    template <typename Number2>
    void
    add(const size_type  n_elements,
        const size_type *indices,
        const Number2 *  values);

    void
    print(std::ostream &     out,
          const unsigned int precision  = 3,
          const bool         scientific = true) const;

    std::size_t
    memory_consumption() const;

  protected:
#ifdef DEAL_II_WITH_TRILINOS

    void
    import(const Epetra_MultiVector &multivector,
           const IndexSet &          locally_owned_elements,
           VectorOperation::values   operation,
           const MPI_Comm &          mpi_comm,
           const std::shared_ptr<const CommunicationPatternBase>
             &communication_pattern);
#endif

    unsigned int
    global_to_local(const types::global_dof_index global_index) const
    {
      // the following will throw an exception if the global_index is not
      // in the remaining_elements
      return static_cast<unsigned int>(
        stored_elements.index_within_set(global_index));
    }

    void
    resize_val(const size_type new_allocated_size);

#if defined(DEAL_II_WITH_TRILINOS) && defined(DEAL_II_WITH_MPI)

    EpetraWrappers::CommunicationPattern
    create_epetra_comm_pattern(const IndexSet &source_index_set,
                               const MPI_Comm &mpi_comm);
#endif

    IndexSet stored_elements;

    IndexSet source_stored_elements;

    std::shared_ptr<CommunicationPatternBase> comm_pattern;

    std::unique_ptr<Number[], decltype(free) *> values;

    mutable std::shared_ptr<::parallel::internal::TBBPartitioner>
      thread_loop_partitioner;

    template <typename Number2>
    friend class ReadWriteVector;

  private:
    template <typename Functor>
    class FunctorTemplate
    {
    public:
      FunctorTemplate(ReadWriteVector<Number> &parent, const Functor &functor);

      virtual void
      operator()(const size_type begin, const size_type end);

    private:
      ReadWriteVector &parent;

      const Functor &functor;
    };
  };

  /*---------------------------- Inline functions ---------------------------*/

#ifndef DOXYGEN

  template <typename Number>
  inline ReadWriteVector<Number>::ReadWriteVector()
    : Subscriptor()
    , values(nullptr, free)
  {
    // virtual functions called in constructors and destructors never use the
    // override in a derived class
    // for clarity be explicit on which function is called
    ReadWriteVector<Number>::reinit(0, true);
  }



  template <typename Number>
  inline ReadWriteVector<Number>::ReadWriteVector(
    const ReadWriteVector<Number> &v)
    : Subscriptor()
    , values(nullptr, free)
  {
    this->operator=(v);
  }



  template <typename Number>
  inline ReadWriteVector<Number>::ReadWriteVector(const size_type size)
    : Subscriptor()
    , values(nullptr, free)
  {
    // virtual functions called in constructors and destructors never use the
    // override in a derived class
    // for clarity be explicit on which function is called
    ReadWriteVector<Number>::reinit(size, false);
  }



  template <typename Number>
  inline ReadWriteVector<Number>::ReadWriteVector(
    const IndexSet &locally_stored_indices)
    : Subscriptor()
    , values(nullptr, free)
  {
    // virtual functions called in constructors and destructors never use the
    // override in a derived class
    // for clarity be explicit on which function is called
    ReadWriteVector<Number>::reinit(locally_stored_indices);
  }



  template <typename Number>
  inline typename ReadWriteVector<Number>::size_type
  ReadWriteVector<Number>::size() const
  {
    return stored_elements.size();
  }



  template <typename Number>
  inline typename ReadWriteVector<Number>::size_type
  ReadWriteVector<Number>::n_elements() const
  {
    return stored_elements.n_elements();
  }



  template <typename Number>
  inline const IndexSet &
  ReadWriteVector<Number>::get_stored_elements() const
  {
    return stored_elements;
  }



  template <typename Number>
  inline typename ReadWriteVector<Number>::iterator
  ReadWriteVector<Number>::begin()
  {
    return values.get();
  }



  template <typename Number>
  inline typename ReadWriteVector<Number>::const_iterator
  ReadWriteVector<Number>::begin() const
  {
    return values.get();
  }



  template <typename Number>
  inline typename ReadWriteVector<Number>::iterator
  ReadWriteVector<Number>::end()
  {
    return values.get() + this->n_elements();
  }



  template <typename Number>
  inline typename ReadWriteVector<Number>::const_iterator
  ReadWriteVector<Number>::end() const
  {
    return values.get() + this->n_elements();
  }



  template <typename Number>
  inline Number
  ReadWriteVector<Number>::operator()(const size_type global_index) const
  {
    return values[global_to_local(global_index)];
  }



  template <typename Number>
  inline Number &
  ReadWriteVector<Number>::operator()(const size_type global_index)
  {
    return values[global_to_local(global_index)];
  }



  template <typename Number>
  inline Number ReadWriteVector<Number>::
                operator[](const size_type global_index) const
  {
    return operator()(global_index);
  }



  template <typename Number>
  inline Number &ReadWriteVector<Number>::
                 operator[](const size_type global_index)
  {
    return operator()(global_index);
  }



  template <typename Number>
  template <typename Number2>
  inline void
  ReadWriteVector<Number>::extract_subvector_to(
    const std::vector<size_type> &indices,
    std::vector<Number2> &        extracted_values) const
  {
    for (size_type i = 0; i < indices.size(); ++i)
      extracted_values[i] = operator()(indices[i]);
  }



  template <typename Number>
  template <typename ForwardIterator, typename OutputIterator>
  inline void
  ReadWriteVector<Number>::extract_subvector_to(
    ForwardIterator       indices_begin,
    const ForwardIterator indices_end,
    OutputIterator        values_begin) const
  {
    while (indices_begin != indices_end)
      {
        *values_begin = operator()(*indices_begin);
        indices_begin++;
        values_begin++;
      }
  }



  template <typename Number>
  inline Number
  ReadWriteVector<Number>::local_element(const size_type local_index) const
  {
    AssertIndexRange(local_index, this->n_elements());

    return values[local_index];
  }



  template <typename Number>
  inline Number &
  ReadWriteVector<Number>::local_element(const size_type local_index)
  {
    AssertIndexRange(local_index, this->n_elements());

    return values[local_index];
  }



  template <typename Number>
  template <typename Number2>
  inline void
  ReadWriteVector<Number>::add(const std::vector<size_type> &indices,
                               const std::vector<Number2> &  values)
  {
    AssertDimension(indices.size(), values.size());
    add(indices.size(), indices.data(), values.data());
  }



  template <typename Number>
  template <typename Number2>
  inline void
  ReadWriteVector<Number>::add(const std::vector<size_type> &  indices,
                               const ReadWriteVector<Number2> &values)
  {
    const size_type size = indices.size();
    for (size_type i = 0; i < size; ++i)
      {
        Assert(
          numbers::is_finite(values[i]),
          ExcMessage(
            "The given value is not finite but either infinite or Not A Number (NaN)"));
        this->operator()(indices[i]) += values[indices[i]];
      }
  }



  template <typename Number>
  template <typename Number2>
  inline void
  ReadWriteVector<Number>::add(const size_type  n_indices,
                               const size_type *indices,
                               const Number2 *  values_to_add)
  {
    for (size_type i = 0; i < n_indices; ++i)
      {
        Assert(
          numbers::is_finite(values[i]),
          ExcMessage(
            "The given value is not finite but either infinite or Not A Number (NaN)"));
        this->operator()(indices[i]) += values_to_add[i];
      }
  }



  template <typename Number>
  template <typename Functor>
  inline ReadWriteVector<Number>::FunctorTemplate<Functor>::FunctorTemplate(
    ReadWriteVector<Number> &parent,
    const Functor &          functor)
    : parent(parent)
    , functor(functor)
  {}



  template <typename Number>
  template <typename Functor>
  void
  ReadWriteVector<Number>::FunctorTemplate<Functor>::
  operator()(const size_type begin, const size_type end)
  {
    for (size_type i = begin; i < end; ++i)
      functor(parent.values[i]);
  }

#endif // ifndef DOXYGEN

} // end of namespace LinearAlgebra



template <typename Number>
inline void
swap(LinearAlgebra::ReadWriteVector<Number> &u,
     LinearAlgebra::ReadWriteVector<Number> &v)
{
  u.swap(v);
}


DEAL_II_NAMESPACE_CLOSE

#endif