Commit 29a0b0b8 authored by Carsten Gräser's avatar Carsten Gräser
Browse files

Add support for regularizing truncated diagonal entries

This adds support for regularization of truncated diagonal
entries in `BoxConstrainedQuadraticFunctionalConstrainedLinearization`
If this is enabled, diagonal entries of truncated rows are set
to 1. In case of nested matrices this is applied recusively,
such that only the truncated 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.
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.
parent bbe3371b
Pipeline #41749 passed with stage
in 11 minutes and 30 seconds
......@@ -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,23 @@ 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 recusively,
* 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.
* 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 +136,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 +165,7 @@ private:
const BitVector& ignore_;
double truncationTolerance_;
bool regularizeDiagonal_;
Vector negativeGradient_;
Matrix hessian_;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment