From c214ae3559b61474a228d22b22e46cd51bc59703 Mon Sep 17 00:00:00 2001
From: Pauli Virtanen <pauli.t.virtanen@jyu.fi>
Date: Tue, 2 Aug 2022 17:31:05 +0300
Subject: [PATCH] array: bugs + eigen conversion

---
 src/array.hpp | 21 +++++++++++++++++++--
 src/main.cpp  | 20 +++++++++++++++-----
 2 files changed, 34 insertions(+), 7 deletions(-)

diff --git a/src/array.hpp b/src/array.hpp
index bdce598..025f2ff 100644
--- a/src/array.hpp
+++ b/src/array.hpp
@@ -67,7 +67,7 @@ public:
             std::array<size_t, NDim-1> stride;
             size_t offset;
 
-            offset = stride[axis] * pos;
+            offset = offset_ + stride_[axis] * pos;
 
             for (size_t i = 0, j = 0; i < NDim; ++i) {
                 if (i == axis)
@@ -89,7 +89,7 @@ public:
             std::array<size_t, NDim> stride;
             size_t offset;
 
-            offset = stride[axis] * begin;
+            offset = offset_ + stride_[axis] * begin;
 
             for (size_t i = 0; i < NDim; ++i) {
                 if (i == axis)
@@ -143,6 +143,23 @@ public:
 #endif
             return idx;
         }
+
+    typedef Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> EigenMatrix;
+    typedef Eigen::Map<EigenMatrix, 0, Eigen::Stride<Eigen::Dynamic,Eigen::Dynamic> > EigenMap;
+
+    operator EigenMap() const
+        {
+            static_assert(NDim == 2, "matrixes must be two-dimensional");
+            return EigenMap((Scalar *)data_.data() + offset_,
+                            shape_[0], shape_[1],
+                            {(long)stride_[1], (long)stride_[0]});
+        }
+
+    operator EigenMatrix() const
+        {
+            return EigenMap(*this);
+        }
 };
 
+
 #endif
diff --git a/src/main.cpp b/src/main.cpp
index 1e02e12..f4aaac3 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -3,7 +3,7 @@
 #include "action.hpp"
 
 
-template <typename Scalar=double>
+template <typename Scalar>
 void dump(Array<Scalar,3> Q)
 {
     for (size_t i = 0; i < Q.template shape<0>(); ++i) {
@@ -21,12 +21,22 @@ void dump(Array<Scalar,3> Q)
 
 int main()
 {
-    std::vector<double> v(2*2*2);
-    Array<double,3> Q(v,{2,2,2},false);
+    std::vector<ADComplex> v(2*2*2);
+    Array<ADComplex,3> Q0(v,{2,2,2},true);
+    auto Q = Q0.slice<0>(1);
 
-    S_2(Q.slice<2>());
+    std::cout << Q0.index(1u,0u,0u) << std::endl;
 
-    dump<>(Q);
+    S_2(Q);
+
+    std::cout << Q.index(0u,0u) << std::endl;
+
+    Array<ADComplex,3>::EigenMap mat(Q);
+
+    mat(1,1) = 123;
+
+    dump(Q0);
+    std::cout << mat << std::endl;
 
     return 0;
 }
-- 
GitLab