Terminal value pinning with identity block in hessian

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));
             }