Add floating support to BLAS/LAPACK backends

benchmarks/small_inv.cpp

@@ -12,7 +12,7 @@ static void inv(benchmark::State &state) {
   for (int i = 0; i < N; ++i)
     for (int j = 0; j < N; ++j) W(i, j) = (i > j ? 0.5 + i + 2.5 * j : i * 0.8 - j - 0.5);
 
-  while (state.KeepRunning()) { benchmark::DoNotOptimize(Wi = inverse(W)); }
+  while (state.KeepRunning()) { benchmark::DoNotOptimize(Wi = nda::linalg::inv(W)); }
 }
 
 BENCHMARK_TEMPLATE(inv, 1);

c++/nda/_impl_basic_array_view_common.hpp

@@ -438,11 +438,13 @@ auto &operator=(R const &rhs) noexcept
 private:
 // Implementation of the assignment from an n-dimensional array type.
 template <typename RHS>
-void assign_from_ndarray(RHS const &rhs) { // FIXME noexcept {
+void assign_from_ndarray(RHS const &rhs) 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();
+  if (this->shape() != rhs.shape()) {
+    std::cerr << "Error in assign_from_ndarray: Size mismatch:"
+              << "\n LHS.shape() = " << this->shape() << "\n RHS.shape() = " << rhs.shape() << std::endl;
+    std::terminate();
+  }
 #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");
@@ -453,6 +455,10 @@ void assign_from_ndarray(RHS const &rhs) { // FIXME noexcept {
   // 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;
@@ -470,7 +476,10 @@ void assign_from_ndarray(RHS const &rhs) { // FIXME noexcept {
         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 (n_bl_dst * bl_size_dst != n_bl_src * bl_size_src) {
+          std::cerr << "Error in assign_from_ndarray: Incompatible block sizes" << std::endl;
+          std::terminate();
+        }
         // 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;
@@ -494,15 +503,15 @@ void assign_from_ndarray(RHS const &rhs) { // FIXME noexcept {
   }
   // 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";
+    std::cerr << "Error in assign_from_ndarray: Elementwise assignment not implemented for arrays/views on the GPU" << std::endl;
+    std::terminate();
   }
   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 (mem::on_host<self_t>) {
     if constexpr (has_layout_strided_1d<self_t>) {
       const long L             = size();
@@ -517,23 +526,23 @@ void fill_with_scalar(Scalar const &scalar) noexcept {
     } else {
       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
+  } 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 {
-      // 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. ";
+        std::cerr << "Error in fill_with_scalar: Only block strided arrays/views are supported on the GPU";
+        std::terminate();
       }
     }
   }

c++/nda/arithmetic.hpp

@@ -12,8 +12,9 @@
 
 #include "./concepts.hpp"
 #include "./declarations.hpp"
+#include "./linalg/inv.hpp"
 #include "./linalg/matmul.hpp"
-#include "./linalg/det_and_inverse.hpp"
+#include "./linalg/matvecmul.hpp"
 #include "./macros.hpp"
 #include "./stdutil/complex.hpp"
 #include "./traits.hpp"
@@ -420,10 +421,10 @@ namespace nda {
       static_assert(r_algebra != 'A', "Error in nda::operator*: Can not multiply a matrix by an array");
       if constexpr (r_algebra == 'M')
         // matrix * matrix
-        return matmul(std::forward<L>(l), std::forward<R>(r));
+        return linalg::matmul(std::forward<L>(l), std::forward<R>(r));
       else
         // matrix * vector
-        return matvecmul(std::forward<L>(l), std::forward<R>(r));
+        return linalg::matvecmul(std::forward<L>(l), std::forward<R>(r));
     }
   }
 
@@ -495,7 +496,7 @@ namespace nda {
     // two matrices: M / M
     if constexpr (l_algebra == 'M') {
       static_assert(r_algebra == 'M', "Error in nda::operator*: Can not divide a matrix by an array/vector");
-      return std::forward<L>(l) * inverse(matrix<get_value_t<R>>{std::forward<R>(r)});
+      return std::forward<L>(l) * linalg::inv(matrix<get_value_t<R>>{std::forward<R>(r)});
     }
   }
 
@@ -532,7 +533,7 @@ namespace nda {
   Array auto operator/(S &&s, A &&a) { // NOLINT (S&& is mandatory for proper concept Array <: typename to work)
     static constexpr char algebra = get_algebra<A>;
     if constexpr (algebra == 'M')
-      return s * inverse(matrix<get_value_t<A>>{std::forward<A>(a)});
+      return s * linalg::inv(matrix<get_value_t<A>>{std::forward<A>(a)});
     else
       return expr<'/', std::decay_t<S>, A>{s, std::forward<A>(a)};
   }

c++/nda/basic_array.hpp

@@ -14,14 +14,16 @@
 #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/slice_static.hpp"
 #include "./layout_transforms.hpp"
 #include "./macros.hpp"
+#include "./matrix_functions.hpp"
 #include "./mem/address_space.hpp"
+#include "./mem/fill.hpp"
 #include "./mem/memcpy.hpp"
 #include "./mem/policies.hpp"
 #include "./stdutil/array.hpp"
@@ -31,17 +33,14 @@
 #include <array>
 #include <complex>
 #include <concepts>
+#include <exception>
 #include <initializer_list>
+#include <iostream>
 #include <random>
 #include <ranges>
 #include <type_traits>
 #include <utility>
 
-#ifdef NDA_ENFORCE_BOUNDCHECK
-#include <exception>
-#include <iostream>
-#endif // NDA_ENFORCE_BOUNDCHECK
-
 namespace nda {
 
   /**

c++/nda/basic_array_view.hpp

@@ -14,18 +14,18 @@
 #include "./clef.hpp"
 #include "./concepts.hpp"
 #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 "./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/fill.hpp"
 #include "./mem/policies.hpp"
 #include "./traits.hpp"
 
@@ -34,16 +34,13 @@
 #include <algorithm>
 #include <array>
 #include <cstring>
+#include <exception>
+#include <iostream>
 #include <memory>
 #include <ranges>
 #include <type_traits>
 #include <utility>
 
-#ifdef NDA_ENFORCE_BOUNDCHECK
-#include <exception>
-#include <iostream>
-#endif // NDA_ENFORCE_BOUNDCHECK
-
 namespace std {
 
   /**