Edit

meson.build

-project('usadelndsoc', 'cpp')
+project('usadelndsoc', 'cpp',
+        default_options : ['cpp_std=c++17'])
 
 cppad_dep = dependency('cppad')
-deps = [cppad_dep]
-executable('main', 'src/main.cpp', dependencies: deps)
+eigen3_dep = dependency('eigen3')
+
+deps = [cppad_dep, eigen3_dep]
+executable('main', 'src/main.cpp', dependencies : deps)

src/action.hpp

+#ifndef ACTION_HPP_
+#define ACTION_HPP_
+
+#include "common.hpp"
+
+
+inline ADComplex S_2(Array<ADComplex,3> Q)
+{
+    size_t i = 2;
+    return Q(0u, i, 2u);
+}
+
+
+
+#endif

src/common.h

-#ifndef COMMON_H_
-#define COMMON_H_
-
-#include <cppad/cppad.hpp>
-
-typedef CppAD::AD<double> ADDouble;
-typedef std::complex<ADDouble> ADComplex;
-
-#endif

src/common.hpp

+#ifndef COMMON_H_
+#define COMMON_H_
+
+#include <array>
+#include <tuple>
+#include <iostream>
+
+#include <cppad/cppad.hpp>
+#include <cppad/example/cppad_eigen.hpp>
+
+#include <Eigen/Core>
+
+typedef CppAD::AD<double> ADDouble;
+typedef std::complex<ADDouble> ADComplex;
+
+template <typename Scalar, int NDim>
+class Array
+{
+private:
+    std::vector<Scalar>& data_;
+    std::array<size_t, NDim> shape_;
+    std::array<size_t, NDim> stride_;
+
+protected:
+    template <typename... Idx>
+    size_t index(Idx... idxs) const
+        {
+            const std::array<size_t,NDim> m{idxs...};
+            size_t idx = 0;
+
+            for (size_t i = 0; i < NDim; ++i)
+                idx += stride_[i] * m[i];
+
+            return idx;
+        }
+    
+public:
+    Array(std::vector<Scalar>& data, std::array<size_t, NDim> shape) : data_(data), shape_(shape)
+        {
+            if (NDim > 0) {
+                size_t i = NDim - 1;
+                stride_[i] = 1;
+                while (i-- > 0)
+                    stride_[i] = shape[i+1] * stride_[i+1];
+            }
+        }
+
+    template <typename... Idx>
+    Scalar& operator()(Idx... idxs)
+        {
+            return data_[index(idxs...)];
+        }
+
+    template <typename... Idx>
+    const Scalar&& operator()(Idx... idxs) const
+        {
+            return data_[index(idxs...)];
+        }
+};
+
+
+#endif

src/main.cpp

-#include "common.h"
+#include "common.hpp"
+#include "action.hpp"
 
 
 int main()
 {
+    std::vector<double> v(3*3*3);
+    Array<double,3> Q(v,{3,3,3});
+
+    for (size_t i = 0; i < 3; ++i) {
+        for (size_t j = 0; i < 3; ++i) {
+            for (size_t k = 0; i < 3; ++i) {
+                Q(i,j,k) = i+j+k;
+            }
+        }
+    }
+
+    for (size_t i = 0; i < 3; ++i) {
+        for (size_t j = 0; i < 3; ++i) {
+            for (size_t k = 0; i < 3; ++i) {
+                std::cout << Q(i,j,k) << std::endl;
+            }
+        }
+    }
 
     return 0;
 }