18
19
20
21[[nodiscard]]
constexpr auto const &indexmap()
const noexcept {
return lay; }
24
25
26
27[[nodiscard]] storage_t
const &storage()
const &
noexcept {
return sto; }
30
31
32
33[[nodiscard]] storage_t &storage() &
noexcept {
return sto; }
36
37
38
39[[nodiscard]] storage_t storage() &&
noexcept {
return std::move(sto); }
42
43
44
45
46
47[[nodiscard]]
constexpr auto stride_order()
const noexcept {
return lay.stride_order; }
50
51
52
53[[nodiscard]] ValueType
const *data()
const noexcept {
return sto.data(); }
56
57
58
59[[nodiscard]] ValueType *data()
noexcept {
return sto.data(); }
62
63
64
65[[nodiscard]]
auto const &shape()
const noexcept {
return lay.lengths(); }
68
69
70
71[[nodiscard]]
auto const &strides()
const noexcept {
return lay.strides(); }
74
75
76
77[[nodiscard]]
long size()
const noexcept {
return lay.size(); }
80
81
82
83[[nodiscard]]
long is_contiguous()
const noexcept {
return lay.is_contiguous(); }
86
87
88
89[[nodiscard]]
bool empty()
const {
return sto.is_null(); }
92[[nodiscard]]
bool is_empty()
const noexcept {
return sto.is_null(); }
95
96
97
98[[nodiscard]]
long extent(
int i)
const noexcept {
99#ifdef NDA_ENFORCE_BOUNDCHECK
100 if (i < 0 || i >= rank) {
101 std::cerr <<
"Error in extent: Dimension " << i <<
" is incompatible with array of rank " << rank << std::endl;
105 return lay.lengths()[i];
109[[nodiscard]]
long shape(
int i)
const noexcept {
return extent(i); }
112
113
114
115[[nodiscard]]
auto indices()
const noexcept {
return itertools::product_range(shape()); }
118
119
120
121static constexpr bool is_stride_order_C()
noexcept {
return layout_t::is_stride_order_C(); }
124
125
126
127static constexpr bool is_stride_order_Fortran()
noexcept {
return layout_t::is_stride_order_Fortran(); }
130
131
132
133
134
135
136
137
138decltype(
auto) operator()(_linear_index_t idx)
const noexcept {
139 if constexpr (layout_t::layout_prop == layout_prop_e::strided_1d)
140 return sto[idx.value * lay.min_stride()];
141 else if constexpr (layout_t::layout_prop == layout_prop_e::contiguous)
142 return sto[idx.value];
144 static_assert(always_false<layout_t>,
"Internal error in array/view: Calling this type with a _linear_index_t is not allowed");
148decltype(
auto) operator()(_linear_index_t idx)
noexcept {
149 if constexpr (layout_t::layout_prop == layout_prop_e::strided_1d)
150 return sto[idx.value * lay.min_stride()];
151 else if constexpr (layout_t::layout_prop == layout_prop_e::contiguous)
152 return sto[idx.value];
154 static_assert(always_false<layout_t>,
"Internal error in array/view: Calling this type with a _linear_index_t is not allowed");
159#ifdef NDA_ENFORCE_BOUNDCHECK
160static constexpr bool has_no_boundcheck =
false;
162static constexpr bool has_no_boundcheck =
true;
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182template <
char ResultAlgebra,
bool SelfIsRvalue,
typename Self,
typename... Ts>
183FORCEINLINE
static decltype(
auto) call(Self &&self, Ts
const &...idxs)
noexcept(has_no_boundcheck) {
185 using r_v_t = std::conditional_t<std::is_const_v<std::remove_reference_t<Self>>, ValueType
const, ValueType>;
188 if constexpr (clef::is_any_lazy<Ts...>) {
190 return clef::make_expr_call(std::forward<Self>(self), idxs...);
191 }
else if constexpr (
sizeof...(Ts) == 0) {
193 return basic_array_view<r_v_t, Rank, LayoutPolicy, Algebra, AccessorPolicy, OwningPolicy>{self.lay, self.sto};
196 static_assert(((layout_t::
template argument_is_allowed_for_call_or_slice<Ts> + ...) > 0),
197 "Error in array/view: Slice arguments must be convertible to range, ellipsis, or long (or string if the layout permits it)");
200 static constexpr int n_args_long = (layout_t::
template argument_is_allowed_for_call<Ts> + ...);
202 if constexpr (n_args_long == rank) {
204 long offset = self.lay(idxs...);
205 if constexpr (is_view
or not SelfIsRvalue) {
207 return AccessorPolicy::
template accessor<ValueType>::access(self.sto.data(), offset);
210 return ValueType{self.sto[offset]};
214 auto const [offset, idxm] = self.lay.slice(idxs...);
215 static constexpr auto res_rank =
decltype(idxm)::rank();
217 static constexpr char newAlgebra = (ResultAlgebra ==
'M' and (res_rank == 1) ?
'V' : ResultAlgebra);
219 using r_layout_p =
typename detail::layout_to_policy<std::decay_t<
decltype(idxm)>>::type;
220 return basic_array_view<ValueType, res_rank, r_layout_p, newAlgebra, AccessorPolicy, OwningPolicy>{std::move(idxm), {self.sto, offset}};
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247template <
typename... Ts>
248FORCEINLINE
decltype(
auto) operator()(Ts
const &...idxs)
const &
noexcept(has_no_boundcheck) {
249 static_assert((rank == -1)
or (
sizeof...(Ts) == rank)
or (
sizeof...(Ts) == 0)
or (ellipsis_is_present<Ts...>
and (
sizeof...(Ts) <= rank + 1)),
250 "Error in array/view: Incorrect number of parameters in call operator");
251 return call<Algebra,
false>(*
this, idxs...);
255template <
typename... Ts>
256FORCEINLINE
decltype(
auto) operator()(Ts
const &...idxs) &
noexcept(has_no_boundcheck) {
257 static_assert((rank == -1)
or (
sizeof...(Ts) == rank)
or (
sizeof...(Ts) == 0)
or (ellipsis_is_present<Ts...>
and (
sizeof...(Ts) <= rank + 1)),
258 "Error in array/view: Incorrect number of parameters in call operator");
259 return call<Algebra,
false>(*
this, idxs...);
263template <
typename... Ts>
264FORCEINLINE
decltype(
auto) operator()(Ts
const &...idxs) &&
noexcept(has_no_boundcheck) {
265 static_assert((rank == -1)
or (
sizeof...(Ts) == rank)
or (
sizeof...(Ts) == 0)
or (ellipsis_is_present<Ts...>
and (
sizeof...(Ts) <= rank + 1)),
266 "Error in array/view: Incorrect number of parameters in call operator");
267 return call<Algebra,
true>(*
this, idxs...);
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
288decltype(
auto) operator[](T
const &idx)
const &
noexcept(has_no_boundcheck) {
289 static_assert((rank == 1),
"Error in array/view: Subscript operator is only available for rank 1 views/arrays in C++17/20");
290 return call<Algebra,
false>(*
this, idx);
295decltype(
auto) operator[](T
const &x) &
noexcept(has_no_boundcheck) {
296 static_assert((rank == 1),
"Error in array/view: Subscript operator is only available for rank 1 views/arrays in C++17/20");
297 return call<Algebra,
false>(*
this, x);
302decltype(
auto) operator[](T
const &x) &&
noexcept(has_no_boundcheck) {
303 static_assert((rank == 1),
"Error in array/view: Subscript operator is only available for rank 1 views/arrays in C++17/20");
304 return call<Algebra,
true>(*
this, x);
308static constexpr int iterator_rank = (has_strided_1d(layout_t::layout_prop) ? 1 : Rank);
311using const_iterator = array_iterator<iterator_rank, ValueType
const,
typename AccessorPolicy::
template accessor<ValueType>::pointer>;
314using iterator = array_iterator<iterator_rank, ValueType,
typename AccessorPolicy::
template accessor<ValueType>::pointer>;
318template <
typename Iterator>
319[[nodiscard]]
auto make_iterator(
bool at_end)
const noexcept {
320 if constexpr (iterator_rank == Rank) {
322 if constexpr (layout_t::is_stride_order_C()) {
324 return Iterator{indexmap().lengths(), indexmap().strides(), sto.data(), at_end};
327 return Iterator{nda::permutations::apply(layout_t::stride_order, indexmap().lengths()),
328 nda::permutations::apply(layout_t::stride_order, indexmap().strides()), sto.data(), at_end};
332 return Iterator{std::array<
long, 1>{size()}, std::array<
long, 1>{indexmap().min_stride()}, sto.data(), at_end};
338[[nodiscard]] const_iterator begin()
const noexcept {
return make_iterator<const_iterator>(
false); }
341[[nodiscard]] const_iterator cbegin()
const noexcept {
return make_iterator<const_iterator>(
false); }
344iterator begin()
noexcept {
return make_iterator<iterator>(
false); }
347[[nodiscard]] const_iterator end()
const noexcept {
return make_iterator<const_iterator>(
true); }
350[[nodiscard]] const_iterator cend()
const noexcept {
return make_iterator<const_iterator>(
true); }
353iterator end()
noexcept {
return make_iterator<iterator>(
true); }
356
357
358
359
360
361
362
363
364
365
366
367template <
typename RHS>
368auto &operator+=(RHS
const &rhs)
noexcept {
369 static_assert(
not is_const,
"Error in array/view: Can not assign to a const view");
370 return operator=(*
this + rhs);
374
375
376
377
378
379
380
381
382
383
384
385template <
typename RHS>
386auto &operator-=(RHS
const &rhs)
noexcept {
387 static_assert(
not is_const,
"Error in array/view: Can not assign to a const view");
388 return operator=(*
this - rhs);
392
393
394
395
396
397
398
399
400
401
402
403template <
typename RHS>
404auto &operator*=(RHS
const &rhs)
noexcept {
405 static_assert(
not is_const,
"Error in array/view: Can not assign to a const view");
406 return operator=((*
this) * rhs);
410
411
412
413
414
415
416
417
418
419
420
421template <
typename RHS>
422auto &operator/=(RHS
const &rhs)
noexcept {
423 static_assert(
not is_const,
"Error in array/view: Can not assign to a const view");
424 return operator=(*
this / rhs);
428
429
430
431
432
433
434
435template <std::ranges::contiguous_range R>
436auto &operator=(R
const &rhs)
noexcept
437 requires(Rank == 1
and not MemoryArray<R>)
439 *
this = basic_array_view{rhs};
445template <
typename RHS>
446void assign_from_ndarray(RHS
const &rhs) {
447#ifdef NDA_ENFORCE_BOUNDCHECK
448 if (
this->shape() != rhs.shape())
449 NDA_RUNTIME_ERROR <<
"Error in assign_from_ndarray: Size mismatch:"
450 <<
"\n LHS.shape() = " <<
this->shape() <<
"\n RHS.shape() = " << rhs.shape();
453 static_assert(std::is_assignable_v<value_type &, get_value_t<RHS>>,
"Error in assign_from_ndarray: Incompatible value types");
456 static constexpr bool both_in_memory = MemoryArray<self_t>
and MemoryArray<RHS>;
459 static constexpr bool same_stride_order = get_layout_info<self_t>.stride_order == get_layout_info<RHS>.stride_order;
462 if constexpr (both_in_memory
and same_stride_order) {
463 if (rhs.empty())
return;
465 static constexpr bool both_1d_strided = has_layout_strided_1d<self_t>
and has_layout_strided_1d<RHS>;
466 if constexpr (mem::on_host<self_t, RHS>
and both_1d_strided) {
468 for (
long i = 0; i < size(); ++i) (*
this)(_linear_index_t{i}) = rhs(_linear_index_t{i});
470 }
else if constexpr (!mem::on_host<self_t, RHS>
and have_same_value_type_v<self_t, RHS>) {
472 auto bl_layout_dst = get_block_layout(*
this);
473 auto bl_layout_src = get_block_layout(rhs);
474 if (bl_layout_dst && bl_layout_src) {
475 auto [n_bl_dst, bl_size_dst, bl_str_dst] = *bl_layout_dst;
476 auto [n_bl_src, bl_size_src, bl_str_src] = *bl_layout_src;
478 if (n_bl_dst * bl_size_dst != n_bl_src * bl_size_src) NDA_RUNTIME_ERROR <<
"Error in assign_from_ndarray: Incompatible block sizes";
480 if (n_bl_dst == 1 && n_bl_src > 1) {
482 bl_size_dst /= n_bl_src;
483 bl_str_dst = bl_size_dst;
485 if (n_bl_src == 1 && n_bl_dst > 1) {
487 bl_size_src /= n_bl_dst;
488 bl_str_src = bl_size_src;
491 if (n_bl_dst == n_bl_src && bl_size_dst == bl_size_src) {
492 mem::memcpy2D<mem::get_addr_space<self_t>, mem::get_addr_space<RHS>>((
void *)data(), bl_str_dst *
sizeof(value_type), (
void *)rhs.data(),
493 bl_str_src *
sizeof(value_type), bl_size_src *
sizeof(value_type),
501 if constexpr (mem::on_device<self_t> || mem::on_device<RHS>) {
502 NDA_RUNTIME_ERROR <<
"Error in assign_from_ndarray: Fallback to elementwise assignment not implemented for arrays/views on the GPU";
504 nda::for_each(shape(), [
this, &rhs](
auto const &...args) { (*
this)(args...) = rhs(args...); });
508template <
typename Scalar>
509void fill_with_scalar(Scalar
const &scalar)
noexcept {
511 if constexpr (has_layout_strided_1d<self_t>) {
512 const long L = size();
513 auto *
__restrict const p = data();
514 if constexpr (has_contiguous_layout<self_t>) {
515 for (
long i = 0; i < L; ++i) p[i] = scalar;
517 const long stri = indexmap().min_stride();
518 const long Lstri = L * stri;
519 for (
long i = 0; i < Lstri; i += stri) p[i] = scalar;
523 for (
auto &x : *
this) x = scalar;
528template <
typename Scalar>
529void assign_from_scalar(Scalar
const &scalar)
noexcept {
530 static_assert(!is_const,
"Error in assign_from_ndarray: Cannot assign to a const view");
531 if constexpr (Algebra !=
'M') {
533 fill_with_scalar(scalar);
538 if constexpr (is_scalar_or_convertible_v<Scalar>)
541 fill_with_scalar(Scalar{0 * scalar});
542 const long imax = std::min(extent(0), extent(1));
543 for (
long i = 0; i < imax; ++i) operator()(i, i) = scalar;