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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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/mpi.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/types.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/lac/vector_operation.h>

#include <cstring>
#include <cstdlib>
#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;
  }
}

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

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

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

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

    ReadWriteVector ();

    ReadWriteVector (const ReadWriteVector<Number> &in_vector);

    explicit ReadWriteVector (const size_type size);

    explicit ReadWriteVector (const IndexSet &locally_stored_indices);

    ~ReadWriteVector () = 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