merge with unstable branch

Fikret Ardal · Fikret Ardal · commit 7747a8c4ee3a · 2025-06-30T00:56:53.000+03:00
diff --git a/benchmarks/constr_heap_vs_sso.cpp b/benchmarks/constr_heap_vs_sso.cpp
@@ -42,7 +42,7 @@ struct toy_mat1_t {
 
 struct toy_mat2_t {
   sso<1000>::handle<double> storage;
-  toy_mat2_t(){}; // Custom constructor -> Not an aggregate type
+  toy_mat2_t() {}; // Custom constructor -> Not an aggregate type
 };
 
 BENCH_EXPR(constr_toy_mat1, toy_mat1_t{});
diff --git a/c++/nda/_impl_basic_array_view_common.hpp b/c++/nda/_impl_basic_array_view_common.hpp
@@ -525,19 +525,39 @@ void assign_from_ndarray(RHS const &rhs) { // FIXME noexcept {
 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             = indexmap().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;
+  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 {
-      const long stri  = indexmap().min_stride();
-      const long Lstri = L * stri;
-      for (long i = 0; i != Lstri; i += stri) p[i] = scalar;
+      for (auto &x : *this) x = scalar;
+    }
+  } else if constexpr (mem::on_device<self_t> or mem::on_unified<self_t>) { // on device
+    if constexpr (has_layout_strided_1d<self_t>) {                          // possibly contiguous
+      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 {
+      // check for 2D layout
+      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 yet for general layout. ";
+      }
     }
-  } else {
-    // no compile-time memory layout guarantees
-    for (auto &x : *this) x = scalar;
   }
 }
 
@@ -556,7 +576,6 @@ void assign_from_scalar(Scalar const &scalar) noexcept {
       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;
+    diagonal(*this).fill_with_scalar(scalar);
   }
 }
diff --git a/c++/nda/basic_array_view.hpp b/c++/nda/basic_array_view.hpp
@@ -16,13 +16,16 @@
 #include "./declarations.hpp"
 #include "./exceptions.hpp"
 #include "./iterators.hpp"
+#include "layout/slice_static.hpp"
 #include "./layout/for_each.hpp"
 #include "./layout/idx_map.hpp"
 #include "./layout/permutation.hpp"
 #include "./layout/range.hpp"
 #include "./macros.hpp"
 #include "./mem/address_space.hpp"
 #include "./mem/memcpy.hpp"
+#include "./mem/memset.hpp"
+#include "./mem/fill.hpp"
 #include "./mem/policies.hpp"
 #include "./traits.hpp"
 #include "./config.hpp"
diff --git a/c++/nda/basic_functions.hpp b/c++/nda/basic_functions.hpp
@@ -605,7 +605,7 @@ namespace nda {
     // slicing helper function
     auto slice_Axis = [](Array auto &a, range r) {
       auto all_or_range = std::make_tuple(range::all, r);
-      return [&]<auto... Is>(std::index_sequence<Is...>) { return a(std::get<Is == Axis>(all_or_range)...); }(std::make_index_sequence<rank>{});
+      return [&]<auto... Is>(std::index_sequence<Is...>) { return a(std::get < Is == Axis > (all_or_range)...); }(std::make_index_sequence<rank>{});
     };
 
     // initialize concatenated array
diff --git a/c++/nda/blas/gemm.hpp b/c++/nda/blas/gemm.hpp
@@ -50,9 +50,15 @@ namespace nda::blas {
     EXPECTS(a.extent(1) == b.extent(0));
     EXPECTS(a.extent(0) == c.extent(0));
     EXPECTS(b.extent(1) == c.extent(1));
+
+    if (beta == 0.0) {
+      c = 0 * alpha;
+    } else {
+      c *= beta;
+    }
+
     for (int i = 0; i < a.extent(0); ++i) {
       for (int j = 0; j < b.extent(1); ++j) {
-        c(i, j) = beta * c(i, j);
         for (int k = 0; k < a.extent(1); ++k) c(i, j) += alpha * a(i, k) * b(k, j);
       }
     }
diff --git a/c++/nda/blas/gemv.hpp b/c++/nda/blas/gemv.hpp
@@ -46,8 +46,14 @@ namespace nda::blas {
   void gemv_generic(get_value_t<A> alpha, A const &a, X const &x, get_value_t<A> beta, Y &&y) { // NOLINT (temporary views are allowed here)
     EXPECTS(a.extent(1) == x.extent(0));
     EXPECTS(a.extent(0) == y.extent(0));
+
+    if (beta == 0.0) {
+      y = 0 * alpha;
+    } else {
+      y *= beta;
+    }
+
     for (int i = 0; i < a.extent(0); ++i) {
-      y(i) = beta * y(i);
       for (int k = 0; k < a.extent(1); ++k) y(i) += alpha * a(i, k) * x(k);
     }
   }
diff --git a/c++/nda/clef/literals.hpp b/c++/nda/clef/literals.hpp
@@ -20,8 +20,7 @@ namespace nda::clef::literals {
    */
 
   // Define literal placeholders starting from the end of the allowed index spectrum.
-#define PH(I)                                                                                                                                        \
-  (placeholder<63 - (I)> {})
+#define PH(I) (placeholder<63 - (I)>{})
 
   /// Generic placeholder #1.
   constexpr auto i_ = PH(0);
@@ -58,6 +57,6 @@ namespace nda::clef::literals {
 
 #undef PH
 
-/** @} */
+  /** @} */
 
 } // namespace nda::clef::literals
diff --git a/c++/nda/clef/make_lazy.hpp b/c++/nda/clef/make_lazy.hpp
@@ -122,11 +122,8 @@ namespace nda::clef {
   }                                                                                                                                                  \
                                                                                                                                                      \
   template <typename... Args>                                                                                                                        \
-  auto operator()(Args &&...args) &&                                                                                                                 \
-    requires(nda::clef::is_any_lazy<Args...>)                                                                                                        \
-  {                                                                                                                                                  \
-    return make_expr_call(std::move(*this), std::forward<Args>(args)...);                                                                            \
-  }
+     auto operator()(Args &&...args)                                                                                                                 \
+     && requires(nda::clef::is_any_lazy<Args...>) { return make_expr_call(std::move(*this), std::forward<Args>(args)...); }
 
   /** @} */
 
diff --git a/c++/nda/gtest_tools.hpp b/c++/nda/gtest_tools.hpp
@@ -126,7 +126,7 @@ ::testing::AssertionResult array_almost_zero(X const &x) {
   nda::array<nda::get_value_t<X>, nda::get_rank<X>> x_reg = x;
 
   constexpr double eps = 1.e-10;
-  const auto max = max_element(abs(x_reg));
+  const auto max       = max_element(abs(x_reg));
   if (x_reg.size() == 0 || max < eps)
     return ::testing::AssertionSuccess();
   else
diff --git a/c++/nda/h5.hpp b/c++/nda/h5.hpp
@@ -51,7 +51,8 @@ namespace nda {
     // Given an array/view, prepare and return the corresponding h5::array_view to be written/read into.
     template <MemoryArray A>
     auto prepare_h5_array_view(const A &a) {
-      auto [parent_shape, h5_strides] = h5::array_interface::get_parent_shape_and_h5_strides(a.indexmap().strides().data(), A::rank, a.shape().data());
+      auto [parent_shape, h5_strides] =
+         h5::array_interface::get_parent_shape_and_h5_strides(a.indexmap().strides().data(), A::rank, a.shape().data());
       auto v = h5::array_interface::array_view{h5::hdf5_type<get_value_t<A>>(), (void *)a.data(), A::rank, is_complex_v<typename A::value_type>};
       for (int u = 0; u < A::rank; ++u) {
         v.slab.count[u]   = a.shape()[u];
diff --git a/c++/nda/layout/slice_static.hpp b/c++/nda/layout/slice_static.hpp
@@ -274,7 +274,7 @@ namespace nda::slice_static {
       static constexpr std::array<int, P> new_stride_order = slice_stride_order(IdxMap::stride_order, n_of_p);
 
       // compile-time layout properties of the resulting index map
-      static constexpr bool has_only_rangeall_and_long = ((std::is_constructible_v<long, Args> or std::is_base_of_v<range::all_t, Args>)and...);
+      static constexpr bool has_only_rangeall_and_long = ((std::is_constructible_v<long, Args> or std::is_base_of_v<range::all_t, Args>) and ...);
       static constexpr layout_prop_e li =
          slice_layout_prop(P, has_only_rangeall_and_long, args_is_rangeall, IdxMap::stride_order, IdxMap::layout_prop, e_pos, e_len);
 
diff --git a/c++/nda/linalg/matmul.hpp b/c++/nda/linalg/matmul.hpp
@@ -46,9 +46,9 @@ namespace nda {
     static constexpr bool is_valid_gemm_triple = []() {
       using blas::has_F_layout;
       if constexpr (has_F_layout<C>) {
-        return !(conj_A and has_F_layout<A>)and!(conj_B and has_F_layout<B>);
+        return !(conj_A and has_F_layout<A>) and !(conj_B and has_F_layout<B>);
       } else {
-        return !(conj_B and !has_F_layout<B>)and!(conj_A and !has_F_layout<A>);
+        return !(conj_B and !has_F_layout<B>) and !(conj_A and !has_F_layout<A>);
       }
     }();
 
diff --git a/c++/nda/mapped_functions.hpp b/c++/nda/mapped_functions.hpp
@@ -40,8 +40,7 @@ namespace nda {
 
     // Get the complex conjugate of a scalar.
     template <nda::Scalar S>
-    auto conj(S x)
-    {
+    auto conj(S x) {
       if constexpr (is_complex_v<S>) {
         return std::conj(x);
       } else {
diff --git a/c++/nda/mem/address_space.hpp b/c++/nda/mem/address_space.hpp
@@ -160,15 +160,15 @@ namespace nda::mem {
 
   /// Constexpr variable that is true if all given types have the same address space.
   template <typename A0, typename... A>
-  static constexpr bool have_same_addr_space = ((get_addr_space<A0> == get_addr_space<A>)and... and true);
+  static constexpr bool have_same_addr_space = ((get_addr_space<A0> == get_addr_space<A>) and ... and true);
 
   /// Constexpr variable that is true if all given types have an address space compatible with `Host`.
   template <typename... Ts>
-  static constexpr bool have_host_compatible_addr_space = ((on_host<Ts> or on_unified<Ts>)and...);
+  static constexpr bool have_host_compatible_addr_space = ((on_host<Ts> or on_unified<Ts>) and ...);
 
   /// Constexpr variable that is true if all given types have an address space compatible with `Device`.
   template <typename... Ts>
-  static constexpr bool have_device_compatible_addr_space = ((on_device<Ts> or on_unified<Ts>)and...);
+  static constexpr bool have_device_compatible_addr_space = ((on_device<Ts> or on_unified<Ts>) and ...);
 
   /// Constexpr variable that is true if all given types have compatible address spaces.
   template <typename... Ts>
diff --git a/c++/nda/mem/allocators.hpp b/c++/nda/mem/allocators.hpp
@@ -14,6 +14,8 @@
 #include "./aligned_alloc.hpp"
 #include "./malloc.hpp"
 #include "./memset.hpp"
+#include "./memcpy.hpp"
+#include "./fill.hpp"
 #include "../macros.hpp"
 
 #include <algorithm>
diff --git a/c++/nda/mem/fill.hpp b/c++/nda/mem/fill.hpp
@@ -0,0 +1,114 @@
+// Copyright (c) 2022-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: Miguel Morales, Nils Wentzell
+
+#pragma once
+
+#include <cstdlib>
+#include <algorithm>
+#include <vector>
+#include <span>
+#include <ranges>
+
+#include "address_space.hpp"
+#include "../traits.hpp"
+
+namespace nda::mem {
+
+  /**
+   * @brief Fills a range of memory with a specified value.
+   *
+   * The behavior depends on the AddressSpace (Host, Device, or Unified).
+   *
+   * @tparam AdrSp The address space (e.g., Host, Device, Unified).
+   * @tparam T The type of the elements to fill.
+   * @param first Pointer to the beginning of the range.
+   * @param count Number of elements to fill.
+   * @param value The value to fill the range with.
+   * @return Pointer to the end of the filled range.
+   */
+  template <AddressSpace AdrSp, typename T>
+    requires(nda::is_scalar_or_convertible_v<T>)
+  T *fill_n(T *first, size_t count, const T &value) {
+    check_adr_sp_valid<AdrSp>();
+    static_assert(nda::have_device == nda::have_cuda, "Adjust function for new device types");
+
+    if constexpr (AdrSp == Host) {
+      return std::fill_n(first, count, value);
+    } else { // Device or Unified
+      auto value_bytes = std::as_bytes(std::span(&value, 1));
+      bool is_zero     = std::ranges::all_of(value_bytes, [](auto b) { return b == std::byte{0}; });
+      if (is_zero) {
+        device_error_check(cudaMemset(first, 0, count * sizeof(T)), "cudaMemset");
+      } else {
+        for (int n = 0; n < sizeof(T); ++n) {
+          const int byte_value [[maybe_unused]] = static_cast<int>(value_bytes[n]);
+          device_error_check(cudaMemset2D((char *)(first) + n, sizeof(T), byte_value, 1, count), "cudaMemset2D");
+        }
+      }
+      return first + count;
+    }
+  }
+
+  /**
+   * @brief Fills a range of memory between two pointers with a specified value.
+   *
+   * Internally calls `fill_n`.
+   *
+   * @tparam AdrSp The address space (e.g., Host, Device, Unified).
+   * @tparam T The type of the elements to fill.
+   * @param first Pointer to the beginning of the range.
+   * @param end Pointer to the end of the range.
+   * @param value The value to fill the range with.
+   * @return Pointer to the end of the filled range.
+   */
+  template <AddressSpace AdrSp, typename T>
+    requires(nda::is_scalar_or_convertible_v<T>)
+  T *fill(T *first, T *end, const T &value) {
+    if (std::distance(first, end) > 0) return fill_n<AdrSp>(first, std::distance(first, end), value);
+    return first;
+  }
+
+  /**
+   * @brief Fills a 2D memory region with a specified value.
+   *
+   * The behavior depends on the AddressSpace (Host, Device, or Unified).
+   *
+   * @tparam AdrSp The address space (e.g., Host, Device, Unified).
+   * @tparam T The type of the elements to fill.
+   * @param first Pointer to the beginning of the 2D memory region.
+   * @param pitch The memory pitch between rows.
+   * @param width The number of elements to fill in each row.
+   * @param height The number of rows to fill.
+   * @param value The value to fill the 2D region with.
+   */
+  template <AddressSpace AdrSp, typename T>
+    requires(nda::is_scalar_or_convertible_v<T>)
+  void fill2D_n(T *first [[maybe_unused]], size_t pitch [[maybe_unused]], size_t width, size_t height, const T &value) {
+    check_adr_sp_valid<AdrSp>();
+    static_assert(nda::have_device == nda::have_cuda, "Adjust function for new device types");
+    static_assert(AdrSp == mem::Device or AdrSp == mem::Unified, "Not implemented for host memory");
+
+    bool is_zero = std::ranges::all_of(std::as_bytes(std::span(&value, 1)), [](auto b) { return b == std::byte{0}; });
+    if (is_zero) {
+      device_error_check(cudaMemset2D(first, pitch * sizeof(T), 0, width * sizeof(T), height), "cudaMemset2D");
+    } else {
+      std::vector<T> v(width * height, value);
+      device_error_check(cudaMemcpy2D(first, pitch * sizeof(T), v.data(), width * sizeof(T), width * sizeof(T), height, cudaMemcpyDefault),
+                         "cudaMemcpy2D");
+    }
+  }
+
+} // namespace nda::mem
diff --git a/c++/nda/traits.hpp b/c++/nda/traits.hpp
@@ -15,6 +15,7 @@
 #include <ranges>
 #include <type_traits>
 #include <utility>
+#include <tuple>
 
 namespace nda {
 
diff --git a/packaging/Dockerfile.msan b/packaging/Dockerfile.msan
diff --git a/test/c++/nda_basic_array_and_view.cpp b/test/c++/nda_basic_array_and_view.cpp
diff --git a/test/c++/nda_cuda.cpp b/test/c++/nda_cuda.cpp
diff --git a/test/c++/nda_math_functions.cpp b/test/c++/nda_math_functions.cpp
