From 691484526704a23ee6f589c81d51f16ef93c80e2 Mon Sep 17 00:00:00 2001 From: vishalkenchan Date: Fri, 10 Nov 2017 10:15:12 +0100 Subject: [PATCH] add a vector adaptor to use Trilinos/ROL library; added cmake variable DEAL_II_TRILINOS_WITH_ROL --- cmake/configure/configure_2_trilinos.cmake | 12 + doc/news/changes/major/20171102VishalBoddu | 7 + include/deal.II/base/config.h.in | 3 + .../deal.II/optimization/rol/vector_adaptor.h | 489 ++++++++++++++++++ source/CMakeLists.txt | 3 +- source/optimization/rol/CMakeLists.txt | 29 ++ 6 files changed, 542 insertions(+), 1 deletion(-) create mode 100644 doc/news/changes/major/20171102VishalBoddu create mode 100644 include/deal.II/optimization/rol/vector_adaptor.h create mode 100644 source/optimization/rol/CMakeLists.txt diff --git a/cmake/configure/configure_2_trilinos.cmake b/cmake/configure/configure_2_trilinos.cmake index deae1c811b..d157125011 100644 --- a/cmake/configure/configure_2_trilinos.cmake +++ b/cmake/configure/configure_2_trilinos.cmake @@ -54,6 +54,18 @@ MACRO(FEATURE_TRILINOS_FIND_EXTERNAL var) ENDIF() ENDFOREACH() + FOREACH(_optional_module + ROL + ) + ITEM_MATCHES(_module_found ${_optional_module} ${Trilinos_PACKAGE_LIST}) + IF(_module_found) + MESSAGE(STATUS "Found ${_optional_module}") + SET(DEAL_II_TRILINOS_WITH_${_optional_module} ON) + ELSE() + MESSAGE(STATUS "Module ${_optional_module} not found!") + ENDIF() + ENDFOREACH() + IF((TRILINOS_VERSION_MAJOR EQUAL 11 AND NOT (TRILINOS_VERSION_MINOR LESS 14)) OR diff --git a/doc/news/changes/major/20171102VishalBoddu b/doc/news/changes/major/20171102VishalBoddu new file mode 100644 index 0000000000..c9c5d22de9 --- /dev/null +++ b/doc/news/changes/major/20171102VishalBoddu @@ -0,0 +1,7 @@ +New: The new namespace Rol contains an adaptor class that provides +the implementation of the ROL::Vector interface. +The Trilinos package, Rapid Optimization Library (ROL), can solve unconstrained +and constrained optimization problems, and optimization problems under +uncertainty. +
+(Vishal Boddu 2017/11/02) diff --git a/include/deal.II/base/config.h.in b/include/deal.II/base/config.h.in index fcf1065a62..a19dffbc91 100644 --- a/include/deal.II/base/config.h.in +++ b/include/deal.II/base/config.h.in @@ -153,6 +153,9 @@ #cmakedefine DEAL_II_USE_MT_POSIX #cmakedefine DEAL_II_USE_MT_POSIX_NO_BARRIERS +/* cmake/configure/configure_2_trilinos.cmake */ +#cmakedefine DEAL_II_TRILINOS_WITH_ROL + /* * Depending on the use of threads, we will have to make some variables * volatile. We do this here in a very old-fashioned C-style, but still diff --git a/include/deal.II/optimization/rol/vector_adaptor.h b/include/deal.II/optimization/rol/vector_adaptor.h new file mode 100644 index 0000000000..cc068a2ba5 --- /dev/null +++ b/include/deal.II/optimization/rol/vector_adaptor.h @@ -0,0 +1,489 @@ +//----------------------------------------------------------- +// +// 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 at +// the top level of the deal.II distribution. +// +//--------------------------------------------------------------- + +#ifndef dealii_optimization_rol_vector_adaptor_h +#define dealii_optimization_rol_vector_adaptor_h + +#include + +#ifdef DEAL_II_TRILINOS_WITH_ROL +#include + +#include +#include + +#include + + +DEAL_II_NAMESPACE_OPEN + + +using namespace dealii; + + +/** + * A namespace that provides an interface to the + * + * Rapid Optimization Library (ROL), a Trilinos package. + */ +namespace Rol +{ + + /** + * An adaptor that provides the implementation of the ROL::Vector interface + * for vectors of type VectorType. + * + * This class supports vectors that satisfy the following requirements: + * + * The VectorType should contain the following types. + * ``` + * VectorType::size_type; // The type for size of the vector. + * VectorType::value_type; // The type for elements stored in the vector. + * VectorType::real_type; // The type for real-valued numbers. + * ``` + * + * However, ROL doesn't distinguish VectorAdaptor::value_type from + * VectorAdaptor::real_type. This is due to ROL's assumption that the + * VectorAdaptor::value_type itself is a type for real-valued numbers. + * Therefore, VectorAdaptor supports vectors whose real_type is + * convertible to value_type in the sense that + * std::is_convertible::value yields + * true. + * + * The VectorType should contain the following methods. + * @code + * // Reinitialize the current vector using a given vector's + * // size (and the parallel distribution) without copying + * // the elements. + * VectorType::reinit(const VectorType &, ...); + * + * // Globally add a given vector to the current. + * VectorType::operator+=(const VectorType &); + * + * // Scale all elements by a given scalar. + * VectorType::operator*=(const VectorType::value_type &); + * + * // Perform dot product with a given vector. + * VectorType::operator*=(const VectorType &); + * + * // Scale all elements of the current vector and globally + * // add a given vector to it. + * VectorType::add(const VectorType::value_type, const VectorType &); + * + * // Copies the data of a given vector to the current. + * // Resize the current vector if necessary (MPI safe). + * VectorType::operation=(const VectorType &); + * + * // Return the global size of the current vector. + * VectorType::size(); + * + * // Return L^2 norm of the current vector + * VectorType::l2_norm(); + * + * // Iterator to the start of the (locally owned) element + * // of the current vector. + * VectorType::begin(); + * + * // Iterator to the one past the last (locally owned) + * // element of the current vector. + * VectorType::end(); + * + * // Compress the vector i.e., flush the buffers of the + * // vector object if it has any. + * VectorType::compress(VectorOperation::insert); + * @endcode + * + * @note The current implementation in ROL doesn't support vector sizes above + * the largest value of int type. Some of the vectors in deal.II (@ref Vector) + * may not satisfy the above requirements. + * + * @author Vishal Boddu, 2017 + */ + template + class VectorAdaptor : public ROL::Vector + { + + /** + * A typedef for size type of VectorType. + */ + using size_type = typename VectorType::size_type; + + /** + * A typedef for element type stored in the VectorType. + */ + using value_type = typename VectorType::value_type; + + /** + * A typedef for real-valued numbers. + */ + using real_type = typename VectorType::real_type; + + static_assert (std::is_convertible::value, + "The real_type of the current VectorType is not " + "convertible to the value_type."); + + private: + + /** + * Teuchos smart reference counting pointer to the underlying vector of type + * VectorType. + */ + Teuchos::RCP vector_ptr; + + public: + + /** + * Constructor. + */ + VectorAdaptor (const Teuchos::RCP &vector_ptr); + + /** + * Returns the Teuchos smart reference counting pointer to + * the wrapper vector, #vector_ptr. + */ + Teuchos::RCP getVector (); + + /** + * Returns the Teuchos smart reference counting pointer to const vector. + */ + Teuchos::RCP getVector () const; + + /** + * Return the dimension (global vector size) of the wrapped vector. + */ + int dimension () const; + + /** + * Set the wrapper vector to a given ROL::Vector @p rol_vector by + * overwriting its contents. + * + * If the current wrapper vector has ghost elements, + * then VectorType::operator=(const VectorType&) should still + * be allowed on it. + */ + void set (const ROL::Vector &rol_vector); + + /** + * Perform addition. + */ + void plus (const ROL::Vector &rol_vector); + + /** + * Scale the wrapper vector by @p alpha and add ROL::Vector @p rol_vector + * to it. + */ + void axpy (const value_type alpha, + const ROL::Vector &rol_vector); + + /** + * Scale the wrapper vector. + */ + void scale (const value_type alpha); + + /** + * Return the dot product with a given ROL::Vector @p rol_vector. + */ + value_type dot (const ROL::Vector &rol_vector) const; + + /** + * Return the \f$ L^{2} \f$ norm of the wrapped vector. + * + * The returned type is of VectorAdaptor::value_type so as to maintain + * consistency with ROL::Vector and + * more importantly to not to create an overloaded version namely, + * VectorAdaptor::real_type norm() const; + * if real_type and value_type are not of the same type. + */ + value_type norm() const; + + /** + * Return a clone of the wrapped vector. + */ + Teuchos::RCP> clone() const; + + /** + * Create and return a Teuchos smart reference counting pointer to the basis + * vector corresponding to the @p i \f${}^{th}\f$ element of + * the wrapper vector. + */ + Teuchos::RCP> basis (const int i) const; + + /** + * Apply unary function @p f to all the elements of the wrapped vector. + */ + void applyUnary (const ROL::Elementwise::UnaryFunction &f); + + /** + * Apply binary function @p f along with ROL::Vector @p rol_vector to all + * the elements of the wrapped vector. + */ + void applyBinary (const ROL::Elementwise::UnaryFunction &f, + const ROL::Vector &rol_vector); + + /** + * Return the accumulated value on applying reduction operation @p r on + * all the elements of the wrapped vector. + */ + value_type reduce (const ROL::Elementwise::ReductionOp &r) const; + + /** + * Print the wrapped vector to the output stream @p outStream. + */ + void print (std::ostream &outStream) const; + + }; + + + /*------------------------------member definitions--------------------------*/ +#ifndef DOXYGEN + + + template + VectorAdaptor:: + VectorAdaptor (const Teuchos::RCP &vector_ptr) + : + vector_ptr (vector_ptr) + {} + + + + template + Teuchos::RCP + VectorAdaptor::getVector () + { + return vector_ptr; + } + + + + template + Teuchos::RCP + VectorAdaptor::getVector () const + { + return vector_ptr; + } + + + + template + void + VectorAdaptor::set (const ROL::Vector &rol_vector) + { + const VectorAdaptor &vector_adaptor = + Teuchos::dyn_cast(rol_vector); + + (*vector_ptr) = *(vector_adaptor.getVector()); + } + + + + template + void + VectorAdaptor::plus (const ROL::Vector &rol_vector) + { + Assert (this->dimension() == rol_vector.dimension(), + ExcDimensionMismatch(this->dimension(), rol_vector.dimension())); + + const VectorAdaptor &vector_adaptor = + Teuchos::dyn_cast(rol_vector); + + *vector_ptr += *(vector_adaptor.getVector()); + } + + + + template + void + VectorAdaptor::axpy (const value_type alpha, + const ROL::Vector &rol_vector) + { + Assert (this->dimension() == rol_vector.dimension(), + ExcDimensionMismatch(this->dimension(), rol_vector.dimension())); + + const VectorAdaptor &vector_adaptor = + Teuchos::dyn_cast(rol_vector); + + vector_ptr->add (alpha, *(vector_adaptor.getVector())); + } + + + + template + int + VectorAdaptor::dimension () const + { + Assert (vector_ptr->size() < std::numeric_limits::max(), + ExcMessage("The size of the vector being used is greater than " + "largest value of type int.")); + return static_cast(vector_ptr->size()); + } + + + + template + void + VectorAdaptor::scale (const value_type alpha) + { + (*vector_ptr) *= alpha; + } + + + + template + typename VectorType::value_type + VectorAdaptor:: + dot (const ROL::Vector &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::value_type + VectorAdaptor::norm() const + { + return vector_ptr->l2_norm(); + } + + + + template + Teuchos::RCP > + VectorAdaptor::clone() const + { + Teuchos::RCP< VectorType> vec_ptr = Teuchos::rcp (new VectorType); + (*vec_ptr) = (*vector_ptr); + + return Teuchos::rcp (new VectorAdaptor (vec_ptr)); + } + + + + template + Teuchos::RCP > + VectorAdaptor::basis (const int i) const + { + Teuchos::RCP 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 e = Teuchos::rcp (new VectorAdaptor(vec_ptr)); + + return e; + } + + + + template + void + VectorAdaptor:: + applyUnary (const ROL::Elementwise::UnaryFunction &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 + void + VectorAdaptor:: + applyBinary (const ROL::Elementwise::UnaryFunction &f, + const ROL::Vector &rol_vector) + { + Assert (this->dimension() == rol_vector.dimension(), + ExcDimensionMismatch(this->dimension(), rol_vector.dimension())); + + const VectorAdaptor &vector_adaptor = + Teuchos::dyn_cast(rol_vector); + + const VectorType &given_rol_vector = *(vector_adaptor.getVector()); + + const typename VectorType::iterator + vend = vector_ptr->end(), + rolend = given_rol_vector.end(); + + for (typename VectorType::iterator + l_iterator = vector_ptr->begin(), r_iterator = given_rol_vector.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::value_type + VectorAdaptor:: + reduce (const ROL::Elementwise::ReductionOp &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 + void + VectorAdaptor::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 diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index d9c364a350..f789a8597f 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -1,6 +1,6 @@ ## --------------------------------------------------------------------- ## -## Copyright (C) 2012 - 2016 by the deal.II authors +## Copyright (C) 2012 - 2017 by the deal.II authors ## ## This file is part of the deal.II library. ## @@ -50,6 +50,7 @@ ADD_SUBDIRECTORY(meshworker) ADD_SUBDIRECTORY(opencascade) ADD_SUBDIRECTORY(particles) ADD_SUBDIRECTORY(physics) +ADD_SUBDIRECTORY(optimization/rol) ADD_SUBDIRECTORY(non_matching) ADD_SUBDIRECTORY(sundials) diff --git a/source/optimization/rol/CMakeLists.txt b/source/optimization/rol/CMakeLists.txt new file mode 100644 index 0000000000..7cbeef9794 --- /dev/null +++ b/source/optimization/rol/CMakeLists.txt @@ -0,0 +1,29 @@ +## --------------------------------------------------------------------- +## +## 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 at +## the top level of the deal.II distribution. +## +## --------------------------------------------------------------------- + +INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) + +SET(_src + ) + +SET(_inst + ) + +FILE(GLOB _header + ${CMAKE_SOURCE_DIR}/include/deal.II/optimization/rol/*.h + ) + +DEAL_II_ADD_LIBRARY(obj_rol OBJECT ${_src} ${_header} ${_inst}) +EXPAND_INSTANTIATIONS(obj_rol "${_inst}") -- 2.39.5