From 48d53f26791bccc88aa4bdb0faf1fabca0ca7be6 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Fri, 6 Jul 2012 11:08:52 +0200
Subject: [PATCH] Handle 0 * 1/x for x close to 0 in addHessian
Since 0 * nan = nan
---
dune/tectonic/localnonlinearity.hh | 74 ++++++++++++++++++++----------
1 file changed, 51 insertions(+), 23 deletions(-)
diff --git a/dune/tectonic/localnonlinearity.hh b/dune/tectonic/localnonlinearity.hh
index 6fcefe4c..ee3d8164 100644
--- a/dune/tectonic/localnonlinearity.hh
+++ b/dune/tectonic/localnonlinearity.hh
@@ -54,39 +54,67 @@ template <int dimension> class LocalNonlinearity {
}
}
- /*
- H''(|x|) x \otimes x / |x|^2
- + H'(|x|) [|x|^2 id - x \otimes x] / |x|^3
-
- For i == j, this is (writing k and using the notation from below)
-
- h2 * x[k] * x[k] / normX2 + h1 [normX2 - x[k]^2]/normX^3
- = h2 * x[k] * x[k] / normX2 + h1ox - h1 x[k]^2/normX^3
- = (h2-h1ox) * x[k] * x[k] / normX2 + h1ox
-
- For i != j, this is (again, using the notation from below)
-
- h2 * x[i] * x[j] / normX2 - h1 * x[i] * x[j]/normX^3
- = (h2 - h1ox) * x[i] * x[j] / normX2
+ /** Formula for the derivative:
+
+ \f{align*}{
+ \frac {d^2}{dz^2} H(|z|)
+ &= \frac d{dz} \left( H'(|z|) \otimes \frac z{|z|} \right)\\
+ &= H''(|z|) \frac z{|z|} \otimes \frac z{|z|}
+ + H'(|z|) \left( \frac {|z| \operatorname{id} - z \otimes z/|z|}{|z|^2}
+ \right)\\
+ &= \frac {H''(|z|)}{|z|^2} z \otimes z
+ + \frac {H'(|z|)}{|z|} \operatorname{id}
+ - \frac {H'(|z|)}{|z|^3} z \otimes z\\
+ &= \left( \frac {H''(|z|)}{|z|^2} - \frac {H'(|z|)}{|z|^3} \right) z
+ \otimes z
+ + \frac {H'(|z|)}{|z|} \operatorname{id}
+ \f}
*/
void addHessian(VectorType const &x, MatrixType &A) const {
- if (x == VectorType(0))
- return;
-
double const normX2 = x.two_norm2();
- double const normX = sqrt(normX2);
- double const h1ox = func->rightDifferential(normX) / normX;
- double const h2 = func->second_deriv(normX);
+ double const normX = std::sqrt(normX2);
+ double const normX3 = normX * normX2;
+
+ double const H1 = func->rightDifferential(normX);
+ double const H2 = func->second_deriv(normX);
+
+ // TODO: potential optimisation: factor out (H1 / normX), get rid of normX3
+ double const weight1 = H2 / normX2;
+ double const weight2 = -H1 / normX3;
+ double const weight3 = H1 / normX;
+
+ // {{{ In what follows, we handle the case 0 * (1/x) = 0 with x
+ // close to 0.
+ double tensorweight = 0;
+
+ if (std::isnan(weight1))
+ assert(H2 == 0);
+ else
+ tensorweight += weight1;
+
+ if (std::isnan(weight2))
+ assert(H1 == 0);
+ else
+ tensorweight += weight2;
+
+ double idweight = 0;
+ if (std::isnan(weight3))
+ assert(H1 == 0);
+ else
+ idweight += weight3;
+ // }}}
for (int i = 0; i < dimension; ++i)
for (int j = 0; j < i; ++j) {
- double const entry = (h2 - h1ox) * x[i] * x[j] / normX2;
+ double const entry = tensorweight * x[i] * x[j];
A[i][j] += entry;
A[j][i] += entry;
}
- for (int k = 0; k < dimension; ++k)
- A[k][k] += (h2 - h1ox) * x[k] * x[k] / normX2 + h1ox;
+ for (int k = 0; k < dimension; ++k) {
+ double const entry = tensorweight * x[k] * x[k];
+ A[k][k] += entry + idweight;
+ }
}
void addGradient(VectorType const &x, VectorType &y) const {
--
GitLab