Add MyBlockProblem class

based on blocknonlineargsproblem

src/myblockproblem.hh

+// Based on dune/tnnmg/problem-classes/blocknonlineargsproblem.hh
+
+// #include <dune/common/parametertree.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/tnnmg/problem-classes/bisection.hh>
+#include <dune/tnnmg/problem-classes/convexproblem.hh>
+#include <dune/tnnmg/problem-classes/nonlinearity.hh>
+#include <dune/tnnmg/problem-classes/onedconvexfunction.hh>
+
+/** \brief Base class for problems where each block can be solved with a scalar
+ * Gauss-Seidel method */
+template <class ConvexProblemTypeTEMPLATE> class MyBlockProblem {
+public:
+  typedef ConvexProblemTypeTEMPLATE ConvexProblemType;
+  typedef typename ConvexProblemType::NonlinearityType NonlinearityType;
+  typedef typename ConvexProblemType::VectorType VectorType;
+  typedef typename ConvexProblemType::MatrixType MatrixType;
+  typedef typename ConvexProblemType::LocalVectorType LocalVectorType;
+  typedef typename ConvexProblemType::LocalMatrixType LocalMatrixType;
+
+  static const int block_size = ConvexProblemType::block_size;
+
+  /** \brief Solves one local system using a scalar Gauss-Seidel method */
+  class IterateObject;
+
+  MyBlockProblem(ConvexProblemType& problem) : problem(problem) {
+    bisection = Bisection(0.0, 1.0, 1e-15, true, 1e-14);
+  };
+
+  /** \brief Constructs and returns an iterate object */
+  IterateObject getIterateObject() { return IterateObject(bisection, problem); }
+
+private:
+  // problem data
+  ConvexProblemType& problem;
+
+  // commonly used minimization stuff
+  Bisection bisection;
+};
+
+/** \brief Solves one local system using a scalar Gauss-Seidel method */
+template <class ConvexProblemTypeTEMPLATE>
+class MyBlockProblem<ConvexProblemTypeTEMPLATE>::IterateObject {
+  friend class MyBlockProblem;
+
+protected:
+  /** \brief Constructor, protected so only friends can instantiate it
+   * \param bisection The class used to do a scalar bisection
+   * \param problem The problem including quadratic part and nonlinear/nonsmooth
+   * part
+   */
+  IterateObject(const Bisection& bisection, ConvexProblemType& problem)
+      : problem(problem),
+        bisection(bisection),
+        local_J(1.0, 0.0, problem.phi, 0, 0) {};
+
+public:
+  /** \brief Set the current iterate */
+  void setIterate(VectorType& u) {
+    // TODO How should the rang-1 term be handled
+    // ????????????????????????????????
+    // s = problem.Am*u;
+    problem.phi.setVector(u);
+
+    this->u = u;
+    return;
+  };
+
+  /** \brief Update the i-th block of the current iterate */
+  void updateIterate(const LocalVectorType& ui, int i) {
+    // TODO How should the rang-1 term be handled
+    // ????????????????????????????????
+    // s += (ui-u[i]) * problem.Am[i];
+
+    for (size_t j = 0; j < block_size; ++j)
+      problem.phi.updateEntry(i, ui[j], j);
+
+    u[i] = ui;
+    return;
+  };
+
+  /** \brief Minimize a local problem using a scalar Gauss-Seidel method
+   * \param[out] ui The solution
+   * \param i Block number
+   * \param ignore Set of degrees of freedom to leave untouched
+   *
+   * \return The minimizer of the local functional in the variable ui
+   */
+  void solveLocalProblem(
+      LocalVectorType& ui, int i,
+      const typename Dune::BitSetVector<block_size>::const_reference ignore) {
+    const LocalMatrixType* Aii;
+    LocalVectorType b(0.0);
+
+    typename MatrixType::row_type::ConstIterator it = problem.A[i].begin();
+    typename MatrixType::row_type::ConstIterator end = problem.A[i].end();
+    for (; it != end; ++it) {
+      int col = it.index();
+      if (col == i)
+        Aii = &(*it);
+      else
+        it->mmv(u[col], b);
+    }
+    b *= problem.a;
+    b += problem.f[i];
+    // TODO How should the rang-1 term be handled
+    // ????????????????????????????????
+
+    ui = u[i];
+    LocalVectorType ui_old = ui;
+    local_J.i = i;
+    for (size_t j = 0; j < block_size; ++j) {
+      if (ignore.test(j))
+        continue;
+
+      local_J.A = problem.a * (*Aii)[j][j];
+      local_J.b = b[j];
+      local_J.j = j;
+
+      typename LocalMatrixType::row_type::ConstIterator it = (*Aii)[j].begin();
+      typename LocalMatrixType::row_type::ConstIterator end = (*Aii)[j].end();
+      for (; it != end; ++it)
+        if (it.index() != j)
+          local_J.b -= (*it) * problem.a * ui[it.index()];
+
+      ui[j] = bisection.minimize(local_J, ui[j], ui[j], bisectionsteps);
+
+      // we need to update the intermediate local values ...
+      problem.phi.updateEntry(i, ui[j], j);
+    }
+    // ... and to restore the old ones after computation
+    for (int j = 0; j < block_size; ++j)
+      problem.phi.updateEntry(i, ui_old[j], j);
+
+    //            local_J.A += problem.am*problem.Am[i]*problem.Am[i];
+    //            local_J.b = problem.f[i] + problem.a*local_J.b;
+    //            local_J.b -= (s-u[i]*problem.Am[i])*problem.am*problem.Am[i];
+
+    //            ui = bisection.minimize<OneDConvexFunction>(local_J, u[i],
+    // u[i], bisectionsteps);
+    return;
+  };
+
+private:
+  // problem data
+  ConvexProblemType& problem;
+
+  // commonly used minimization stuff
+  Bisection bisection;
+  OneDConvexFunction<NonlinearityType> local_J;
+
+  // state data for smoothing procedure used by:
+  // setIterate, updateIterate, solveLocalProblem
+  VectorType u;
+
+  /** \brief Keeps track of the total number of bisection steps that were
+   * performed */
+  int bisectionsteps;
+};