rol_adaptor.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) 2017 - 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_trilinos_rol_adaptor_h
#define dealii_trilinos_rol_adaptor_h
 
#include <deal.II/base/config.h>
 
#ifdef DEAL_II_TRILINOS_WITH_ROL
#  include <deal.II/base/exceptions.h>
#  include <deal.II/base/index_set.h>
#  include <deal.II/base/mpi.h>
#  include <deal.II/base/types.h>
 
#  include <ROL_Vector.hpp>
 
#  include <limits>
#  include <tuple>
#  include <type_traits>
 
#endif // DEAL_II_TRILINOS_WITH_ROL
 
DEAL_II_NAMESPACE_OPEN
 
#ifdef DEAL_II_TRILINOS_WITH_ROL
namespace TrilinosWrappers
{
  template <typename VectorType>
  class ROLAdaptor : public ROL::Vector<typename VectorType::value_type>
  {
    using size_type = typename VectorType::size_type;
 
    using value_type = typename VectorType::value_type;
 
    using real_type = typename VectorType::real_type;
 
    static_assert(std::is_convertible_v<real_type, value_type>,
                  "The real_type of the current VectorType is not "
                  "convertible to the value_type.");
 
  private:
    ROL::Ptr<VectorType> vector_ptr;
 
    IndexSet optimization_space;
 
    ::types::global_dof_index global_opt_dimension;
 
    ::types::global_dof_index local_opt_start_index;
 
  public:
    ROLAdaptor(const ROL::Ptr<VectorType> &vector_ptr);
 
    ROLAdaptor(const ROL::Ptr<VectorType> &vector_ptr,
               const IndexSet             &optimization_space);
 
    ROL::Ptr<VectorType>
    getVector();
 
    ROL::Ptr<const VectorType>
    getVector() const;
 
    int
    dimension() const override;
 
    void
    set(const ROL::Vector<value_type> &rol_vector) override;
 
    void
    plus(const ROL::Vector<value_type> &rol_vector) override;
 
    void
    axpy(const value_type               alpha,
         const ROL::Vector<value_type> &rol_vector) override;
 
    void
    scale(const value_type alpha) override;
 
    value_type
    dot(const ROL::Vector<value_type> &rol_vector) const override;
 
    value_type
    norm() const override;
 
    ROL::Ptr<ROL::Vector<value_type>>
    clone() const override;
 
    ROL::Ptr<ROL::Vector<value_type>>
    basis(const int i) const override;
 
    void
    applyUnary(const ROL::Elementwise::UnaryFunction<value_type> &f) override;
 
    void
    applyBinary(const ROL::Elementwise::BinaryFunction<value_type> &f,
                const ROL::Vector<value_type> &rol_vector) override;
 
    value_type
    reduce(const ROL::Elementwise::ReductionOp<value_type> &r) const override;
 
    void
    print(std::ostream &outStream) const override;
  };
 
 
  /*------------------------------member definitions--------------------------*/
#  ifndef DOXYGEN
 
 
  template <typename VectorType>
  ROLAdaptor<VectorType>::ROLAdaptor(const ROL::Ptr<VectorType> &vector_ptr)
    : ROLAdaptor(vector_ptr, vector_ptr->locally_owned_elements())
  {}
 
 
 
  template <typename VectorType>
  ROLAdaptor<VectorType>::ROLAdaptor(const ROL::Ptr<VectorType> &vector_ptr,
                                     const IndexSet &optimization_space)
    : vector_ptr(vector_ptr)
    , optimization_space(optimization_space)
  {
    Assert(
      optimization_space.is_subset_of(vector_ptr->locally_owned_elements()),
      ExcMessage(
        "Provided IndexSet needs to be a subset of locally owned indices."));
 
    std::tie(local_opt_start_index, global_opt_dimension) =
      Utilities::MPI::partial_and_total_sum(optimization_space.n_elements(),
                                            vector_ptr->get_mpi_communicator());
 
    // The return type of ROL::Vector::dimension() has to be int,
    // so we check this requirement of ROL here.
    Assert(global_opt_dimension < static_cast<::types::global_dof_index>(
                                    std::numeric_limits<int>::max()),
           ExcMessage("The number of elements to optimize is greater than the "
                      "largest value of type int."));
  }
 
 
 
  template <typename VectorType>
  ROL::Ptr<VectorType>
  ROLAdaptor<VectorType>::getVector()
  {
    return vector_ptr;
  }
 
 
 
  template <typename VectorType>
  ROL::Ptr<const VectorType>
  ROLAdaptor<VectorType>::getVector() const
  {
    return vector_ptr;
  }
 
 
 
  template <typename VectorType>
  void
  ROLAdaptor<VectorType>::set(const ROL::Vector<value_type> &other_)
  {
    Assert(dynamic_cast<const ROLAdaptor *>(&other_) != nullptr,
           ExcInternalError());
    const ROLAdaptor &other = dynamic_cast<const ROLAdaptor &>(other_);
 
    // Perform a deep copy of the vector pointed to by vector_ptr.
    (*vector_ptr) = *(other.getVector());
 
    optimization_space    = other.optimization_space;
    global_opt_dimension  = other.global_opt_dimension;
    local_opt_start_index = other.local_opt_start_index;
  }
 
 
 
  template <typename VectorType>
  void
  ROLAdaptor<VectorType>::plus(const ROL::Vector<value_type> &other_)
  {
    Assert(dynamic_cast<const ROLAdaptor *>(&other_) != nullptr,
           ExcInternalError());
    const ROLAdaptor &other = dynamic_cast<const ROLAdaptor &>(other_);
 
    Assert(vector_ptr->has_ghost_elements() == false, ExcGhostsPresent());
    Assert(optimization_space.size() == vector_ptr->size(),
           ExcMessage("Optimization space is out-of-sync. "
                      "Please create a new wrapper."));
    Assert(optimization_space == other.optimization_space,
           ExcMessage("Optimization spaces of vectors do not match."));
 
    for (const auto i : optimization_space)
      (*vector_ptr)[i] += (*other.getVector())[i];
 
    vector_ptr->compress(VectorOperation::add);
  }
 
 
 
  template <typename VectorType>
  void
  ROLAdaptor<VectorType>::axpy(const value_type               alpha,
                               const ROL::Vector<value_type> &other_)
  {
    Assert(dynamic_cast<const ROLAdaptor *>(&other_) != nullptr,
           ExcInternalError());
    const ROLAdaptor &other = dynamic_cast<const ROLAdaptor &>(other_);
 
    Assert(vector_ptr->has_ghost_elements() == false, ExcGhostsPresent());
    Assert(optimization_space.size() == vector_ptr->size(),
           ExcMessage("Optimization space is out-of-sync. "
                      "Please create a new wrapper."));
    Assert(optimization_space == other.optimization_space,
           ExcMessage("Optimization spaces of vectors do not match."));
 
    for (const auto i : optimization_space)
      (*vector_ptr)[i] += alpha * (*other.getVector())[i];
 
    vector_ptr->compress(VectorOperation::add);
  }
 
 
 
  template <typename VectorType>
  int
  ROLAdaptor<VectorType>::dimension() const
  {
    Assert(optimization_space.size() == vector_ptr->size(),
           ExcMessage("Optimization space is out-of-sync. "
                      "Please create a new wrapper."));
 
    return static_cast<int>(global_opt_dimension);
  }
 
 
 
  template <typename VectorType>
  void
  ROLAdaptor<VectorType>::scale(const value_type alpha)
  {
    Assert(vector_ptr->has_ghost_elements() == false, ExcGhostsPresent());
    Assert(optimization_space.size() == vector_ptr->size(),
           ExcMessage("Optimization space is out-of-sync. "
                      "Please create a new wrapper."));
 
    for (const auto i : optimization_space)
      (*vector_ptr)[i] *= alpha;
 
    vector_ptr->compress(VectorOperation::insert);
  }
 
 
 
  template <typename VectorType>
  typename VectorType::value_type
  ROLAdaptor<VectorType>::dot(const ROL::Vector<value_type> &other_) const
  {
    Assert(dynamic_cast<const ROLAdaptor *>(&other_) != nullptr,
           ExcInternalError());
    const ROLAdaptor &other = dynamic_cast<const ROLAdaptor &>(other_);
 
    Assert(optimization_space.size() == vector_ptr->size(),
           ExcMessage("Optimization space is out-of-sync. "
                      "Please create a new wrapper."));
    Assert(optimization_space == other.optimization_space,
           ExcMessage("Optimization spaces of vectors do not match."));
 
    value_type dot(0);
    for (const auto i : optimization_space)
      dot += (*vector_ptr)[i] * (*other.getVector())[i];
 
    return Utilities::MPI::sum<value_type>(dot,
                                           vector_ptr->get_mpi_communicator());
  }
 
 
 
  template <typename VectorType>
  typename VectorType::value_type
  ROLAdaptor<VectorType>::norm() const
  {
    return std::sqrt(this->dot(*this));
  }
 
 
 
  template <typename VectorType>
  ROL::Ptr<ROL::Vector<typename VectorType::value_type>>
  ROLAdaptor<VectorType>::clone() const
  {
    // create new vector with same size as wrapped one
    ROL::Ptr<VectorType> clone_ptr = ROL::makePtr<VectorType>();
    clone_ptr->reinit(*vector_ptr, false);
 
    return ROL::makePtr<ROLAdaptor>(clone_ptr, optimization_space);
    // TODO: also somehow pass the partial sums etc?
  }
 
 
 
  template <typename VectorType>
  ROL::Ptr<ROL::Vector<typename VectorType::value_type>>
  ROLAdaptor<VectorType>::basis(const int global_opt_index) const
  {
    Assert(optimization_space.size() == vector_ptr->size(),
           ExcMessage("Optimization space is out-of-sync. "
                      "Please create a new wrapper."));
    AssertIndexRange(global_opt_index, global_opt_dimension);
 
    // clone the internal vector as basis
    ROL::Ptr<VectorType> basis_ptr = ROL::makePtr<VectorType>();
    basis_ptr->reinit(*vector_ptr, false);
 
    // check whether the i-th basis corresponds to one of the
    // locally owned elements of the wrapped vector
    if ((global_opt_index >= local_opt_start_index) &&
        (global_opt_index <
         local_opt_start_index + optimization_space.n_elements()))
      {
        const ::types::global_dof_index local_opt_index =
          global_opt_index - local_opt_start_index;
 
        // global dof index corresponding to i-th basis in optimization
        const ::types::global_dof_index global_dof_index =
          optimization_space.nth_index_in_set(local_opt_index);
 
        (*basis_ptr)[global_dof_index] = 1.;
      }
 
    basis_ptr->compress(VectorOperation::insert);
 
    return ROL::makePtr<ROLAdaptor>(basis_ptr, optimization_space);
    // TODO: also somehow pass the partial sums etc?
  }
 
 
 
  template <typename VectorType>
  void
  ROLAdaptor<VectorType>::applyUnary(
    const ROL::Elementwise::UnaryFunction<value_type> &f)
  {
    Assert(vector_ptr->has_ghost_elements() == false, ExcGhostsPresent());
    Assert(optimization_space.size() == vector_ptr->size(),
           ExcMessage("Optimization space is out-of-sync. "
                      "Please create a new wrapper."));
 
    for (const auto i : optimization_space)
      (*vector_ptr)[i] = f.apply((*vector_ptr)[i]);
 
    vector_ptr->compress(VectorOperation::insert);
  }
 
 
 
  template <typename VectorType>
  void
  ROLAdaptor<VectorType>::applyBinary(
    const ROL::Elementwise::BinaryFunction<value_type> &f,
    const ROL::Vector<value_type>                      &other_)
  {
    Assert(dynamic_cast<const ROLAdaptor *>(&other_) != nullptr,
           ExcInternalError());
    const ROLAdaptor &other = dynamic_cast<const ROLAdaptor &>(other_);
 
    Assert(vector_ptr->has_ghost_elements() == false, ExcGhostsPresent());
    Assert(optimization_space.size() == vector_ptr->size(),
           ExcMessage("Optimization space is out-of-sync. "
                      "Please create a new wrapper."));
    Assert(optimization_space == other.optimization_space,
           ExcMessage("Optimization spaces of vectors do not match."));
 
    for (const auto i : optimization_space)
      (*vector_ptr)[i] = f.apply((*vector_ptr)[i], (*other.getVector())[i]);
 
    vector_ptr->compress(VectorOperation::insert);
  }
 
 
 
  template <typename VectorType>
  typename VectorType::value_type
  ROLAdaptor<VectorType>::reduce(
    const ROL::Elementwise::ReductionOp<value_type> &r) const
  {
    Assert(optimization_space.size() == vector_ptr->size(),
           ExcMessage("Optimization space is out-of-sync. "
                      "Please create a new wrapper."));
 
    value_type result = r.initialValue();
 
    // local reduction
    for (const auto i : optimization_space)
      r.reduce((*vector_ptr)[i], result);
 
    // global reduction
    const auto combiner = [&r](const value_type a,
                               const value_type b) -> value_type {
      auto out = b;
      r.reduce(a, out);
      return out;
    };
 
    return Utilities::MPI::all_reduce<value_type>(
      result, vector_ptr->get_mpi_communicator(), combiner);
  }
 
 
 
  template <typename VectorType>
  void
  ROLAdaptor<VectorType>::print(std::ostream &outStream) const
  {
    vector_ptr->print(outStream);
  }
 
 
#  endif // DOXYGEN
 
 
} // namespace TrilinosWrappers
 
#endif // DEAL_II_TRILINOS_WITH_ROL
 
DEAL_II_NAMESPACE_CLOSE
 
#endif // dealii_trilinos_rol_adaptor_h