diff --git a/meson.build b/meson.build
index 928852aac311eedb009061edc1b981838a93a929..b83c57dfb36deb459ee2487a8a86db8e6a28e351 100644
--- a/meson.build
+++ b/meson.build
@@ -1,5 +1,7 @@
-project('usadelndsoc', 'cpp',
- default_options : ['cpp_std=c++17'])
+project('usadelndsoc', ['cpp'],
+ version : '0.1',
+ meson_version: '>= 0.55.0',
+ default_options : ['cpp_std=c++17', 'build.cpp_std=c++17'])
cppad_dep = dependency('cppad')
eigen3_dep = dependency('eigen3')
diff --git a/src/array.hpp b/src/array.hpp
index b8faf19170c042cf282d748638a240b84bce9f19..bbdec4f37f03210b2624b2d909e82847eb810823 100644
--- a/src/array.hpp
+++ b/src/array.hpp
@@ -19,6 +19,7 @@ struct Shape<>
{
static const size_t head = 1;
static const size_t ndim = 0;
+ static const size_t size = 1;
static const bool fixed = true;
static constexpr size_t dim(size_t i) { return 0; }
static constexpr size_t stride(size_t i) { return 1; }
@@ -31,6 +32,7 @@ struct Shape<Dim0,Dims...>
typedef Shape<Dims...> Tail;
static const size_t ndim = 1 + sizeof...(Dims);
static const bool fixed = (Dim0 != Dynamic) && Tail::fixed;
+ static const size_t size = fixed ? Dim0*Tail::size : Dynamic;
static constexpr size_t dim(size_t i)
{
@@ -51,12 +53,22 @@ struct Shape<Dim0,Dims...>
}
}
- static constexpr std::array<size_t, ndim> dims() { return {Dim0, Dims...}; }
+ static constexpr std::array<size_t, ndim> shape() { return {Dim0, Dims...}; }
+
+ static constexpr std::array<size_t, ndim> strides()
+ {
+ std::array<size_t, ndim> m;
+
+ for (size_t i = 0; i < ndim; ++i)
+ m[i] = stride(i);
+ return m;
+ }
};
namespace detail
{
- /*! Extract Nth tail shape */
+ /*! Extract Nth tail shape
+ */
template <size_t I, typename Shape>
struct TailNth;
@@ -71,47 +83,116 @@ namespace detail
{
using type = typename TailNth<I-1, typename Shape::Tail>::type;
};
-}
-/*! Simple data-by-reference strided array class a la Fortran
- */
-template <typename Scalar, typename Shape, typename Storage = std::vector<Scalar> >
-class Array
-{
-private:
- Storage& data_;
- std::array<size_t, Shape::ndim> shape_;
- std::array<size_t, Shape::ndim> stride_;
- size_t offset_ = 0;
+ /*! Base class for compile-time fixed size arrays
+ */
+ template <typename Shape>
+ class fixed_base
+ {
+ protected:
+ size_t offset_;
- void data_bounds_check() const
- {
+ template <typename Storage>
+ void data_bounds_check(Storage& data) const
+ {
#ifndef NO_BOUNDS_CHECK
- size_t total = 1;
+ if (Shape::size > data.size())
+ throw std::out_of_range("data array too small");
+#endif
+ }
- for (size_t i = 0; i < Shape::ndim; ++i) {
- total *= shape_[i];
- if ((Shape::dim(i) != Dynamic && Shape::dim(i) != shape_[i])
- || shape_[i] == Dynamic)
- throw std::out_of_range("mismatch with fixed shape");
+ public:
+ fixed_base(size_t offset=0) : offset_(offset)
+ {
+ static_assert(Shape::fixed, "array shape must be compile-time fixed");
}
- total += offset_;
- if (total > data_.size())
- throw std::out_of_range("data array too small");
+
+ constexpr size_t dim(size_t i) const { return Shape::dim(i); }
+ constexpr size_t stride(size_t i) const { return Shape::stride(i); }
+
+ constexpr std::array<size_t, Shape::ndim> shape() const { return Shape::shape(); }
+ constexpr std::array<size_t, Shape::ndim> strides() const { return Shape::strides(); }
+ constexpr size_t offset() const { return offset_; }
+ };
+
+ /*! Base class for compile-time dynamic size arrays
+ */
+ template <typename Shape>
+ class dynamic_base
+ {
+ protected:
+ std::array<size_t, Shape::ndim> shape_;
+ std::array<size_t, Shape::ndim> strides_;
+ size_t offset_;
+
+ template <typename Storage>
+ void data_bounds_check(Storage& data) const
+ {
+#ifndef NO_BOUNDS_CHECK
+ size_t total = 1;
+
+ for (size_t i = 0; i < Shape::ndim; ++i) {
+ total *= dim(i);
+ if (shape_[i] != dim(i))
+ throw std::out_of_range("mismatch with fixed shape");
+ }
+ total += offset_;
+ if (total > data.size())
+ throw std::out_of_range("data array too small");
#endif
}
- void init_strides()
- {
- for (size_t i = Shape::ndim; i > 0; --i) {
- if (Shape::stride(i-1) != Dynamic)
- stride_[i-1] = Shape::stride(i-1);
- else if (i == Shape::ndim)
- stride_[i-1] = 1;
- else
- stride_[i-1] = shape_[i] * stride_[i];
+ public:
+ dynamic_base(const std::array<size_t, Shape::ndim> shape, size_t offset)
+ : shape_(shape), offset_(offset)
+ {
+ static_assert(!Shape::fixed, "array shape must not be compile-time fixed");
+
+ for (size_t i = Shape::ndim; i > 0; --i) {
+ if (Shape::stride(i-1) != Dynamic)
+ strides_[i-1] = Shape::stride(i-1);
+ else if (i == Shape::ndim)
+ strides_[i-1] = 1;
+ else
+ strides_[i-1] = shape_[i] * strides_[i];
+ }
}
- }
+
+ size_t dim(size_t axis) const
+ {
+#ifndef NO_BOUNDS_CHECK
+ if (axis >= Shape::ndim)
+ throw std::out_of_range("array axis must be < ndim");
+#endif
+ size_t n = Shape::dim(axis);
+ return n == Dynamic ? shape_[axis] : n;
+ }
+
+ size_t stride(size_t axis) const
+ {
+#ifndef NO_BOUNDS_CHECK
+ if (axis >= Shape::ndim)
+ throw std::out_of_range("array axis must be < ndim");
+#endif
+ size_t n = Shape::stride(axis);
+ return n == Dynamic ? strides_[axis] : n;
+ }
+
+ constexpr std::array<size_t, Shape::ndim>& shape() const { return shape_; }
+ constexpr std::array<size_t, Shape::ndim>& strides() const { return strides_; }
+ constexpr size_t offset() const { return offset_; }
+ };
+}
+
+/*! Simple data-by-reference strided array class a la Fortran
+ */
+template <typename Scalar, typename Shape, typename Storage = std::vector<Scalar> >
+class Array : public std::conditional_t<Shape::fixed, detail::fixed_base<Shape>, detail::dynamic_base<Shape> >
+{
+private:
+ typedef std::conditional_t<Shape::fixed, detail::fixed_base<Shape>, detail::dynamic_base<Shape> > Base;
+
+ Storage& data_;
template <size_t axis>
static constexpr Eigen::Index eigen_shape()
@@ -122,67 +203,48 @@ private:
public:
Array(Storage& data, const std::array<size_t, Shape::ndim> shape, const size_t offset=0)
- : data_(data), shape_(shape), offset_(offset)
+ : Base(shape, offset), data_(data)
{
- init_strides();
- data_bounds_check();
+ Base::data_bounds_check(data);
}
- Array(Storage& data)
- : data_(data), offset_(0)
+ Array(Storage& data, const size_t offset=0)
+ : Base(offset), data_(data)
{
- static_assert(Shape::fixed, "array shape is not compile-time fixed");
- shape_ = Shape::dims();
- init_strides();
- data_bounds_check();
+ Base::data_bounds_check(data);
}
template <typename... Idx>
- Scalar& operator()(Idx... idxs) { return data_[index(idxs...)]; }
-
- template <typename... Idx>
- const Scalar&& operator()(Idx... idxs) const { return data_[index(idxs...)]; }
-
- Storage& data() { return data_; }
-
- template <size_t axis>
- const size_t shape() const
+ Scalar& operator()(Idx... idxs)
{
- static_assert(axis < Shape::ndim, "invalid axis");
- constexpr size_t n = Shape::dim(axis);
- if (n != Dynamic)
- return n;
- return shape_[axis];
+ static_assert(sizeof...(Idx) == Shape::ndim,
+ "number of indices must be equal to the number of dimensions");
+ return data_[index(idxs...)];
}
- template <size_t axis>
- const size_t stride() const
+ template <typename... Idx>
+ const Scalar&& operator()(Idx... idxs) const
{
- static_assert(axis < Shape::ndim, "invalid axis");
- if (Shape::stride(axis) != Dynamic)
- return Shape::stride(axis);
- return stride_[axis];
+ static_assert(sizeof...(Idx) == Shape::ndim,
+ "number of indices must be equal to the number of dimensions");
+ return data_[index(idxs...)];
}
- size_t offset() const { return offset_; }
+ Storage& data() { return data_; }
template <typename... Idx>
size_t index(Idx... idxs) const
{
- static_assert(sizeof...(idxs) == Shape::ndim,
- "number of indices must equal the number of dimensions");
+ static_assert(sizeof...(idxs) <= Shape::ndim,
+ "number of indices must be smaller than the number of dimensions");
const std::array<size_t, Shape::ndim> m{idxs...};
- size_t idx = offset_;
+ size_t idx = Base::offset();
- for (size_t i = 0; i < Shape::ndim; ++i) {
- if (Shape::stride(i) != Dynamic) {
- idx += Shape::stride(i) * m[i];
- } else {
- idx += stride_[i] * m[i];
- }
+ for (size_t i = 0; i < sizeof...(idxs); ++i) {
+ idx += Base::stride(i) * m[i];
#ifndef NO_BOUNDS_CHECK
- if (m[i] >= shape_[i])
+ if (m[i] >= Base::dim(i))
throw std::out_of_range("index out of bounds");
#endif
}
@@ -194,31 +256,34 @@ public:
return idx;
}
+ /*! Extract partial matrix (compile-time fixed size) */
template <typename... Idx>
- Array<Scalar, typename detail::TailNth<sizeof...(Idx), Shape>::type > part(Idx... idxs) const
+ typename std::enable_if<detail::TailNth<sizeof...(Idx), Shape>::type::fixed,
+ Array<Scalar, typename detail::TailNth<sizeof...(Idx), Shape>::type > >::type
+ part(Idx... idxs) const
{
- constexpr size_t nidxs = sizeof...(idxs);
+ constexpr size_t nidxs = sizeof...(Idx);
static_assert(nidxs < Shape::ndim,
"number of indices must be small than the number of dimensions");
- const std::array<size_t, Shape::ndim> m{idxs...};
- std::array<size_t, Shape::ndim - nidxs> shape;
- size_t offset = offset_;
+ size_t offset = index(idxs...);
+ return {data_, offset};
+ }
- for (size_t i = 0; i < nidxs; ++i) {
- if (Shape::stride(i) != Dynamic) {
- offset += Shape::stride(i) * m[i];
- } else {
- offset += stride_[i] * m[i];
- }
-#ifndef NO_BOUNDS_CHECK
- if (m[i] >= shape_[i])
- throw std::out_of_range("index out of bounds");
-#endif
- }
+ /*! Extract partial matrix (compile-time dynamic size) */
+ template <typename... Idx>
+ typename std::enable_if<!detail::TailNth<sizeof...(Idx), Shape>::type::fixed,
+ Array<Scalar, typename detail::TailNth<sizeof...(Idx), Shape>::type > >::type
+ part(Idx... idxs) const
+ {
+ constexpr size_t nidxs = sizeof...(Idx);
+ static_assert(nidxs < Shape::ndim,
+ "number of indices must be small than the number of dimensions");
+ size_t offset = index(idxs...);
+ std::array<size_t, Shape::ndim - nidxs> shape;
for (size_t i = nidxs; i < Shape::ndim; ++i)
- shape[i - nidxs] = shape_[i];
+ shape[i - nidxs] = Base::dim(i);
return {data_, shape, offset};
}
@@ -228,8 +293,8 @@ public:
EigenMap matrix() const
{
static_assert(Shape::ndim == 2, "matrix must be two-dimensional");
- return {(Scalar *)data_.data() + offset_,
- static_cast<Eigen::Index>(shape<0>()), static_cast<Eigen::Index>(shape<1>())};
+ return {(Scalar *)data_.data() + Base::offset(),
+ static_cast<Eigen::Index>(Base::dim(0)), static_cast<Eigen::Index>(Base::dim(1))};
}
operator EigenMap() const { return matrix(); }
diff --git a/src/main.cpp b/src/main.cpp
index 867f7f7d95800b2ac7a38993f1f71e0e490ac21e..8efa8a7dc404e5725eaf14191ca3ddfa75ba7270 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -6,10 +6,11 @@ using namespace std::complex_literals;
template <typename Scalar, typename Shape>
void dump(Array<Scalar,Shape> Q)
+#if 1
{
- for (size_t i = 0; i < Q.template shape<0>(); ++i) {
- for (size_t j = 0; j < Q.template shape<1>(); ++j) {
- for (size_t k = 0; k < Q.template shape<2>(); ++k) {
+ for (size_t i = 0; i < Q.dim(0); ++i) {
+ for (size_t j = 0; j < Q.dim(1); ++j) {
+ for (size_t k = 0; k < Q.dim(2); ++k) {
std::cout << i << "," << j << "," << k << " = " << Q(i,j,k) << std::endl;
}
}
@@ -19,28 +20,28 @@ void dump(Array<Scalar,Shape> Q)
std::cout << Q.data()[i] << std::endl;
}
}
+#else
+;
+#endif
int main()
{
- std::vector<ADComplex> v(2*2*2);
- Array<ADComplex,Shape<Dynamic,2,2> > Q0(v, {2,2,2});
- auto Q = Q0.part(1u);
-
- std::cout << Q0.index(1u,0u,0u) << std::endl;
+ std::vector<double> v(2*2*2);
+ Array<double,Shape<Dynamic,2,2> > Q0(v,{2,2,2});
+ auto Q = Q0.part(0u);
S_2(Q);
- std::cout << Q.index(0u,0u) << std::endl;
+ dump(Q0);
auto mat = Q.matrix();
- mat = mat * mat * mat * adcomplex(1j);
- mat(0,0) = mat(1,1) = mat.trace();
+ mat = mat * mat * mat;
- dump(Q0);
- std::cout << mat << std::endl;
+ if constexpr(true) {
+ }
- static_assert(std::is_same<decltype(mat), Eigen::Map<Eigen::Matrix<ADComplex, 2, 2, Eigen::RowMajor> > >::value, "wrong type");
+ dump(Q0);
return 0;
}