From 26317b1df71a3717063d0a9af36bcee8614c532e Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Wed, 26 Oct 2011 15:19:10 +0200
Subject: [PATCH] Add MyBlockProblem class
based on blocknonlineargsproblem
---
src/myblockproblem.hh | 160 ++++++++++++++++++++++++++++++++++++++++++
1 file changed, 160 insertions(+)
create mode 100644 src/myblockproblem.hh
diff --git a/src/myblockproblem.hh b/src/myblockproblem.hh
new file mode 100644
index 00000000..a7c535f3
--- /dev/null
+++ b/src/myblockproblem.hh
@@ -0,0 +1,160 @@
+// 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;
+};
--
GitLab