diff options
| -rw-r--r-- | langevin.cpp | 5 | ||||
| -rw-r--r-- | stereographic.hpp | 20 |
2 files changed, 14 insertions, 11 deletions
diff --git a/langevin.cpp b/langevin.cpp index 3bc12ed..05b436e 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -7,6 +7,8 @@ #include "pcg-cpp/include/pcg_random.hpp" #include "randutils/randutils.hpp" +#include <eigen3/Eigen/LU> + #include "complex_normal.hpp" #include "p-spin.hpp" #include "stereographic.hpp" @@ -57,7 +59,8 @@ Vector findSaddle(const Tensor& J, const Vector& z0, double γ0, double δ, doub while (W > ε) { double γ = pow(r.variate<double, std::normal_distribution>(0, γ0), 2); - Vector ζNew = ζ - ddH.ldlt().solve(dH); + Vector dζ = ddH.partialPivLu().solve(dH); + Vector ζNew = ζ - dζ; double WNew; Vector dWNew; 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}; } |