diff options
Diffstat (limited to 'stereographic.hpp')
-rw-r--r-- | stereographic.hpp | 20 |
1 files changed, 10 insertions, 10 deletions
diff --git a/stereographic.hpp b/stereographic.hpp index 59432b7..4c3c831 100644 --- a/stereographic.hpp +++ b/stereographic.hpp @@ -58,25 +58,25 @@ Eigen::Tensor<Scalar, 3> stereographicChristoffel(const Vector& ζ, const Matrix unsigned N = ζ.size() + 1; Eigen::Tensor<Scalar, 3> dJ(N - 1, N - 1, N); - Scalar b = ζ.transpose() * ζ; + Scalar b = 1.0 + (Scalar)(ζ.transpose() * ζ); for (unsigned i = 0; i < N - 1; i++) { for (unsigned j = 0; j < N - 1; j++) { for (unsigned k = 0; k < N - 1; k++) { - dJ(i, j, k) = 16 * sqrt(N) * ζ(i) * ζ(j) * ζ(k) / pow(1.0 + b, 3); + dJ(i, j, k) = 16 * sqrt(N) * ζ(i) * ζ(j) * ζ(k) / pow(b, 3); if (i == j) { - dJ(i, j, k) -= 4 * sqrt(N) * (1.0 + b) * ζ(k) / pow(1.0 + b, 3); + dJ(i, j, k) -= 4 * sqrt(N) * ζ(k) / pow(b, 2); } if (i == k) { - dJ(i, j, k) -= 4 * sqrt(N) * (1.0 + b) * ζ(j) / pow(1.0 + b, 3); + dJ(i, j, k) -= 4 * sqrt(N) * ζ(j) / pow(b, 2); } if (j == k) { - dJ(i, j, k) -= 4 * sqrt(N) * (1.0 + b) * ζ(i) / pow(1.0 + b, 3); + dJ(i, j, k) -= 4 * sqrt(N) * ζ(i) / pow(b, 2); } } - dJ(i, j, N - 1) = - 16 * sqrt(N) * ζ(i) * ζ(j) / pow(1.0 + b, 3); + dJ(i, j, N - 1) = - 16 * sqrt(N) * ζ(i) * ζ(j) / pow(b, 3); if (i == j) { - dJ(i, j, N - 1) += 4 * sqrt(N) * (1.0 + b) * ζ(i) / pow(1.0 + b, 3); + dJ(i, j, N - 1) += 4 * sqrt(N) * ζ(i) / pow(b, 2); } } } @@ -103,13 +103,13 @@ std::tuple<Scalar, Vector, Matrix> stereographicHamGradHess(const Tensor& J, con grad = gInvJac * gradZ; - Eigen::Tensor<Scalar, 3> Γ = stereographicChristoffel(ζ, gInvJac.conjugate()); + 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>(0, 0)}; + std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(1, 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}; } |