basic_array.hpp

Go to the documentation of this file.

// Copyright (c) 2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
// Copyright (c) 2018 Centre national de la recherche scientifique (CNRS)
// Copyright (c) 2018-2023 Simons Foundation
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0.txt
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// Authors: Olivier Parcollet, Nils Wentzell
 
/**
 * @file
 * @brief Provides the generic class for arrays.
 */
 
#pragma once
 
#include "./accessors.hpp"
#include "./basic_array_view.hpp"
#include "./basic_functions.hpp"
#include "./concepts.hpp"
#include "./exceptions.hpp"
#include "./iterators.hpp"
#include "./layout/for_each.hpp"
#include "./layout/permutation.hpp"
#include "./layout/range.hpp"
#include "./layout_transforms.hpp"
#include "./macros.hpp"
#include "./mem/address_space.hpp"
#include "./mem/memcpy.hpp"
#include "./mem/policies.hpp"
#include "./stdutil/array.hpp"
#include "./traits.hpp"
 
#include <algorithm>
#include <array>
#include <complex>
#include <concepts>
#include <initializer_list>
#include <random>
#include <ranges>
#include <type_traits>
#include <utility>
 
#ifdef NDA_ENFORCE_BOUNDCHECK
#include <exception>
#include <iostream>
#endif // NDA_ENFORCE_BOUNDCHECK
 
namespace nda {
 
  /**
   * @ingroup arrays_views
   * @brief A generic multi-dimensional array.
   *
   * @details Together with nda::basic_array_view, this class forms the backbone of the **nda** library. It is templated
   * with the following parameters:
   *
   * - `ValueType`: This is the type of the elements stored in the array. Most of the time, this will be a scalar type
   * like an int, double or std::complex<double>, but it can also be a more complex type like a custom class or a
   * another nda::basic_array.
   * - `Rank`: Integer specifying the number of dimensions of the array. This is a compile-time constant.
   * - `LayoutPolicy`: The layout policy specifies how the array views the memory it uses and how it accesses its
   * elements. It provides a mapping from multi-dimensional to linear indices and vice versa (see @ref layout_pols).
   * - `Algebra`: The algebra specifies how an array behaves when it is used in an expression. Possible values are 'A'
   * (array), 'M' (matrix) and 'V' (vector) (see nda::get_algebra).
   * - `ContainerPolicy`: The container policy specifies how and where the data is stored. It is responsible for
   * allocating/deallocating the memory (see @ref mem_handles).
   *
   * In contrast to views (see nda::basic_array_view), regular arrays own the memory they use for data storage.
   *
   * @code{.cpp}
   * // create a regular 3x2 array of ones
   * auto arr = nda::ones<int>(3, 2);
   * std::cout << arr << std::endl;
   *
   * // assign the value 42 to the first row
   * arr(0, nda::ellipsis{}) = 42;
   * std::cout << arr << std::endl;
   * @endcode
   *
   * Output:
   *
   * @code{bash}
   *
   * [[1,1]
   *  [1,1]
   *  [1,1]]
   *
   * [[42,42]
   *  [1,1]
   *  [1,1]]
   * @endcode
   *
   * Arrays and views share a lot of the same operations and functionalities. To turn a view into a regular array, use
   * nda::make_regular.
   *
   * @tparam ValueType Type stored in the array.
   * @tparam Rank Number of dimensions of the array.
   * @tparam LayoutPolicy Policy determining the memory layout.
   * @tparam Algebra Algebra of the array.
   * @tparam ContainerPolicy Policy determining how and where the data is stored.
   */
  template <typename ValueType, int Rank, typename LayoutPolicy, char Algebra, typename ContainerPolicy>
  class basic_array {
    // Compile-time checks.
    static_assert(!std::is_const_v<ValueType>, "Error in nda::basic_array: ValueType cannot be const");
    static_assert((Algebra != 'N'), "Internal error in nda::basic_array: Algebra 'N' not supported");
    static_assert((Algebra != 'M') or (Rank == 2), "Internal error in nda::basic_array: Algebra 'M' requires a rank 2 array");
    static_assert((Algebra != 'V') or (Rank == 1), "Internal error in nda::basic_array: Algebra 'V' requires a rank 1 array");
 
    public:
    /// Type of the values in the array (can not be const).
    using value_type = ValueType;
 
    /// Type of the memory layout policy (see @ref layout_pols).
    using layout_policy_t = LayoutPolicy;
 
    /// Type of the memory layout (an nda::idx_map).
    using layout_t = typename LayoutPolicy::template mapping<Rank>;
 
    /// Type of the container policy (see @ref mem_pols).
    using container_policy_t = ContainerPolicy;
 
    /// Type of the memory handle (see @ref mem_handles).
    using storage_t = typename ContainerPolicy::template handle<ValueType>;
 
    /// The associated regular type.
    using regular_type = basic_array;
 
    /// Number of dimensions of the array.
    static constexpr int rank = Rank;
 
    // Compile-time check.
    static_assert(has_contiguous(layout_t::layout_prop), "Error in nda::basic_array: Memory layout has to be contiguous");
 
    private:
    // Type of the array itself.
    using self_t = basic_array;
 
    // Type of the accessor policy for views (no no_alias_accessor).
    using AccessorPolicy = default_accessor;
 
    // Type of the owning policy for views.
    using OwningPolicy = borrowed<storage_t::address_space>;
 
    // Constexpr variable that is true if the value_type is const (never for basic_array).
    static constexpr bool is_const = false;
 
    // Constexpr variable that is true if the array is a view (never for basic_array).
    static constexpr bool is_view = false;
 
    // Memory layout of the array, i.e. the nda::idx_map.
    layout_t lay;
 
    // Memory handle of the array.
    storage_t sto;
 
    // Construct an array with a given shape and initialize the memory with zeros.
    template <std::integral Int = long>
    basic_array(std::array<Int, Rank> const &shape, mem::init_zero_t) : lay{shape}, sto{lay.size(), mem::init_zero} {}
 
    public:
    /**
     * @brief Convert the current array to a view with an 'A' (array) algebra.
     * @return An nda::basic_array_view of the current array.
     */
    auto as_array_view() { return basic_array_view<ValueType, Rank, LayoutPolicy, 'A', AccessorPolicy, OwningPolicy>{*this}; };
 
    /**
     * @brief Convert the current array to a view with an 'A' (array) algebra.
     * @return An nda::basic_array_view of the current array with const value type.
     */
    auto as_array_view() const { return basic_array_view<const ValueType, Rank, LayoutPolicy, 'A', AccessorPolicy, OwningPolicy>{*this}; };
 
    /// @deprecated Create the transpose of a 2-dimensional array. Use nda::transpose instead.
    [[deprecated]] auto transpose()
      requires(Rank == 2)
    {
      return permuted_indices_view<encode(std::array<int, 2>{1, 0})>(*this);
    }
 
    /// @deprecated Create the transpose of a 2-dimensional array. Use nda::transpose instead.
    [[deprecated]] auto transpose() const
      requires(Rank == 2)
    {
      return permuted_indices_view<encode(std::array<int, 2>{1, 0})>(*this);
    }
 
    /// Default constructor constructs an empty array with a default constructed memory handle and layout.
    basic_array() {}; // NOLINT (user-defined constructor to avoid value initialization of the sso buffer)
 
    /// Default move constructor moves the memory handle and layout.
    basic_array(basic_array &&) = default;
 
    /// Default copy constructor copies the memory handle and layout.
    explicit basic_array(basic_array const &a) = default;
 
    /**
     * @brief Construct an array from another array with a different algebra and/or container policy.
     *
     * @details It takes the memory layout of the other array and copies the data from the memory handle.
     *
     * @tparam A Algebra of the other array.
     * @tparam CP Container policy of the other array.
     * @param a Other array.
     */
    template <char A, typename CP>
    explicit basic_array(basic_array<ValueType, Rank, LayoutPolicy, A, CP> a) noexcept : lay(a.indexmap()), sto(std::move(a.storage())) {}
 
    /**
     * @brief Construct an array with the given dimensions.
     *
     * @details The integer type must be convertible to long and there must be exactly `Rank` arguments. It depends on
     * the value type and the container policy whether the data is initialized with zeros or not.
     *
     * @tparam Ints Integer types.
     * @param is Extent (number of elements) along each dimension.
     */
    template <std::integral... Ints>
      requires(sizeof...(Ints) == Rank)
    explicit basic_array(Ints... is) {
      // setting the layout and storage in the constructor body improves error messages for wrong # of args
      lay = layout_t{std::array{long(is)...}}; // NOLINT (for better error messages)
      sto = storage_t{lay.size()};             // NOLINT (for better error messages)
    }
 
    /**
     * @brief Construct a 1-dimensional array with the given size and initialize each element to the given scalar value.
     *
     * @tparam Int Integer type.
     * @tparam RHS Type of the scalar value to initialize the array with.
     * @param sz Size of the array.
     * @param val Value to initialize the array with.
     */
    template <std::integral Int, typename RHS>
    explicit basic_array(Int sz, RHS const &val)
      requires((Rank == 1 and is_scalar_for_v<RHS, basic_array>))
       : lay(layout_t{std::array{long(sz)}}), sto{lay.size()} {
      assign_from_scalar(val);
    }
 
    /**
     * @brief Construct an array with the given shape.
     *
     * @details It depends on the value type and the container policy whether the data is initialized with zeros or not.
     *
     * @tparam Int Integer type.
     * @param shape Shape of the array.
     */
    template <std::integral Int = long>
    explicit basic_array(std::array<Int, Rank> const &shape)
      requires(std::is_default_constructible_v<ValueType>)
       : lay(shape), sto(lay.size()) {}
 
    /**
     * @brief Construct an array with the given memory layout.
     *
     * @details It depends on the value type and the container policy whether the data is initialized with zeros or not.
     *
     * @param layout Memory layout.
     */
    explicit basic_array(layout_t const &layout)
      requires(std::is_default_constructible_v<ValueType>)
       : lay{layout}, sto{lay.size()} {}
 
    /**
     * @brief Construct an array with the given memory layout and with an existing memory handle/storage.
     *
     * @details The memory handle is moved and the layout is copied into the array.
     *
     * @param layout Memory layout.
     * @param storage Memory handle/storage.
     */
    explicit basic_array(layout_t const &layout, storage_t &&storage) noexcept : lay{layout}, sto{std::move(storage)} {}
 
    /**
     * @brief Construct an array from an nda::ArrayOfRank object with the same rank by copying each element.
     *
     * @tparam A nda::ArrayOfRank type.
     * @param a nda::ArrayOfRank object.
     */
    template <ArrayOfRank<Rank> A>
      requires(HasValueTypeConstructibleFrom<A, value_type>)
    basic_array(A const &a) : lay(a.shape()), sto{lay.size(), mem::do_not_initialize} {
      static_assert(std::is_constructible_v<value_type, get_value_t<A>>, "Error in nda::basic_array: Incompatible value types in constructor");
      if constexpr (std::is_trivial_v<ValueType> or is_complex_v<ValueType>) {
        // trivial and complex value types can use the optimized assign_from_ndarray
        assign_from_ndarray(a);
      } else {
        // general value types may not be default constructible -> use placement new
        nda::for_each(lay.lengths(), [&](auto const &...is) { new (sto.data() + lay(is...)) ValueType{a(is...)}; });
      }
    }
 
    /**
     * @brief Construct an array from an nda::ArrayInitializer object.
     *
     * @details An nda::ArrayInitializer is typically returned by delayed operations (see e.g. nda::mpi_gather).
     * The constructor can then be used to create the resulting array.
     *
     * @tparam Initializer nda::ArrayInitializer type.
     * @param initializer nda::ArrayInitializer object.
     */
    template <ArrayInitializer<basic_array> Initializer> // can not be explicit
    basic_array(Initializer const &initializer) : basic_array{initializer.shape()} {
      initializer.invoke(*this);
    }
 
    private:
    // Get the corresponding shape from an initializer list.
    static std::array<long, 1> shape_from_init_list(std::initializer_list<ValueType> const &l) noexcept { return {long(l.size())}; }
 
    // Get the corresponding shape from a nested initializer list.
    template <typename L>
    static auto shape_from_init_list(std::initializer_list<L> const &l) noexcept {
      const auto [min, max] = std::minmax_element(std::begin(l), std::end(l), [](auto &&x, auto &&y) { return x.size() < y.size(); });
      EXPECTS_WITH_MESSAGE(min->size() == max->size(),
                           "Error in nda::basic_array: Arrays can only be initialized with rectangular initializer lists");
      return stdutil::front_append(shape_from_init_list(*max), long(l.size()));
    }
 
    public:
    /**
     * @brief Construct a 1-dimensional array from an initializer list.
     * @param l Initializer list.
     */
    basic_array(std::initializer_list<ValueType> const &l)
      requires(Rank == 1)
       : lay(std::array<long, 1>{long(l.size())}), sto{lay.size(), mem::do_not_initialize} {
      long i = 0;
      for (auto const &x : l) { new (sto.data() + lay(i++)) ValueType{x}; }
    }
 
    /**
     * @brief Construct a 2-dimensional array from a double nested initializer list.
     * @param l2 Initializer list.
     */
    basic_array(std::initializer_list<std::initializer_list<ValueType>> const &l2)
      requires(Rank == 2)
       : lay(shape_from_init_list(l2)), sto{lay.size(), mem::do_not_initialize} {
      long i = 0, j = 0;
      for (auto const &l1 : l2) {
        for (auto const &x : l1) { new (sto.data() + lay(i, j++)) ValueType{x}; }
        j = 0;
        ++i;
      }
    }
 
    /**
     * @brief Construct a 3-dimensional array from a triple nested initializer list.
     * @param l3 Initializer list.
     */
    basic_array(std::initializer_list<std::initializer_list<std::initializer_list<ValueType>>> const &l3)
      requires(Rank == 3)
       : lay(shape_from_init_list(l3)), sto{lay.size(), mem::do_not_initialize} {
      long i = 0, j = 0, k = 0;
      for (auto const &l2 : l3) {
        for (auto const &l1 : l2) {
          for (auto const &x : l1) { new (sto.data() + lay(i, j, k++)) ValueType{x}; }
          k = 0;
          ++j;
        }
        j = 0;
        ++i;
      }
    }
 
    /**
     * @brief Construct a 2-dimensional array from another 2-dimensional array with a different algebra.
     *
     * @details The given array is moved into the constructed array. Note that for nda::stack_array or other array
     * types, this might still involve a copy of the data.
     *
     * @tparam A2 Algebra of the given array.
     * @param a Other array.
     */
    template <char A2>
    explicit basic_array(basic_array<ValueType, 2, LayoutPolicy, A2, ContainerPolicy> &&a) noexcept
      requires(Rank == 2)
       : basic_array{a.indexmap(), std::move(a).storage()} {}
 
    /**
     * @brief Make a zero-initialized array with the given shape.
     *
     * @tparam Int Integer type.
     * @param shape Shape of the array.
     * @return Zero-initialized array.
     */
    template <std::integral Int = long>
    static basic_array zeros(std::array<Int, Rank> const &shape)
      requires(std::is_standard_layout_v<ValueType> && std::is_trivially_copyable_v<ValueType>)
    {
      return basic_array{stdutil::make_std_array<long>(shape), mem::init_zero};
    }
 
    /**
     * @brief Make a zero-initialized array with the given dimensions.
     *
     * @tparam Ints Integer types.
     * @param is Extent (number of elements) along each dimension.
     * @return Zero-initialized array.
     */
    template <std::integral... Ints>
    static basic_array zeros(Ints... is)
      requires(sizeof...(Ints) == Rank)
    {
      return zeros(std::array<long, Rank>{is...});
    }
 
    /**
     * @brief Make a one-initialized array with the given shape.
     *
     * @tparam Int Integer type.
     * @param shape Shape of the array.
     * @return One-initialized array.
     */
    template <std::integral Int = long>
    static basic_array ones(std::array<Int, Rank> const &shape)
      requires(nda::is_scalar_v<ValueType>)
    {
      auto res = basic_array{stdutil::make_std_array<long>(shape)};
      res()    = ValueType{1};
      return res;
    }
 
    /**
     * @brief Make a one-initialized array with the given dimensions.
     *
     * @tparam Ints Integer types.
     * @param is Extent (number of elements) along each dimension.
     * @return One-initialized array.
     */
    template <std::integral... Ints>
    static basic_array ones(Ints... is)
      requires(sizeof...(Ints) == Rank)
    {
      return ones(std::array<long, Rank>{is...});
    }
 
    /**
     * @brief Make a random-initialized array with the given shape.
     *
     * @details The random values are take from a uniform distribution over [0, 1). For a complex array, both real and
     * imaginary parts are initialized with random values.
     *
     * @tparam Int Integer type.
     * @param shape Shape of the array.
     * @return Random-initialized array.
     */
    template <std::integral Int = long>
    static basic_array rand(std::array<Int, Rank> const &shape)
      requires(std::is_floating_point_v<ValueType> or nda::is_complex_v<ValueType>)
    {
      using namespace std::complex_literals;
      auto static gen  = std::mt19937(std::random_device{}());
      auto static dist = std::uniform_real_distribution<>(0.0, 1.0);
      auto res         = basic_array{shape};
      if constexpr (nda::is_complex_v<ValueType>)
        for (auto &x : res) x = dist(gen) + 1i * dist(gen);
      else
        for (auto &x : res) x = dist(gen);
      return res;
    }
 
    /**
     * @brief Make a random-initialized array with the given dimensions.
     *
     * @details The random values are take from a uniform distribution over [0, 1). For a complex array, both real and
     * imaginary parts are initialized with random values.
     *
     * @tparam Ints Integer types.
     * @param is Extent (number of elements) along each dimension.
     * @return Random-initialized array.
     */
    template <std::integral... Ints>
    static basic_array rand(Ints... is)
      requires(sizeof...(Ints) == Rank)
    {
      return rand(std::array<long, Rank>{is...});
    }
 
    /// Default move assignment moves the memory handle and layout from the right hand side array.
    basic_array &operator=(basic_array &&) = default;
 
    /// Default copy assignment copies the memory handle and layout from the right hand side array.
    basic_array &operator=(basic_array const &) = default;
 
    /**
     * @brief Assignment operator makes a deep copy of another array with a different algebra and/or container policy.
     *
     * @details A new array object is constructed from the right hand side array and then moved into the current array.
     * This will invalidate all references/views to the existing storage.
     *
     * @tparam A Algebra of the other array.
     * @tparam CP Container policy of the other array.
     * @param rhs Right hand side of the assignment operation.
     */
    template <char A, typename CP>
    basic_array &operator=(basic_array<ValueType, Rank, LayoutPolicy, A, CP> const &rhs) {
      *this = basic_array{rhs};
      return *this;
    }
 
    /**
     * @brief Assignment operator makes a deep copy of an nda::ArrayOfRank object.
     *
     * @details The array is first resized to the shape of the right hand side and the elements are copied. This might
     * invalidate all references/views to the existing storage.
     *
     * @tparam RHS nda::ArrayOfRank type.
     * @param rhs Right hand side of the assignment operation.
     */
    template <ArrayOfRank<Rank> RHS>
    basic_array &operator=(RHS const &rhs) {
      resize(rhs.shape());
      assign_from_ndarray(rhs);
      return *this;
    }
 
    /**
     * @brief Assignment operator assigns a scalar to the array.
     *
     * @details The behavior depends on the algebra of the array:
     * - 'A' (array) and 'V' (vector): The scalar is assigned to all elements of the array.
     * - 'M' (matrix): The scalar is assigned to the diagonal elements of the shorter dimension.
     *
     * @tparam RHS Type of the scalar.
     * @param rhs Right hand side of the assignment operation.
     */
    template <typename RHS>
    basic_array &operator=(RHS const &rhs) noexcept
      requires(is_scalar_for_v<RHS, basic_array>)
    {
      assign_from_scalar(rhs);
      return *this;
    }
 
    /**
     * @brief Assignment operator uses an nda::ArrayInitializer to assign to the array.
     *
     * @details The array is resized to the shape of the initializer. This might invalidate all references/views to the
     * existing storage.
     *
     * @tparam Initializer nda::ArrayInitializer type.
     * @param initializer Initializer object.
     */
    template <ArrayInitializer<basic_array> Initializer>
    basic_array &operator=(Initializer const &initializer) {
      resize(initializer.shape());
      initializer.invoke(*this);
      return *this;
    }
 
    /**
     * @brief Resize the array to a new shape.
     *
     * @details Resizing is only performed if the storage is not null and if the new size is different from the previous
     * size. If resizing is performed, the content of the resulting array is undefined since it makes no copy of the
     * previous data and all references/views to the existing storage will be invalidated.
     *
     * @tparam Ints Integer types.
     * @param is New extent (number of elements) along each dimension.
     */
    template <std::integral... Ints>
    void resize(Ints const &...is) {
      static_assert(std::is_copy_constructible_v<ValueType>, "Error in nda::basic_array: Resizing requires the value_type to be copy constructible");
      static_assert(sizeof...(is) == Rank, "Error in nda::basic_array: Resizing requires exactly Rank arguments");
      resize(std::array<long, Rank>{long(is)...});
    }
 
    /**
     * @brief Resize the array to a new shape.
     *
     * @details Resizing is only performed if the storage is not null and if the new size is different from the previous
     * size. If resizing is performed, the content of the resulting array is undefined since it makes no copy of the
     * previous data and all references/views to the existing storage will be invalidated.
     *
     * @param shape New shape of the array.
     */
    [[gnu::noinline]] void resize(std::array<long, Rank> const &shape) {
      lay = layout_t(shape);
      if (sto.is_null() or (sto.size() != lay.size())) sto = storage_t{lay.size()};
    }
 
// include common functionality of arrays and views
// Copyright (c) 2019-2024 Simons Foundation
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0.txt
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// Authors: Thomas Hahn, Miguel Morales, Olivier Parcollet, Nils Wentzell
 
/**
 * @brief Get the memory layout of the view/array.
 * @return nda::idx_map specifying the layout of the view/array.
 */
[[nodiscard]] constexpr auto const &indexmap() const noexcept { return lay; }
 
/**
 * @brief Get the data storage of the view/array.
 * @return A const reference to the memory handle of the view/array.
 */
[[nodiscard]] storage_t const &storage() const & noexcept { return sto; }
 
/**
 * @brief Get the data storage of the view/array.
 * @return A reference to the memory handle of the view/array.
 */
[[nodiscard]] storage_t &storage() & noexcept { return sto; }
 
/**
 * @brief Get the data storage of the view/array.
 * @return A copy of the memory handle of the view/array.
 */
[[nodiscard]] storage_t storage() && noexcept { return std::move(sto); }
 
/**
 * @brief Get the stride order of the memory layout of the view/array (see nda::idx_map for more details on how we
 * define stride orders).
 *
 * @return `std::array<int, Rank>` object specifying the stride order.
 */
[[nodiscard]] constexpr auto stride_order() const noexcept { return lay.stride_order; }
 
/**
 * @brief Get a pointer to the actual data (in general this is not the beginning of the memory block for a view).
 * @return Const pointer to the first element of the view/array.
 */
[[nodiscard]] ValueType const *data() const noexcept { return sto.data(); }
 
/**
 * @brief Get a pointer to the actual data (in general this is not the beginning of thr memory block for a view).
 * @return Pointer to the first element of the view/array.
 */
[[nodiscard]] ValueType *data() noexcept { return sto.data(); }
 
/**
 * @brief Get the shape of the view/array.
 * @return `std::array<long, Rank>` object specifying the shape of the view/array.
 */
[[nodiscard]] auto const &shape() const noexcept { return lay.lengths(); }
 
/**
 * @brief Get the strides of the view/array (see nda::idx_map for more details on how we define strides).
 * @return `std::array<long, Rank>` object specifying the strides of the view/array.
 */
[[nodiscard]] auto const &strides() const noexcept { return lay.strides(); }
 
/**
 * @brief Get the total size of the view/array.
 * @return Number of elements contained in the view/array.
 */
[[nodiscard]] long size() const noexcept { return lay.size(); }
 
/**
 * @brief Is the memory layout of the view/array contiguous?
 * @return True if the nda::idx_map is contiguous, false otherwise.
 */
[[nodiscard]] long is_contiguous() const noexcept { return lay.is_contiguous(); }
 
/**
 * @brief Is the view/array empty?
 * @return True if the view/array does not contain any elements.
 */
[[nodiscard]] bool empty() const { return sto.is_null(); }
 
/// @deprecated Use empty() instead.
[[nodiscard]] bool is_empty() const noexcept { return sto.is_null(); }
 
/**
 * @brief Get the extent of the i<sup>th</sup> dimension.
 * @return Number of elements along the i<sup>th</sup> dimension.
 */
[[nodiscard]] long extent(int i) const noexcept {
#ifdef NDA_ENFORCE_BOUNDCHECK
  if (i < 0 || i >= rank) {
    std::cerr << "Error in extent: Dimension " << i << " is incompatible with array of rank " << rank << std::endl;
    std::terminate();
  }
#endif
  return lay.lengths()[i];
}
 
/// @deprecated Use `extent(i)` or `shape()[i]` instead.
[[nodiscard]] long shape(int i) const noexcept { return extent(i); }
 
/**
 * @brief Get a range that generates all valid index tuples.
 * @return An `itertools::multiplied` range that can be used to iterate over all valid index tuples.
 */
[[nodiscard]] auto indices() const noexcept { return itertools::product_range(shape()); }
 
/**
 * @brief Is the stride order of the view/array in C-order?
 * @return True if the stride order of the nda::idx_map is C-order, false otherwise.
 */
static constexpr bool is_stride_order_C() noexcept { return layout_t::is_stride_order_C(); }
 
/**
 * @brief Is the stride order of the view/array in Fortran-order?
 * @return True if the stride order of the nda::idx_map is Fortran-order, false otherwise.
 */
static constexpr bool is_stride_order_Fortran() noexcept { return layout_t::is_stride_order_Fortran(); }
 
/**
 * @brief Access the element of the view/array at the given nda::_linear_index_t.
 *
 * @details The linear index specifies the position of the element in the view/array and not the position of the
 * element w.r.t. to the data pointer (i.e. any possible strides should not be taken into account).
 *
 * @param idx nda::_linear_index_t object.
 * @return Const reference to the element at the given linear index.
 */
decltype(auto) operator()(_linear_index_t idx) const noexcept {
  if constexpr (layout_t::layout_prop == layout_prop_e::strided_1d)
    return sto[idx.value * lay.min_stride()];
  else if constexpr (layout_t::layout_prop == layout_prop_e::contiguous)
    return sto[idx.value];
  else
    static_assert(always_false<layout_t>, "Internal error in array/view: Calling this type with a _linear_index_t is not allowed");
}
 
/// Non-const overload of @ref nda::basic_array_view::operator()(_linear_index_t) const.
decltype(auto) operator()(_linear_index_t idx) noexcept {
  if constexpr (layout_t::layout_prop == layout_prop_e::strided_1d)
    return sto[idx.value * lay.min_stride()];
  else if constexpr (layout_t::layout_prop == layout_prop_e::contiguous)
    return sto[idx.value];
  else
    static_assert(always_false<layout_t>, "Internal error in array/view: Calling this type with a _linear_index_t is not allowed");
}
 
private:
// Constexpr variable that is true if bounds checking is disabled.
#ifdef NDA_ENFORCE_BOUNDCHECK
static constexpr bool has_no_boundcheck = false;
#else
static constexpr bool has_no_boundcheck = true;
#endif
 
public:
/**
 * @brief Implementation of the function call operator.
 *
 * @details This function is an implementation detail an should be private. Since the Green's function library in
 * TRIQS uses this function, it is kept public (for now).
 *
 * @tparam ResultAlgebra Algebra of the resulting view/array.
 * @tparam SelfIsRvalue True if the view/array is an rvalue.
 * @tparam Self Type of the calling view/array.
 * @tparam T Types of the arguments.
 *
 * @param self Calling view.
 * @param idxs Multi-dimensional index consisting of `long`, `nda::range`, `nda::range::all_t`, nda::ellipsis or lazy
 * arguments.
 * @return Result of the function call depending on the given arguments and type of the view/array.
 */
template <char ResultAlgebra, bool SelfIsRvalue, typename Self, typename... Ts>
FORCEINLINE static decltype(auto) call(Self &&self, Ts const &...idxs) noexcept(has_no_boundcheck) {
  // resulting value type
  using r_v_t = std::conditional_t<std::is_const_v<std::remove_reference_t<Self>>, ValueType const, ValueType>;
 
  // behavior depends on the given arguments
  if constexpr (clef::is_any_lazy<Ts...>) {
    // if there are lazy arguments, e.g. as in A(i_) << i_, a lazy expression is returned
    return clef::make_expr_call(std::forward<Self>(self), idxs...);
  } else if constexpr (sizeof...(Ts) == 0) {
    // if no arguments are given, a full view is returned
    return basic_array_view<r_v_t, Rank, LayoutPolicy, Algebra, AccessorPolicy, OwningPolicy>{self.lay, self.sto};
  } else {
    // otherwise we check the arguments and either access a single element or make a slice
    static_assert(((layout_t::template argument_is_allowed_for_call_or_slice<Ts> + ...) > 0),
                  "Error in array/view: Slice arguments must be convertible to range, ellipsis, or long (or string if the layout permits it)");
 
    // number of arguments convertible to long
    static constexpr int n_args_long = (layout_t::template argument_is_allowed_for_call<Ts> + ...);
 
    if constexpr (n_args_long == rank) {
      // access a single element
      long offset = self.lay(idxs...);
      if constexpr (is_view or not SelfIsRvalue) {
        // if the calling object is a view or an lvalue, we return a reference
        return AccessorPolicy::template accessor<ValueType>::access(self.sto.data(), offset);
      } else {
        // otherwise, we return a copy of the value
        return ValueType{self.sto[offset]};
      }
    } else {
      // access a slice of the view/array
      auto const [offset, idxm]      = self.lay.slice(idxs...);
      static constexpr auto res_rank = decltype(idxm)::rank();
      // resulting algebra
      static constexpr char newAlgebra = (ResultAlgebra == 'M' and (res_rank == 1) ? 'V' : ResultAlgebra);
      // resulting layout policy
      using r_layout_p = typename detail::layout_to_policy<std::decay_t<decltype(idxm)>>::type;
      return basic_array_view<ValueType, res_rank, r_layout_p, newAlgebra, AccessorPolicy, OwningPolicy>{std::move(idxm), {self.sto, offset}};
    }
  }
}
 
public:
/**
 * @brief Function call operator to access the view/array.
 *
 * @details Depending on the type of the calling object and the given arguments, this function call does the following:
 * - If any of the arguments is lazy, an nda::clef::expr with the nda::clef::tags::function tag is returned.
 * - If no arguments are given, a full view of the calling object is returned:
 *   - If the calling object itself or its value type is const, a view with a const value type is returned.
 *   - Otherwise, a view with a non-const value type is returned.
 * - If the number of arguments is equal to the rank of the calling object and all arguments are convertible to `long`,
 * a single element is accessed:
 *   - If the calling object is a view or an lvalue, a (const) reference to the element is returned.
 *   - Otherwise, a copy of the element is returned.
 * - Otherwise a slice of the calling object is returned with the same value type and accessor and owning policies as
 * the calling object. The algebra of the slice is the same as well, except if a 1-dimensional slice of a matrix is
 * taken. In this case, the algebra is changed to 'V'.
 *
 * @tparam Ts Types of the function arguments.
 * @param idxs Multi-dimensional index consisting of `long`, `nda::range`, `nda::range::all_t`, nda::ellipsis or lazy
 * arguments.
 * @return Result of the function call depending on the given arguments and type of the view/array.
 */
template <typename... Ts>
FORCEINLINE decltype(auto) operator()(Ts const &...idxs) const & noexcept(has_no_boundcheck) {
  static_assert((rank == -1) or (sizeof...(Ts) == rank) or (sizeof...(Ts) == 0) or (ellipsis_is_present<Ts...> and (sizeof...(Ts) <= rank + 1)),
                "Error in array/view: Incorrect number of parameters in call operator");
  return call<Algebra, false>(*this, idxs...);
}
 
/// Non-const overload of `nda::basic_array_view::operator()(Ts const &...) const &`.
template <typename... Ts>
FORCEINLINE decltype(auto) operator()(Ts const &...idxs) & noexcept(has_no_boundcheck) {
  static_assert((rank == -1) or (sizeof...(Ts) == rank) or (sizeof...(Ts) == 0) or (ellipsis_is_present<Ts...> and (sizeof...(Ts) <= rank + 1)),
                "Error in array/view: Incorrect number of parameters in call operator");
  return call<Algebra, false>(*this, idxs...);
}
 
/// Rvalue overload of `nda::basic_array_view::operator()(Ts const &...) const &`.
template <typename... Ts>
FORCEINLINE decltype(auto) operator()(Ts const &...idxs) && noexcept(has_no_boundcheck) {
  static_assert((rank == -1) or (sizeof...(Ts) == rank) or (sizeof...(Ts) == 0) or (ellipsis_is_present<Ts...> and (sizeof...(Ts) <= rank + 1)),
                "Error in array/view: Incorrect number of parameters in call operator");
  return call<Algebra, true>(*this, idxs...);
}
 
/**
 * @brief Subscript operator to access the 1-dimensional view/array.
 *
 * @details Depending on the type of the calling object and the given argument, this subscript operation does the
 * following:
 * - If the argument is lazy, an nda::clef::expr with the nda::clef::tags::function tag is returned.
 * - If the argument is convertible to `long`, a single element is accessed:
 *   - If the calling object is a view or an lvalue, a (const) reference to the element is returned.
 *   - Otherwise, a copy of the element is returned.
 * - Otherwise a slice of the calling object is returned with the same value type, algebra and accessor and owning
 * policies as the calling object.
 *
 * @tparam T Type of the argument.
 * @param idx 1-dimensional index that is either a `long`, `nda::range`, `nda::range::all_t`, nda::ellipsis or a lazy
 * argument.
 * @return Result of the subscript operation depending on the given argument and type of the view/array.
 */
template <typename T>
decltype(auto) operator[](T const &idx) const & noexcept(has_no_boundcheck) {
  static_assert((rank == 1), "Error in array/view: Subscript operator is only available for rank 1 views/arrays in C++17/20");
  return call<Algebra, false>(*this, idx);
}
 
/// Non-const overload of `nda::basic_array_view::operator[](T const &) const &`.
template <typename T>
decltype(auto) operator[](T const &x) & noexcept(has_no_boundcheck) {
  static_assert((rank == 1), "Error in array/view: Subscript operator is only available for rank 1 views/arrays in C++17/20");
  return call<Algebra, false>(*this, x);
}
 
/// Rvalue overload of `nda::basic_array_view::operator[](T const &) const &`.
template <typename T>
decltype(auto) operator[](T const &x) && noexcept(has_no_boundcheck) {
  static_assert((rank == 1), "Error in array/view: Subscript operator is only available for rank 1 views/arrays in C++17/20");
  return call<Algebra, true>(*this, x);
}
 
/// Rank of the nda::array_iterator for the view/array.
static constexpr int iterator_rank = (has_strided_1d(layout_t::layout_prop) ? 1 : Rank);
 
/// Const iterator type of the view/array.
using const_iterator = array_iterator<iterator_rank, ValueType const, typename AccessorPolicy::template accessor<ValueType>::pointer>;
 
/// Iterator type of the view/array.
using iterator = array_iterator<iterator_rank, ValueType, typename AccessorPolicy::template accessor<ValueType>::pointer>;
 
private:
// Make an iterator for the view/array depending on its type.
template <typename Iterator>
[[nodiscard]] auto make_iterator(bool at_end) const noexcept {
  if constexpr (iterator_rank == Rank) {
    // multi-dimensional iterator
    if constexpr (layout_t::is_stride_order_C()) {
      // C-order case (array_iterator already traverses the data in C-order)
      return Iterator{indexmap().lengths(), indexmap().strides(), sto.data(), at_end};
    } else {
      // general case (we need to permute the shape and the strides according to the stride order of the layout)
      return Iterator{nda::permutations::apply(layout_t::stride_order, indexmap().lengths()),
                      nda::permutations::apply(layout_t::stride_order, indexmap().strides()), sto.data(), at_end};
    }
  } else {
    // 1-dimensional iterator
    return Iterator{std::array<long, 1>{size()}, std::array<long, 1>{indexmap().min_stride()}, sto.data(), at_end};
  }
}
 
public:
/// Get a const iterator to the beginning of the view/array.
[[nodiscard]] const_iterator begin() const noexcept { return make_iterator<const_iterator>(false); }
 
/// Get a const iterator to the beginning of the view/array.
[[nodiscard]] const_iterator cbegin() const noexcept { return make_iterator<const_iterator>(false); }
 
/// Get an iterator to the beginning of the view/array.
iterator begin() noexcept { return make_iterator<iterator>(false); }
 
/// Get a const iterator to the end of the view/array.
[[nodiscard]] const_iterator end() const noexcept { return make_iterator<const_iterator>(true); }
 
/// Get a const iterator to the end of the view/array.
[[nodiscard]] const_iterator cend() const noexcept { return make_iterator<const_iterator>(true); }
 
/// Get an iterator to the end of the view/array.
iterator end() noexcept { return make_iterator<iterator>(true); }
 
/**
 * @brief Addition assignment operator.
 *
 * @details It first performs the (lazy) addition with the right hand side operand and then assigns the result to the
 * left hand side view/array.
 *
 * See nda::operator+(L &&, R &&) and nda::operator+(A &&, S &&) for more details.
 *
 * @tparam RHS nda::Scalar or nda::Array type.
 * @param rhs Right hand side operand of the addition assignment operation.
 * @return Reference to this object.
 */
template <typename RHS>
auto &operator+=(RHS const &rhs) noexcept {
  static_assert(not is_const, "Error in array/view: Can not assign to a const view");
  return operator=(*this + rhs);
}
 
/**
 * @brief Subtraction assignment operator.
 *
 * @details It first performs the (lazy) subtraction with the right hand side operand and then assigns the result to
 * the left hand side view/array.
 *
 * See nda::operator-(L &&, R &&) and nda::operator-(A &&, S &&) for more details.
 *
 * @tparam RHS nda::Scalar or nda::Array type.
 * @param rhs Right hand side operand of the subtraction assignment operation.
 * @return Reference to this object.
 */
template <typename RHS>
auto &operator-=(RHS const &rhs) noexcept {
  static_assert(not is_const, "Error in array/view: Can not assign to a const view");
  return operator=(*this - rhs);
}
 
/**
 * @brief Multiplication assignment operator.
 *
 * @details It first performs the (lazy) multiplication with the right hand side operand and then assigns the result
 * to the left hand side view/array.
 *
 * See nda::operator*(L &&, R &&) and nda::operator*(A &&, S &&) for more details.
 *
 * @tparam RHS nda::Scalar or nda::Array type.
 * @param rhs Right hand side operand of the multiplication assignment operation.
 * @return Reference to this object.
 */
template <typename RHS>
auto &operator*=(RHS const &rhs) noexcept {
  static_assert(not is_const, "Error in array/view: Can not assign to a const view");
  return operator=((*this) * rhs);
}
 
/**
 * @brief Division assignment operator.
 *
 * @details It first performs the (lazy) division with the right hand side operand and then assigns the result to the
 * left hand side view/array.
 *
 * See nda::operator/(L &&, R &&) and nda::operator/(A &&, S &&) for more details.
 *
 * @tparam RHS nda::Scalar or nda::Array type.
 * @param rhs Right hand side operand of the division assignment operation.
 * @return Reference to this object.
 */
template <typename RHS>
auto &operator/=(RHS const &rhs) noexcept {
  static_assert(not is_const, "Error in array/view: Can not assign to a const view");
  return operator=(*this / rhs);
}
 
/**
 * @brief Assignment operator makes a deep copy of a general contiguous range and assigns it to the 1-dimensional
 * view/array.
 *
 * @tparam R Range type.
 * @param rhs Right hand side range object.
 * @return Reference to this object.
 */
template <std::ranges::contiguous_range R>
auto &operator=(R const &rhs) noexcept
  requires(Rank == 1 and not MemoryArray<R>)
{
  *this = basic_array_view{rhs};
  return *this;
}
 
private:
// Implementation of the assignment from an n-dimensional array type.
template <typename RHS>
void assign_from_ndarray(RHS const &rhs) { // FIXME noexcept {
#ifdef NDA_ENFORCE_BOUNDCHECK
  if (this->shape() != rhs.shape())
    NDA_RUNTIME_ERROR << "Error in assign_from_ndarray: Size mismatch:"
                      << "\n LHS.shape() = " << this->shape() << "\n RHS.shape() = " << rhs.shape();
#endif
  // compile-time check if assignment is possible
  static_assert(std::is_assignable_v<value_type &, get_value_t<RHS>>, "Error in assign_from_ndarray: Incompatible value types");
 
  // are both operands nda::MemoryArray types?
  static constexpr bool both_in_memory = MemoryArray<self_t> and MemoryArray<RHS>;
 
  // do both operands have the same stride order?
  static constexpr bool same_stride_order = get_layout_info<self_t>.stride_order == get_layout_info<RHS>.stride_order;
 
  // prefer optimized options if possible
  if constexpr (both_in_memory and same_stride_order) {
    if (rhs.empty()) return;
    // are both operands strided in 1d?
    static constexpr bool both_1d_strided = has_layout_strided_1d<self_t> and has_layout_strided_1d<RHS>;
    if constexpr (mem::on_host<self_t, RHS> and both_1d_strided) {
      // vectorizable copy on host
      for (long i = 0; i < size(); ++i) (*this)(_linear_index_t{i}) = rhs(_linear_index_t{i});
      return;
    } else if constexpr (!mem::on_host<self_t, RHS> and have_same_value_type_v<self_t, RHS>) {
      // check for block-layout and use mem::memcpy2D if possible
      auto bl_layout_dst = get_block_layout(*this);
      auto bl_layout_src = get_block_layout(rhs);
      if (bl_layout_dst && bl_layout_src) {
        auto [n_bl_dst, bl_size_dst, bl_str_dst] = *bl_layout_dst;
        auto [n_bl_src, bl_size_src, bl_str_src] = *bl_layout_src;
        // check that the total memory size is the same
        if (n_bl_dst * bl_size_dst != n_bl_src * bl_size_src) NDA_RUNTIME_ERROR << "Error in assign_from_ndarray: Incompatible block sizes";
        // if either destination or source consists of a single block, we can chunk it up to make the layouts compatible
        if (n_bl_dst == 1 && n_bl_src > 1) {
          n_bl_dst = n_bl_src;
          bl_size_dst /= n_bl_src;
          bl_str_dst = bl_size_dst;
        }
        if (n_bl_src == 1 && n_bl_dst > 1) {
          n_bl_src = n_bl_dst;
          bl_size_src /= n_bl_dst;
          bl_str_src = bl_size_src;
        }
        // copy only if block-layouts are compatible, otherwise continue to fallback
        if (n_bl_dst == n_bl_src && bl_size_dst == bl_size_src) {
          mem::memcpy2D<mem::get_addr_space<self_t>, mem::get_addr_space<RHS>>((void *)data(), bl_str_dst * sizeof(value_type), (void *)rhs.data(),
                                                                               bl_str_src * sizeof(value_type), bl_size_src * sizeof(value_type),
                                                                               n_bl_src);
          return;
        }
      }
    }
  }
  // otherwise fallback to elementwise assignment
  if constexpr (mem::on_device<self_t> || mem::on_device<RHS>) {
    NDA_RUNTIME_ERROR << "Error in assign_from_ndarray: Fallback to elementwise assignment not implemented for arrays/views on the GPU";
  }
  nda::for_each(shape(), [this, &rhs](auto const &...args) { (*this)(args...) = rhs(args...); });
}
 
// Implementation to fill a view/array with a constant scalar value.
template <typename Scalar>
void fill_with_scalar(Scalar const &scalar) noexcept {
  // we make a special implementation if the array is strided in 1d or contiguous
  if constexpr (has_layout_strided_1d<self_t>) {
    const long L             = size();
    auto *__restrict const p = data(); // no alias possible here!
    if constexpr (has_contiguous_layout<self_t>) {
      for (long i = 0; i < L; ++i) p[i] = scalar;
    } else {
      const long stri  = indexmap().min_stride();
      const long Lstri = L * stri;
      for (long i = 0; i < Lstri; i += stri) p[i] = scalar;
    }
  } else {
    // no compile-time memory layout guarantees
    for (auto &x : *this) x = scalar;
  }
}
 
// Implementation of the assignment from a scalar value.
template <typename Scalar>
void assign_from_scalar(Scalar const &scalar) noexcept {
  static_assert(!is_const, "Error in assign_from_ndarray: Cannot assign to a const view");
  if constexpr (Algebra != 'M') {
    // element-wise assignment for non-matrix algebras
    fill_with_scalar(scalar);
  } else {
    // a scalar has to be interpreted as a unit matrix for matrix algebras (the scalar in the shortest diagonal)
    // FIXME : A priori faster to put 0 everywhere and then change the diag to avoid the if.
    // FIXME : Benchmark and confirm.
    if constexpr (is_scalar_or_convertible_v<Scalar>)
      fill_with_scalar(0);
    else
      fill_with_scalar(Scalar{0 * scalar}); // FIXME : improve this
    const long imax = std::min(extent(0), extent(1));
    for (long i = 0; i < imax; ++i) operator()(i, i) = scalar;
  }
}
  };
 
  // Class template argument deduction guides.
  template <MemoryArray A>
  basic_array(A &&a) -> basic_array<get_value_t<A>, get_rank<A>, get_contiguous_layout_policy<get_rank<A>, get_layout_info<A>.stride_order>,
                                    get_algebra<A>, heap<mem::get_addr_space<A>>>;
 
  template <Array A>
  basic_array(A &&a) -> basic_array<get_value_t<A>, get_rank<A>, C_layout, get_algebra<A>, heap<>>;
 
} // namespace nda