vector_adaptor.h

//-----------------------------------------------------------
//
//    Copyright (C) 2017 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_optimization_rol_vector_adaptor_h
#define dealii_optimization_rol_vector_adaptor_h

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

#ifdef DEAL_II_TRILINOS_WITH_ROL
#  include <deal.II/base/exceptions.h>

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

#  include <ROL_Vector.hpp>

#  include <type_traits>


DEAL_II_NAMESPACE_OPEN


using namespace dealii;


namespace Rol
{
  template <typename VectorType>
  class VectorAdaptor : 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<real_type, value_type>::value,
                  "The real_type of the current VectorType is not "
                  "convertible to the value_type.");

  private:
    Teuchos::RCP<VectorType> vector_ptr;

  public:
    VectorAdaptor(const Teuchos::RCP<VectorType> &vector_ptr);

    Teuchos::RCP<VectorType>
    getVector();

    Teuchos::RCP<const VectorType>
    getVector() const;

    int
    dimension() const;

    void
    set(const ROL::Vector<value_type> &rol_vector);

    void
    plus(const ROL::Vector<value_type> &rol_vector);

    void
    axpy(const value_type alpha, const ROL::Vector<value_type> &rol_vector);

    void
    scale(const value_type alpha);

    value_type
    dot(const ROL::Vector<value_type> &rol_vector) const;

    value_type
    norm() const;

    Teuchos::RCP<ROL::Vector<value_type>>
    clone() const;

    Teuchos::RCP<ROL::Vector<value_type>>
    basis(const int i) const;

    void
    applyUnary(const ROL::Elementwise::UnaryFunction<value_type> &f);

    void
    applyBinary(const ROL::Elementwise::BinaryFunction<value_type> &f,
                const ROL::Vector<value_type> &                     rol_vector);

    value_type
    reduce(const ROL::Elementwise::ReductionOp<value_type> &r) const;

    void
    print(std::ostream &outStream) const;
  };


  /*------------------------------member definitions--------------------------*/
#  ifndef DOXYGEN


  template <typename VectorType>
  VectorAdaptor<VectorType>::VectorAdaptor(
    const Teuchos::RCP<VectorType> &vector_ptr)
    : vector_ptr(vector_ptr)
  {}



  template <typename VectorType>
  Teuchos::RCP<VectorType>
  VectorAdaptor<VectorType>::getVector()
  {
    return vector_ptr;
  }



  template <typename VectorType>
  Teuchos::RCP<const VectorType>
  VectorAdaptor<VectorType>::getVector() const
  {
    return vector_ptr;
  }



  template <typename VectorType>
  void
  VectorAdaptor<VectorType>::set(const ROL::Vector<value_type> &rol_vector)
  {
    const VectorAdaptor &vector_adaptor =
      Teuchos::dyn_cast<const VectorAdaptor>(rol_vector);

    (*vector_ptr) = *(vector_adaptor.getVector());
  }



  template <typename VectorType>
  void
  VectorAdaptor<VectorType>::plus(const ROL::Vector<value_type> &rol_vector)
  {
    Assert(this->dimension() == rol_vector.dimension(),
           ExcDimensionMismatch(this->dimension(), rol_vector.dimension()));

    const VectorAdaptor &vector_adaptor =
      Teuchos::dyn_cast<const VectorAdaptor>(rol_vector);

    *vector_ptr += *(vector_adaptor.getVector());
  }



  template <typename VectorType>
  void
  VectorAdaptor<VectorType>::axpy(const value_type               alpha,
                                  const ROL::Vector<value_type> &rol_vector)
  {
    Assert(this->dimension() == rol_vector.dimension(),
           ExcDimensionMismatch(this->dimension(), rol_vector.dimension()));

    const VectorAdaptor &vector_adaptor =
      Teuchos::dyn_cast<const VectorAdaptor>(rol_vector);

    vector_ptr->add(alpha, *(vector_adaptor.getVector()));
  }



  template <typename VectorType>
  int
  VectorAdaptor<VectorType>::dimension() const
  {
    Assert(vector_ptr->size() < std::numeric_limits<int>::max(),
           ExcMessage("The size of the vector being used is greater than "
                      "largest value of type int."));
    return static_cast<int>(vector_ptr->size());
  }



  template <typename VectorType>
  void
  VectorAdaptor<VectorType>::scale(const value_type alpha)
  {
    (*vector_ptr) *= alpha;
  }



  template <typename VectorType>
  typename VectorType::value_type
  VectorAdaptor<VectorType>::dot(
    const ROL::Vector<value_type> &rol_vector) const
  {
    Assert(this->dimension() == rol_vector.dimension(),
           ExcDimensionMismatch(this->dimension(), rol_vector.dimension()));

    const VectorAdaptor &vector_adaptor =
      Teuchos::dyn_cast<const VectorAdaptor>(rol_vector);

    return (*vector_ptr) * (*vector_adaptor.getVector());
  }



  template <typename VectorType>
  typename VectorType::value_type
  VectorAdaptor<VectorType>::norm() const
  {
    return vector_ptr->l2_norm();
  }



  template <typename VectorType>
  Teuchos::RCP<ROL::Vector<typename VectorType::value_type>>
  VectorAdaptor<VectorType>::clone() const
  {
    Teuchos::RCP<VectorType> vec_ptr = Teuchos::rcp(new VectorType);
    (*vec_ptr)                       = (*vector_ptr);

    return Teuchos::rcp(new VectorAdaptor(vec_ptr));
  }



  template <typename VectorType>
  Teuchos::RCP<ROL::Vector<typename VectorType::value_type>>
  VectorAdaptor<VectorType>::basis(const int i) const
  {
    Teuchos::RCP<VectorType> vec_ptr = Teuchos::rcp(new VectorType);

    // Zero all the entries in dealii vector.
    vec_ptr->reinit(*vector_ptr, false);

    if (vector_ptr->locally_owned_elements().is_element(i))
      vec_ptr->operator[](i) = 1.;

    vec_ptr->compress(VectorOperation::insert);

    Teuchos::RCP<VectorAdaptor> e = Teuchos::rcp(new VectorAdaptor(vec_ptr));

    return e;
  }



  template <typename VectorType>
  void
  VectorAdaptor<VectorType>::applyUnary(
    const ROL::Elementwise::UnaryFunction<value_type> &f)
  {
    const typename VectorType::iterator vend = vector_ptr->end();

    for (typename VectorType::iterator iterator = vector_ptr->begin();
         iterator != vend;
         iterator++)
      *iterator = f.apply(*iterator);

    vector_ptr->compress(VectorOperation::insert);
  }



  template <typename VectorType>
  void
  VectorAdaptor<VectorType>::applyBinary(
    const ROL::Elementwise::BinaryFunction<value_type> &f,
    const ROL::Vector<value_type> &                     rol_vector)
  {
    Assert(this->dimension() == rol_vector.dimension(),
           ExcDimensionMismatch(this->dimension(), rol_vector.dimension()));

    const VectorAdaptor &vector_adaptor =
      Teuchos::dyn_cast<const VectorAdaptor>(rol_vector);

    const VectorType &given_rol_vector = *(vector_adaptor.getVector());

    const typename VectorType::iterator       vend   = vector_ptr->end();
    const typename VectorType::const_iterator rolend = given_rol_vector.end();

    typename VectorType::const_iterator r_iterator = given_rol_vector.begin();
    for (typename VectorType::iterator l_iterator = vector_ptr->begin();
         l_iterator != vend && r_iterator != rolend;
         l_iterator++, r_iterator++)
      *l_iterator = f.apply(*l_iterator, *r_iterator);

    vector_ptr->compress(VectorOperation::insert);
  }



  template <typename VectorType>
  typename VectorType::value_type
  VectorAdaptor<VectorType>::reduce(
    const ROL::Elementwise::ReductionOp<value_type> &r) const
  {
    typename VectorType::value_type result = r.initialValue();

    const typename VectorType::iterator vend = vector_ptr->end();

    for (typename VectorType::iterator iterator = vector_ptr->begin();
         iterator != vend;
         iterator++)
      r.reduce(*iterator, result);
    // Parallel reduction among processes is handled internally.

    return result;
  }



  template <typename VectorType>
  void
  VectorAdaptor<VectorType>::print(std::ostream &outStream) const
  {
    vector_ptr->print(outStream);
  }


#  endif // DOXYGEN


} // namespace Rol


DEAL_II_NAMESPACE_CLOSE


#endif // DEAL_II_TRILINOS_WITH_ROL

#endif // dealii_optimization_rol_vector_adaptor_h