action: initial implementation fo S_2

src/action.hpp

@@ -8,33 +8,55 @@
 #ifndef ACTION_HPP_
 #define ACTION_HPP_
 
+#include <Eigen/Core>
+#include <Eigen/LU>
+
 #include "common.hpp"
 #include "array.hpp"
 
-template <typename Scalar, typename Shape>
-inline auto get_Q_matrix(array::Array<Scalar, Shape> Q)
+
+template <typename Scalar>
+inline Eigen::Matrix<Scalar,4,4>
+get_Q_matrix(array::Array<Scalar, array::Shape<2,2,2> > Q)
 {
-    return 0;
+    auto g = Q.part(0u).matrix();
+    auto gt = Q.part(1u).matrix();
+    Eigen::Matrix<Scalar,2,2> I, N, Nt;
+    Eigen::Matrix<Scalar,4,4> Qm;
+
+    I << 1, 0,
+        0, 1;
+
+    N = (I + g*gt).inverse();
+    Nt = (I + gt*g).inverse();
+
+    Qm.template block<2,2>(0,0) = N * (I - g * gt);
+    Qm.template block<2,2>(0,2) = N * 2 * g * gt;
+    Qm.template block<2,2>(2,0) = Nt * 2 * gt * g;
+    Qm.template block<2,2>(2,2) = Nt * (I + gt * g);
+
+    return Qm;
 }
 
 template <typename Scalar, typename Shape, typename Shape2>
-inline Scalar S_2(array::Array<Scalar, Shape> Q, array::Array<Complex, Shape2> U, double h)
+inline Complex S_2(array::Array<Scalar, Shape> Q, array::Array<Complex, Shape2> U, double h)
 {
     static const int N_i[4] = {0, 1, -1, 0};
     static const int N_j[4] = {1, 0, 0, -1};
-    Scalar S = Scalar(0);
+    static const size_t invk[4] = {3, 2, 1, 0};
+    Complex S = Scalar(0);
 
     size_t nx = Q.dim(0);
     size_t ny = Q.dim(1);
 
     for (size_t i = 0; i < nx; ++i) {
         for (size_t j = 0; j < ny; ++j) {
-            auto Q1 = Q.part(i,j).matrix();
+            auto Q1 = get_Q_matrix(Q.part(i,j));
 
             /* Sum over neighbors */
             for (size_t k = 0; k < 4; ++k) {
-                size_t di = N_i[k];
-                size_t dj = N_j[k];
+                int di = N_i[k];
+                int dj = N_j[k];
                 
                 if (i == 0 && di < 0)
                         continue;
@@ -48,11 +70,12 @@ inline Scalar S_2(array::Array<Scalar, Shape> Q, array::Array<Complex, Shape2> U
                 size_t i2 = i + di;
                 size_t j2 = j + dj;
 
-                auto Q2 = Q.part(i2,j2).matrix();
-                auto U12 = U.part(i,j,k);
-                auto U21 = U.part(j,i,k);
+                auto Q2 = get_Q_matrix(Q.part(i,j));
+                auto U12 = U.part(i,j,k).matrix();
+                auto U21 = U.part(i2,j2,invk[k]).matrix();
+                auto M = (Q1 - U12*Q2*U21);
 
-                S += 2/(4*h) * (Q1 - U12*Q2*U21).trace();
+                S += 2/(4*h) * (M*M).trace();
             }
         }
     }

src/main.cpp

@@ -25,22 +25,12 @@ void dump(Array<Scalar,Shape> Q)
 
 int main()
 {
-    std::vector<double> v(10*10*2*2);
-    Array<double,Shape<Dynamic,Dynamic,2,2> > Q0(v.data(), v.size(), {10,10,2,2});
-    auto Q = Q0.part(0u,1u);
+    StoredArray<double,Shape<Dynamic,Dynamic,2,2,2> > Q({10,10,2,2,2});
+    StoredArray<Complex,Shape<Dynamic,Dynamic,4,4,4> > U({10,10,4,4,4});
 
-    S_2(Q0);
+    Complex res = S_2(Q, U, 0.1);
 
-    std::cout << "size = " << Q.size() << std::endl;
+    std::cout << res << std::endl;
 
-    dump(Q0);
-
-    auto mat = Q.matrix();
-
-    mat = mat * mat * mat;
-
-    static_assert(std::is_same<decltype(mat), Eigen::Map<Eigen::Matrix<double, 2, 2, Eigen::RowMajor>, Eigen::Aligned8 > >::value, "wrong type");
-
-    dump(Q0);
     return 0;
 }