diff options
Diffstat (limited to 'stereographic.hpp')
| -rw-r--r-- | stereographic.hpp | 7 |
1 files changed, 4 insertions, 3 deletions
diff --git a/stereographic.hpp b/stereographic.hpp index 4c3c831..515890e 100644 --- a/stereographic.hpp +++ b/stereographic.hpp @@ -99,17 +99,18 @@ std::tuple<Scalar, Vector, Matrix> stereographicHamGradHess(const Tensor& J, con std::tie(hamiltonian, gradZ, hessZ) = hamGradHess(J, z); Matrix g = stereographicMetric(jacobian); - Matrix gInvJac = g.ldlt().solve(jacobian); + // g is Hermitian, which is positive definite. + Matrix gInvJac = g.llt().solve(jacobian); grad = gInvJac * gradZ; Eigen::Tensor<Scalar, 3> Γ = stereographicChristoffel(ζ, gInvJac); Vector df = jacobian * gradZ; Eigen::Tensor<Scalar, 1> dfT = Eigen::TensorMap<Eigen::Tensor<Scalar, 1>>(df.data(), {df.size()}); - std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(1, 0)}; + std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(0, 0)}; Eigen::Tensor<Scalar, 2> H2 = Γ.contract(dfT, ip); - hess = jacobian * hessZ * jacobian.transpose() + Eigen::Map<Matrix>(H2.data(), ζ.size(), ζ.size()); + hess = jacobian * hessZ * jacobian.transpose() - Eigen::Map<Matrix>(H2.data(), ζ.size(), ζ.size()).transpose(); return {hamiltonian, grad, hess}; } |