basic_array_view.hpp

Go to the documentation of this file.

// Copyright (c) 2019--present, The Simons Foundation
// This file is part of TRIQS/nda and is licensed under the Apache License, Version 2.0.
// SPDX-License-Identifier: Apache-2.0
// See LICENSE in the root of this distribution for details.

 
#pragma once
 
#include "./basic_functions.hpp"
#include "./clef.hpp"
#include "./concepts.hpp"
#include "./declarations.hpp"
#include "./device.hpp"
#include "./iterators.hpp"
#include "./layout/for_each.hpp"
#include "./layout/idx_map.hpp"
#include "./layout/permutation.hpp"
#include "./layout/range.hpp"
#include "./layout/slice_static.hpp"
#include "./macros.hpp"
#include "./matrix_functions.hpp"
#include "./mem/address_space.hpp"
#include "./mem/fill.hpp"
#include "./mem/memcpy.hpp"
#include "./mem/memset.hpp"
#include "./mem/policies.hpp"
#include "./traits.hpp"
 
#include <itertools/itertools.hpp>
 
#include <algorithm>
#include <array>
#include <cstring>
#include <exception>
#include <iostream>
#include <memory>
#include <ranges>
#include <type_traits>
#include <utility>
 
namespace std {

  template <typename V1, int R1, typename LP1, char A1, typename AP1, typename OP1, typename V2, int R2, typename LP2, char A2, typename AP2,
            typename OP2>
  void swap(nda::basic_array_view<V1, R1, LP1, A1, AP1, OP1> &a, nda::basic_array_view<V2, R2, LP2, A2, AP2, OP2> &b) = delete;
 
} // namespace std
 
namespace nda {

  template <typename ValueType, int Rank, typename LayoutPolicy, char Algebra, typename AccessorPolicy, typename OwningPolicy>
  class basic_array_view {
    // Compile-time checks.
    static_assert((Algebra != 'N'), "Internal error in nda::basic_array_view: Algebra 'N' not supported");
    static_assert((Algebra != 'M') or (Rank == 2), "Internal error in nda::basic_array_view: Algebra 'M' requires a rank 2 view");
    static_assert((Algebra != 'V') or (Rank == 1), "Internal error in nda::basic_array_view: Algebra 'V' requires a rank 1 view");
 
    public:
    using value_type = ValueType;

    using layout_policy_t = LayoutPolicy;

    using layout_t = typename LayoutPolicy::template mapping<Rank>;

    using accessor_policy_t = AccessorPolicy;

    using owning_policy_t = OwningPolicy;

    using storage_t = typename OwningPolicy::template handle<ValueType>;

    using regular_type = basic_array<std::remove_const_t<ValueType>, Rank, C_layout, Algebra, heap<mem::get_addr_space<storage_t>>>;

    static constexpr int rank = Rank;
 
    // Compile-time check.
    static_assert(not mem::have_device_compatible_addr_space<storage_t> or nda::have_device,
                  "Error in nda::basic_array_view: Device compatible address space requires compiling with device support.");
 
    private:
    // Type of the view itself.
    using self_t = basic_array_view;
 
    // Constexpr variable that is true if the view is a view (always for basic_array_view).
    static constexpr bool is_view = true;
 
    // Constexpr variable that is true if the value_type is const.
    static constexpr bool is_const = std::is_const_v<ValueType>;
 
    // Memory layout of the view, i.e. the nda::idx_map.
    layout_t lay;
 
    // Memory handle of the view.
    storage_t sto;
 
    // Declare any nda::basic_array as a friend.
    template <typename T, int R, typename L, char A, typename CP>
    friend class basic_array;
 
    // Declare any other nda::basic_array_view as a friend.
    template <typename T, int R, typename L, char A, typename AP, typename OP>
    friend class basic_array_view;
 
    // Check the layout compatibility of this view with another layout.
    template <typename L>
    static constexpr bool requires_runtime_check = not layout_property_compatible(L::template mapping<Rank>::layout_prop, layout_t::layout_prop);
 
    public:
    // FIXME : TRIQS PORTING
    // private constructor for the previous friend
    basic_array_view(layout_t const &idxm, storage_t st) : lay(idxm), sto(std::move(st)) {}
 
    public:
    // backward : FIXME : temporary to be removed
    [[deprecated]] auto as_array_view() { return basic_array_view<ValueType, Rank, LayoutPolicy, 'A', AccessorPolicy, OwningPolicy>{*this}; };

    [[deprecated]] auto as_array_view() const {
      return basic_array_view<const ValueType, Rank, LayoutPolicy, 'A', AccessorPolicy, OwningPolicy>{*this};
    };

    basic_array_view() = default;

    basic_array_view(basic_array_view &&) = default;

    basic_array_view(basic_array_view const &) = default;

    template <MemoryArrayOfRank<Rank> A>
      requires((get_layout_info<A>.stride_order == layout_t::stride_order_encoded)
               and (std::is_same_v<std::remove_const_t<ValueType>, get_value_t<A>>)
               and (std::is_const_v<ValueType> or !std::is_const_v<typename std::decay_t<A>::value_type>))
    explicit(requires_runtime_check<typename std::decay_t<A>::layout_policy_t>)
       basic_array_view(A &&a) noexcept // NOLINT (should we forward the reference?)
       : lay(a.indexmap()), sto(a.storage()) {}

    basic_array_view(std::array<long, Rank> const &shape, ValueType *p) noexcept : basic_array_view(layout_t{shape}, p) {}

    basic_array_view(layout_t const &idxm, ValueType *p) noexcept : lay(idxm), sto{p} {}

    template <size_t N>
      requires(Rank == 1)
    explicit basic_array_view(std::array<ValueType, N> &a) noexcept : basic_array_view{{long(N)}, a.data()} {}

    template <size_t N>
      requires(Rank == 1 and std::is_const_v<ValueType>)
    explicit basic_array_view(std::array<std::remove_const_t<ValueType>, N> const &a) noexcept : basic_array_view{{long(N)}, a.data()} {}

    template <std::ranges::contiguous_range R>
      requires(Rank == 1 and not MemoryArray<R>
               and (std::is_same_v<std::ranges::range_value_t<R>, ValueType> or std::is_same_v<const std::ranges::range_value_t<R>, ValueType>))
    explicit basic_array_view(R &rg) noexcept : basic_array_view{{long(std::ranges::size(rg))}, std::to_address(std::begin(rg))} {}

    basic_array_view &operator=(basic_array_view const &rhs) noexcept {
      assign_from_ndarray(rhs);
      return *this;
    }

    template <ArrayOfRank<Rank> RHS>
    basic_array_view &operator=(RHS const &rhs) noexcept {
      static_assert(!is_const, "Cannot assign to an nda::basic_array_view with const value_type");
      assign_from_ndarray(rhs);
      return *this;
    }

    template <typename RHS>
    basic_array_view &operator=(RHS const &rhs) noexcept
      requires(is_scalar_for_v<RHS, basic_array_view>)
    {
      static_assert(!is_const, "Cannot assign to an nda::basic_array_view with const value_type");
      assign_from_scalar(rhs);
      return *this;
    }

    template <ArrayInitializer<basic_array_view> Initializer>
    basic_array_view &operator=(Initializer const &initializer) noexcept {
      EXPECTS(shape() == initializer.shape());
      initializer.invoke(*this);
      return *this;
    }

    template <typename T, int R, typename LP, char A, typename AP, typename OP>
    void rebind(basic_array_view<T, R, LP, A, AP, OP> v) noexcept {
      static_assert(R == Rank, "Cannot rebind a view to another view of a different rank");
      static_assert(std::is_same_v<std::remove_const_t<T>, std::remove_const_t<ValueType>>,
                    "Cannot rebind a view to another view with a different value_type (except for const)");
      static constexpr bool same_type = std::is_same_v<T, ValueType>;
      static_assert(same_type or is_const, "Cannot rebind a view with non-const value_type to another view with const value_type");
      if constexpr (same_type) {
        // FIXME Error message in layout error !
        lay = v.lay;
        sto = v.sto;
      } else if constexpr (is_const) {
        // the last if is always trivially true but in case of an error in the static_assert above,
        // it improves the error message by not compiling the = afterwards
        lay = layout_t{v.indexmap()};
        sto = storage_t{v.storage()};
      }
    }

    friend void swap(basic_array_view &a, basic_array_view &b) noexcept {
      std::swap(a.lay, b.lay);
      std::swap(a.sto, b.sto);
    }

    friend void deep_swap(basic_array_view a, basic_array_view b) noexcept {
      auto tmp = make_regular(a);
      a        = b;
      b        = tmp;
    }
 
// include common functionality of arrays and views
// Copyright (c) 2019--present, The Simons Foundation
// This file is part of TRIQS/nda and is licensed under the Apache License, Version 2.0.
// SPDX-License-Identifier: Apache-2.0
// See LICENSE in the root of this distribution for details.

[[nodiscard]] constexpr auto const &indexmap() const noexcept { return lay; }

[[nodiscard]] storage_t const &storage() const & noexcept { return sto; }

[[nodiscard]] storage_t &storage() & noexcept { return sto; }

[[nodiscard]] storage_t storage() && noexcept { return std::move(sto); }

[[nodiscard]] constexpr auto stride_order() const noexcept { return lay.stride_order; }

[[nodiscard]] ValueType const *data() const noexcept { return sto.data(); }

[[nodiscard]] ValueType *data() noexcept { return sto.data(); }

[[nodiscard]] auto const &shape() const noexcept { return lay.lengths(); }

[[nodiscard]] auto const &strides() const noexcept { return lay.strides(); }

[[nodiscard]] long size() const noexcept { return lay.size(); }

[[nodiscard]] long is_contiguous() const noexcept { return lay.is_contiguous(); }

[[nodiscard]] long has_positive_strides() const noexcept { return lay.has_positive_strides(); }

[[nodiscard]] bool empty() const { return sto.is_null(); }

[[nodiscard]] bool is_empty() const noexcept { return sto.is_null(); }

[[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];
}

[[nodiscard]] long shape(int i) const noexcept { return extent(i); }

[[nodiscard]] auto indices() const noexcept { return itertools::product_range(shape()); }

static constexpr bool is_stride_order_C() noexcept { return layout_t::is_stride_order_C(); }

static constexpr bool is_stride_order_Fortran() noexcept { return layout_t::is_stride_order_Fortran(); }

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");
}

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:
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<r_v_t>::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<r_v_t, res_rank, r_layout_p, newAlgebra, AccessorPolicy, OwningPolicy>{std::move(idxm), {self.sto, offset}};
    }
  }
}
 
public:
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...);
}

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...);
}

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...);
}

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);
}

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);
}

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);
}

static constexpr int iterator_rank = (has_strided_1d(layout_t::layout_prop) ? 1 : Rank);

using const_iterator = array_iterator<iterator_rank, ValueType const, typename AccessorPolicy::template accessor<ValueType>::pointer>;

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:
[[nodiscard]] const_iterator begin() const noexcept { return make_iterator<const_iterator>(false); }

[[nodiscard]] const_iterator cbegin() const noexcept { return make_iterator<const_iterator>(false); }

iterator begin() noexcept { return make_iterator<iterator>(false); }

[[nodiscard]] const_iterator end() const noexcept { return make_iterator<const_iterator>(true); }

[[nodiscard]] const_iterator cend() const noexcept { return make_iterator<const_iterator>(true); }

iterator end() noexcept { return make_iterator<iterator>(true); }

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);
}

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);
}

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);
}

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);
}

template <std::ranges::contiguous_range R>
auto &operator=(R const &rhs) noexcept
  requires(Rank == 1 and not MemoryArray<R> and not is_scalar_for_v<R, self_t>)
{
  *this = array_const_view<std::ranges::range_value_t<R>, 1>{rhs};
  return *this;
}
 
private:
// Implementation of the assignment from an n-dimensional array type.
template <typename RHS>
void assign_from_ndarray(RHS const &rhs) {
#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;
 
  // compile-time check for device arrays to avoid runtime errors
  static_assert(!(mem::on_device<self_t> or mem::on_device<RHS>) or (both_in_memory and same_stride_order and have_same_value_type_v<self_t, RHS>),
                "Error in assign_from_ndarray: Assignment to/from device arrays is not supported for the given types.");
 
  // 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) {
  if constexpr (mem::on_host<self_t>) {
    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 {
      for (auto &x : *this) x = scalar;
    }
  } else if constexpr (mem::on_device<self_t> or mem::on_unified<self_t>) {
    if constexpr (has_layout_strided_1d<self_t>) {
      if constexpr (has_contiguous_layout<self_t>) {
        mem::fill_n<mem::get_addr_space<self_t>>(data(), size(), value_type(scalar));
      } else {
        const long stri = indexmap().min_stride();
        mem::fill2D_n<mem::get_addr_space<self_t>>(data(), stri, 1, size(), value_type(scalar));
      }
    } else {
      auto bl_layout = get_block_layout(*this);
      if (bl_layout) {
        auto [n_bl, bl_size, bl_str] = *bl_layout;
        mem::fill2D_n<mem::get_addr_space<self_t>>(data(), bl_str, bl_size, n_bl, value_type(scalar));
      } else {
        // MAM: implement recursive call to fill_with_scalar on (i,nda::ellipsis{})
        NDA_RUNTIME_ERROR << "fill_with_scalar: Not implemented on device for generic (non-blocked) layout.";
      }
    }
  }
}
 
// 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
    diagonal(*this).fill_with_scalar(scalar);
  }
}
  };
 
  // Class template argument deduction guides.
  template <MemoryArray A>
  basic_array_view(A &&a)
     -> basic_array_view<std::conditional_t<std::is_const_v<std::remove_reference_t<A>>, const typename std::decay_t<A>::value_type,
                                            typename std::decay_t<A>::value_type>,
                         get_rank<A>, typename std::decay_t<A>::layout_policy_t, get_algebra<A>, default_accessor, borrowed<mem::get_addr_space<A>>>;
 
  template <typename V, size_t R>
  basic_array_view(std::array<V, R> &a)
     -> basic_array_view<V, 1, nda::basic_layout<nda::static_extents(R), nda::C_stride_order<1>, nda::layout_prop_e::contiguous>, 'V',
                         default_accessor, borrowed<>>;
 
  template <typename V, size_t R>
  basic_array_view(std::array<V, R> const &a)
     -> basic_array_view<const V, 1, nda::basic_layout<nda::static_extents(R), nda::C_stride_order<1>, nda::layout_prop_e::contiguous>, 'V',
                         default_accessor, borrowed<>>;
 
  template <std::ranges::contiguous_range R>
  basic_array_view(R &r) -> basic_array_view<std::conditional_t<std::is_const_v<R>, const typename R::value_type, typename R::value_type>, 1,
                                             C_layout, 'V', default_accessor, borrowed<>>;
 
} // namespace nda