diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 79fd5bb6d6..0b2043a723 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -32,6 +32,7 @@ TiledArray/bitset.h
TiledArray/block_range.h
TiledArray/dense_shape.h
TiledArray/dist_array.h
+TiledArray/dist_array_vector.h
TiledArray/distributed_storage.h
TiledArray/error.h
TiledArray/external/madness.h
diff --git a/src/TiledArray/dist_array_vector.h b/src/TiledArray/dist_array_vector.h
new file mode 100644
index 0000000000..f5fe9fb1f5
--- /dev/null
+++ b/src/TiledArray/dist_array_vector.h
@@ -0,0 +1,823 @@
+/*
+ * This file is a part of TiledArray.
+ * Copyright (C) 2013 Virginia Tech
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see .
+ *
+ */
+
+#ifndef TILEDARRAY_DIST_ARRAY_VECTOR_H__INCLUDED
+#define TILEDARRAY_DIST_ARRAY_VECTOR_H__INCLUDED
+
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+namespace TiledArray {
+
+namespace detail {
+
+/// A vector of objects that, if there is only one object, can be implicitly
+/// converted to it
+
+/// @tparam T the type of contained objects; must be copy constructable and
+/// move constructable, does not have to be copy or move assignable
+template
+class VectorBase : public container::svector {
+ public:
+ using base_type = container::svector;
+
+ VectorBase() = default;
+ /// copy ctor delegates to ctor that uses emplace_back
+ VectorBase(const VectorBase &other)
+ : VectorBase(static_cast(other)) {}
+ /// move ctor delegates to ctor that uses emplace_back
+ VectorBase(VectorBase &&other)
+ : VectorBase(static_cast(other)) {}
+
+ /// constructs a vector of default-constructed objects
+ /// @param[in] sz the target size of the vector
+ explicit VectorBase(std::size_t sz) : base_type(sz) {}
+
+ /// copies objects from a std::vector
+ /// @param v a std::vector from which the objects will be copied
+ explicit VectorBase(const std::vector &v) {
+ this->reserve(v.size());
+ for (auto &&e : v) {
+ this->emplace_back(e);
+ }
+ }
+
+ /// moves objects from a std::vector
+ /// @param v a std::vector from which the objects will be moved
+ explicit VectorBase(std::vector &&v) {
+ this->reserve(v.size());
+ for (auto &&e : v) {
+ this->emplace_back(std::move(e));
+ }
+ }
+
+ /// copies objects from a container::svector
+ /// @param v a container::svector from which the objects will be copied
+ template
+ explicit VectorBase(const container::svector &v) {
+ this->reserve(v.size());
+ for (auto &&e : v) {
+ this->emplace_back(e);
+ }
+ }
+
+ /// moves objects from a container::svector
+ /// @param v a container::svector from which the objects will be moved
+ template
+ explicit VectorBase(container::svector &&v) {
+ this->reserve(v.size());
+ for (auto &&e : v) {
+ this->emplace_back(std::move(e));
+ }
+ }
+
+ template
+ VectorBase(B &&begin, E &&end)
+ : container::svector{std::forward(begin), std::forward(end)} {}
+
+ /// copy assignment
+ VectorBase &operator=(const VectorBase &other) {
+ this->clear();
+ if (this->size() < other.size()) this->reserve(other.size());
+ for (auto &&e : other) {
+ this->emplace_back(e);
+ }
+ return *this;
+ }
+
+ /// move assignment
+ VectorBase &operator=(VectorBase &&other) {
+ this->clear();
+ if (this->size() < other.size()) this->reserve(other.size());
+ for (auto &&e : other) {
+ this->emplace_back(std::move(e));
+ }
+ return *this;
+ }
+
+ /// converts to a nonconst lvalue ref to the first object is there is only one
+ /// @throw TiledArray::Exception if the number of objects is not equal to 1
+ operator T &() & {
+ if (this->size() == 1)
+ return this->operator[](0);
+ else
+ throw TiledArray::Exception(
+ "casting Vector to T only possible if it contains 1 "
+ "T");
+ }
+
+ /// converts to a const lvalue ref to the first object is there is only one
+ /// @throw TiledArray::Exception if the number of objects is not equal to 1
+ operator const T &() const & {
+ if (this->size() == 1)
+ return this->operator[](0);
+ else
+ throw TiledArray::Exception(
+ "casting Vector to T only possible if it contains 1 "
+ "T");
+ }
+
+ /// converts to a nonconst rvalue ref to the first object is there is only one
+ /// @throw TiledArray::Exception if the number of objects is not equal to 1
+ operator T &&() && {
+ if (this->size() == 1)
+ return std::move(this->operator[](0));
+ else
+ throw TiledArray::Exception(
+ "casting Vector to T only possible if it contains 1 "
+ "T");
+ }
+
+ /// @return true if nonempty
+ explicit operator bool() const { return !this->empty(); }
+}; // VectorBase
+
+} // namespace detail
+
+/// a vector of Shape objects
+template
+class ShapeVector : public detail::VectorBase {
+ public:
+ using base_type = detail::VectorBase;
+ using T = typename base_type::value_type;
+ using value_type = typename Shape_::value_type;
+
+ ShapeVector() = default;
+
+ ShapeVector(std::size_t i) : base_type(i) {}
+
+ template , T> && ...)>>
+ ShapeVector(Args &&..._t) : base_type{std::forward(_t)...} {}
+
+ ShapeVector(const std::vector &v) : base_type(v) {}
+
+ ShapeVector(std::vector &&v) : base_type(std::move(v)) {}
+
+ template
+ ShapeVector(const container::svector &v) : base_type(v) {}
+
+ template
+ ShapeVector(container::svector &&v) : base_type(std::move(v)) {}
+
+ template
+ ShapeVector(B &&begin, E &&end)
+ : base_type(std::forward(begin), std::forward(end)) {}
+
+ template
+ ShapeVector &operator=(const S &shape) {
+ static_cast(*this) = shape;
+ return *this;
+ }
+
+ using base_type::operator T &;
+ using base_type::operator const T &;
+ using base_type::operator T &&;
+
+ /// @return true if nonempty and each element is nonnull (i.e., is
+ /// initialized)
+ using base_type::operator bool;
+
+ template
+ ShapeVector transform(const Op &op) {
+ ShapeVector result;
+ result.reserve(this->size());
+ for (auto &&v : *this) result.emplace_back(v.transform(op));
+ return result;
+ }
+
+ /// @return minimum value of threshold of elements
+ value_type threshold() {
+ value_type thresh = (*this)[0].threshold();
+ for (auto &&v : *this) {
+ value_type thresh_val = v.threshold();
+ if (thresh_val < thresh) thresh = thresh_val;
+ }
+ return thresh;
+ }
+};
+
+namespace expressions {
+
+template
+using enable_if_expression = std::void_t::engine_type>;
+
+template
+struct is_tsr_expression : public std::false_type {};
+template
+struct is_tsr_expression> : public std::true_type {};
+template
+constexpr const bool is_tsr_expression_v = is_tsr_expression::value;
+
+/// a vector of Expr objects
+template >
+class ExprVector : public TiledArray::detail::VectorBase {
+ public:
+ using base_type = TiledArray::detail::VectorBase;
+ using T = typename base_type::value_type;
+
+ ExprVector() = default;
+ ExprVector(const ExprVector &) = default;
+ ExprVector(ExprVector &&) = default;
+
+ ExprVector(std::size_t i) : base_type(i) {}
+
+ template , T> && ...)>>
+ ExprVector(Args &&..._t) : base_type{std::forward(_t)...} {}
+
+ ExprVector(const std::vector &v) : base_type(v) {}
+
+ ExprVector(std::vector &&v) : base_type(std::move(v)) {}
+
+ template
+ ExprVector(const container::svector &v) : base_type(v) {}
+
+ template
+ ExprVector(container::svector &&v) : base_type(std::move(v)) {}
+
+ template
+ ExprVector(B &&begin, E &&end)
+ : base_type(std::forward(begin), std::forward(end)) {}
+
+ // Apply unary operator- to each member of ExprVector
+ ExprVector operator-() const {
+ auto expr_sz = (*this).size();
+ ExprVector result;
+ for (auto i = 0; i < expr_sz; ++i) {
+ result[i] = -(*this)[i];
+ }
+ return result;
+ }
+
+ template
+ ExprVector &operator=(const ExprVector &expr) {
+ const auto expr_sz = expr.size();
+ if (this->size() ==
+ 0) { // if null, and this is TsrExpr, resize the underlying
+ // DistArrayVector to match the size of expr and rebuild this
+ if constexpr (is_tsr_expression_v) {
+ TA_ASSERT(arrayvec_ptr_);
+ TA_ASSERT(!arrayvec_annotation_.empty());
+ arrayvec_ptr_->resize(expr_sz);
+ static_cast(*this) =
+ static_cast((*arrayvec_ptr_)(arrayvec_annotation_));
+ }
+ } else { // if this is nonull must match the size of expr
+ TA_ASSERT(this->size() == expr_sz);
+ }
+ for (size_t i = 0; i != expr_sz; ++i) {
+ (*this)[i] = expr[i];
+ }
+ return *this;
+ }
+
+ template
+ ExprVector &operator=(ExprVector &&expr) {
+ return this->operator=(static_cast &>(expr));
+ }
+
+ ExprVector &operator=(ExprVector &&expr) {
+ return this->operator=(std::move(expr));
+ }
+ ExprVector &operator=(const ExprVector &expr) {
+ return this->operator=(expr);
+ }
+
+ template >>
+ ExprVector &operator=(E &&expr) {
+ container::svector> vec;
+ vec.reserve(1);
+ vec.emplace_back(std::forward(expr));
+ return this->operator=(
+ ExprVector>(std::move(vec)));
+ }
+
+ using base_type::operator T &;
+ using base_type::operator const T &;
+ using base_type::operator T &&;
+
+ /// @return true if nonempty and each element is nonnull (i.e., is
+ /// initialized)
+ using base_type::operator bool;
+
+#define TA_EXPRVEC_UNARYREDUCTION_DEF(op) \
+ auto op(World &world) const { \
+ container::svector result; \
+ result.reserve(this->size()); \
+ for (auto &&v : (*this)) result.emplace_back(v.op(world)); \
+ return result; \
+ } \
+ auto op() const { return this->op(this->default_world()); }
+
+ TA_EXPRVEC_UNARYREDUCTION_DEF(norm)
+
+#undef TA_EXPRVEC_UNARYREDUCTION_DEF
+
+#define TA_EXPRVEC_BINARYREDUCTION_DEF(op) \
+ template > \
+ auto op(const E &right_expr, World &world) const { \
+ container::svector result; \
+ result.reserve(this->size()); \
+ for (auto &&v : (*this)) result.emplace_back(v.op(right_expr, world)); \
+ return result; \
+ } \
+ template > \
+ auto op(const E &right_expr) const { \
+ return this->op(right_expr, this->default_world()); \
+ } \
+ template > \
+ auto op(const ExprVector &right_expr, World &world) const { \
+ container::svector result; \
+ const auto sz = this->size(); \
+ TA_ASSERT(sz == right_expr.size()); \
+ result.reserve(sz); \
+ for (size_t i = 0; i != sz; ++i) \
+ result.emplace_back((*this)[i].op(right_expr[i], world)); \
+ return result; \
+ } \
+ template > \
+ auto op(const ExprVector &right_expr) const { \
+ return this->op(right_expr, this->default_world()); \
+ }
+
+ TA_EXPRVEC_BINARYREDUCTION_DEF(dot)
+
+#undef TA_EXPRVEC_BINARYREDUCTION_DEF
+
+ World &default_world() const {
+ if (*this) {
+ if constexpr (has_array::value)
+ return (*this)[0].array().world();
+ else
+ return TiledArray::get_default_world();
+ } else
+ throw TiledArray::Exception(
+ "ExprVector::default_world() called on null object");
+ }
+
+ /// calls \c set_shape(shape) on each member
+ /// \param[in] shape the shape object to use for each member of this
+ /// \return reference to this
+ ExprVector &set_shape(typename EngineParamOverride<
+ typename Expr_::engine_type>::shape_type const &shape) {
+ for (auto &&v : (*this)) {
+ v.set_shape(shape);
+ }
+ return *this;
+ }
+
+ /// calls \c set_shape(shape[i]) on member \c i
+ /// \param[in] shape the vector of shapes
+ /// \return reference to this
+ ExprVector &set_shape(
+ ShapeVector::shape_type> const &shapes) {
+ auto sz = this->size();
+ TA_ASSERT(sz == shapes.size());
+ size_t i = 0;
+ for (auto &&v : (*this)) {
+ v.set_shape(shapes[i]);
+ ++i;
+ }
+ return *this;
+ }
+
+ /// Expression plus-assignment operator
+
+ /// \tparam D The derived expression type
+ /// \param other The expression that will be added to this array
+ template
+ ExprVector &operator+=(const ExprVector &other) {
+ const auto sz = this->size();
+ TA_ASSERT(sz == other.size());
+ size_t i = 0;
+ for (auto &&v : *this) {
+ v += other[i];
+ ++i;
+ }
+ return *this;
+ }
+
+ /// Expression conjugate applied to each element
+ ExprVector &conj() {
+ size_t i = 0;
+ for (auto &&v : *this) {
+ v.conj();
+ ++i;
+ }
+ return *this;
+ }
+
+ private:
+ // for TsrExpr type need to keep ptr to the bound DistArrayVector
+ template
+ constexpr static decltype(auto) arrayvec_ptr_init() {
+ if constexpr (is_tsr_expression_v) {
+ TiledArray::DistArrayVector<
+ typename ExprTrait::array_type::value_type,
+ typename ExprTrait::array_type::policy_type> *result = nullptr;
+ return result;
+ } else
+ return nullptr;
+ }
+ decltype(arrayvec_ptr_init()) arrayvec_ptr_ =
+ arrayvec_ptr_init();
+ std::string arrayvec_annotation_;
+
+ template
+ friend class TiledArray::DistArrayVector;
+
+ void set_arrayvec_ptr(void *ptr, const std::string &annotation) {
+ if constexpr (is_tsr_expression_v) {
+ arrayvec_ptr_ = static_cast())>(ptr);
+ arrayvec_annotation_ = annotation;
+ }
+ }
+};
+
+template
+auto operator*(const double &factor, const ExprVector &expr) {
+ const auto expr_sz = expr.size();
+ ExprVector result;
+ result.reserve(expr_sz);
+ for (size_t i = 0; i != expr_sz; ++i) {
+ result.emplace_back(factor * expr[i]);
+ }
+ return result;
+}
+
+template >
+const Expr &to_base_expr(const ExprVector &e) {
+ return static_cast &>(static_cast(e));
+}
+
+// template
+// ExprVector> operator*(const ExprVector& left,
+// const Expr& right);
+
+#define TA_EXPRVEC_BINARYOP_DEF(op, result_type) \
+ template \
+ ExprVector> operator op( \
+ const ExprVector &left, const Expr &right) { \
+ const auto sz = left.size(); \
+ ExprVector> result; \
+ result.reserve(sz); \
+ for (size_t i = 0; i != sz; ++i) result.emplace_back(left[i] op right); \
+ return result; \
+ } \
+ template \
+ ExprVector> operator op( \
+ const Expr &left, const ExprVector &right) { \
+ const auto sz = right.size(); \
+ ExprVector> result; \
+ result.reserve(sz); \
+ for (size_t i = 0; i != sz; ++i) result.emplace_back(left op right[i]); \
+ return result; \
+ } \
+ template \
+ ExprVector> operator op( \
+ const ExprVector &left, const ExprVector &right) { \
+ if (left.size() == 1 && right.size() > 1) { \
+ const auto sz = right.size(); \
+ ExprVector> result; \
+ result.reserve(sz); \
+ for (size_t i = 0; i != sz; ++i) \
+ result.emplace_back(left[0] op right[i]); \
+ return result; \
+ } else if (left.size() > 1 && right.size() == 1) { \
+ const auto sz = left.size(); \
+ ExprVector> result; \
+ result.reserve(sz); \
+ for (size_t i = 0; i != sz; ++i) \
+ result.emplace_back(left[i] op right[0]); \
+ return result; \
+ } else { \
+ TA_ASSERT(left.size() == right.size()); \
+ const auto sz = left.size(); \
+ ExprVector> result; \
+ result.reserve(sz); \
+ for (size_t i = 0; i != sz; ++i) \
+ result.emplace_back(left[i] op right[i]); \
+ return result; \
+ } \
+ }
+
+// fwd declare expressions
+template
+class AddExpr;
+template
+class SubtExpr;
+template
+class MultExpr;
+
+TA_EXPRVEC_BINARYOP_DEF(+, AddExpr)
+TA_EXPRVEC_BINARYOP_DEF(-, SubtExpr)
+TA_EXPRVEC_BINARYOP_DEF(*, MultExpr)
+#undef TA_EXPRVEC_BINARYOP_DEF
+
+} // namespace expressions
+
+/// A vector of DistArray objects
+
+/// Each DistArray object is distributed over the same World, i.e. the vector is
+/// not distributed. If you need a vector distributed across a world, with each
+/// DistArray object belonging to a local world, see
+/// TiledArray/vector_of_arrays.h .
+/// @internal Rationale: sometimes it is more convenient to represent a sequence
+/// (vector) of arrays
+/// as a vector, rather than a "fused" array with an extra dimension.
+/// The latter makes sense when operations are homogeneous across the
+/// extra dimension. If extents of modes will depend on the vector
+/// index then fusing the arrays becomes cumbersome (N.B.: if the
+/// "inner" arrays do not need to be distributed then you should use a
+/// tensor-of-tensor, i.e. DistArray>>)
+/// @tparam T a DistArray type
+template
+class DistArrayVector : public detail::VectorBase> {
+ public:
+ using base_type = detail::VectorBase>;
+ using array_type = TiledArray::DistArray;
+ using eval_type = typename array_type::eval_type;
+ using policy_type = Policy;
+ using value_type = array_type;
+ using element_type = typename array_type::element_type;
+ using scalar_type = typename array_type::scalar_type;
+ using T = array_type;
+
+ DistArrayVector() = default;
+
+ DistArrayVector(std::size_t i) : base_type(i) {}
+
+ template , T> && ...)>>
+ DistArrayVector(Args &&..._t) : base_type(std::forward(_t)...) {}
+
+ DistArrayVector(const std::vector &t_pack) : base_type(t_pack) {}
+
+ DistArrayVector(std::vector &&t_pack) : base_type(std::move(t_pack)) {}
+
+ template
+ DistArrayVector(const container::svector &v) : base_type(v) {}
+
+ template
+ DistArrayVector(container::svector &&v) : base_type(std::move(v)) {}
+
+ template
+ DistArrayVector(B &&begin, E &&end)
+ : base_type(std::forward(begin), std::forward(end)) {}
+
+ using base_type::operator T &;
+ using base_type::operator const T &;
+ using base_type::operator T &&;
+
+ /// @return true is nonempty and each element is nonnull (i.e., is
+ /// initialized)
+ explicit operator bool() const {
+ bool result = false;
+ if (this->size()) {
+ result = true;
+ for (auto &&v : (*this)) {
+ result = result && static_cast(v);
+ }
+ }
+ return result;
+ }
+
+ /// @return a World to which all DistArrays belong
+ /// @throw TiledArray::Exception if this is null
+ World &world() const {
+ if (!(*this))
+ throw TiledArray::Exception(
+ "DistArrayVector::world() called on null object");
+ auto &w = (*this)[0].world();
+#ifndef NDEBUG
+ for (size_t v = 1; v < this->size(); ++v)
+ TA_ASSERT(w.id() == (*this)[v].world().id());
+#endif
+ return w;
+ }
+
+ /// @return a TiledRange object describing all objects
+ /// @throw TiledArray::Exception if this is null
+ const TiledRange &trange() const {
+ if (!(*this))
+ throw TiledArray::Exception(
+ "DistArrayVector::trange() called on null object");
+ auto &tr = (*this)[0].trange();
+#ifndef NDEBUG
+ for (size_t v = 1; v < this->size(); ++v)
+ TA_ASSERT(tr == (*this)[v].trange());
+#endif
+ return tr;
+ }
+
+ /// Calls truncate() on each element
+ /// @throw TiledArray::Exception if this is null
+ void truncate() {
+ if (!(*this))
+ throw TiledArray::Exception(
+ "DistArrayVector::truncate() called on null object");
+ for (size_t v = 1; v < this->size(); ++v) (*this)[v].truncate();
+ }
+
+ /// Calls is_initalized() on each element and return true if all
+ /// @return a bool true if all elemments are initialized
+ /// @throw TiledArray::Exception if this is null
+ bool is_initialized() {
+ if (!(*this))
+ throw TiledArray::Exception(
+ "DistArrayVector::is_initalized() called on null object");
+ for (auto &&v : *this) {
+ if (!v.is_initalized()) return false;
+ };
+ return true;
+ }
+
+ /// returns the shapes of the member arrays
+ /// @return vector of shapes
+ /// @throw TiledArray::Exception if this is null
+ ShapeVector shape() const {
+ if (!(*this))
+ throw TiledArray::Exception(
+ "DistArrayVector::shape() called on null object");
+ ShapeVector result;
+ result.reserve(this->size());
+ for (auto &&v : *this) result.emplace_back(v.shape());
+ return result;
+ }
+
+ /// Create a tensor expression
+
+ /// \param annotation An annotation string
+ /// \return A const tensor expression object
+ TiledArray::expressions::ExprVector<
+ TiledArray::expressions::TsrExpr>
+ operator()(const std::string &annotation) const {
+ TA_ASSERT(*this);
+ TiledArray::expressions::ExprVector<
+ TiledArray::expressions::TsrExpr>
+ result;
+ for (auto &&elem : (*this)) {
+ result.emplace_back(elem(annotation));
+ }
+ return result;
+ }
+
+ /// Create a tensor expression
+
+ /// \param annotation An annotation string
+ /// \return A non-const tensor expression object
+ TiledArray::expressions::ExprVector>
+ operator()(const std::string &annotation) {
+ TiledArray::expressions::ExprVector<
+ TiledArray::expressions::TsrExpr>
+ result;
+ if (!*this) result.set_arrayvec_ptr(this, annotation);
+ for (auto &&elem : (*this)) {
+ result.emplace_back(elem(annotation));
+ }
+ return result;
+ }
+
+ /// \return A DistArray that is the sum of elements in DistArrayVector
+ DistArray vector_sum() const {
+ DistArray result = (*this)[0];
+ if ((*this).size() > 1) {
+ for (auto i = 1; i < (*this).size(); ++i) {
+ result("i,j") += (*this)[i]("i,j");
+ }
+ }
+ return result;
+ }
+
+}; // DistArrayVector
+
+template
+auto rank(const DistArrayVector &a) {
+ MPQC_ASSERT(a.size() > 0);
+ size_t result = rank(a[0]);
+ for (size_t v = 1; v < a.size(); ++v) MPQC_ASSERT(result == rank(a[v]));
+ return result;
+}
+
+template
+size_t volume(const DistArrayVector &a) {
+ MPQC_ASSERT(a.size() > 0);
+ size_t result = volume(a[0]);
+ for (size_t v = 1; v < a.size(); ++v) MPQC_ASSERT(result == volume(a[v]));
+ return result;
+}
+
+template ::type>::value>::type>
+inline auto foreach_inplace(DistArrayVector &arg, Op &&op,
+ bool fence = true) {
+ const auto array_sz = arg.size();
+ for (size_t i = 0; i != array_sz; ++i) {
+ foreach_inplace(arg[i], op);
+ }
+}
+
+template
+inline auto norm2(const DistArrayVector &a) -> decltype(
+ norm2(std::declval::array_type>())) {
+ using scalar_type = typename DistArrayVector::scalar_type;
+ scalar_type norm2_squared = 0;
+ for (auto &t : a) {
+ norm2_squared += squared_norm(t);
+ }
+ return double(std::sqrt(norm2_squared));
+}
+
+template
+inline void zero(DistArrayVector &a) {
+ for (auto &t : a) {
+ TiledArray::math::linalg::zero(t);
+ }
+}
+
+template
+inline auto dot(const DistArrayVector &a,
+ const DistArrayVector &b) {
+ TA_ASSERT(a.size() == b.size());
+ typename DistArrayVector::scalar_type result = 0;
+ for (decltype(a.size()) i = 0; i != a.size(); ++i) {
+ result += dot(a[i], b[i]);
+ // it seems to have some random issue (hang queue error) when calling
+ // dot_product adding fence() seems to help prevent it
+ TiledArray::get_default_world().gop.fence();
+ }
+ return result;
+}
+
+template
+inline void axpy(DistArrayVector &y, Scalar a,
+ const DistArrayVector &x) {
+ TA_ASSERT(x.size() == y.size());
+ for (decltype(x.size()) i = 0; i != x.size(); ++i) {
+ TiledArray::math::linalg::axpy(y[i], a, x[i]);
+ }
+}
+
+// template
+// inline DistArrayVector copy(DistArrayVector &a) {
+// return a;
+// }
+
+template
+inline void scale(DistArrayVector &y, Scalar a) {
+ for (auto &t : y) {
+ TiledArray::math::linalg::scale(t, a);
+ }
+}
+
+} // namespace TiledArray
+
+namespace madness {
+
+template
+auto trace(const std::vector> &vec) {
+ T result = 0;
+ for (auto &&v : vec) result += v.get();
+ return result;
+}
+
+template
+auto trace(
+ const TiledArray::container::svector, Size> &vec) {
+ T result = 0;
+ for (auto &&v : vec) result += v.get();
+ return result;
+}
+} // namespace madness
+
+#endif // TILEDARRAY_DIST_ARRAY_VECTOR_H__INCLUDED
diff --git a/src/TiledArray/type_traits.h b/src/TiledArray/type_traits.h
index b6ce5ec3b4..2cdda7c4c0 100644
--- a/src/TiledArray/type_traits.h
+++ b/src/TiledArray/type_traits.h
@@ -90,6 +90,8 @@ class DensePolicy;
struct ZeroTensor;
template
class DistArray;
+template
+class DistArrayVector;
namespace detail {
template
class LazyArrayTile;
@@ -1192,6 +1194,15 @@ struct is_array> : public std::true_type {};
template
static constexpr bool is_array_v = is_array::value;
+template
+struct is_arrayvec : public std::false_type {};
+
+template
+struct is_arrayvec> : public std::true_type {};
+
+template
+static constexpr bool is_arrayvec_v = is_arrayvec::value;
+
template
using trange_t = typename T::trange_type;
diff --git a/src/tiledarray_fwd.h b/src/tiledarray_fwd.h
index af41f63c6f..79acaf0b19 100644
--- a/src/tiledarray_fwd.h
+++ b/src/tiledarray_fwd.h
@@ -47,15 +47,15 @@ class SparsePolicy;
template
class Tensor;
-typedef Tensor > TensorD;
-typedef Tensor > TensorI;
-typedef Tensor > TensorF;
-typedef Tensor > TensorL;
+typedef Tensor> TensorD;
+typedef Tensor> TensorI;
+typedef Tensor> TensorF;
+typedef Tensor> TensorL;
typedef Tensor,
- Eigen::aligned_allocator > >
+ Eigen::aligned_allocator>>
TensorZ;
typedef Tensor,
- Eigen::aligned_allocator > >
+ Eigen::aligned_allocator>>
TensorC;
// CUDA tensor
@@ -90,32 +90,36 @@ class DistArray;
// Dense Array Typedefs
template
-using TArray = DistArray >, DensePolicy>;
+using TArray = DistArray>, DensePolicy>;
typedef TArray TArrayD;
typedef TArray TArrayI;
typedef TArray TArrayF;
typedef TArray TArrayL;
-typedef TArray > TArrayZ;
-typedef TArray > TArrayC;
+typedef TArray> TArrayZ;
+typedef TArray> TArrayC;
// Sparse Array Typedefs
template
using TSpArray =
- DistArray >, SparsePolicy>;
+ DistArray>, SparsePolicy>;
typedef TSpArray TSpArrayD;
typedef TSpArray TSpArrayI;
typedef TSpArray TSpArrayF;
typedef TSpArray TSpArrayL;
-typedef TSpArray > TSpArrayZ;
-typedef TSpArray > TSpArrayC;
+typedef TSpArray> TSpArrayZ;
+typedef TSpArray> TSpArrayC;
// type alias for backward compatibility: the old Array has static type,
// DistArray is rank-polymorphic
template >,
+ typename Tile = Tensor>,
typename Policy = DensePolicy>
using Array = DistArray;
+// vector pack of DistArray's
+template
+class DistArrayVector;
+
} // namespace TiledArray
#ifndef TILEDARRAY_DISABLE_NAMESPACE_TA
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index d75c0f8e5f..260d8b1795 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -68,6 +68,7 @@ set(ta_test_src_files ta_test.cpp
index_list.cpp
bipartite_index_list.cpp
dist_array.cpp
+ dist_array_vector.cpp
conversions.cpp
eigen.cpp
dist_op_dist_cache.cpp
diff --git a/tests/array_fixture.h b/tests/array_fixture.h
index 259bd09852..44b96850ff 100644
--- a/tests/array_fixture.h
+++ b/tests/array_fixture.h
@@ -42,6 +42,17 @@ struct ArrayFixture : public TiledRangeFixture {
TiledArray::World& world;
ArrayN a;
SpArrayN b;
+
+ template
+ const auto& array() const {
+ if constexpr (std::is_same_v)
+ return b;
+ else if constexpr (std::is_same_v)
+ return a;
+ else
+ abort();
+ }
+
}; // struct ArrayFixture
#endif // TILEDARRAY_TEST_ARRAY_FIXTURE_H__INCLUDED
diff --git a/tests/dist_array_vector.cpp b/tests/dist_array_vector.cpp
new file mode 100644
index 0000000000..ded2f60895
--- /dev/null
+++ b/tests/dist_array_vector.cpp
@@ -0,0 +1,92 @@
+/*
+ * This file is a part of TiledArray.
+ * Copyright (C) 2021 Virginia Tech
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see .
+ *
+ */
+
+#include "../src/TiledArray/dist_array_vector.h"
+#include "unit_test_config.h"
+
+#include
+
+using namespace TiledArray;
+
+using ArrayVecN = DistArrayVector, DensePolicy>;
+using SpArrayVecN = DistArrayVector, SparsePolicy>;
+
+// These are all of the template parameters we are going to test over
+using tparams = boost::mpl::list;
+
+BOOST_FIXTURE_TEST_SUITE(arrayvec_suite, ArrayFixture)
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(constructors, ArrayVec, tparams) {
+ const auto& arr = array();
+
+ // default ctor
+ {
+ BOOST_REQUIRE_NO_THROW(ArrayVec{});
+ ArrayVec av;
+ BOOST_CHECK(!av); // empty vec is null
+ BOOST_CHECK_EQUAL(av.size(), 0);
+ }
+
+ // vector with null arrays
+ {
+ const auto sz = 2;
+ BOOST_REQUIRE_NO_THROW(ArrayVec{sz});
+ ArrayVec av{sz};
+ BOOST_CHECK(!av); // nonempty vec of null arrays is null
+ BOOST_CHECK_EQUAL(av.size(), sz);
+ }
+
+ // vector with nonnull arrays
+ {
+ const auto sz = 2;
+ ArrayVec av{sz};
+ BOOST_CHECK_EQUAL(av.size(), sz);
+ BOOST_REQUIRE_NO_THROW(av[0] = arr);
+ BOOST_REQUIRE_NO_THROW(av.at(1) = arr);
+ BOOST_CHECK(av); // nonempty vec of non-null arrays is nonnull
+ }
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(expressions, ArrayVec, tparams) {
+ const auto& arr = array();
+
+ // array -> arrayvec
+ {
+ ArrayVec av;
+ BOOST_REQUIRE_NO_THROW(av("i,j,k") = arr("i,j,k"));
+ BOOST_CHECK_EQUAL(av.size(), 1);
+ }
+
+ // arrayvec op array -> arrayvec
+ {
+ ArrayVec av0, av1, av2, av3;
+ av0("i,j,k") = arr("i,k,j");
+ // these don't compile
+ BOOST_REQUIRE_NO_THROW(av1("i,j,k") = av0("i,j,k") + arr("i,k,j"));
+ BOOST_CHECK_EQUAL(av1.size(), 1);
+ BOOST_REQUIRE_NO_THROW(av2("i,j,k") = arr("i,k,j") + av1("i,j,k"));
+ BOOST_CHECK_EQUAL(av2.size(), 1);
+ BOOST_REQUIRE_NO_THROW(av3("i,j,k") =
+ 3.0 * (1 * av2("i,j,k") - 2 * arr("i,k,j")) *
+ av1("i,j,k"));
+ BOOST_CHECK_EQUAL(av3.size(), 1);
+ }
+}
+
+BOOST_AUTO_TEST_SUITE_END()