diff options
Diffstat (limited to 'stereographic.hpp')
-rw-r--r-- | stereographic.hpp | 19 |
1 files changed, 10 insertions, 9 deletions
diff --git a/stereographic.hpp b/stereographic.hpp index 515890e..a3f2cc6 100644 --- a/stereographic.hpp +++ b/stereographic.hpp @@ -54,7 +54,8 @@ Matrix stereographicMetric(const Matrix& J) { return J * J.adjoint(); } -Eigen::Tensor<Scalar, 3> stereographicChristoffel(const Vector& ζ, const Matrix& gInvJac) { +// Gives the Christoffel symbol, with Γ_(i1, i2)^(i3). +Eigen::Tensor<Scalar, 3> stereographicChristoffel(const Vector& ζ, const Matrix& gInvJacConj) { unsigned N = ζ.size() + 1; Eigen::Tensor<Scalar, 3> dJ(N - 1, N - 1, N); @@ -83,7 +84,7 @@ Eigen::Tensor<Scalar, 3> stereographicChristoffel(const Vector& ζ, const Matrix std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(2, 1)}; - return dJ.contract(Eigen::TensorMap<Eigen::Tensor<const Scalar, 2>>(gInvJac.data(), N - 1, N), ip); + return dJ.contract(Eigen::TensorMap<Eigen::Tensor<const Scalar, 2>>(gInvJacConj.data(), N - 1, N), ip); } std::tuple<Scalar, Vector, Matrix> stereographicHamGradHess(const Tensor& J, const Vector& ζ) { @@ -98,19 +99,19 @@ std::tuple<Scalar, Vector, Matrix> stereographicHamGradHess(const Tensor& J, con Matrix hessZ; std::tie(hamiltonian, gradZ, hessZ) = hamGradHess(J, z); - Matrix g = stereographicMetric(jacobian); - // g is Hermitian, which is positive definite. - Matrix gInvJac = g.llt().solve(jacobian); + Matrix metric = stereographicMetric(jacobian); - grad = gInvJac * gradZ; + grad = metric.llt().solve(jacobian) * gradZ; - Eigen::Tensor<Scalar, 3> Γ = stereographicChristoffel(ζ, gInvJac); + /* This is much slower to calculate than the marginal speedup it offers... + Eigen::Tensor<Scalar, 3> Γ = stereographicChristoffel(ζ, gInvJac.conjugate()); 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>(0, 0)}; + std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(2, 0)}; Eigen::Tensor<Scalar, 2> H2 = Γ.contract(dfT, ip); + */ - hess = jacobian * hessZ * jacobian.transpose() - Eigen::Map<Matrix>(H2.data(), ζ.size(), ζ.size()).transpose(); + hess = jacobian * hessZ * jacobian.transpose(); // - Eigen::Map<Matrix>(H2.data(), ζ.size(), ζ.size()); return {hamiltonian, grad, hess}; } |