Merge branch 'feature/regularize-truncated-diagonal' into 'master'

Add support for regularizing truncated diagonal entries

See merge request !26

dune/tnnmg/functionals/bcqfconstrainedlinearization.hh

@@ -61,21 +61,26 @@ class BoxConstrainedQuadraticFunctionalConstrainedLinearization
   }
 
   template<class NM, class RBV, class CBV>
-  static void truncateMatrix(NM& A, const RBV& rowTruncationFlags, const CBV& colTruncationFlags)
+  static void truncateMatrix(NM& A, const RBV& rowTruncationFlags, const CBV& colTruncationFlags, bool regularizeDiagonal)
   {
     namespace H = Dune::Hybrid;
     using namespace Dune::MatrixVector;
     if constexpr (IsNumber<NM>())
     {
       if (rowTruncationFlags or colTruncationFlags)
-        A = 0;
+      {
+        if (regularizeDiagonal)
+          A = 1;
+        else
+          A = 0;
+      }
     }
     else
     {
       H::forEach(H::integralRange(H::size(rowTruncationFlags)), [&](auto&& i) {
         auto&& Ai = A[i];
         sparseRangeFor(Ai, [&](auto&& Aij, auto&& j) {
-            This::truncateMatrix(Aij, rowTruncationFlags[i], colTruncationFlags[j]);
+            This::truncateMatrix(Aij, rowTruncationFlags[i], colTruncationFlags[j], regularizeDiagonal and (i==j));
         });
       });
     }
@@ -93,7 +98,8 @@ public:
   BoxConstrainedQuadraticFunctionalConstrainedLinearization(const F& f, const BitVector& ignore) :
     f_(f),
     ignore_(ignore),
-    truncationTolerance_(1e-10)
+    truncationTolerance_(1e-10),
+    regularizeDiagonal_(false)
   {}
 
   void setTruncationTolerance(double tolerance)
@@ -101,6 +107,24 @@ public:
     truncationTolerance_ = tolerance;
   }
 
+  /**
+   * \brief Enable/disable regularization of truncated diagonal entries
+   *
+   * If this is enabled, diagonal entries of truncated rows are set
+   * to 1. In case of nested matrices this is applied recursively,
+   * such that only the scalar diagonal entries of nontrivial
+   * diagonal blocks are modified. This can be handy if the truncated
+   * linearized problem should be treated by a linear solver
+   * that is not capable of handling truncated rows and columns.
+   * By default this is disabled, such that linear solvers for the
+   * linearized problem have to be robust with respect to singular
+   * problems due to truncation.
+   */
+  void enableDiagonalRegularization(bool regularizeDiagonal)
+  {
+    regularizeDiagonal_ = regularizeDiagonal;
+  }
+
   void bind(const Vector& x)
   {
     negativeGradient_ = derivative(f_)(x);
@@ -113,7 +137,7 @@ public:
 
     // truncate gradient and hessian
     truncateVector(negativeGradient_, truncationFlags_);
-    truncateMatrix(hessian_, truncationFlags_, truncationFlags_);
+    truncateMatrix(hessian_, truncationFlags_, truncationFlags_, regularizeDiagonal_);
   }
 
   void extendCorrection(ConstrainedVector& cv, Vector& v) const
@@ -142,6 +166,7 @@ private:
   const BitVector& ignore_;
 
   double truncationTolerance_;
+  bool regularizeDiagonal_;
 
   Vector negativeGradient_;
   Matrix hessian_;