From 55dd443a4c5fc527b03609b000e92782916e60e2 Mon Sep 17 00:00:00 2001
From: Pauli Virtanen <pauli.t.virtanen@jyu.fi>
Date: Tue, 16 Aug 2022 18:34:11 +0300
Subject: [PATCH] Terminal value pinning with identity block in hessian
---
src/action.hpp | 23 +++++++++++++----------
1 file changed, 13 insertions(+), 10 deletions(-)
diff --git a/src/action.hpp b/src/action.hpp
index f60aa70..8e35dfe 100644
--- a/src/action.hpp
+++ b/src/action.hpp
@@ -162,7 +162,7 @@ inline Scalar S_H(const array::ArrayView<const Mask, Shape0>& mask,
Scalar S{0.0};
- if (alpha == 0)
+ if (alpha == Complex(0))
return S;
for (int i = 0; i < nx-1; ++i) {
@@ -263,16 +263,19 @@ inline Scalar S(const array::ArrayView<const Mask, Shape0>& mask,
if (mask(i,j) == Mask::VACUUM) {
continue;
} else if (mask(i,j) == Mask::TERMINAL) {
- auto Qv = get_Q_matrix(Q.part(i,j));
- for (size_t p = 0; p < 4; ++p) {
- for (size_t q = 0; q < 4; ++q) {
- auto w = flat_value(Qv(p,q));
- Qmat(i,j,p,q) = w;
-
- /* Pin value to the flat value */
- S = S + (Qv(p,q) - w) * (Qv(p,q) - w);
- }
+ const Scalar *Qp = Q.part(i,j).data();
+ Array<Scalar,Shape<2,2,2> > Qpart;
+
+ for (size_t p = 0; p < Qpart.size(); ++p) {
+ auto w = flat_value(Qp[p]);
+
+ Qpart.data()[p] = w;
+
+ /* Pin value to the flat value */
+ S = S + 0.5 * (Qp[p] - w) * (Qp[p] - w);
}
+
+ Qmat.part(i,j).matrix() = get_Q_matrix(Qpart);
} else {
Qmat.part(i,j).matrix() = get_Q_matrix(Q.part(i,j));
}
--
GitLab