Merge branch 'feature/new_tnnmg'

dune/tnnmg/CMakeLists.txt

 add_subdirectory("functionals")
 add_subdirectory("iterationsteps")
+add_subdirectory("linearsolvers")
+add_subdirectory("localsolvers")
 add_subdirectory("nonlinearities")
 add_subdirectory("problem-classes")
+add_subdirectory("projections")
 add_subdirectory("solvers")
+add_subdirectory("test")
+
+
+install(FILES
+    concepts.hh
+    typetraits.hh
+    DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tnnmg)

dune/tnnmg/concepts.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set ts=8 sw=2 et sts=2:
+#ifndef DUNE_TNNMG_CONCEPTS_HH
+#define DUNE_TNNMG_CONCEPTS_HH
+
+
+
+namespace Dune {
+namespace TNNMG {
+namespace Concept {
+
+
+
+struct IsLessThanComparable
+{
+
+  struct RawLessThanComparable
+  {
+    template<class T1, class T2>
+    auto require(const T1& t1, const T2& t2) -> decltype(
+      t1<t2
+    );
+  };
+
+  struct LessThanComparableHelper
+  {
+
+    // if one of both types is not a tuple we can use the raw comparable check
+    template<class T1, class T2,
+      std::enable_if_t<(not IsTupleOrDerived<T1>::value) or (not IsTupleOrDerived<T2>::value), int> = 0>
+    static constexpr auto check(const T1&, const T2&)
+    {
+      return std::integral_constant<bool, models<RawLessThanComparable, T1, T2>()>();
+    }
+
+    // tuples of different length are not comparable
+    template<class... T, class... U,
+      std::enable_if_t<not (std::tuple_size<std::tuple<T...>>::value == std::tuple_size<std::tuple<U...>>::value), int> = 0>
+    static constexpr auto check(const std::tuple<T...>&, const std::tuple<U...>&)
+    {
+      return std::false_type();
+    }
+
+    // empty tuples are always comparable
+    static constexpr auto check(const std::tuple<>&, const std::tuple<>&)
+    {
+      return std::true_type();
+    }
+
+    // check if tuples are comparable
+    template<class T0, class... T, class U0, class... U,
+      std::enable_if_t<(std::tuple_size<std::tuple<T...>>::value == std::tuple_size<std::tuple<U...>>::value), int> = 0>
+    static constexpr auto check(const std::tuple<T0, T...>&, const std::tuple<U0, U...>&)
+    {
+      using HeadComparable = decltype(check(std::declval<T0>(), std::declval<U0>()));
+      using TailComparable = decltype(check(std::declval<std::tuple<T...>>(), std::declval<std::tuple<U...>>()));
+      return std::integral_constant<bool, HeadComparable::value and TailComparable::value>();
+    }
+  };
+
+  template<class T1, class T2>
+  auto require(const T1& t1, const T2& t2) -> decltype(
+    Dune::Concept::requireTrue<decltype(LessThanComparableHelper::check(t1, t2))::value>()
+  );
+};
+
+
+} // namespace Concept
+} // namespace TNNMG
+} // namespace Dune
+
+
+
+#endif // DUNE_TNNMG_CONCEPTS_HH

dune/tnnmg/functionals/CMakeLists.txt

 add_subdirectory("test")
 
 install(FILES
+    affineoperator.hh
+    bcqfconstrainedlinearization.hh
+    boxconstrainedquadraticfunctional.hh
     convexfunctional.hh
     cubeindicatorfunctional.hh
     nonsmoothconvexfunctional.hh
     normfunctional.hh
+    quadraticfunctional.hh
+    quadraticfunctionalconstrainedlinearization.hh
     stronglyconvexfunctional.hh
     sumfunctional.hh
     zerofunctional.hh

dune/tnnmg/functionals/affineoperator.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_TNNMG_FUNCTIONALS_AFFINEOPERATOR_HH
+#define DUNE_TNNMG_FUNCTIONALS_AFFINEOPERATOR_HH
+
+
+
+namespace Dune {
+namespace TNNMG {
+
+namespace Imp {
+
+  template<class MM, class X, class Y>
+  static auto mmv(const MM& m, const X& x, Y&& y, PriorityTag<1>) -> void_t<decltype(m.mmv(x, y))>
+  {
+    m.mmv(x, y);
+  }
+
+  template<class MM, class X, class Y>
+  static void mmv(const MM& m, const X& x, Y&& y, PriorityTag<0>)
+  {
+    y -= m*x;
+  }
+} // namespace Imp
+
+
+/** \brief The first derivative of a functional
+ *
+ * The functional is thought to be smooth here: truncation and other ways to deal
+ * with nonsmoothness happens elsewhere.
+ *
+ * \tparam M Matrix type
+ * \tparam V Vector type
+ */
+template<class M, class V>
+class AffineOperator
+{
+  /** \brief The second derivative of a functional: A fixed matrix wrapped as a function */
+  class ConstantMatrixFunction
+  {
+  public:
+    ConstantMatrixFunction(const M& value) :
+      value_(&value)
+    {}
+
+    const M& operator()(const V& d)
+    { return *value_; }
+
+  private:
+    const M* value_;
+  };
+
+public:
+
+  using Matrix = M;
+  using Vector = V;
+
+  /** \brief Constructor */
+  AffineOperator(const Matrix& matrix, const Vector& linearTerm) :
+    matrix_(&matrix),
+    linearTerm_(&linearTerm)
+  {}
+
+  /** \brief Evaluate the first derivative of the functional at a given point */
+  Vector operator()(const Vector& v) const
+  {
+    Vector result;
+    Solvers::resizeInitializeZero(result,v);
+    matrix_->umv(v, result);
+    result -= *linearTerm_;
+    return result;
+  }
+
+  /** \brief Get the second derivative of the functional
+   * \todo This could be a lambda but gcc-4.9.2 (shipped with Debian
+   * Jessie) does not accept inline friends with auto return type
+   * due to bug-59766. This is fixed in gcc-4.9.3.
+   */
+  friend ConstantMatrixFunction derivative(const AffineOperator& f)
+  {
+    return ConstantMatrixFunction(*f.matrix_);
+  }
+
+protected:
+  const Matrix* matrix_;
+  const Vector* linearTerm_;
+};
+
+
+}   // TNNMG
+}   // Dune
+
+#endif   // DUNE_TNNMG_FUNCTIONALS_AFFINEOPERATOR_HH

dune/tnnmg/functionals/bcqfconstrainedlinearization.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_TNNMG_FUNCTIONAL_BCQFCONSTRAINEDLINEARIZATION_HH
+#define DUNE_TNNMG_FUNCTIONAL_BCQFCONSTRAINEDLINEARIZATION_HH
+
+
+#include <cstddef>
+
+#include <dune/common/typetraits.hh>
+#include <dune/common/hybridutilities.hh>
+
+#include <dune/matrix-vector/algorithm.hh>
+
+
+
+namespace Dune {
+namespace TNNMG {
+
+
+
+/**
+ * \brief A constrained linearization for BoxConstrainedQuadraticFunctional
+ */
+template<class F, class BV>
+class BoxConstrainedQuadraticFunctionalConstrainedLinearization
+{
+  using This = BoxConstrainedQuadraticFunctionalConstrainedLinearization<F,BV>;
+
+  template<class NV, class NBV, class T>
+  static void determineTruncation(const NV& x, const NV& lower, const NV& upper, NBV&& truncationFlags, const T& truncationTolerance)
+  {
+    namespace H = Dune::Hybrid;
+    H::ifElse(IsNumber<NV>(), [&](auto id){
+      if ((id(x) <= id(lower)+truncationTolerance) || (id(x) >= id(upper) - truncationTolerance))
+        id(truncationFlags) = true;
+    }, [&](auto id){
+      H::forEach(H::integralRange(H::size(id(x))), [&](auto&& i) {
+        This::determineTruncation(x[i], lower[i], upper[i], truncationFlags[i], truncationTolerance);
+      });
+    });
+  }
+
+  template<class NV, class NBV>
+  static void truncateVector(NV& x, const NBV& truncationFlags)
+  {
+    namespace H = Dune::Hybrid;
+    H::ifElse(IsNumber<NV>(), [&](auto id){
+      if (id(truncationFlags))
+        id(x) = 0;
+    }, [&](auto id){
+      H::forEach(H::integralRange(H::size(id(x))), [&](auto&& i) {
+        This::truncateVector(x[i], truncationFlags[i]);
+      });
+    });
+  }
+
+  template<class NM, class RBV, class CBV>
+  static void truncateMatrix(NM& A, const RBV& rowTruncationFlags, const CBV& colTruncationFlags)
+  {
+    namespace H = Dune::Hybrid;
+    using namespace Dune::MatrixVector;
+    H::ifElse(IsNumber<NM>(), [&](auto id){
+      if(id(rowTruncationFlags) or id(colTruncationFlags))
+        A = 0;
+    }, [&](auto id){
+      H::forEach(H::integralRange(H::size(id(rowTruncationFlags))), [&](auto&& i) {
+        auto&& Ai = A[i];
+        sparseRangeFor(Ai, [&](auto&& Aij, auto&& j) {
+            This::truncateMatrix(Aij, rowTruncationFlags[i], colTruncationFlags[j]);
+        });
+      });
+    });
+  }
+
+public:
+  using Matrix = typename F::Matrix;
+  using Vector = typename F::Vector;
+  using BitVector = BV;
+  using ConstrainedMatrix = Matrix;
+  using ConstrainedVector = Vector;
+  using ConstrainedBitVector = BitVector;
+
+
+  BoxConstrainedQuadraticFunctionalConstrainedLinearization(const F& f, const BitVector& ignore) :
+    f_(f),
+    ignore_(ignore),
+    truncationTolerance_(1e-10)
+  {}
+
+  void bind(const Vector& x)
+  {
+    negativeGradient_ = derivative(f_)(x);
+    negativeGradient_ *= -1;
+    hessian_ = derivative(derivative(f_))(x);
+    truncationFlags_ = ignore_;
+
+    // determine which components to truncate
+    determineTruncation(x, f_.lowerObstacle(), f_.upperObstacle(), truncationFlags_, truncationTolerance_);
+
+    // truncate gradient and hessian
+    truncateVector(negativeGradient_, truncationFlags_);
+    truncateMatrix(hessian_, truncationFlags_, truncationFlags_);
+  }
+
+  void extendCorrection(ConstrainedVector& cv, Vector& v) const
+  {
+    v = cv;
+    truncateVector(v, truncationFlags_);
+  }
+
+  const BitVector& truncated() const
+  {
+    return truncationFlags_;
+  }
+
+  const auto& negativeGradient() const
+  {
+    return negativeGradient_;
+  }
+
+  const auto& hessian() const
+  {
+    return hessian_;
+  }
+
+private:
+  const F& f_;
+  const BitVector& ignore_;
+
+  double truncationTolerance_;
+
+  Vector negativeGradient_;
+  Matrix hessian_;
+  BitVector truncationFlags_;
+};
+
+
+
+} // end namespace TNNMG
+} // end namespace Dune
+
+
+
+
+#endif // DUNE_TNNMG_FUNCTIONAL_BCQFCONSTRAINEDLINEARIZATION

dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set ts=8 sw=2 et sts=2:
+#ifndef DUNE_TNNMG_FUNCTIONALS_BOXCONSTRAINEDQUADRATICFUNCTIONAL_HH
+#define DUNE_TNNMG_FUNCTIONALS_BOXCONSTRAINEDQUADRATICFUNCTIONAL_HH
+
+#include <cstddef>
+#include <type_traits>
+
+#include <dune/common/concept.hh>
+#include <dune/common/typetraits.hh>
+#include <dune/common/hybridutilities.hh>
+
+#include <dune/matrix-vector/algorithm.hh>
+
+#include <dune/solvers/common/copyorreference.hh>
+
+#include <dune/tnnmg/concepts.hh>
+
+#include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
+
+#include <dune/tnnmg/functionals/quadraticfunctional.hh>
+
+
+
+
+
+namespace Dune {
+namespace TNNMG {
+
+
+
+template<class V, class L, class U,
+  std::enable_if_t<models<Concept::IsLessThanComparable, V, L>(), int> = 0>
+bool checkInsideBox(const V& v, const L& lower, const U& upper)
+{
+  return (lower <= v) && (v <= upper);
+}
+
+template<class V, class L, class U,
+  std::enable_if_t<not models<Concept::IsLessThanComparable, V, L>(), int> = 0>
+bool checkInsideBox(const V& v, const L& lower, const U& upper)
+{
+  namespace H = Dune::Hybrid;
+  bool result = true;
+  H::forEach(H::integralRange(H::size(v)), [&](auto&& i) {
+    result &= checkInsideBox(v[i], lower[i], upper[i]);
+  });
+  return result;
+}
+
+
+
+template<class V, class L, class U, class T,
+  std::enable_if_t<std::is_arithmetic<V>::value, int> = 0>
+void directionalObstacles(const V& x, const V& v, const L& lower, const U& upper, T& defectLower, T& defectUpper)
+{
+  if (v>0)
+  {
+    defectLower = std::max(defectLower, (lower - x)/v);
+    defectUpper = std::min(defectUpper, (upper - x)/v);
+  }
+  else if (v<0)
+  {
+    defectLower = std::max(defectLower, (upper - x)/v);
+    defectUpper = std::min(defectUpper, (lower - x)/v);
+  }
+}
+
+template<class V, class L, class U, class T,
+  std::enable_if_t<not std::is_arithmetic<V>::value, int> = 0>
+void directionalObstacles(const V& x, const V& v, const L& lower, const U& upper, T& defectLower, T& defectUpper)
+{
+  namespace H = Dune::Hybrid;
+  H::forEach(H::integralRange(H::size(v)), [&](auto&& i) {
+    directionalObstacles(x[i], v[i], lower[i], upper[i], defectLower, defectUpper);
+  });
+}
+
+
+
+template<class M, class V, class R>
+class ShiftedBoxConstrainedQuadraticFunctional;
+
+template<class M, class V, class L, class U, class R>
+class BoxConstrainedQuadraticFunctional;
+
+
+
+
+/** \brief Coordinate restriction of a quadratic functional
+ *
+ *  \tparam M Global matrix type
+ *  \tparam V global vector type
+ *  \tparam R Range type
+ */
+template<class M, class V, class R>
+class ShiftedBoxConstrainedQuadraticFunctional : public
+  ShiftedQuadraticFunctional<M,V,R>
+{
+  using Base = ShiftedQuadraticFunctional<M,V,R>;
+public:
+
+  using Matrix = M;
+  using Vector = V;
+  using LowerObstacle = Vector;
+  using UpperObstacle = Vector;
+  using Range = R;
+
+  ShiftedBoxConstrainedQuadraticFunctional(const Matrix& quadraticPart, const Vector& linearPart, const Vector& lowerObstacle, const Vector& upperObstacle, const Vector& origin) :
+    Base(quadraticPart, linearPart, origin),
+    originalLowerObstacle_(&lowerObstacle),
+    originalUpperObstacle_(&upperObstacle)
+  {}
+
+  Range operator()(const Vector& v) const
+  {
+    auto temp = Base::origin();
+    temp += v;
+    if (checkInsideBox(temp, originalLowerObstacle(), originalUpperObstacle()))
+      return Base::operator()(v);
+    return std::numeric_limits<Range>::max();
+  }
+
+  const Vector& originalLowerObstacle() const
+  {
+    return *originalLowerObstacle_;
+  }
+
+  const Vector& originalUpperObstacle() const
+  {
+    return *originalUpperObstacle_;
+  }
+
+protected:
+  const Vector* originalLowerObstacle_;
+  const Vector* originalUpperObstacle_;
+};
+
+// \ToDo This should be an inline friend of ShiftedBoxConstrainedQuadraticFunctional
+// but gcc-4.9.2 shipped with debian jessie does not accept
+// inline friends with auto return type due to bug-59766.
+// Notice, that this is fixed in gcc-4.9.3.
+template<class Index, class M, class V, class R>
+auto coordinateRestriction(const ShiftedBoxConstrainedQuadraticFunctional<M, V, R>& f, const Index& i)
+{
+  using Range = R;
+  using LocalMatrix = std::decay_t<decltype(f.quadraticPart()[i][i])>;
+  using LocalVector = std::decay_t<decltype(f.originalLinearPart()[i])>;
+  using LocalLowerObstacle = std::decay_t<decltype(f.originalLowerObstacle()[i])>;
+  using LocalUpperObstacle = std::decay_t<decltype(f.originalUpperObstacle()[i])>;
+
+  using namespace Dune::MatrixVector;
+  namespace H = Dune::Hybrid;
+
+  LocalVector ri = f.originalLinearPart()[i];
+  const LocalMatrix* Aii_p = nullptr;
+
+  const auto& Ai = f.quadraticPart()[i];
+  sparseRangeFor(Ai, [&](auto&& Aij, auto&& j) {
+    // TODO Here we must implement a wrapper to guarantee that this will work with proxy matrices!
+    H::ifElse(H::equals(j, i), [&](auto&& id){
+      Aii_p = id(&Aij);
+    });
+    Imp::mmv(Aij, f.origin()[j], ri, PriorityTag<1>());
+  });
+
+  LocalLowerObstacle dli = f.originalLowerObstacle()[i];
+  LocalUpperObstacle dui = f.originalUpperObstacle()[i];
+  dli -= f.origin()[i];
+  dui -= f.origin()[i];
+
+  return BoxConstrainedQuadraticFunctional<LocalMatrix&, LocalVector, LocalLowerObstacle, LocalUpperObstacle, Range>(*Aii_p, std::move(ri), std::move(dli), std::move(dui));
+}
+
+
+
+/** \brief Coordinate restriction of a quadratic functional
+ *
+ *  \tparam M Global matrix type
+ *  \tparam V Global vector type
+ *  \tparam L Global lower obstacle type
+ *  \tparam U Global upper obstacle type
+ *  \tparam R Range type
+ */
+template<class M, class V, class L, class U, class R=double>
+class BoxConstrainedQuadraticDirectionalRestriction :
+  public QuadraticDirectionalRestriction<M,V,R>
+{
+  using Base = QuadraticDirectionalRestriction<M,V,R>;
+public:
+
+  using GlobalMatrix = typename Base::GlobalMatrix;
+  using GlobalVector = typename Base::GlobalVector;
+  using GlobalLowerObstacle = L;
+  using GlobalUpperObstacle = U;
+
+  using Matrix = typename Base::Matrix;
+  using Vector = typename Base::Vector;
+
+  using LowerObstacle = Vector;
+  using UpperObstacle = Vector;
+
+  using Range = typename Base::Range;
+
+  BoxConstrainedQuadraticDirectionalRestriction(const GlobalMatrix& matrix, const GlobalVector& linearTerm, const GlobalLowerObstacle& lower, const GlobalLowerObstacle& upper, const GlobalVector& origin, const GlobalVector& direction) :
+    Base(matrix, linearTerm, origin, direction)
+  {
+    defectLower_ = -std::numeric_limits<Vector>::max();
+    defectUpper_ = std::numeric_limits<Vector>::max();
+    directionalObstacles(origin, direction, lower, upper, defectLower_, defectUpper_);
+  }
+
+  Range operator()(const Vector& v) const
+  {
+    DUNE_THROW(Dune::NotImplemented, "Evaluation of QuadraticCoordinateRestriction not implemented");
+  }
+
+  const LowerObstacle& lowerObstacle() const
+  {
+    return defectLower_;
+  }
+
+  const UpperObstacle& upperObstacle() const
+  {
+    return defectUpper_;
+  }
+
+protected:
+  LowerObstacle defectLower_;
+  UpperObstacle defectUpper_;
+};
+
+
+
+/** \brief Box constrained quadratic functional
+ *
+ *  \tparam M Matrix type
+ *  \tparam V Vector type
+ *  \tparam L Lower obstacle type
+ *  \tparam U Upper obstacle type
+ *  \tparam R Range type