Merge branch 'develop' into fix/from-jax-array-issues

CHANGELOG

@@ -1,3 +1,28 @@
+Version 0.0.7:
+  - Overhaul the (internal) ScalarType API:
+    . the overloads of convert_copy have been removed in favour of the variant that takes two ScalarPointers;
+    . the single versions of add(_inplace) and friends have been replaced with more flexible add_into;
+      batch compute methods and friends;
+    . replaced single value uminus with uminus into with similar signature to to add_into and friends;
+    . removed single value copy method;
+  - Added constructor for ScalarPointer from type_id and pointer.
+  - Implementations of ScalarType methods that are essentially the same for all types are implemented
+    in a common implementation layer.
+  - Added threading support in platform
+  - add_into and friends have threading support if available and enabled.
+  - Added default implementation of type_id_of so that non-specialized types
+  look for a ScalarType object.
+  - Greatly simplified the design of ScalarMatrix - it now only supports full,
+  dense matrices.
+  - Redesigned the interface between the Scalar linear algebra and
+  MKL/BLAS+LAPACK.
+  - Added function to query ring characteristics of a ScalarType - currently
+  unused.
+
+
+
+
+
 Version 0.0.6:
   - Externally sourced streams (sound-file streams) now support setting channel
   types/schema in factory function.

algebra/src/libalgebra_lite_internal/sparse_mutable_ref_scalar_trait.h

@@ -99,12 +99,15 @@ class SparseMutableRefScalarImpl : public ScalarInterface
     void assign(const void* data, const string& type_id) override
     {
         value_type tmp = static_cast<const value_type&>(m_data);
-        type()->convert_copy({type(), &tmp}, data, 1, type_id);
+        type()->convert_copy(
+                {type(), &tmp}, {scalars::get_type(type_id), data}, 1
+        );
         m_data = tmp;
     }
     ScalarPointer to_pointer() override
     {
-        RPY_THROW(std::runtime_error,
+        RPY_THROW(
+                std::runtime_error,
                 "cannot get non-const pointer to proxy reference type"
         );
     }

algebra/src/lite_context.h

@@ -431,7 +431,7 @@ OutType LiteContext<Coefficients>::construct_impl(
     std::vector<scalar_type> tmp;
     if (data.data.type() != ctype()) {
         tmp.resize(data.data.size());
-        ctype()->convert_copy(tmp.data(), data.data, size);
+        ctype()->convert_copy({ctype(), tmp.data()}, data.data, size);
         data_ptr = tmp.data();
     } else {
         data_ptr = data.data.raw_cast<const scalar_type>();

cmake/roughpy_helpers.cmake

@@ -199,7 +199,9 @@ function(add_roughpy_component _name)
     _get_component_name(_component_name ${_component})
 
     if (NOT _lib_type STREQUAL INTERFACE)
-        set(_private_include_dirs "${CMAKE_CURRENT_LIST_DIR}/include/roughpy/${_component}/")
+        set(_private_include_dirs
+                "${CMAKE_CURRENT_LIST_DIR}/include/roughpy/${_component}/"
+                "${CMAKE_CURRENT_LIST_DIR}/src")
         if (ARG_PVT_INCLUDE_DIRS)
             foreach (_pth IN LISTS ARG_PVT_INCLUDE_DIRS)
                 list(APPEND _private_include_dirs ${CMAKE_CURRENT_LIST_DIR}/${_pth})

core/include/roughpy/core/helpers.h

@@ -75,6 +75,11 @@ count_bits(T val) noexcept
     return std::bitset<CHAR_BIT * sizeof(T)>(val).count();
 }
 
+RPY_NO_DISCARD constexpr size_t static_log2p1(size_t value) noexcept
+{
+    return (value == 0) ? 0 : 1 + static_log2p1(value >> 1);
+}
+
 /**
  * @brief
  * @tparam T
@@ -107,11 +112,9 @@ class MaybeOwned
         return *this;
     }
 
-    RPY_NO_DISCARD
-    operator T*() const noexcept { return p_data; }
+    RPY_NO_DISCARD operator T*() const noexcept { return p_data; }
 
-    RPY_NO_DISCARD
-    operator bool() const noexcept { return p_data != nullptr; }
+    RPY_NO_DISCARD operator bool() const noexcept { return p_data != nullptr; }
 };
 
 }// namespace rpy

platform/CMakeLists.txt

@@ -6,17 +6,19 @@ project(Roughpy_Platform VERSION 0.0.1)
 add_roughpy_component(Platform
         SOURCES
         src/configuration.cpp
+        src/threading/openmp_threading.cpp
         PUBLIC_HEADERS
         include/roughpy/platform.h
         include/roughpy/platform/filesystem.h
         include/roughpy/platform/configuration.h
         include/roughpy/platform/serialization.h
-
+        include/roughpy/platform/threads.h
         PUBLIC_DEPS
         Boost::boost
         Boost::url
         Boost::disable_autolinking
         cereal::cereal
+        OpenMP::OpenMP_CXX
         NEEDS
         RoughPy::Core
         )

platform/include/roughpy/platform/threads.h

@@ -0,0 +1,59 @@
+//
+// Created by sam on 28/07/23.
+//
+
+#ifndef ROUGHPY_PLATFORM_THREADS_H
+#define ROUGHPY_PLATFORM_THREADS_H
+
+#include <roughpy/core/alloc.h>
+#include <roughpy/core/macros.h>
+#include <roughpy/core/types.h>
+
+#define RPY_ALLOW_THREADING
+#if defined(_OPENMP) && defined(RPY_ALLOW_THREADING)
+#  define RPY_THREADING_OPENMP 1
+#endif
+
+namespace rpy {
+namespace platform {
+
+enum struct ThreadBackend
+{
+    Disabled = 0,
+    OpenMP = 1
+};
+
+struct ThreadState {
+    /// The threading library responsible for managing multi-threaded
+    /// computation.
+    ThreadBackend backend;
+
+    /// The maximum number of threads that can be spawned by a single operation.
+    int max_threads;
+
+    /// The maximum possible number of threads that can be spawned
+    int max_available_threads;
+
+    /// Is multithreading enabled
+    bool is_enabled;
+};
+
+RPY_NO_DISCARD RPY_EXPORT ThreadState get_thread_state();
+
+RPY_NO_DISCARD RPY_EXPORT constexpr bool threading_available()
+{
+#ifndef RPY_ALLOW_THREADING
+    return false;
+#else
+    return true;
+#endif
+}
+
+RPY_EXPORT void set_num_threads(int num_threads);
+
+RPY_EXPORT void set_threading_enabled(bool state);
+
+}// namespace platform
+}// namespace rpy
+
+#endif// ROUGHPY_PLATFORMS_THREADS_H

platform/src/threading/openmp_threading.cpp

@@ -0,0 +1,37 @@
+//
+// Created by sam on 28/07/23.
+//
+
+#include <roughpy/platform/threads.h>
+
+#include <atomic>
+
+#include <omp.h>
+
+#if !defined(RPY_THREADING_OPENMP) || !RPY_THREADING_OPENMP
+#  error "OpenMP backend cannot be used unless OpenMP is avaiable."
+#endif
+
+using namespace rpy;
+using namespace rpy::platform;
+
+static const string s_omp_backend_name = "OpenMP";
+
+static std::atomic_bool s_enable_omp_threading = true;
+
+ThreadState rpy::platform::get_thread_state()
+{
+    return {ThreadBackend::OpenMP, omp_get_num_threads(), omp_get_max_threads(),
+            s_enable_omp_threading.load(std::memory_order_relaxed)};
+}
+
+
+void rpy::platform::set_num_threads(int num_threads)
+{
+    omp_set_num_threads(num_threads);
+}
+
+void rpy::platform::set_threading_enabled(bool state)
+{
+    s_enable_omp_threading.store(state, std::memory_order_relaxed);
+}

roughpy/src/scalars/scalars.cpp

@@ -55,7 +55,7 @@ void python::assign_py_object_to_scalar(
     if (py::isinstance<py::float_>(object)) {
         *ptr = object.cast<double>();
     } else if (py::isinstance<py::int_>(object)) {
-        *ptr = object.cast<long long>();
+        *ptr = object.cast<int64_t>();
     } else if (RPyPolynomial_Check(object.ptr())) {
         *ptr = RPyPolynomial_cast(object.ptr());
     } else {
@@ -168,7 +168,7 @@ static inline void dl_copy_strided(
 {
     if (ndim == 1) {
         if (strides[0] == 1) {
-            dst.type()->convert_copy(dst.ptr(), src, shape[0]);
+            dst.type()->convert_copy(dst, src, shape[0]);
         } else {
             for (std::int64_t i = 0; i < shape[0]; ++i) {
                 dst[i] = src[i * strides[0]];
@@ -574,7 +574,8 @@ scalars::KeyScalarArray python::py_to_buffer(
         // The only way type can still be null is if there are no elements.
         if (options.type != nullptr) {
             options.type->convert_copy(
-                    result, info.ptr, static_cast<dimn_t>(info.size), type_id
+                    result, {scalars::get_type(type_id), info.ptr},
+                    static_cast<dimn_t>(info.size)
             );
             options.shape.assign(info.shape.begin(), info.shape.end());
         }

roughpy/src/streams/brownian_stream.cpp

@@ -81,7 +81,7 @@ Brownian_from_generator(const py::args& args, const py::kwargs& kwargs)
 
     // TODO: Finish this
 
-    BrownianStream inner(pmd.scalar_type->get_rng(), std::move(md));
+    BrownianStream inner(pmd.scalar_type->get_rng("", {}), std::move(md));
     Stream stream(std::move(inner));
 
     return py::reinterpret_steal<py::object>(

scalars/CMakeLists.txt

@@ -16,6 +16,16 @@ set(RPY_USE_MKL ${MKL_FOUND})
 configure_file(scalar_blas_defs.h.in scalar_blas_defs.h @ONLY)
 add_roughpy_component(Scalars
         SOURCES
+        src/linear_algebra/blas.h
+        src/linear_algebra/blas_complex_double.cpp
+        src/linear_algebra/blas_complex_float.cpp
+        src/linear_algebra/blas_double.cpp
+        src/linear_algebra/blas_float.cpp
+        src/linear_algebra/lapack.h
+        src/linear_algebra/lapack_complex_double.cpp
+        src/linear_algebra/lapack_complex_float.cpp
+        src/linear_algebra/lapack_double.cpp
+        src/linear_algebra/lapack_float.cpp
         src/scalar_type.cpp
         src/scalar_pointer.cpp
         src/scalar_interface.cpp
@@ -28,31 +38,38 @@ add_roughpy_component(Scalars
         src/key_scalar_array.cpp
         src/scalar_stream.cpp
         src/standard_scalar_type.h
-        src/half_type.h
-        src/half_type.cpp
-        src/float_type.cpp
-        src/float_type.h
-        src/double_type.cpp
-        src/double_type.h
-        src/RationalType.cpp
-        src/RationalType.h
+        src/standard_linalg.h
         src/standard_random_generator.cpp
         src/standard_random_generator.h
         src/scalar_matrix.cpp
         src/scalar_blas.cpp
-        src/float_blas.cpp
-        src/float_blas.h
-        src/scalar_blas_impl.h
-        src/half_random_generator.cpp
-        src/half_random_generator.h
-        src/b_float_16_type.cpp
-        src/b_float_16_type.h
-        src/bfloat16_random_generator.cpp
-        src/bfloat16_random_generator.h
         src/random_impl.cpp
         src/random_impl.h
-        src/rational_poly_scalar_type.cpp
-        src/rational_poly_scalar_type.h
+        src/scalar_blas_impl.h
+        src/scalar_implementations/half/half_type.h
+        src/scalar_implementations/half/half_type.cpp
+        src/scalar_implementations/float/float_type.cpp
+        src/scalar_implementations/float/float_type.h
+        src/scalar_implementations/double/double_type.cpp
+        src/scalar_implementations/double/double_type.h
+        src/scalar_implementations/rational/RationalType.cpp
+        src/scalar_implementations/rational/RationalType.h
+        src/scalar_implementations/float/float_blas.cpp
+        src/scalar_implementations/float/float_blas.h
+        src/scalar_implementations/half/half_random_generator.cpp
+        src/scalar_implementations/half/half_random_generator.h
+        src/scalar_implementations/bfloat16/b_float_16_type.cpp
+        src/scalar_implementations/bfloat16/b_float_16_type.h
+        src/scalar_implementations/bfloat16/bfloat16_random_generator.cpp
+        src/scalar_implementations/bfloat16/bfloat16_random_generator.h
+        src/scalar_implementations/rational_poly/rational_poly_scalar_type.cpp
+        src/scalar_implementations/rational_poly/rational_poly_scalar_type.h
+        src/scalar_implementations/double/double_lin_alg.cpp
+        src/scalar_implementations/double/double_lin_alg.h
+        #src/scalar_implementations/complex_double/complex_double_type.cpp
+        #src/scalar_implementations/complex_double/complex_double_type.h
+        #src/scalar_implementations/complex_float/complex_float_type.cpp
+        #src/scalar_implementations/complex_float/complex_float_type.h
         ${CMAKE_CURRENT_BINARY_DIR}/scalar_blas_defs.h
         PUBLIC_HEADERS
         include/roughpy/scalars/scalars_fwd.h

scalars/include/roughpy/scalars/scalar.h

@@ -87,7 +87,7 @@ class Scalar : private ScalarPointer
             const auto& id = type_id_of<ScalarArg>();
             if (type == nullptr) { type = ScalarType::for_id(id); }
             ScalarPointer::operator=(type->allocate(1));
-            type->convert_copy(to_mut_pointer(), &arg, 1, id);
+            type->convert_copy(to_mut_pointer(), {id, &arg}, 1);
         }
     }
 
@@ -115,13 +115,14 @@ class Scalar : private ScalarPointer
             ScalarPointer::operator=(p_type->allocate(1));
         }
 
-        const auto& type_id = type_id_of<T>();
         if ((m_flags & interface_flag) != 0) {
             static_cast<ScalarInterface*>(const_cast<void*>(p_data))
-                    ->assign(std::addressof(arg), type_id);
+                    ->assign(std::addressof(arg), type_id_of<T>());
         } else {
-            p_type->convert_copy(static_cast<const ScalarPointer&>(*this),
-                                 std::addressof(arg), 1, type_id);
+            const auto& type_id = type_id_of<T>();
+            p_type->convert_copy(static_cast<ScalarPointer&>(*this),
+                                 {type_id, std::addressof(arg)}, 1);
+
         }
 
         return *this;
@@ -240,7 +241,7 @@ struct type_of_T_not_defined {
     static T cast(ScalarPointer scalar)
     {
         T result;
-        scalar.type()->convert_copy({nullptr, &result}, scalar.ptr(), 1,
+        scalar.type()->convert_copy({nullptr, &result}, scalar, 1,
                                     type_id_of<T>());
         return result;
     }