psblas_vector.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2019 - 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_psblas_vector_h
#define dealii_psblas_vector_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/index_set.h>
#include <deal.II/base/types.h>
 
#include <deal.II/lac/read_vector.h>
#include <deal.II/lac/vector_operation.h>
#include <deal.II/lac/vector_type_traits.h>
 
#include <cstddef>
#include <memory>
 
#ifdef DEAL_II_WITH_PSBLAS
 
#  include <deal.II/lac/psblas_common.h>
 
#endif // DEAL_II_WITH_PSBLAS
 
DEAL_II_NAMESPACE_OPEN
 
#ifdef DEAL_II_WITH_PSBLAS
namespace PSCToolkitWrappers
{

  class Vector : public ReadVector<double>
  {
  private:
    class VectorReference
    {
    private:
      using size_type = types::global_dof_index;
 
      using value_type = double;

      VectorReference(Vector &vector, const size_type index)
        : vector(vector)
        , index(index)
      {}
 
    public:
      const VectorReference &
      operator=(const value_type &s) const
      {
        Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
        Assert(
          vector.owned_elements.is_element(index),
          ExcMessage(
            "You are trying to write to an element of the vector that is not "
            "locally owned. This is not allowed for the current interface to"
            " PSBLAS vectors."));
 
        // Make sure the operation is consistent with the last one
        Assert(vector.last_action == VectorOperation::insert ||
                 vector.last_action == VectorOperation::unknown,
               ExcWrongMode(VectorOperation::insert, vector.last_action));
 
        std::vector<size_type>  idx{index};
        std::vector<value_type> value{s};
        vector.set(idx, value);
        vector.last_action = VectorOperation::insert;
        return *this;
      }

      const VectorReference &
      operator+=(const value_type &s) const
      {
        Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
        Assert(vector.last_action == VectorOperation::add ||
                 vector.last_action == VectorOperation::unknown,
               ExcWrongMode(VectorOperation::add, vector.last_action));
 
        vector.last_action = VectorOperation::add;
 
        // First check for early return
        if (s == 0.)
          return *this;
 
        std::vector<size_type>  idx{index};
        std::vector<value_type> value{s};
        vector.add(idx, value);
        return *this;
      }

      const VectorReference &
      operator-=(const value_type &s) const
      {
        Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
        Assert(vector.last_action == VectorOperation::add ||
                 vector.last_action == VectorOperation::unknown,
               ExcWrongMode(VectorOperation::add, vector.last_action));
 
        vector.last_action = VectorOperation::add;
 
        // First check for early return
        if (s == 0.)
          return *this;
 
        std::vector<size_type>  idx{index};
        std::vector<value_type> value{-s};
        vector.add(idx, value);
        return *this;
      }

      const VectorReference &
      operator*=(const value_type &s) const
      {
        Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
        Assert((vector.last_action == VectorOperation::insert) ||
                 (vector.last_action == VectorOperation::unknown),
               ExcWrongMode(VectorOperation::insert, vector.last_action));
 
        vector.last_action = VectorOperation::insert;
        if (s == 1.)
          return *this;
 
        std::vector<size_type>  idx{index};
        value_type              new_value = static_cast<value_type>(*this) * s;
        std::vector<value_type> value{new_value};
        vector.set(idx, value);
 
        return *this;
      }

      const VectorReference &
      operator/=(const value_type &s) const
      {
        Assert(!vector.has_ghost_elements(), ExcGhostsPresent());
 
        Assert((vector.last_action == VectorOperation::insert) ||
                 (vector.last_action == VectorOperation::unknown),
               ExcWrongMode(VectorOperation::insert, vector.last_action));
        std::vector<size_type>  idx{index};
        value_type              new_value = static_cast<value_type>(*this) / s;
        std::vector<value_type> value{new_value};
        vector.set(idx, value);
        vector.last_action = VectorOperation::insert;
        return *this;
      }
 
      /*
       * Convert the reference to an actual value, i.e. return the value of
       * the referenced element of the vector.
       */
      operator value_type() const
      {
        AssertIndexRange(index, vector.size());
        if (vector.ghosted)
          {
            AssertThrow(vector.ghost_indices.is_element(index) ||
                          vector.owned_elements.is_element(index),
                        ExcMessage(
                          "You are trying to access an element of a vector "
                          "that is neither a locally owned element nor a "
                          "ghost element of the vector."));
          }
        else
          {
            AssertThrow(
              vector.owned_elements.is_element(index),
              ExcAccessToNonlocalElement(index,
                                         *vector.owned_elements.begin(),
                                         (*vector.owned_elements.begin() +
                                          vector.locally_owned_size())));
          }
        return psb_c_dgetelem(vector.psblas_vector,
                              index,
                              vector.psblas_descriptor.get());
      };
 
    private:
      Vector &vector;
 
      const size_type index;
 
      friend class Vector;
    };
 
  public:
    using size_type = ::types::global_dof_index;
 
    using value_type = double;

    Vector();

    Vector(const Vector &);

    explicit Vector(const IndexSet &local_partitioning,
                    const MPI_Comm  communicator);

    Vector(const IndexSet &local_partitioning,
           const IndexSet &ghost_indices,
           const MPI_Comm  communicator);

    ~Vector();

    void
    reinit(const IndexSet &local_partitioning,
           const MPI_Comm  communicator,
           const bool      omit_zeroing_entries = false);

    void
    reinit(const Vector &v, const bool omit_zeroing_entries = false);

    void
    reinit(const IndexSet &local_partitioning,
           const IndexSet &ghost_indices,
           const MPI_Comm  communicator);
 

    Vector &
    operator=(const Vector &v);

    size_type
    size() const override;

    virtual void
    extract_subvector_to(
      const ArrayView<const types::global_dof_index> &indices,
      const ArrayView<value_type>                    &elements) const override;

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

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

    size_type
    locally_owned_size() const;

    void
    set(const std::vector<size_type>  &indices,
        const std::vector<value_type> &values);

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

    void
    add(const value_type s, const Vector &V);

    void
    add(const value_type s);

    void
    scale(const Vector &v);

    value_type
    add_and_dot(const value_type a, const Vector &v, const Vector &W);
 
    /*
     * Scaling and vector addition, i.e.  <tt>*this = s*(*this)+V</tt>.
     */
    void
    sadd(const value_type s, const Vector &V);
 
    /*
     * Scaling and vector addition, i.e.  <tt>*this = s*(*this)+a*V</tt>.
     */
    void
    sadd(const value_type s, const value_type a, const Vector &V);
 
    /*
     * Assignment *this = a*V.
     */
    void
    equ(const value_type a, const Vector &v);

    value_type
    operator()(const size_type index) const;

    VectorReference
    operator()(const size_type index);

    value_type
    operator[](const size_type index) const;

    VectorReference
    operator[](const size_type index);

    value_type
    operator*(const Vector &v) const;

    Vector &
    operator-=(const Vector &v);

    Vector &
    operator+=(const Vector &v);

    Vector &
    operator=(const value_type s);

    const IndexSet &
    locally_owned_elements() const;

    const IndexSet &
    ghost_elements() const;

    bool
    has_ghost_elements() const;

    void
    update_ghost_values() const;

    void
    swap(Vector &v);

    void
    compress(const VectorOperation::values operation);

    value_type *
    begin();

    const value_type *
    begin() const;

    value_type *
    end();

    const value_type *
    end() const;
 
    MPI_Comm
    get_mpi_communicator() const;

    psb_c_descriptor *
    get_psblas_descriptor() const;

    psb_c_dvector *
    get_psblas_vector() const;

    void
    clear();

    Vector::value_type
    linfty_norm() const;

    Vector::value_type
    l1_norm() const;

    Vector::value_type
    l2_norm() const;

    Vector::value_type
    mean_value() const;

    bool
    all_zero() const;

    std::size_t
    memory_consumption() const;
 
  private:
    /*
     * Pointer to the underlying PSBLAS vector.
     */
    psb_c_dvector *psblas_vector;
 
    /*
     * Pointer to PSBLAS context.
     */
    psb_c_ctxt *psblas_context;
 
    /*
     * Shared pointer to the PSBLAS descriptor.
     */
    std::shared_ptr<psb_c_descriptor> psblas_descriptor;

    MPI_Comm communicator;

    IndexSet owned_elements;

    IndexSet ghost_indices;

    bool ghosted;

    internal::State state;

    VectorOperation::values last_action;
 
    friend class SparseMatrix;
  };
 
 
  /* ----------------------------- Inline functions ---------------- */
 
 
  inline PSCToolkitWrappers::Vector::value_type *
  PSCToolkitWrappers::Vector::begin()
  {
    return psb_c_dvect_f_get_pnt(psblas_vector);
  }
 
 
 
  inline const PSCToolkitWrappers::Vector::value_type *
  PSCToolkitWrappers::Vector::begin() const
  {
    return psb_c_dvect_f_get_pnt(psblas_vector);
  }
 
 
 
  inline PSCToolkitWrappers::Vector::value_type *
  PSCToolkitWrappers::Vector::end()
  {
    return psb_c_dvect_f_get_pnt(psblas_vector) + locally_owned_size();
  }
 
 
 
  inline const PSCToolkitWrappers::Vector::value_type *
  PSCToolkitWrappers::Vector::end() const
  {
    return psb_c_dvect_f_get_pnt(psblas_vector) + locally_owned_size();
  }
 
 
 
  inline Vector::size_type
  Vector::size() const
  {
    return owned_elements.size();
  }
 
 
 
  inline const IndexSet &
  Vector::locally_owned_elements() const
  {
    return owned_elements;
  }
 
 
 
  inline const IndexSet &
  Vector::ghost_elements() const
  {
    return ghost_indices;
  }
 
 
 
  inline bool
  Vector::has_ghost_elements() const
  {
    return ghosted;
  }
 
 
 
  inline Vector::value_type
  Vector::operator()(const Vector::size_type index) const
  {
    return psb_c_dgetelem(psblas_vector, index, psblas_descriptor.get());
  }
 
 
 
  inline Vector::VectorReference
  Vector::operator()(const size_type index)
  {
    return VectorReference(*this, index);
  }
 
 
 
  inline Vector::value_type
  Vector::operator[](const Vector::size_type index) const
  {
    return operator()(index);
  }
 
 
 
  inline Vector::VectorReference
  Vector::operator[](const size_type index)
  {
    return operator()(index);
  }
 
 
 
  inline void
  Vector::extract_subvector_to(
    const ArrayView<const types::global_dof_index> &indices,
    const ArrayView<double>                        &elements) const
  {
    AssertDimension(indices.size(), elements.size());
    extract_subvector_to(indices.begin(), indices.end(), elements.begin());
  }
 
 
 
  inline void
  Vector::extract_subvector_to(const std::vector<size_type>    &indices,
                               std::vector<Vector::value_type> &values) const
  {
    AssertDimension(indices.size(), values.size());
    extract_subvector_to(indices.begin(), indices.end(), values.begin());
  }
 
 
 
  template <typename ForwardIterator, typename OutputIterator>
  inline void
  Vector::extract_subvector_to(ForwardIterator indices_begin,
                               ForwardIterator indices_end,
                               OutputIterator  output) const
  {
    if (indices_begin == indices_end)
      return;
 
    if (ghosted)
      {
        types::global_dof_index begin = *owned_elements.begin();
        types::global_dof_index end   = begin + owned_elements.n_elements();
 
        auto input = indices_begin;
        while (input != indices_end)
          {
            const auto index = static_cast<psb_l_t>(*input);
            AssertThrow((index >= begin && index < end) ||
                          ghost_indices.is_element(index),
                        ExcMessage(
                          "You are trying to access an element of a vector "
                          "that is neither a locally owned element nor a "
                          "ghost element of the vector."));
            *output =
              psb_c_dgetelem(psblas_vector, index, psblas_descriptor.get());
 
 
            ++input;
            ++output;
          }
      }
    else
      {
        // no ghost elements, so we can
        // just access the local
        // elements directly
        while (indices_begin != indices_end)
          {
            const size_type index = *indices_begin;
            Assert(owned_elements.is_element(index),
                   ExcMessage("You are accessing elements of a vector without "
                              "ghost elements that are not actually owned by "
                              "this vector. A typical case where this may "
                              "happen is if you are passing a non-ghosted "
                              "(completely distributed) vector to a function "
                              "that expects a vector that stores ghost "
                              "elements for all locally relevant or locally "
                              "active vector entries."));
 
            *output =
              psb_c_dgetelem(psblas_vector, index, psblas_descriptor.get());
 
            ++indices_begin;
            ++output;
          }
      }
  }
 
 
} // namespace PSCToolkitWrappers

template <>
struct is_serial_vector<PSCToolkitWrappers::Vector> : std::false_type
{};
 
#endif // DEAL_II_WITH_PSBLAS
 
DEAL_II_NAMESPACE_CLOSE
#endif