Merge pull request #232 from LibRapid/test

Updates to benchmarks and some small performance improvements

Author: Pencilcaseman

CMakeLists.txt

@@ -37,6 +37,8 @@ option(LIBRAPID_CODE_COV "Compile LibRapid C++ with Coverage" OFF)
 option(LIBRAPID_STRICT "Force all warnings into errors (use with caution)" OFF)
 option(LIBRAPID_QUIET "Hide warnings generated WITHIN LibRapid's source" OFF)
 
+option(LIBRAPID_USE_PRECOMPILED_HEADER "Use precompiled headers to speed up compilation" OFF)
+
 option(LIBRAPID_GET_FFTW "Clone and use FFTW -- WARNING: SEE DOCUMENTATION" OFF)
 option(LIBRAPID_GET_BLAS "Download pre-built OpenBLAS binaries and use them" OFF)
 option(LIBRAPID_GET_MULTIPREC "Download generic multiprecision libraries, as opposed to trying to find one on the system" OFF)
@@ -106,7 +108,9 @@ target_compile_definitions(${module_name} PUBLIC LIBRAPID_SOURCE="${LIBRAPID_SOU
 # clang-format on
 
 # Precompiled Headers
-target_precompile_headers(${module_name} PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/librapid/include/librapid/core/librapidPch.hpp")
+if (LIBRAPID_USE_PRECOMPILED_HEADER)
+    target_precompile_headers(${module_name} PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/librapid/include/librapid/core/librapidPch.hpp")
+endif ()
 
 # Extract system information
 set(IS_LINUX OFF)

docs/source/cmakeIntegration.md

@@ -61,6 +61,25 @@ DEFAULT: OFF
 Disable all warnings from LibRapid. This is useful if you are using LibRapid as a dependency and want a cleaner
 compilation output. Warnings should be minimal in the first place, but this option is provided just in case.
 
+### ``LIBRAPID_USE_PRECOMPILED_HEADER``
+
+```
+DEFAULT: OFF
+```
+
+Enable the use of precompiled headers within LibRapid's compilation. This can be useful to accelerate compilation, but
+can lead to some strange build errors, which is why it is disabled by default.
+
+:::{warning}
+One such build error occurs on some macOS systems with GCC. The resulting error is something along the lines of:
+
+```
+Unknown flag -Xarch_amd64
+```
+
+If you encounter this error, try disabling ``LIBRAPID_USE_PRECOMPILED_HEADER``.
+:::
+
 ### ``LIBRAPID_GET_FFTW``
 
 ```

docs/source/index.md

@@ -68,7 +68,8 @@ Getting Started <gettingStarted>
 CMake Integration <cmakeIntegration>
 API Reference <apiReference>
 Tutorials <tutorials>
-Performance and Benchmarks <performance/performance>
+Performance <performance/performance>
+Benchmark Results <BenchmarkResults/BenchmarkResults>
 Caution <caution>
 ```
 

examples/templateCMakeLists.txt

@@ -30,6 +30,10 @@ set(LIBRAPID_STRICT OFF) # Enable all warnings and treat them as errors
 
 set(LIBRAPID_QUIET OFF) # Silence all warnings
 
+set(LIBRAPID_USE_PRECOMPILED_HEADER OFF) # Use precompiled headers to accelerate
+                                         # compilation at the cost of stability
+                                         # on some systems.
+
 set(LIBRAPID_USE_BLAS ON) # Attempt to use a BLAS library -- if not found,
                           # LibRapid falls back to less optimized routines
 

librapid/include/librapid/array/arrayContainer.hpp

@@ -73,7 +73,8 @@ namespace librapid {
 		struct IsArrayContainer : std::false_type {};
 
 		template<typename ShapeType, typename StorageScalar>
-		struct IsArrayContainer<array::ArrayContainer<ShapeType, StorageScalar>> : std::true_type {};
+		struct IsArrayContainer<array::ArrayContainer<ShapeType, StorageScalar>> : std::true_type {
+		};
 
 		LIBRAPID_DEFINE_AS_TYPE(typename StorageScalar,
 								array::ArrayContainer<Shape COMMA StorageScalar>);
@@ -93,7 +94,6 @@ namespace librapid {
 			using Scalar	  = typename StorageType::Scalar;
 			using Packet	  = typename typetraits::TypeInfo<Scalar>::Packet;
 			using Backend	  = typename typetraits::TypeInfo<ArrayContainer>::Backend;
-			using Iterator	  = detail::ArrayIterator<GeneralArrayView<ArrayContainer, ShapeType>>;
 
 			using DirectSubscriptType	 = typename detail::SubscriptType<StorageType>::Direct;
 			using DirectRefSubscriptType = typename detail::SubscriptType<StorageType>::Ref;
@@ -330,19 +330,19 @@ namespace librapid {
 
 			/// \brief Return an iterator to the beginning of the array container
 			/// \return Iterator
-			LIBRAPID_ALWAYS_INLINE Iterator begin() const noexcept;
+			LIBRAPID_ALWAYS_INLINE auto begin() const noexcept;
 
 			/// \brief Return an iterator to the end of the array container
 			/// \return Iterator
-			LIBRAPID_ALWAYS_INLINE Iterator end() const noexcept;
+			LIBRAPID_ALWAYS_INLINE auto end() const noexcept;
 
 			/// \brief Return an iterator to the beginning of the array container
 			/// \return Iterator
-			LIBRAPID_ALWAYS_INLINE Iterator begin();
+			LIBRAPID_ALWAYS_INLINE auto begin();
 
 			/// \brief Return an iterator to the end of the array container
 			/// \return Iterator
-			LIBRAPID_ALWAYS_INLINE Iterator end();
+			LIBRAPID_ALWAYS_INLINE auto end();
 
 			template<typename T, typename Char, typename Ctx>
 			void str(const fmt::formatter<T, Char> &format, char bracket, char separator,
@@ -765,12 +765,11 @@ namespace librapid {
 		ArrayContainer<ShapeType_, StorageType_>::packet(size_t index) const -> Packet {
 			auto ptr = LIBRAPID_ASSUME_ALIGNED(m_storage.begin());
 
-#if defined(LIBRAPID_NATIVE_ARCH) && !defined(LIBRAPID_OSX)
-			// On MacOS (and other platforms??) we cannot use aligned loads in arrays due to one
-			// annoying edge case. Normally, all SIMD loads will be aligned to a 64-byte boundary.
-			// Say, however, this array is a sub-array of a larger array. If the outer dimension
-			// of the larger array does not result in a 64-byte alignment, the data of *this* array
-			// will not be correctly aligned, hence causing a segfault.
+#if defined(LIBRAPID_NATIVE_ARCH)
+			LIBRAPID_ASSERT(
+			  reinterpret_cast<uintptr_t>(ptr) % typetraits::TypeInfo<Scalar>::packetWidth == 0,
+			  "ArrayContainer::packet called on unaligned storage");
+
 			return xsimd::load_aligned(ptr + index);
 #else
 			return xsimd::load_unaligned(ptr + index);
@@ -786,10 +785,15 @@ namespace librapid {
 		template<typename ShapeType_, typename StorageType_>
 		LIBRAPID_ALWAYS_INLINE void
 		ArrayContainer<ShapeType_, StorageType_>::writePacket(size_t index, const Packet &value) {
+			auto ptr = LIBRAPID_ASSUME_ALIGNED(m_storage.begin());
+
 #if defined(LIBRAPID_NATIVE_ARCH)
-			value.store_aligned(m_storage.begin() + index);
+			LIBRAPID_ASSERT(
+			  reinterpret_cast<uintptr_t>(ptr) % typetraits::TypeInfo<Scalar>::packetWidth == 0,
+			  "ArrayContainer::packet called on unaligned storage");
+			value.store_aligned(ptr + index);
 #else
-			value.store_unaligned(m_storage.begin() + index);
+			value.store_unaligned(ptr + index);
 #endif
 		}
 
@@ -881,24 +885,24 @@ namespace librapid {
 
 		template<typename ShapeType_, typename StorageType_>
 		LIBRAPID_ALWAYS_INLINE auto ArrayContainer<ShapeType_, StorageType_>::begin() const noexcept
-		  -> Iterator {
-			return Iterator(GeneralArrayView(*this), 0);
+		  -> auto {
+			return detail::ArrayIterator(createGeneralArrayView(*this), 0);
 		}
 
 		template<typename ShapeType_, typename StorageType_>
 		LIBRAPID_ALWAYS_INLINE auto ArrayContainer<ShapeType_, StorageType_>::end() const noexcept
-		  -> Iterator {
-			return Iterator(GeneralArrayView(*this), m_shape[0]);
+		  -> auto {
+			return detail::ArrayIterator(createGeneralArrayView(*this), m_shape[0]);
 		}
 
 		template<typename ShapeType_, typename StorageType_>
-		LIBRAPID_ALWAYS_INLINE auto ArrayContainer<ShapeType_, StorageType_>::begin() -> Iterator {
-			return Iterator(GeneralArrayView(*this), 0);
+		LIBRAPID_ALWAYS_INLINE auto ArrayContainer<ShapeType_, StorageType_>::begin() -> auto {
+			return detail::ArrayIterator(createGeneralArrayView(*this), 0);
 		}
 
 		template<typename ShapeType_, typename StorageType_>
-		LIBRAPID_ALWAYS_INLINE auto ArrayContainer<ShapeType_, StorageType_>::end() -> Iterator {
-			return Iterator(GeneralArrayView(*this), m_shape[0]);
+		LIBRAPID_ALWAYS_INLINE auto ArrayContainer<ShapeType_, StorageType_>::end() -> auto {
+			return detail::ArrayIterator(createGeneralArrayView(*this), m_shape[0]);
 		}
 
 		template<typename ShapeType_, typename StorageType_>

librapid/include/librapid/array/arrayIterator.hpp

@@ -5,7 +5,7 @@ namespace librapid::detail {
 	template<typename T>
 	class ArrayIterator {
 	public:
-		using IndexType = int64_t;
+		using IndexType = size_t;
 
 		/// Default constructor should never be used
 		ArrayIterator() = delete;
@@ -14,6 +14,8 @@ namespace librapid::detail {
 
 		explicit LIBRAPID_ALWAYS_INLINE ArrayIterator(const T &array, IndexType index);
 
+		explicit LIBRAPID_ALWAYS_INLINE ArrayIterator(T &&array, IndexType index);
+
 		/// Copy an ArrayIterator object (const)
 		/// \param other The array to copy
 		LIBRAPID_ALWAYS_INLINE ArrayIterator(const ArrayIterator &other) = default;
@@ -53,6 +55,10 @@ namespace librapid::detail {
 	LIBRAPID_ALWAYS_INLINE ArrayIterator<T>::ArrayIterator(const T &array, IndexType index) :
 			m_array(array), m_index(index) {}
 
+	template<typename T>
+	LIBRAPID_ALWAYS_INLINE ArrayIterator<T>::ArrayIterator(T &&array, IndexType index) :
+			m_array(std::move(array)), m_index(index) {}
+
 	template<typename T>
 	LIBRAPID_ALWAYS_INLINE ArrayIterator<T> &ArrayIterator<T>::operator++() {
 		++m_index;

librapid/include/librapid/array/assignOps.hpp

@@ -328,7 +328,7 @@ namespace librapid {
                                                          ::librapid::detail::LibRapidType::Scalar,
                                                        int> = 0>
         LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto dataSourceExtractor(const T &obj) {
-            return obj.storage().begin().get();
+            return obj.storage().begin();
         }
 
         template<typename T, typename std::enable_if_t<typetraits::TypeInfo<T>::type ==
@@ -518,7 +518,7 @@ namespace librapid {
               std::make_index_sequence<argSize>(),
               filename,
               kernelName,
-              reinterpret_cast<Scalar *>(lhs.storage().begin().get()),
+              reinterpret_cast<Scalar *>(lhs.storage().begin()),
               function);
         }
     }  // namespace detail

librapid/include/librapid/array/linalg/arrayMultiply.hpp

@@ -106,10 +106,14 @@ namespace librapid {
 			/// \return Class of the array multiplication
 			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE MatmulClass matmulClass() const;
 
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE ShapeType calculateShape() const;
+
 			/// \brief Determine the shape of the result
 			/// \return Shape of the result
 			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE ShapeType shape() const;
 
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE size_t size() const;
+
 			/// \brief Determine the number of dimensions of the result
 			/// \return Number of dimensions of the result
 			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE int64_t ndim() const;
@@ -173,6 +177,9 @@ namespace librapid {
 			ScalarA m_alpha; // Scaling factor for A
 			TypeB m_b;		 // Second array
 			ScalarB m_beta;	 // Scaling factor for B
+
+			ShapeType m_shape;
+			size_t m_size;
 		};
 
 		template<typename ShapeTypeA, typename StorageTypeA, typename ShapeTypeB,
@@ -182,23 +189,26 @@ namespace librapid {
 										   TypeB &&b, Beta beta) :
 				m_transA(transA),
 				m_transB(transB), m_a(std::forward<TypeA>(a)), m_alpha(static_cast<ScalarA>(alpha)),
-				m_b(std::forward<TypeB>(b)), m_beta(static_cast<ScalarB>(beta)) {}
+				m_b(std::forward<TypeB>(b)), m_beta(static_cast<ScalarB>(beta)),
+				m_shape(calculateShape()), m_size(m_shape.size()) {}
 
 		template<typename ShapeTypeA, typename StorageTypeA, typename ShapeTypeB,
 				 typename StorageTypeB, typename Alpha, typename Beta>
 		ArrayMultiply<ShapeTypeA, StorageTypeA, ShapeTypeB, StorageTypeB, Alpha,
 					  Beta>::ArrayMultiply(TypeA &&a, TypeB &&b) :
 				m_transA(false),
 				m_transB(false), m_a(std::forward<TypeA>(a)), m_alpha(1),
-				m_b(std::forward<TypeB>(b)), m_beta(0) {}
+				m_b(std::forward<TypeB>(b)), m_beta(0), m_shape(calculateShape()),
+				m_size(m_shape.size()) {}
 
 		template<typename ShapeTypeA, typename StorageTypeA, typename ShapeTypeB,
 				 typename StorageTypeB, typename Alpha, typename Beta>
 		ArrayMultiply<ShapeTypeA, StorageTypeA, ShapeTypeB, StorageTypeB, Alpha,
 					  Beta>::ArrayMultiply(bool transA, bool transB, TypeA &&a, TypeB &&b) :
 				m_transA(transA),
 				m_transB(transB), m_a(std::forward<TypeA>(a)), m_alpha(1),
-				m_b(std::forward<TypeB>(b)), m_beta(0) {}
+				m_b(std::forward<TypeB>(b)), m_beta(0), m_shape(calculateShape()),
+				m_size(m_shape.size()) {}
 
 		template<typename ShapeTypeA, typename StorageTypeA, typename ShapeTypeB,
 				 typename StorageTypeB, typename Alpha, typename Beta>
@@ -268,8 +278,8 @@ namespace librapid {
 
 		template<typename ShapeTypeA, typename StorageTypeA, typename ShapeTypeB,
 				 typename StorageTypeB, typename Alpha, typename Beta>
-		auto ArrayMultiply<ShapeTypeA, StorageTypeA, ShapeTypeB, StorageTypeB, Alpha, Beta>::shape()
-		  const -> ShapeType {
+		auto ArrayMultiply<ShapeTypeA, StorageTypeA, ShapeTypeB, StorageTypeB, Alpha,
+						   Beta>::calculateShape() const -> ShapeType {
 			const auto &shapeA		= m_a.shape();
 			const auto &shapeB		= m_b.shape();
 			MatmulClass matmulClass = this->matmulClass();
@@ -294,6 +304,21 @@ namespace librapid {
 			return {1};
 		}
 
+		template<typename ShapeTypeA, typename StorageTypeA, typename ShapeTypeB,
+				 typename StorageTypeB, typename Alpha, typename Beta>
+		auto ArrayMultiply<ShapeTypeA, StorageTypeA, ShapeTypeB, StorageTypeB, Alpha, Beta>::shape()
+		  const -> ShapeType {
+			return m_shape;
+		}
+
+		template<typename ShapeTypeA, typename StorageTypeA, typename ShapeTypeB,
+				 typename StorageTypeB, typename Alpha, typename Beta>
+		auto
+		ArrayMultiply<ShapeTypeA, StorageTypeA, ShapeTypeB, StorageTypeB, Alpha, Beta>::size() const
+		  -> size_t {
+			return m_size;
+		}
+
 		template<typename ShapeTypeA, typename StorageTypeA, typename ShapeTypeB,
 				 typename StorageTypeB, typename Alpha, typename Beta>
 		auto

librapid/include/librapid/array/linalg/transpose.hpp

@@ -441,8 +441,6 @@ namespace librapid {
 			using ArrayType		 = TransposeType;
 			using BaseType		 = typename std::decay_t<TransposeType>;
 			using Scalar		 = typename typetraits::TypeInfo<BaseType>::Scalar;
-			using Reference		 = BaseType &;
-			using ConstReference = const BaseType &;
 			using ShapeType		 = typename BaseType::ShapeType;
 			using Backend		 = typename typetraits::TypeInfo<BaseType>::Backend;
 
@@ -624,8 +622,8 @@ namespace librapid {
 				else {
 					if (m_inputShape.ndim() == 2) {
 						int64_t blockSize		= global::cacheLineSize / sizeof(Scalar);
-						auto *__restrict outPtr = out.storage().begin().get();
-						auto *__restrict inPtr	= m_array.storage().begin().get();
+						auto *__restrict outPtr = out.storage().begin();
+						auto *__restrict inPtr	= m_array.storage().begin();
 						detail::cuda::transposeImpl(
 						  outPtr, inPtr, m_inputShape[0], m_inputShape[1], m_alpha, blockSize);
 					} else {