diff --git a/CMakeLists.txt b/CMakeLists.txt index 2cabe7d6..fa53f39d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -84,7 +84,17 @@ cmaize_find_or_build_optional_dependency( CMAKE_ARGS BUILD_TESTING=OFF ENABLE_EIGEN_SUPPORT=ON ) -set(DEPENDENCIES utilities parallelzone Boost::boost eigen sigma) + +cmaize_find_or_build_dependency( + WeaklyTypedFloat + NAME WeaklyTypedFloat + URL https://www.github.com/NWChemEx/weaklytypedfloat + VERSION master + BUILD_TARGET wtf + FIND_TARGET nwx::wtf +) + +set(DEPENDENCIES utilities parallelzone Boost::boost eigen sigma wtf) if("${ENABLE_CUTENSOR}") include(cmake/FindcuTENSOR.cmake) diff --git a/include/tensorwrapper/buffer/mdbuffer.hpp b/include/tensorwrapper/buffer/mdbuffer.hpp new file mode 100644 index 00000000..72f5c765 --- /dev/null +++ b/include/tensorwrapper/buffer/mdbuffer.hpp @@ -0,0 +1,64 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include + +namespace tensorwrapper::buffer { + +/** @brief A multidimensional (MD) buffer. + * + * This class is a dense multidimensional buffer of floating-point values. + */ +class MDBuffer { +private: + using traits_type = types::ClassTraits; + +public: + /// Add types to public API + ///@{ + using buffer_type = typename traits_type::buffer_type; + using pimpl_type = typename traits_type::pimpl_type; + using pimpl_pointer = typename traits_type::pimpl_pointer; + using rank_type = typename traits_type::rank_type; + using shape_type = typename traits_type::shape_type; + ///@} + + MDBuffer() noexcept; + + template + MDBuffer(shape_type shape, std::vector elements) { + MDBuffer(std::move(shape), buffer_type(std::move(elements))); + } + + MDBuffer(shape_type shape, buffer_type buffer); + + rank_type rank() const; + +private: + explicit MDBuffer(pimpl_pointer pimpl) noexcept; + + bool has_pimpl_() const noexcept; + + void assert_pimpl_() const; + + pimpl_type& pimpl_(); + const pimpl_type& pimpl_() const; + + pimpl_pointer m_pimpl_; +}; + +} // namespace tensorwrapper::buffer diff --git a/include/tensorwrapper/forward_declarations.hpp b/include/tensorwrapper/forward_declarations.hpp new file mode 100644 index 00000000..16c51064 --- /dev/null +++ b/include/tensorwrapper/forward_declarations.hpp @@ -0,0 +1,36 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once + +namespace tensorwrapper { + +namespace buffer { +namespace detail_ { +class MDBufferPIMPL; +} + +class MDBuffer; + +} // namespace buffer + +namespace shape { + +class Smooth; + +} // namespace shape + +} // namespace tensorwrapper diff --git a/include/tensorwrapper/symmetry/permutation.hpp b/include/tensorwrapper/symmetry/permutation.hpp index e749402b..ca5b349f 100644 --- a/include/tensorwrapper/symmetry/permutation.hpp +++ b/include/tensorwrapper/symmetry/permutation.hpp @@ -161,6 +161,36 @@ class Permutation : public Operation { */ mode_index_type size() const noexcept { return m_cycles_.size(); } + /** @brief Permutes the objects in @p input according to *this. + * + * @tparam T The type of a container-like object. It must support size(), + * and operator[]. + * + * @param[in] input The object to permute. + * + * @return A copy of @p input with its elements permuted according to + * *this. + * + * @throw std::runtime_error if the size of @p input does not match the + * rank of *this. Strong throw guarantee. + */ + template + T apply(T input) const { + if(input.size() != m_rank_) + throw std::runtime_error( + "Input size does not match permutation rank"); + for(const auto& cycle : m_cycles_) { + if(cycle.size() < 2) continue; + T buffer = input; + for(std::size_t i = 0; i < cycle.size(); ++i) { + auto from = cycle[i]; + auto to = cycle[(i + 1) % cycle.size()]; + input[to] = buffer[from]; + } + } + return input; + } + // ------------------------------------------------------------------------- // -- Utility methods // ------------------------------------------------------------------------- diff --git a/include/tensorwrapper/types/floating_point.hpp b/include/tensorwrapper/types/floating_point.hpp index d34ba1f1..fb37346a 100644 --- a/include/tensorwrapper/types/floating_point.hpp +++ b/include/tensorwrapper/types/floating_point.hpp @@ -15,6 +15,7 @@ */ #pragma once +#include #include #ifdef ENABLE_SIGMA #include diff --git a/include/tensorwrapper/types/mdbuffer_traits.hpp b/include/tensorwrapper/types/mdbuffer_traits.hpp new file mode 100644 index 00000000..27c74421 --- /dev/null +++ b/include/tensorwrapper/types/mdbuffer_traits.hpp @@ -0,0 +1,52 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include +#include +#include +#include + +namespace tensorwrapper::types { + +struct MDBufferTraitsCommon { + using value_type = wtf::fp::Float; + using const_reference = wtf::fp::FloatView; + using buffer_type = wtf::buffer::FloatBuffer; + using const_buffer_view = wtf::buffer::BufferView; + using shape_type = shape::Smooth; + using rank_type = typename shape_type::rank_type; + using pimpl_type = tensorwrapper::buffer::detail_::MDBufferPIMPL; + using pimpl_pointer = std::unique_ptr; +}; + +template<> +struct ClassTraits + : public MDBufferTraitsCommon { + using reference = wtf::fp::FloatView; + + using buffer_view = wtf::buffer::BufferView; + using const_buffer_view = wtf::buffer::BufferView; +}; + +template<> +struct ClassTraits + : public MDBufferTraitsCommon { + using reference = wtf::fp::FloatView; + using buffer_view = wtf::buffer::BufferView; +}; + +} // namespace tensorwrapper::types diff --git a/include/tensorwrapper/types/shape_traits.hpp b/include/tensorwrapper/types/shape_traits.hpp new file mode 100644 index 00000000..00d7c6dd --- /dev/null +++ b/include/tensorwrapper/types/shape_traits.hpp @@ -0,0 +1,36 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include +#include +#include + +namespace tensorwrapper::types { + +struct ShapeTraitsCommon { + using size_type = std::size_t; + using rank_type = unsigned short; +}; + +template<> +struct ClassTraits : public ShapeTraitsCommon {}; + +template<> +struct ClassTraits + : public ShapeTraitsCommon {}; + +} // namespace tensorwrapper::types diff --git a/src/tensorwrapper/backends/cutensor/cuda_tensor.cpp b/src/tensorwrapper/backends/cutensor/cuda_tensor.cpp new file mode 100644 index 00000000..ee3dca0f --- /dev/null +++ b/src/tensorwrapper/backends/cutensor/cuda_tensor.cpp @@ -0,0 +1,49 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#include "cuda_tensor.hpp" + +#ifdef ENABLE_CUTENSOR +#include "eigen_tensor.cuh" +#endif + +namespace tensorwrapper::backends::cutensor { + +#define TPARAMS template +#define CUDA_TENSOR CUDATensor + +TPARAMS +void CUDA_TENSOR::contraction_assignment(label_type this_label, + label_type lhs_label, + label_type rhs_label, + const_my_reference lhs, + const_my_reference rhs) { +#ifdef ENABLE_CUTENSOR + cutensor_contraction(this_label, lhs_label, rhs_label, lhs, rhs, + *this); +#else + throw std::runtime_error( + "cuTENSOR backend not enabled. Recompile with -DENABLE_CUTENSOR."); +#endif +} + +#undef CUDA_TENSOR +#undef TPARAMS + +template class CUDATensor; +template class CUDATensor; + +} // namespace tensorwrapper::backends::cutensor diff --git a/src/tensorwrapper/buffer/detail_/eigen_tensor.cu b/src/tensorwrapper/backends/cutensor/cuda_tensor.cu similarity index 74% rename from src/tensorwrapper/buffer/detail_/eigen_tensor.cu rename to src/tensorwrapper/backends/cutensor/cuda_tensor.cu index dd6efc1a..b3e12329 100644 --- a/src/tensorwrapper/buffer/detail_/eigen_tensor.cu +++ b/src/tensorwrapper/backends/cutensor/cuda_tensor.cu @@ -14,12 +14,16 @@ * limitations under the License. */ #ifdef ENABLE_CUTENSOR +#include "cuda_tensor.cuh" #include "cutensor_traits.cuh" -#include "eigen_tensor.cuh" #include #include -namespace tensorwrapper::buffer::detail_ { +namespace tensorwrapper::backends::cutensor { + +// Some common typedefs +using mode_vector_t = std::vector; +using int64_vector_t = std::vector; // Handle cuda errors #define HANDLE_CUDA_ERROR(x) \ @@ -38,12 +42,7 @@ namespace tensorwrapper::buffer::detail_ { if(err != CUTENSOR_STATUS_SUCCESS) { \ printf("Error: %s\n", cutensorGetErrorString(err)); \ exit(-1); \ - } \ - }; - -// Some common typedefs -using mode_vector_t = std::vector; -using int64_vector_t = std::vector; + } // Convert a label into a vector of modes template @@ -53,16 +52,6 @@ mode_vector_t label_to_modes(const LabelType& label) { return mode; } -// Query extent information from an input -template -int64_vector_t get_extents(const InfoType& info) { - int64_vector_t extent; - for(std::size_t i = 0; i < info.rank(); ++i) { - extent.push_back((int64_t)info.extent(i)); - } - return extent; -} - // Compute strides in row major int64_vector_t get_strides(std::size_t N, const int64_vector_t& extent) { int64_vector_t strides; @@ -74,41 +63,56 @@ int64_vector_t get_strides(std::size_t N, const int64_vector_t& extent) { return strides; } +// Query extent information from an input +template +int64_vector_t get_extents(const InfoType& info) { + int64_vector_t extent; + for(std::size_t i = 0; i < info.rank(); ++i) { + extent.push_back((int64_t)info.extent(i)); + } + return extent; +} + // Perform tensor contraction with cuTENSOR template void cutensor_contraction(typename TensorType::label_type c_label, typename TensorType::label_type a_label, typename TensorType::label_type b_label, - typename TensorType::const_shape_reference c_shape, - typename TensorType::const_pimpl_reference A, - typename TensorType::const_pimpl_reference B, - typename TensorType::eigen_reference C) { - using element_t = typename TensorType::element_type; - using eigen_data_t = typename TensorType::eigen_data_type; + const TensorType& A, const TensorType& B, + TensorType& C) { + using element_t = typename TensorType::value_type; + + const auto a_rank = A.rank(); + const auto b_rank = B.rank(); + const auto c_rank = C.rank(); + + const auto& a_shape = A.shape(); + const auto& b_shape = B.shape(); + const auto& c_shape = C.shape(); // GEMM alpha and beta (hardcoded for now) element_t alpha = 1.0; element_t beta = 0.0; + // The extents of each tensor + int64_vector_t a_extents = get_extents(a_shape); + int64_vector_t b_extents = get_extents(b_shape); + int64_vector_t c_extents = get_extents(c_shape); + // The modes of the tensors mode_vector_t a_modes = label_to_modes(a_label); mode_vector_t b_modes = label_to_modes(b_label); mode_vector_t c_modes = label_to_modes(c_label); - // The extents of each tensor - int64_vector_t a_extents = get_extents(A); - int64_vector_t b_extents = get_extents(B); - int64_vector_t c_extents = get_extents(c_shape.as_smooth()); - // The strides of each tensor - int64_vector_t a_strides = get_strides(A.rank(), a_extents); - int64_vector_t b_strides = get_strides(B.rank(), b_extents); - int64_vector_t c_strides = get_strides(c_shape.rank(), c_extents); + int64_vector_t a_strides = get_strides(a_rank, a_extents); + int64_vector_t b_strides = get_strides(b_rank, b_extents); + int64_vector_t c_strides = get_strides(c_rank, c_extents); // The size of each tensor std::size_t a_size = sizeof(element_t) * A.size(); std::size_t b_size = sizeof(element_t) * B.size(); - std::size_t c_size = sizeof(element_t) * c_shape.size(); + std::size_t c_size = sizeof(element_t) * C.size(); // Allocate on device void *A_d, *B_d, *C_d; @@ -118,9 +122,9 @@ void cutensor_contraction(typename TensorType::label_type c_label, // Copy to data to device HANDLE_CUDA_ERROR( - cudaMemcpy(A_d, A.get_immutable_data(), a_size, cudaMemcpyHostToDevice)); + cudaMemcpy(A_d, A.data(), a_size, cudaMemcpyHostToDevice)); HANDLE_CUDA_ERROR( - cudaMemcpy(B_d, B.get_immutable_data(), b_size, cudaMemcpyHostToDevice)); + cudaMemcpy(B_d, B.data(), b_size, cudaMemcpyHostToDevice)); HANDLE_CUDA_ERROR( cudaMemcpy(C_d, C.data(), c_size, cudaMemcpyHostToDevice)); @@ -141,17 +145,17 @@ void cutensor_contraction(typename TensorType::label_type c_label, // Create Tensor Descriptors cutensorTensorDescriptor_t descA; HANDLE_CUTENSOR_ERROR(cutensorCreateTensorDescriptor( - handle, &descA, A.rank(), a_extents.data(), a_strides.data(), + handle, &descA, a_rank, a_extents.data(), a_strides.data(), traits.cutensorDataType, kAlignment)); cutensorTensorDescriptor_t descB; HANDLE_CUTENSOR_ERROR(cutensorCreateTensorDescriptor( - handle, &descB, B.rank(), b_extents.data(), b_strides.data(), + handle, &descB, b_rank, b_extents.data(), b_strides.data(), traits.cutensorDataType, kAlignment)); cutensorTensorDescriptor_t descC; HANDLE_CUTENSOR_ERROR(cutensorCreateTensorDescriptor( - handle, &descC, c_shape.rank(), c_extents.data(), c_strides.data(), + handle, &descC, c_rank, c_extents.data(), c_strides.data(), traits.cutensorDataType, kAlignment)); // Create Contraction Descriptor @@ -232,34 +236,17 @@ void cutensor_contraction(typename TensorType::label_type c_label, #undef HANDLE_CUDA_ERROR // Template instantiations -#define FUNCTION_INSTANTIATE(TYPE, RANK) \ - template void cutensor_contraction>( \ - typename EigenTensor::label_type, \ - typename EigenTensor::label_type, \ - typename EigenTensor::label_type, \ - typename EigenTensor::const_shape_reference, \ - typename EigenTensor::const_pimpl_reference, \ - typename EigenTensor::const_pimpl_reference, \ - typename EigenTensor::eigen_reference) - -#define DEFINE_CUTENSOR_CONTRACTION(TYPE) \ - FUNCTION_INSTANTIATE(TYPE, 0); \ - FUNCTION_INSTANTIATE(TYPE, 1); \ - FUNCTION_INSTANTIATE(TYPE, 2); \ - FUNCTION_INSTANTIATE(TYPE, 3); \ - FUNCTION_INSTANTIATE(TYPE, 4); \ - FUNCTION_INSTANTIATE(TYPE, 5); \ - FUNCTION_INSTANTIATE(TYPE, 6); \ - FUNCTION_INSTANTIATE(TYPE, 7); \ - FUNCTION_INSTANTIATE(TYPE, 8); \ - FUNCTION_INSTANTIATE(TYPE, 9); \ - FUNCTION_INSTANTIATE(TYPE, 10) - -TW_APPLY_FLOATING_POINT_TYPES(DEFINE_CUTENSOR_CONTRACTION); - -#undef DEFINE_CUTENSOR_CONTRACTION +#define FUNCTION_INSTANTIATE(TYPE) \ + template void cutensor_contraction>( \ + typename CUDATensor::label_type, \ + typename CUDATensor::label_type, \ + typename CUDATensor::label_type, const CUDATensor&, \ + const CUDATensor&, CUDATensor&) + +TW_APPLY_FLOATING_POINT_TYPES(FUNCTION_INSTANTIATE); + #undef FUNCTION_INSTANTIATE -} // namespace tensorwrapper::buffer::detail_ +} // namespace tensorwrapper::backends::cutensor #endif diff --git a/src/tensorwrapper/buffer/detail_/eigen_tensor.cuh b/src/tensorwrapper/backends/cutensor/cuda_tensor.cuh similarity index 78% rename from src/tensorwrapper/buffer/detail_/eigen_tensor.cuh rename to src/tensorwrapper/backends/cutensor/cuda_tensor.cuh index bc7d4e0b..b820afca 100644 --- a/src/tensorwrapper/buffer/detail_/eigen_tensor.cuh +++ b/src/tensorwrapper/backends/cutensor/cuda_tensor.cuh @@ -15,9 +15,9 @@ */ #pragma once #ifdef ENABLE_CUTENSOR -#include "eigen_tensor.hpp" +#include "cuda_tensor.hpp" -namespace tensorwrapper::buffer::detail_ { +namespace tensorwrapper::backends::cutensor { /** @brief Performs a tensor contraction on GPU * @@ -36,11 +36,9 @@ template void cutensor_contraction(typename TensorType::label_type c_label, typename TensorType::label_type a_label, typename TensorType::label_type b_label, - typename TensorType::const_shape_reference c_shape, - typename TensorType::const_pimpl_reference A, - typename TensorType::const_pimpl_reference B, - typename TensorType::eigen_reference C); + const TensorType& A, const TensorType& B, + TensorType& C); -} // namespace tensorwrapper::buffer::detail_ +} // namespace tensorwrapper::backends::cutensor #endif diff --git a/src/tensorwrapper/backends/cutensor/cuda_tensor.hpp b/src/tensorwrapper/backends/cutensor/cuda_tensor.hpp new file mode 100644 index 00000000..d0ef8eb9 --- /dev/null +++ b/src/tensorwrapper/backends/cutensor/cuda_tensor.hpp @@ -0,0 +1,73 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include +#include +#include + +namespace tensorwrapper::backends::cutensor { + +/** @brief Wraps using cuTENSOR + * + * @tparam FloatType Floating point type used for the tensor's elements + * + * N.b. The name of this class is chosen to avoid conflict with cuTENSOR. + */ +template +class CUDATensor { +private: + /// Type of *this + using my_type = CUDATensor; + + /// Read-only reference to an object of type my_type + using const_my_reference = const my_type&; + +public: + using value_type = FloatType; + using span_type = std::span; + using shape_type = shape::Smooth; + using const_shape_view = shape::SmoothView; + using label_type = dsl::DummyIndices; + using size_type = std::size_t; + + CUDATensor(span_type data, const_shape_view shape) : + m_data_(data), m_shape_(shape) {} + + void contraction_assignment(label_type this_labels, label_type lhs_labels, + label_type rhs_labels, const_my_reference lhs, + const_my_reference rhs); + + size_type rank() const noexcept { return m_shape_.rank(); } + + size_type size() const noexcept { return m_shape_.size(); } + + auto shape() const noexcept { return m_shape_; } + + auto data() noexcept { return m_data_.data(); } + + auto data() const noexcept { return m_data_.data(); } + +private: + span_type m_data_; + + const_shape_view m_shape_; +}; + +extern template class CUDATensor; +extern template class CUDATensor; + +} // namespace tensorwrapper::backends::cutensor diff --git a/src/tensorwrapper/buffer/detail_/cutensor_traits.cuh b/src/tensorwrapper/backends/cutensor/cutensor_traits.cuh similarity index 92% rename from src/tensorwrapper/buffer/detail_/cutensor_traits.cuh rename to src/tensorwrapper/backends/cutensor/cutensor_traits.cuh index c098c1c5..4846bc54 100644 --- a/src/tensorwrapper/buffer/detail_/cutensor_traits.cuh +++ b/src/tensorwrapper/backends/cutensor/cutensor_traits.cuh @@ -18,7 +18,7 @@ #include #include -namespace tensorwrapper::buffer::detail_ { +namespace tensorwrapper::backends::cutensor { // Traits for cuTENSOR based on the floating point type template @@ -36,6 +36,6 @@ struct cutensor_traits { cutensorComputeDescriptor_t descCompute = CUTENSOR_COMPUTE_DESC_64F; }; -} // namespace tensorwrapper::buffer::detail_ +} // namespace tensorwrapper::backends::cutensor #endif diff --git a/src/tensorwrapper/backends/eigen/eigen_tensor.hpp b/src/tensorwrapper/backends/eigen/eigen_tensor.hpp new file mode 100644 index 00000000..8c7c2760 --- /dev/null +++ b/src/tensorwrapper/backends/eigen/eigen_tensor.hpp @@ -0,0 +1,223 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include +#include +#include +#include +#include +#include +#include +#include + +namespace tensorwrapper::backends::eigen { + +/** @brief API for interacting with Eigen's tensor object. + * + * @tparam FloatType The floating-point type stored in the tensor. + * + * This class defines the API for interacting with Eigen's tensor objects. + * Unfortunately, Eigen's tensor objects are templated on the rank of the + * tensor (and some other stuff) which makes it hard to deal with them + * generically. This class gets the templating down to just the floating-point + * type. + * + * N.b. this class wraps Eigen::TensorMap objects, not Eigen::Tensor objects + * so as to avoid needing to transfer the data to Eigen. Idea is these classes + * can be made on-demand since they just wrap pointers. + */ +template +class EigenTensor { +private: + /// Type of *this + using my_type = EigenTensor; + +public: + /// Pointer to an object of my_type + using eigen_tensor_pointer = std::unique_ptr; + + /// Type of an element in *this + using value_type = FloatType; + + /// Type of a reference to an element in *this + using reference = value_type&; + + /// Type of a read-only reference to an element in *this + using const_reference = const value_type&; + + /// Type of a span to raw memory + using span_type = std::span; + + /// Type of a read-only span to raw memory + using const_span_type = std::span; + + /// Type used to express the shape of *this + using shape_type = shape::Smooth; + + /// Type of a view acting like a read-only shape_type + using const_shape_reference = shape::SmoothView; + + /// Type Eigen uses to express tensor rank + using eigen_rank_type = unsigned int; + + /// Type used to express sizes and extents + using size_type = std::size_t; + + /// Type used toe express multi-dimensional indices + using index_vector = std::vector; + + /// Type used to express strings + using string_type = std::string; + + /// Type of a label + using label_type = dsl::DummyIndices; + + /// Type returned by permuted_copy + using permuted_copy_return_type = + std::pair, eigen_tensor_pointer>; + + virtual ~EigenTensor() noexcept = default; + + permuted_copy_return_type permuted_copy(label_type perm, + label_type this_label) const { + return permuted_copy_(perm, this_label); + } + + /** @brief Retrieves the rank of the wrapped tensor. + * + * @return The rank of the wrapped tensor. + */ + eigen_rank_type rank() const noexcept { return rank_(); } + + /** @brief The total number of elements in *this. + * + * @return The total number of elements in *this. + * + * @throw None No throw guarantee. + */ + size_type size() const noexcept { return size_(); } + + size_type extent(eigen_rank_type i) const { + assert(i < rank()); + return extent_(i); + } + + const_reference get_elem(index_vector index) const { + assert(index.size() == rank()); + return get_elem_(std::move(index)); + } + + void set_elem(index_vector index, value_type new_value) { + assert(index.size() == rank()); + set_elem_(index, new_value); + } + + span_type data() noexcept { return data_(); } + + const_span_type data() const noexcept { return data_(); } + + void fill(value_type value) { fill_(std::move(value)); } + + string_type to_string() const { return to_string_(); } + + std::ostream& add_to_stream(std::ostream& os) const { + return add_to_stream_(os); + } + + void addition_assignment(label_type this_label, label_type lhs_label, + label_type rhs_label, const EigenTensor& lhs, + const EigenTensor& rhs) { + return addition_assignment_(this_label, lhs_label, rhs_label, lhs, rhs); + } + + void subtraction_assignment(label_type this_label, label_type lhs_label, + label_type rhs_label, const EigenTensor& lhs, + const EigenTensor& rhs) { + return subtraction_assignment_(this_label, lhs_label, rhs_label, lhs, + rhs); + } + + void hadamard_assignment(label_type this_label, label_type lhs_label, + label_type rhs_label, const EigenTensor& lhs, + const EigenTensor& rhs) { + return hadamard_assignment_(this_label, lhs_label, rhs_label, lhs, rhs); + } + + void contraction_assignment(label_type this_label, label_type lhs_label, + label_type rhs_label, const EigenTensor& lhs, + const EigenTensor& rhs) { + contraction_assignment_(this_label, lhs_label, rhs_label, lhs, rhs); + } + + void permute_assignment(label_type this_label, label_type rhs_label, + const EigenTensor& rhs) { + return permute_assignment_(this_label, rhs_label, rhs); + } + + void scalar_multiplication(label_type this_label, label_type rhs_label, + FloatType scalar, const EigenTensor& rhs) { + return scalar_multiplication_(this_label, rhs_label, scalar, rhs); + } + +protected: + explicit EigenTensor() noexcept = default; + virtual permuted_copy_return_type permuted_copy_( + label_type perm, label_type this_label) const = 0; + virtual eigen_rank_type rank_() const noexcept = 0; + virtual size_type size_() const = 0; + virtual size_type extent_(eigen_rank_type i) const = 0; + virtual const_reference get_elem_(index_vector index) const = 0; + virtual void set_elem_(index_vector index, value_type new_value) = 0; + virtual void fill_(value_type value) = 0; + virtual string_type to_string_() const = 0; + virtual std::ostream& add_to_stream_(std::ostream& os) const = 0; + virtual span_type data_() noexcept = 0; + virtual const_span_type data_() const noexcept = 0; + virtual void addition_assignment_(label_type this_label, + label_type lhs_label, + label_type rhs_label, + const EigenTensor& lhs, + const EigenTensor& rhs) = 0; + + virtual void subtraction_assignment_(label_type this_label, + label_type lhs_label, + label_type rhs_label, + const EigenTensor& lhs, + const EigenTensor& rhs) = 0; + + virtual void hadamard_assignment_(label_type this_label, + label_type lhs_label, + label_type rhs_label, + const EigenTensor& lhs, + const EigenTensor& rhs) = 0; + + virtual void permute_assignment_(label_type this_label, + label_type rhs_label, + const EigenTensor& rhs) = 0; + + virtual void scalar_multiplication_(label_type this_label, + label_type rhs_label, FloatType scalar, + const EigenTensor& rhs) = 0; + + virtual void contraction_assignment_(label_type this_label, + label_type lhs_label, + label_type rhs_label, + const EigenTensor& lhs, + const EigenTensor& rhs) = 0; +}; + +} // namespace tensorwrapper::backends::eigen diff --git a/src/tensorwrapper/backends/eigen/eigen_tensor_impl.cpp b/src/tensorwrapper/backends/eigen/eigen_tensor_impl.cpp new file mode 100644 index 00000000..53eb677e --- /dev/null +++ b/src/tensorwrapper/backends/eigen/eigen_tensor_impl.cpp @@ -0,0 +1,314 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#include "../../buffer/contraction_planner.hpp" +#include "eigen_tensor_impl.hpp" +#include +#include + +namespace tensorwrapper::backends::eigen { + +#define TPARAMS template +#define EIGEN_TENSOR EigenTensorImpl + +TPARAMS +EIGEN_TENSOR::EigenTensorImpl(std::span data, + const_shape_reference shape) : + m_tensor_(make_from_shape_(data, shape, std::make_index_sequence())) {} + +TPARAMS +auto EIGEN_TENSOR::permuted_copy_(label_type out, label_type in) const + -> permuted_copy_return_type { + using value_type = FloatType; + std::vector buffer(this->size(), value_type{0}); + std::span buffer_span(buffer.data(), buffer.size()); + + // Make a shape::Smooth object for tensor + std::vector old_shape_vec(this->rank()); + for(std::size_t i = 0; i < old_shape_vec.size(); ++i) { + old_shape_vec[i] = this->extent(i); + } + shape_type old_shape(old_shape_vec.begin(), old_shape_vec.end()); + shape_type new_shape(old_shape); + new_shape(out) = old_shape(in); + auto pnew_tensor = + std::make_unique(buffer_span, new_shape); + pnew_tensor->permute_assignment(out, in, *this); + return std::make_pair(std::move(buffer), std::move(pnew_tensor)); +} + +TPARAMS +auto EIGEN_TENSOR::get_elem_(index_vector index) const -> const_reference { + return unwrap_vector_(std::move(index), std::make_index_sequence()); +} + +TPARAMS +void EIGEN_TENSOR::set_elem_(index_vector index, value_type new_value) { + unwrap_vector_(std::move(index), std::make_index_sequence()) = + new_value; +} + +TPARAMS +void EIGEN_TENSOR::fill_(value_type value) { + std::fill(m_tensor_.data(), m_tensor_.data() + m_tensor_.size(), value); +} + +TPARAMS +auto EIGEN_TENSOR::to_string_() const -> string_type { + std::stringstream ss; + add_to_stream_(ss); + return ss.str(); +} + +TPARAMS +std::ostream& EIGEN_TENSOR::add_to_stream_(std::ostream& os) const { + os << std::fixed << std::setprecision(16); + return os << m_tensor_.format(Eigen::TensorIOFormat::Numpy()); +} + +TPARAMS +void EIGEN_TENSOR::addition_assignment_(label_type this_label, + label_type lhs_label, + label_type rhs_label, + const base_type& lhs, + const base_type& rhs) { + auto lambda = [](auto&& lhs, auto&& rhs) { return lhs + rhs; }; + element_wise_op_(lambda, this_label, lhs_label, rhs_label, lhs, rhs); +} + +TPARAMS +void EIGEN_TENSOR::subtraction_assignment_(label_type this_label, + label_type lhs_label, + label_type rhs_label, + const base_type& lhs, + const base_type& rhs) { + auto lambda = [](auto&& lhs, auto&& rhs) { return lhs - rhs; }; + element_wise_op_(lambda, this_label, lhs_label, rhs_label, lhs, rhs); +} + +TPARAMS +void EIGEN_TENSOR::hadamard_assignment_(label_type this_label, + label_type lhs_label, + label_type rhs_label, + const base_type& lhs, + const base_type& rhs) { + auto lambda = [](auto&& lhs, auto&& rhs) { return lhs * rhs; }; + element_wise_op_(lambda, this_label, lhs_label, rhs_label, lhs, rhs); +} + +TPARAMS +void EIGEN_TENSOR::permute_assignment_(label_type this_label, + label_type rhs_label, + const base_type& rhs) { + const auto* rhs_down = dynamic_cast(&rhs); + + if constexpr(Rank <= 1) { + m_tensor_ = rhs_down->m_tensor_; + return; + } else { + if(this_label != rhs_label) { // We need to permute rhs first + // Eigen adopts the opposite definition of permutation from us. + auto r_to_l = this_label.permutation(rhs_label); + // Eigen wants int objects + std::vector r_to_l2(r_to_l.begin(), r_to_l.end()); + m_tensor_ = rhs_down->m_tensor_.shuffle(r_to_l2); + } else { + m_tensor_ = rhs_down->m_tensor_; + } + } +} + +TPARAMS +void EIGEN_TENSOR::scalar_multiplication_(label_type this_label, + label_type rhs_label, + FloatType scalar, + const base_type& rhs) { + const auto* rhs_down = dynamic_cast(&rhs); + + if constexpr(Rank <= 1) { + m_tensor_ = rhs_down->m_tensor_ * scalar; + return; + } else { + if(this_label != rhs_label) { // We need to permute rhs first + auto r_to_l = rhs_label.permutation(this_label); + // Eigen wants int objects + std::vector r_to_l2(r_to_l.begin(), r_to_l.end()); + m_tensor_ = rhs_down->m_tensor_.shuffle(r_to_l2) * scalar; + } else { + m_tensor_ = rhs_down->m_tensor_ * scalar; + } + } +} + +TPARAMS +template +void EIGEN_TENSOR::element_wise_op_(OperationType op, label_type this_label, + label_type lhs_label, label_type rhs_label, + const base_type& lhs, + const base_type& rhs) { + const auto* lhs_down = dynamic_cast(&lhs); + const auto* rhs_down = dynamic_cast(&rhs); + + // Whose indices match whose? + bool this_matches_lhs = (this_label == lhs_label); + bool this_matches_rhs = (this_label == rhs_label); + bool lhs_matches_rhs = (lhs_label == rhs_label); + + // The three possible permutations we may need to apply + auto get_permutation = [](auto&& lhs_, auto&& rhs_) { + auto l_to_r = lhs_.permutation(rhs_); + return std::vector(l_to_r.begin(), l_to_r.end()); + }; + auto r_to_l = get_permutation(rhs_label, lhs_label); + auto l_to_r = get_permutation(lhs_label, rhs_label); + auto this_to_r = get_permutation(this_label, rhs_label); + + auto& lhs_eigen = lhs_down->m_tensor_; + auto& rhs_eigen = rhs_down->m_tensor_; + + if constexpr(Rank <= 1) { + m_tensor_ = op(lhs_eigen, rhs_eigen); + return; + } else { + if(this_matches_lhs && this_matches_rhs) { // No permutations + m_tensor_ = op(lhs_eigen, rhs_eigen); + } else if(this_matches_lhs) { // RHS needs permuted + m_tensor_ = op(lhs_eigen, rhs_eigen.shuffle(r_to_l)); + } else if(this_matches_rhs) { // LHS needs permuted + m_tensor_ = op(lhs_eigen.shuffle(l_to_r), rhs_eigen); + } else if(lhs_matches_rhs) { // This needs permuted + m_tensor_ = op(lhs_eigen, rhs_eigen).shuffle(this_to_r); + } else { // Everything needs permuted + auto lhs_shuffled = lhs_eigen.shuffle(l_to_r); + m_tensor_ = op(lhs_shuffled, rhs_eigen).shuffle(this_to_r); + } + } +} + +template +auto matrix_size(TensorType&& t, std::size_t row_ranks) { + std::size_t nrows = 1; + for(std::size_t i = 0; i < row_ranks; ++i) nrows *= t.extent(i); + + std::size_t ncols = 1; + const auto rank = t.rank(); + for(std::size_t i = row_ranks; i < rank; ++i) ncols *= t.extent(i); + return std::make_pair(nrows, ncols); +} + +TPARAMS +void EIGEN_TENSOR::contraction_assignment_(label_type this_label, + label_type lhs_label, + label_type rhs_label, + const base_type& lhs, + const base_type& rhs) { + buffer::ContractionPlanner plan(this_label, lhs_label, rhs_label); + + // Transpose, Transpose part of TTGT + auto&& [new_lhs_buffer, pnew_lhs_tensor] = + lhs.permuted_copy(plan.lhs_permutation(), lhs_label); + + auto&& [new_rhs_buffer, pnew_rhs_tensor] = + rhs.permuted_copy(plan.rhs_permutation(), rhs_label); + + // Gemm part of TTGT + auto olabels = plan.result_matrix_labels(); + + auto&& [out_buffer, pout_tensor] = this->permuted_copy(olabels, this_label); + + const auto [lrows, lcols] = + matrix_size(*pnew_lhs_tensor, plan.lhs_free().size()); + const auto [rrows, rcols] = + matrix_size(*pnew_rhs_tensor, plan.rhs_dummy().size()); + + // Work out the types of the matrix amd a map + constexpr auto e_dyn = ::Eigen::Dynamic; + constexpr auto e_row_major = ::Eigen::RowMajor; + using matrix_t = ::Eigen::Matrix; + using map_t = ::Eigen::Map; + + map_t lmatrix(new_lhs_buffer.data(), lrows, lcols); + map_t rmatrix(new_rhs_buffer.data(), rrows, rcols); + map_t omatrix(out_buffer.data(), lrows, rcols); + + omatrix = lmatrix * rmatrix; + + // The last transpose part of TTGT + this->permute_assignment(this_label, olabels, *pout_tensor); +} + +#undef EIGEN_TENSOR +#undef TPARAMS + +template +std::unique_ptr> make_eigen_tensor( + std::span data, shape::SmoothView shape) { + switch(shape.rank()) { + case 0: + return std::make_unique>(data, shape); + case 1: + return std::make_unique>(data, shape); + case 2: + return std::make_unique>(data, shape); + case 3: + return std::make_unique>(data, shape); + case 4: + return std::make_unique>(data, shape); + case 5: + return std::make_unique>(data, shape); + case 6: + return std::make_unique>(data, shape); + case 7: + return std::make_unique>(data, shape); + case 8: + return std::make_unique>(data, shape); + case 9: + return std::make_unique>(data, shape); + case 10: + return std::make_unique>(data, + shape); + default: + throw std::runtime_error( + "EigenTensor backend only supports ranks 0 through 10."); + } +} + +#define DEFINE_MAKE_EIGEN_TENSOR(TYPE) \ + template std::unique_ptr> make_eigen_tensor( \ + std::span data, shape::SmoothView shape); + +TW_APPLY_FLOATING_POINT_TYPES(DEFINE_MAKE_EIGEN_TENSOR); + +#undef DEFINE_MAKE_EIGEN_TENSOR + +#define DEFINE_EIGEN_TENSOR(TYPE) \ + template class EigenTensorImpl; \ + template class EigenTensorImpl; \ + template class EigenTensorImpl; \ + template class EigenTensorImpl; \ + template class EigenTensorImpl; \ + template class EigenTensorImpl; \ + template class EigenTensorImpl; \ + template class EigenTensorImpl; \ + template class EigenTensorImpl; \ + template class EigenTensorImpl; \ + template class EigenTensorImpl + +TW_APPLY_FLOATING_POINT_TYPES(DEFINE_EIGEN_TENSOR); + +#undef DEFINE_EIGEN_TENSOR + +} // namespace tensorwrapper::backends::eigen diff --git a/src/tensorwrapper/backends/eigen/eigen_tensor_impl.hpp b/src/tensorwrapper/backends/eigen/eigen_tensor_impl.hpp new file mode 100644 index 00000000..f3d0b60c --- /dev/null +++ b/src/tensorwrapper/backends/eigen/eigen_tensor_impl.hpp @@ -0,0 +1,175 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include "eigen_tensor.hpp" +#include +#include +#include +#include +#include +#include + +namespace tensorwrapper::backends::eigen { + +/// Implements EigenTensor by wrapping eigen::TensorMap +template +class EigenTensorImpl : public EigenTensor { +private: + /// Type of *this + using my_type = EigenTensorImpl; + + /// Type *this inherits from + using base_type = EigenTensor; + +public: + using eigen_tensor_type = Eigen::Tensor; + using eigen_data_type = Eigen::TensorMap; + using eigen_reference = eigen_data_type&; + using const_eigen_reference = const eigen_data_type&; + + ///@{ + using typename base_type::const_reference; + using typename base_type::const_shape_reference; + using typename base_type::const_span_type; + using typename base_type::eigen_rank_type; + using typename base_type::index_vector; + using typename base_type::label_type; + using typename base_type::permuted_copy_return_type; + using typename base_type::reference; + using typename base_type::shape_type; + using typename base_type::size_type; + using typename base_type::span_type; + using typename base_type::string_type; + using typename base_type::value_type; + ///@} + + EigenTensorImpl(std::span data, const_shape_reference shape); + +protected: + permuted_copy_return_type permuted_copy_( + label_type perm, label_type this_label) const override; + + /// Implement rank by returning template parameter + eigen_rank_type rank_() const noexcept override { return Rank; } + + /// Calls Eigen's size() method to implement size() + size_type size_() const noexcept override { return m_tensor_.size(); } + + /// Calls Eigen's dimension(i) method to implement extent(i) + size_type extent_(eigen_rank_type i) const override { + return m_tensor_.dimension(i); + } + + /// Unwraps index vector into Eigen's operator() to get element + const_reference get_elem_(index_vector index) const override; + + /// Unwraps index vector into Eigen's operator() to set element + void set_elem_(index_vector index, value_type new_value) override; + + virtual span_type data_() noexcept override { + return {m_tensor_.data(), size_type(m_tensor_.size())}; + } + virtual const_span_type data_() const noexcept override { + return {m_tensor_.data(), size_type(m_tensor_.size())}; + } + + /// Calls std::fill to set the values + void fill_(value_type value) override; + + /// Calls add_to_stream_ on a stringstream to implement to_string + string_type to_string_() const override; + + /// Relies on Eigen's operator<< to add to stream + std::ostream& add_to_stream_(std::ostream& os) const override; + + void addition_assignment_(label_type this_label, label_type lhs_label, + label_type rhs_label, const base_type& lhs, + const base_type& rhs) override; + + void subtraction_assignment_(label_type this_label, label_type lhs_label, + label_type rhs_label, const base_type& lhs, + const base_type& rhs) override; + + void hadamard_assignment_(label_type this_label, label_type lhs_label, + label_type rhs_label, const base_type& lhs, + const base_type& rhs) override; + + void permute_assignment_(label_type this_label, label_type rhs_label, + const base_type& rhs) override; + + void scalar_multiplication_(label_type this_label, label_type rhs_label, + FloatType scalar, + const base_type& rhs) override; + + void contraction_assignment_(label_type this_labels, label_type lhs_labels, + label_type rhs_labels, const base_type& lhs, + const base_type& rhs) override; + +private: + // Code factorization for implementing element-wise operations + template + void element_wise_op_(OperationType op, label_type this_label, + label_type lhs_label, label_type rhs_label, + const base_type& lhs, const base_type& rhs); + + // Handles TMP needed to create an Eigen TensorMap from a Smooth object + template + auto make_from_shape_(std::span data, + const_shape_reference shape, + std::index_sequence) { + return eigen_data_type(data.data(), shape.extent(I)...); + } + + // Gets an element from the Eigen Tensor by unwrapping a std::vector + template + reference unwrap_vector_(index_vector index, std::index_sequence) { + return m_tensor_(tensorwrapper::detail_::to_long(index.at(I))...); + } + + // Same as mutable version, but result is read-only + template + const_reference unwrap_vector_(index_vector index, + std::index_sequence) const { + return m_tensor_(tensorwrapper::detail_::to_long(index.at(I))...); + } + + // The Eigen tensor *this wraps + eigen_data_type m_tensor_; +}; + +template +std::unique_ptr> make_eigen_tensor( + std::span data, shape::SmoothView shape); + +#define DECLARE_EIGEN_TENSOR(TYPE) \ + extern template class EigenTensorImpl; \ + extern template class EigenTensorImpl; \ + extern template class EigenTensorImpl; \ + extern template class EigenTensorImpl; \ + extern template class EigenTensorImpl; \ + extern template class EigenTensorImpl; \ + extern template class EigenTensorImpl; \ + extern template class EigenTensorImpl; \ + extern template class EigenTensorImpl; \ + extern template class EigenTensorImpl; \ + extern template class EigenTensorImpl + +TW_APPLY_FLOATING_POINT_TYPES(DECLARE_EIGEN_TENSOR); + +#undef DECLARE_EIGEN_TENSOR + +} // namespace tensorwrapper::backends::eigen diff --git a/src/tensorwrapper/buffer/detail_/addition_visitor.hpp b/src/tensorwrapper/buffer/detail_/addition_visitor.hpp new file mode 100644 index 00000000..4e021e8a --- /dev/null +++ b/src/tensorwrapper/buffer/detail_/addition_visitor.hpp @@ -0,0 +1,38 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include + +namespace tensorwrapper::buffer::detail_ { + +/** @brief Dispatches to the appropriate backend based on the FP type. + * + * + * + */ +class AdditionVisitor { +public: + // AdditionVisitor(shape, permutation, shape, permutation) + template + void operator()(std::span lhs, std::span rhs) { + // auto lhs_wrapped = backends::eigen::wrap_span(lhs); + // auto rhs_wrapped = backends::eigen::wrap_span(rhs); + for(std::size_t i = 0; i < lhs.size(); ++i) lhs[i] += rhs[i]; + } +}; + +} // namespace tensorwrapper::buffer::detail_ diff --git a/src/tensorwrapper/buffer/detail_/eigen_dispatch.hpp b/src/tensorwrapper/buffer/detail_/eigen_dispatch.hpp new file mode 100644 index 00000000..90c26509 --- /dev/null +++ b/src/tensorwrapper/buffer/detail_/eigen_dispatch.hpp @@ -0,0 +1,60 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include +#include + +namespace tensorwrapper::buffer::detail_ { + +constexpr std::size_t MaxEigenRank = 8; + +template +using eigen_tensor_type = eigen::Tensor; + +template +using eigen_tensor_map = eigen::TensorMap>; + +template +auto wrap_tensor(std::span s, const shape::Smooth& shape) { + using tensor_type = eigen::Tensor; + using map_type = eigen::TensorMap; + + if constexpr(Rank > MaxEigenRank) { + static_assert( + Rank <= MaxEigenRank, + "Eigen tensors of rank > MaxEigenRank are not supported."); + } else { + if(shape.rank() == Rank) return variant_type(map_type(s)); + } +} + +template +auto eigen_dispatch_impl(VisitorType&& visitor, + eigen::TensorMap>& A, + Args&&... args) { + return visitor(A, std::forward(args)...); +} + +template +auto eigen_tensor_dispatch(std::span s, shape::Smooth shape, + Args&&... args) { + using tensor_type = eigen::Tensor; +} + +} // namespace tensorwrapper::buffer::detail_ diff --git a/src/tensorwrapper/buffer/detail_/eigen_tensor.cpp b/src/tensorwrapper/buffer/detail_/eigen_tensor.cpp index dde5a7e0..84096fb5 100644 --- a/src/tensorwrapper/buffer/detail_/eigen_tensor.cpp +++ b/src/tensorwrapper/buffer/detail_/eigen_tensor.cpp @@ -18,10 +18,6 @@ #include "../contraction_planner.hpp" #include "eigen_tensor.hpp" -#ifdef ENABLE_CUTENSOR -#include "eigen_tensor.cuh" -#endif - namespace tensorwrapper::buffer::detail_ { #define TPARAMS template @@ -100,16 +96,6 @@ void EIGEN_TENSOR::contraction_assignment_(label_type olabels, const_pimpl_reference rhs) { ContractionPlanner plan(olabels, llabels, rlabels); -#ifdef ENABLE_CUTENSOR - // Prepare m_tensor_ - m_tensor_ = allocate_from_shape_(result_shape.as_smooth(), - std::make_index_sequence()); - m_tensor_.setZero(); - - // Dispatch to cuTENSOR - cutensor_contraction(olabels, llabels, rlabels, result_shape, lhs, - rhs, m_tensor_); -#else auto lt = lhs.clone(); auto rt = rhs.clone(); lt->permute_assignment(plan.lhs_permutation(), llabels, lhs); @@ -154,7 +140,7 @@ void EIGEN_TENSOR::contraction_assignment_(label_type olabels, } else { m_tensor_ = tensor; } -#endif + mark_for_rehash_(); } diff --git a/src/tensorwrapper/buffer/detail_/mdbuffer_pimpl.hpp b/src/tensorwrapper/buffer/detail_/mdbuffer_pimpl.hpp new file mode 100644 index 00000000..6f410098 --- /dev/null +++ b/src/tensorwrapper/buffer/detail_/mdbuffer_pimpl.hpp @@ -0,0 +1,53 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include +#include + +namespace tensorwrapper::buffer::detail_ { + +class MDBufferPIMPL { +public: + using parent_type = tensorwrapper::buffer::MDBuffer; + using traits_type = tensorwrapper::types::ClassTraits; + + /// Add types to public API + ///@{ + using value_type = typename traits_type::value_type; + using rank_type = typename traits_type::rank_type; + using buffer_type = typename traits_type::buffer_type; + using shape_type = typename traits_type::shape_type; + ///@} + + MDBufferPIMPL(shape_type shape, buffer_type buffer) noexcept : + m_shape_(std::move(shape)), m_buffer_(std::move(buffer)) {} + + auto& shape() noexcept { return m_shape_; } + + const auto& shape() const noexcept { return m_shape_; } + + auto& buffer() noexcept { return m_buffer_; } + + const auto& buffer() const noexcept { return m_buffer_; } + +private: + shape_type m_shape_; + + buffer_type m_buffer_; +}; + +} // namespace tensorwrapper::buffer::detail_ diff --git a/src/tensorwrapper/buffer/mdbuffer.cpp b/src/tensorwrapper/buffer/mdbuffer.cpp new file mode 100644 index 00000000..fe92be9c --- /dev/null +++ b/src/tensorwrapper/buffer/mdbuffer.cpp @@ -0,0 +1,56 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#include "detail_/addition_visitor.hpp" +#include "detail_/mdbuffer_pimpl.hpp" +#include +#include + +namespace tensorwrapper::buffer { + +MDBuffer::MDBuffer() noexcept : m_pimpl_(nullptr) {} + +MDBuffer::MDBuffer(shape_type shape, buffer_type buffer) : + MDBuffer(std::make_unique(std::move(shape), + std::move(buffer))) {} + +MDBuffer::MDBuffer(pimpl_pointer pimpl) noexcept : m_pimpl_(std::move(pimpl)) {} + +auto MDBuffer::rank() const -> rank_type { + assert_pimpl_(); + return m_pimpl_->shape().rank(); +} + +bool MDBuffer::has_pimpl_() const noexcept { return m_pimpl_ != nullptr; } + +void MDBuffer::assert_pimpl_() const { + if(!has_pimpl_()) { + throw std::runtime_error( + "MDBuffer has no PIMPL. Was it default constructed?"); + } +} + +auto MDBuffer::pimpl_() -> pimpl_type& { + assert_pimpl_(); + return *m_pimpl_; +} + +auto MDBuffer::pimpl_() const -> const pimpl_type& { + assert_pimpl_(); + return *m_pimpl_; +} + +} // namespace tensorwrapper::buffer diff --git a/tests/cxx/unit_tests/tensorwrapper/backends/cutensor/cuda_tensor.cpp b/tests/cxx/unit_tests/tensorwrapper/backends/cutensor/cuda_tensor.cpp new file mode 100644 index 00000000..a71352ca --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/backends/cutensor/cuda_tensor.cpp @@ -0,0 +1,99 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#include "../../testing/testing.hpp" +#include "../testing/contraction_assignment.hpp" +#include + +using namespace tensorwrapper; +using namespace tensorwrapper::backends::cutensor; + +using supported_fp_types = std::tuple; + +TEMPLATE_LIST_TEST_CASE("CUDATensor", "", supported_fp_types) { + using tensor_type = CUDATensor; + using shape_type = typename tensor_type::shape_type; + using label_type = typename tensor_type::label_type; + + std::vector data(16); + for(std::size_t i = 0; i < data.size(); ++i) + data[i] = static_cast(i); + + std::span data_span(data.data(), data.size()); + + shape_type scalar_shape({}); + shape_type vector_shape({16}); + shape_type matrix_shape({4, 4}); + shape_type tensor3_shape({2, 2, 4}); + shape_type tensor4_shape({2, 2, 2, 2}); + + tensor_type scalar(data_span, scalar_shape); + tensor_type vector(data_span, vector_shape); + tensor_type matrix(data_span, matrix_shape); + tensor_type tensor3(data_span, tensor3_shape); + tensor_type tensor4(data_span, tensor4_shape); + + SECTION("rank") { + REQUIRE(scalar.rank() == 0); + REQUIRE(vector.rank() == 1); + REQUIRE(matrix.rank() == 2); + REQUIRE(tensor3.rank() == 3); + REQUIRE(tensor4.rank() == 4); + } + + SECTION("size") { + REQUIRE(scalar.size() == 1); + REQUIRE(vector.size() == 16); + REQUIRE(matrix.size() == 16); + REQUIRE(tensor3.size() == 16); + REQUIRE(tensor4.size() == 16); + } + + SECTION("shape") { + REQUIRE(scalar.shape() == scalar_shape); + REQUIRE(vector.shape() == vector_shape); + REQUIRE(matrix.shape() == matrix_shape); + REQUIRE(tensor3.shape() == tensor3_shape); + REQUIRE(tensor4.shape() == tensor4_shape); + } + + SECTION("data()") { + REQUIRE(scalar.data() == data.data()); + REQUIRE(vector.data() == data.data()); + REQUIRE(matrix.data() == data.data()); + REQUIRE(tensor3.data() == data.data()); + REQUIRE(tensor4.data() == data.data()); + } + + SECTION("data() const") { + REQUIRE(std::as_const(scalar).data() == data.data()); + REQUIRE(std::as_const(vector).data() == data.data()); + REQUIRE(std::as_const(matrix).data() == data.data()); + REQUIRE(std::as_const(tensor3).data() == data.data()); + REQUIRE(std::as_const(tensor4).data() == data.data()); + } + + SECTION("contraction_assignment") { +#ifdef ENABLE_CUTESNSOR + testing::contraction_assignment(); +#else + label_type label(""); + REQUIRE_THROWS_AS( + scalar.contraction_assignment(label, label, label, scalar, scalar), + std::runtime_error); +#endif + } +} diff --git a/tests/cxx/unit_tests/tensorwrapper/backends/eigen/eigen_tensor_impl.cpp b/tests/cxx/unit_tests/tensorwrapper/backends/eigen/eigen_tensor_impl.cpp new file mode 100644 index 00000000..7d1ec863 --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/backends/eigen/eigen_tensor_impl.cpp @@ -0,0 +1,393 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#include "../../testing/testing.hpp" +#include "../testing/addition_assignment.hpp" +#include "../testing/contraction_assignment.hpp" +#include "../testing/hadamard_assignment.hpp" +#include "../testing/permute_assignment.hpp" +#include "../testing/scalar_multiplication.hpp" +#include "../testing/subtraction_assignment.hpp" +#include + +using namespace tensorwrapper; +using namespace tensorwrapper::backends::eigen; + +TEMPLATE_LIST_TEST_CASE("EigenTensorImpl", "", types::floating_point_types) { + using scalar_type = EigenTensorImpl; + using vector_type = EigenTensorImpl; + using matrix_type = EigenTensorImpl; + using tensor3_type = EigenTensorImpl; + using tensor4_type = EigenTensorImpl; + + std::vector data(16); + for(std::size_t i = 0; i < data.size(); ++i) + data[i] = static_cast(i); + + std::span data_span(data.data(), data.size()); + + using shape_type = scalar_type::shape_type; + + shape_type scalar_shape({}); + shape_type vector_shape({16}); + shape_type matrix_shape({4, 4}); + shape_type tensor3_shape({2, 2, 4}); + shape_type tensor4_shape({2, 2, 2, 2}); + + scalar_type scalar(data_span, scalar_shape); + vector_type vector(data_span, vector_shape); + matrix_type matrix(data_span, matrix_shape); + tensor3_type tensor3(data_span, tensor3_shape); + tensor4_type tensor4(data_span, tensor4_shape); + + SECTION("permuted_copy") { + using label_type = typename scalar_type::label_type; + SECTION("scalar") { + label_type label(""); + auto&& [scalar_buffer, pscalar] = + scalar.permuted_copy(label, label); + REQUIRE(scalar_buffer.size() == 1); + REQUIRE(pscalar->get_elem({}) == data[0]); + REQUIRE(scalar_buffer.data() != data.data()); // Ensure a copy + } + + SECTION("vector") { + label_type i("i"); + auto&& [vector_buffer, pvector] = vector.permuted_copy(i, i); + REQUIRE(vector_buffer.size() == 16); + for(std::size_t idx = 0; idx < 16; ++idx) { + REQUIRE(pvector->get_elem({idx}) == data[idx]); + } + REQUIRE(vector_buffer.data() != data.data()); // Ensure a copy + } + + SECTION("matrix") { + label_type ij("i,j"); + label_type ji("j,i"); + SECTION("no permutation") { + auto&& [matrix_buffer, pmatrix] = matrix.permuted_copy(ij, ij); + REQUIRE(matrix_buffer.size() == 16); + for(std::size_t idx = 0; idx < 4; ++idx) { + for(std::size_t jdx = 0; jdx < 4; ++jdx) { + REQUIRE(pmatrix->get_elem({idx, jdx}) == + data[idx * 4 + jdx]); + } + } + REQUIRE(matrix_buffer.data() != data.data()); // Ensure a copy + } + SECTION("Permutation") { + auto&& [matrix_buffer, pmatrix] = matrix.permuted_copy(ji, ij); + REQUIRE(matrix_buffer.size() == 16); + for(std::size_t idx = 0; idx < 4; ++idx) { + for(std::size_t jdx = 0; jdx < 4; ++jdx) { + REQUIRE(pmatrix->get_elem({jdx, idx}) == + data[idx * 4 + jdx]); + } + } + REQUIRE(matrix_buffer.data() != data.data()); // Ensure a copy + } + } + + SECTION("Rank 3 tensor") { + label_type ijk("i,j,k"); + label_type jik("j,i,k"); + SECTION("no permutation") { + auto&& [t_buffer, pt] = tensor3.permuted_copy(ijk, ijk); + REQUIRE(t_buffer.size() == 16); + for(std::size_t idx = 0; idx < 2; ++idx) { + for(std::size_t jdx = 0; jdx < 2; ++jdx) { + for(std::size_t kdx = 0; kdx < 4; ++kdx) { + REQUIRE(pt->get_elem({idx, jdx, kdx}) == + data[idx * 8 + jdx * 4 + kdx]); + } + } + } + REQUIRE(t_buffer.data() != data.data()); // Ensure a copy + } + SECTION("Permutation") { + auto&& [t_buffer, pt] = tensor3.permuted_copy(jik, ijk); + REQUIRE(t_buffer.size() == 16); + for(std::size_t idx = 0; idx < 2; ++idx) { + for(std::size_t jdx = 0; jdx < 2; ++jdx) { + for(std::size_t kdx = 0; kdx < 4; ++kdx) { + REQUIRE(pt->get_elem({jdx, idx, kdx}) == + data[idx * 8 + jdx * 4 + kdx]); + } + } + } + REQUIRE(t_buffer.data() != data.data()); // Ensure a copy + } + } + + SECTION("Rank 4 tensor") { + label_type ijkl("i,j,k,l"); + label_type jikl("j,i,k,l"); + SECTION("no permutation") { + auto&& [t_buffer, pt] = tensor4.permuted_copy(ijkl, ijkl); + REQUIRE(t_buffer.size() == 16); + for(std::size_t idx = 0; idx < 2; ++idx) { + for(std::size_t jdx = 0; jdx < 2; ++jdx) { + for(std::size_t kdx = 0; kdx < 2; ++kdx) { + for(std::size_t ldx = 0; ldx < 2; ++ldx) { + REQUIRE( + pt->get_elem({idx, jdx, kdx, ldx}) == + data[idx * 8 + jdx * 4 + kdx * 2 + ldx]); + } + } + } + } + REQUIRE(t_buffer.data() != data.data()); // Ensure a copy + } + SECTION("Permutation") { + auto&& [t_buffer, pt] = tensor4.permuted_copy(jikl, ijkl); + REQUIRE(t_buffer.size() == 16); + for(std::size_t idx = 0; idx < 2; ++idx) { + for(std::size_t jdx = 0; jdx < 2; ++jdx) { + for(std::size_t kdx = 0; kdx < 2; ++kdx) { + for(std::size_t ldx = 0; ldx < 2; ++ldx) { + REQUIRE( + pt->get_elem({jdx, idx, kdx, ldx}) == + data[idx * 8 + jdx * 4 + kdx * 2 + ldx]); + } + } + } + } + REQUIRE(t_buffer.data() != data.data()); // Ensure a copy + } + } + } + + SECTION("rank") { + REQUIRE(scalar.rank() == 0); + REQUIRE(vector.rank() == 1); + REQUIRE(matrix.rank() == 2); + REQUIRE(tensor3.rank() == 3); + REQUIRE(tensor4.rank() == 4); + } + + SECTION("size") { + REQUIRE(scalar.size() == 1); + REQUIRE(vector.size() == 16); + REQUIRE(matrix.size() == 16); + REQUIRE(tensor3.size() == 16); + REQUIRE(tensor4.size() == 16); + } + + SECTION("extent") { + REQUIRE(vector.extent(0) == 16); + + REQUIRE(matrix.extent(0) == 4); + REQUIRE(matrix.extent(1) == 4); + + REQUIRE(tensor3.extent(0) == 2); + REQUIRE(tensor3.extent(1) == 2); + REQUIRE(tensor3.extent(2) == 4); + + REQUIRE(tensor4.extent(0) == 2); + REQUIRE(tensor4.extent(1) == 2); + REQUIRE(tensor4.extent(2) == 2); + REQUIRE(tensor4.extent(3) == 2); + } + + SECTION("get_elem") { + REQUIRE(scalar.get_elem({}) == data[0]); + + REQUIRE(vector.get_elem({0}) == data[0]); + REQUIRE(vector.get_elem({15}) == data[15]); + + REQUIRE(matrix.get_elem({0, 0}) == data[0]); + REQUIRE(matrix.get_elem({3, 3}) == data[15]); + + REQUIRE(tensor3.get_elem({0, 0, 0}) == data[0]); + REQUIRE(tensor3.get_elem({1, 1, 3}) == data[15]); + + REQUIRE(tensor4.get_elem({0, 0, 0, 0}) == data[0]); + REQUIRE(tensor4.get_elem({1, 1, 1, 1}) == data[15]); + } + + SECTION("set_elem") { + TestType corr(42); + scalar.set_elem({}, 42); + REQUIRE(scalar.get_elem({}) == corr); + + vector.set_elem({5}, 42); + REQUIRE(vector.get_elem({5}) == corr); + + matrix.set_elem({2, 2}, 42); + REQUIRE(matrix.get_elem({2, 2}) == corr); + + tensor3.set_elem({1, 0, 3}, 42); + REQUIRE(tensor3.get_elem({1, 0, 3}) == corr); + + tensor4.set_elem({0, 1, 1, 0}, 42); + REQUIRE(tensor4.get_elem({0, 1, 1, 0}) == corr); + } + + SECTION("data()") { + REQUIRE(scalar.data().data() == data.data()); + REQUIRE(vector.data().data() == data.data()); + REQUIRE(matrix.data().data() == data.data()); + REQUIRE(tensor3.data().data() == data.data()); + REQUIRE(tensor4.data().data() == data.data()); + } + + SECTION("data() const") { + REQUIRE(std::as_const(scalar).data().data() == data.data()); + REQUIRE(std::as_const(vector).data().data() == data.data()); + REQUIRE(std::as_const(matrix).data().data() == data.data()); + REQUIRE(std::as_const(tensor3).data().data() == data.data()); + REQUIRE(std::as_const(tensor4).data().data() == data.data()); + } + + SECTION("fill") { + TestType corr(7); + SECTION("scalar") { + scalar.fill(corr); + REQUIRE(scalar.get_elem({}) == corr); + } + + SECTION("vector") { + vector.fill(corr); + for(std::size_t i = 0; i < vector.size(); ++i) + REQUIRE(vector.get_elem({i}) == corr); + } + + SECTION("matrix") { + matrix.fill(corr); + for(std::size_t i = 0; i < matrix.extent(0); ++i) + for(std::size_t j = 0; j < matrix.extent(1); ++j) + REQUIRE(matrix.get_elem({i, j}) == corr); + } + + SECTION("rank 3 tensor") { + tensor3.fill(corr); + for(std::size_t i = 0; i < tensor3.extent(0); ++i) + for(std::size_t j = 0; j < tensor3.extent(1); ++j) + for(std::size_t k = 0; k < tensor3.extent(2); ++k) + REQUIRE(tensor3.get_elem({i, j, k}) == corr); + } + + SECTION("rank 4 tensor") { + tensor4.fill(corr); + for(std::size_t i = 0; i < tensor4.extent(0); ++i) + for(std::size_t j = 0; j < tensor4.extent(1); ++j) + for(std::size_t k = 0; k < tensor4.extent(2); ++k) + for(std::size_t l = 0; l < tensor4.extent(3); ++l) + REQUIRE(tensor4.get_elem({i, j, k, l}) == corr); + } + } + + SECTION("addition_assignment") { + SECTION("scalar") { + testing::scalar_addition_assignment(); + } + + SECTION("vector") { + testing::vector_addition_assignment(); + } + + SECTION("matrix") { + testing::matrix_addition_assignment(); + } + + SECTION("rank 3 tensor") { + testing::tensor3_addition_assignment(); + } + + SECTION("rank 4 tensor") { + testing::tensor4_addition_assignment(); + } + } + + SECTION("subtraction_assignment") { + SECTION("scalar") { + testing::scalar_subtraction_assignment(); + } + + SECTION("vector") { + testing::vector_subtraction_assignment(); + } + + SECTION("matrix") { + testing::matrix_subtraction_assignment(); + } + + SECTION("rank 3 tensor") { + testing::tensor3_subtraction_assignment(); + } + + SECTION("rank 4 tensor") { + testing::tensor4_subtraction_assignment(); + } + } + + SECTION("hadamard_assignment") { + SECTION("scalar") { + testing::scalar_hadamard_assignment(); + } + + SECTION("vector") { + testing::vector_hadamard_assignment(); + } + + SECTION("matrix") { + testing::matrix_hadamard_assignment(); + } + + SECTION("rank 3 tensor") { + testing::tensor3_hadamard_assignment(); + } + + SECTION("rank 4 tensor") { + testing::tensor4_hadamard_assignment(); + } + } + + SECTION("permute_assignment") { + SECTION("scalar") { testing::scalar_permute_assignment(); } + SECTION("vector") { testing::vector_permute_assignment(); } + SECTION("matrix") { testing::matrix_permute_assignment(); } + SECTION("rank 3 tensor") { + testing::tensor3_permute_assignment(); + } + SECTION("rank 4 tensor") { + testing::tensor4_permute_assignment(); + } + } + + SECTION("scalar_multiplication") { + SECTION("scalar") { + testing::scalar_scalar_multiplication(); + } + SECTION("vector") { + testing::vector_scalar_multiplication(); + } + SECTION("matrix") { + testing::matrix_scalar_multiplication(); + } + SECTION("rank 3 tensor") { + testing::tensor3_scalar_multiplication(); + } + SECTION("rank 4 tensor") { + testing::tensor4_scalar_multiplication(); + } + } + + SECTION("contraction_assignment") { + testing::contraction_assignment_tests< + scalar_type, vector_type, matrix_type, tensor3_type, tensor4_type>(); + } +} diff --git a/tests/cxx/unit_tests/tensorwrapper/backends/testing/addition_assignment.hpp b/tests/cxx/unit_tests/tensorwrapper/backends/testing/addition_assignment.hpp new file mode 100644 index 00000000..9e17d587 --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/backends/testing/addition_assignment.hpp @@ -0,0 +1,72 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include "elementwise_op.hpp" + +namespace tensorwrapper::testing { + +template +void scalar_addition_assignment() { + auto the_op = [](auto&& out_idx, auto&& lhs_idx, auto&& rhs_idx, + auto&& result, auto&& t0, auto&& t1) { + result.addition_assignment(out_idx, lhs_idx, rhs_idx, t0, t1); + }; + auto corr_op = [](auto a, auto b) { return a + b; }; + scalar_binary_assignment(the_op, corr_op); +} + +template +void vector_addition_assignment() { + auto the_op = [](auto&& out_idx, auto&& lhs_idx, auto&& rhs_idx, + auto&& result, auto&& t0, auto&& t1) { + result.addition_assignment(out_idx, lhs_idx, rhs_idx, t0, t1); + }; + auto corr_op = [](auto a, auto b) { return a + b; }; + vector_binary_assignment(the_op, corr_op); +} + +template +void matrix_addition_assignment() { + auto the_op = [](auto&& out_idx, auto&& lhs_idx, auto&& rhs_idx, + auto&& result, auto&& t0, auto&& t1) { + result.addition_assignment(out_idx, lhs_idx, rhs_idx, t0, t1); + }; + auto corr_op = [](auto a, auto b) { return a + b; }; + matrix_binary_assignment(the_op, corr_op); +} + +template +void tensor3_addition_assignment() { + auto the_op = [](auto&& out_idx, auto&& lhs_idx, auto&& rhs_idx, + auto&& result, auto&& t0, auto&& t1) { + result.addition_assignment(out_idx, lhs_idx, rhs_idx, t0, t1); + }; + auto corr_op = [](auto a, auto b) { return a + b; }; + tensor3_binary_assignment(the_op, corr_op); +} + +template +void tensor4_addition_assignment() { + auto the_op = [](auto&& out_idx, auto&& lhs_idx, auto&& rhs_idx, + auto&& result, auto&& t0, auto&& t1) { + result.addition_assignment(out_idx, lhs_idx, rhs_idx, t0, t1); + }; + auto corr_op = [](auto a, auto b) { return a + b; }; + tensor4_binary_assignment(the_op, corr_op); +} + +} // namespace tensorwrapper::testing diff --git a/tests/cxx/unit_tests/tensorwrapper/backends/testing/contraction_assignment.hpp b/tests/cxx/unit_tests/tensorwrapper/backends/testing/contraction_assignment.hpp new file mode 100644 index 00000000..b6837d0e --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/backends/testing/contraction_assignment.hpp @@ -0,0 +1,217 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include +#include + +namespace tensorwrapper::testing { + +template +void contraction_assignment_tests() { + using scalar_value_type = typename ScalarType::value_type; + using vector_value_type = typename VectorType::value_type; + using matrix_value_type = typename MatrixType::value_type; + using tensor3_value_type = typename Tensor3Type::value_type; + using tensor4_value_type = typename Tensor4Type::value_type; + + using shape_type = typename ScalarType::shape_type; + using label_type = typename ScalarType::label_type; + + std::vector scalar_data(1, scalar_value_type(42.0)); + std::vector vector_data(2, vector_value_type(0.0)); + std::vector matrix_data(4, matrix_value_type(0.0)); + std::vector tensor3_data(8, tensor3_value_type(0.0)); + std::vector tensor4_data(16, tensor4_value_type(0.0)); + + for(std::size_t i = 0; i < vector_data.size(); ++i) + vector_data[i] = scalar_value_type(i + 1.0); + + for(std::size_t i = 0; i < matrix_data.size(); ++i) + matrix_data[i] = scalar_value_type(i + 1.0); + + for(std::size_t i = 0; i < tensor3_data.size(); ++i) + tensor3_data[i] = scalar_value_type(i + 1.0); + + std::span vector_data_span(vector_data.data(), + vector_data.size()); + std::span matrix_data_span(matrix_data.data(), + matrix_data.size()); + std::span tensor3_data_span(tensor3_data.data(), + tensor3_data.size()); + std::span tensor4_data_span(tensor4_data.data(), + tensor4_data.size()); + + shape_type scalar_shape{}; + shape_type vector_shape{2}; + shape_type matrix_shape{2, 2}; + shape_type tensor3_shape{2, 2, 2}; + shape_type tensor4_shape{2, 2, 2, 2}; + + ScalarType scalar(scalar_data, scalar_shape); + VectorType vector(vector_data_span, vector_shape); + MatrixType matrix(matrix_data_span, matrix_shape); + Tensor3Type tensor3(tensor3_data_span, tensor3_shape); + Tensor4Type tensor4(tensor4_data, shape_type{2, 2, 2, 2}); + + SECTION("scalar,scalar->") { + label_type o(""); + label_type l(""); + label_type r(""); + scalar.contraction_assignment(o, l, r, scalar, scalar); + + REQUIRE(scalar.get_elem({}) == scalar_value_type(42.0 * 42.0)); + } + + SECTION("i,i->") { + label_type o(""); + label_type l("i"); + label_type r("i"); + scalar.contraction_assignment(o, l, r, vector, vector); + REQUIRE(scalar.get_elem({}) == vector_value_type(5.0)); + } + + SECTION("i,ij->j") { + label_type o("j"); + label_type l("i"); + label_type r("i,j"); + vector.contraction_assignment(o, l, r, vector, matrix); + REQUIRE(vector.get_elem({0}) == vector_value_type(7.0)); + REQUIRE(vector.get_elem({1}) == vector_value_type(10.0)); + } + + SECTION("ij,ji->") { + label_type o(""); + label_type l("i,j"); + label_type r("j,i"); + scalar.contraction_assignment(o, l, r, matrix, matrix); + + REQUIRE(scalar.get_elem({}) == matrix_value_type(29.0)); + } + + SECTION("ij,jk->ik") { + label_type o("i,k"); + label_type l("i,j"); + label_type r("j,k"); + matrix.contraction_assignment(o, l, r, matrix, matrix); + + REQUIRE(matrix.get_elem({0, 0}) == matrix_value_type(7.0)); + REQUIRE(matrix.get_elem({0, 1}) == matrix_value_type(10.0)); + REQUIRE(matrix.get_elem({1, 0}) == matrix_value_type(15.0)); + REQUIRE(matrix.get_elem({1, 1}) == matrix_value_type(22.0)); + } + + SECTION("ijk,ijk->") { + label_type o(""); + label_type l("i,j,k"); + label_type r("i,j,k"); + scalar.contraction_assignment(o, l, r, tensor3, tensor3); + + REQUIRE(scalar.get_elem({}) == scalar_value_type(204.0)); + } + + SECTION("ijk,jik->") { + label_type o(""); + label_type l("i,j,k"); + label_type r("j,i,k"); + scalar.contraction_assignment(o, l, r, tensor3, tensor3); + + REQUIRE(scalar.get_elem({}) == scalar_value_type(196.0)); + } + + SECTION("ijk,jkl->il") { + label_type o("i,l"); + label_type l("i,j,k"); + label_type r("j,k,l"); + matrix.contraction_assignment(o, l, r, tensor3, tensor3); + + REQUIRE(matrix.get_elem({0, 0}) == matrix_value_type(50.0)); + REQUIRE(matrix.get_elem({0, 1}) == matrix_value_type(60.0)); + REQUIRE(matrix.get_elem({1, 0}) == matrix_value_type(114.0)); + REQUIRE(matrix.get_elem({1, 1}) == matrix_value_type(140.0)); + } + + SECTION("ijk,jlk->il") { + label_type o("i,l"); + label_type l("i,j,k"); + label_type r("j,l,k"); + matrix.contraction_assignment(o, l, r, tensor3, tensor3); + + REQUIRE(matrix.get_elem({0, 0}) == matrix_value_type(44.0)); + REQUIRE(matrix.get_elem({0, 1}) == matrix_value_type(64.0)); + REQUIRE(matrix.get_elem({1, 0}) == matrix_value_type(100.0)); + REQUIRE(matrix.get_elem({1, 1}) == matrix_value_type(152.0)); + } + + SECTION("ijk,jlk->li") { + label_type o("l,i"); + label_type l("i,j,k"); + label_type r("j,l,k"); + matrix.contraction_assignment(o, l, r, tensor3, tensor3); + + REQUIRE(matrix.get_elem({0, 0}) == matrix_value_type(44.0)); + REQUIRE(matrix.get_elem({0, 1}) == matrix_value_type(100.0)); + REQUIRE(matrix.get_elem({1, 0}) == matrix_value_type(64.0)); + REQUIRE(matrix.get_elem({1, 1}) == matrix_value_type(152.0)); + } + + // SECTION("ijk,ljm->iklm") { + + SECTION("ijk,ljm->iklm") { + label_type o("i,k,l,m"); + label_type l("i,j,k"); + label_type r("l,j,m"); + tensor4.contraction_assignment(o, l, r, tensor3, tensor3); + + REQUIRE(tensor4.get_elem({0, 0, 0, 0}) == tensor4_value_type(10.0)); + REQUIRE(tensor4.get_elem({0, 0, 0, 1}) == tensor4_value_type(14.0)); + REQUIRE(tensor4.get_elem({0, 0, 1, 0}) == tensor4_value_type(26.0)); + REQUIRE(tensor4.get_elem({0, 0, 1, 1}) == tensor4_value_type(30.0)); + REQUIRE(tensor4.get_elem({0, 1, 0, 0}) == tensor4_value_type(14.0)); + REQUIRE(tensor4.get_elem({0, 1, 0, 1}) == tensor4_value_type(20.0)); + REQUIRE(tensor4.get_elem({0, 1, 1, 0}) == tensor4_value_type(38.0)); + REQUIRE(tensor4.get_elem({0, 1, 1, 1}) == tensor4_value_type(44.0)); + REQUIRE(tensor4.get_elem({1, 0, 0, 0}) == tensor4_value_type(26.0)); + REQUIRE(tensor4.get_elem({1, 0, 0, 1}) == tensor4_value_type(38.0)); + REQUIRE(tensor4.get_elem({1, 0, 1, 0}) == tensor4_value_type(74.0)); + REQUIRE(tensor4.get_elem({1, 0, 1, 1}) == tensor4_value_type(86.0)); + REQUIRE(tensor4.get_elem({1, 1, 0, 0}) == tensor4_value_type(30.0)); + REQUIRE(tensor4.get_elem({1, 1, 0, 1}) == tensor4_value_type(44.0)); + REQUIRE(tensor4.get_elem({1, 1, 1, 0}) == tensor4_value_type(86.0)); + REQUIRE(tensor4.get_elem({1, 1, 1, 1}) == tensor4_value_type(100.0)); + } + + // SECTION("ij,jkl->ikl") { + + SECTION("ij,jkl->ikl") { + label_type o("i,k,l"); + label_type l("i,j"); + label_type r("j,k,l"); + tensor3.contraction_assignment(o, l, r, matrix, tensor3); + + REQUIRE(tensor3.get_elem({0, 0, 0}) == tensor3_value_type(11.0)); + REQUIRE(tensor3.get_elem({0, 0, 1}) == tensor3_value_type(14.0)); + REQUIRE(tensor3.get_elem({0, 1, 0}) == tensor3_value_type(17.0)); + REQUIRE(tensor3.get_elem({0, 1, 1}) == tensor3_value_type(20.0)); + REQUIRE(tensor3.get_elem({1, 0, 0}) == tensor3_value_type(23.0)); + REQUIRE(tensor3.get_elem({1, 0, 1}) == tensor3_value_type(30.0)); + REQUIRE(tensor3.get_elem({1, 1, 0}) == tensor3_value_type(37.0)); + REQUIRE(tensor3.get_elem({1, 1, 1}) == tensor3_value_type(44.0)); + } +} +} // namespace tensorwrapper::testing diff --git a/tests/cxx/unit_tests/tensorwrapper/backends/testing/elementwise_op.hpp b/tests/cxx/unit_tests/tensorwrapper/backends/testing/elementwise_op.hpp new file mode 100644 index 00000000..b19ee125 --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/backends/testing/elementwise_op.hpp @@ -0,0 +1,315 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include +#include + +namespace tensorwrapper::testing { + +template +void scalar_binary_assignment(Fxn1&& the_op, Fxn2&& corr_op) { + using value_type = typename TestType::value_type; + using shape_type = typename TestType::shape_type; + using label_type = typename TestType::label_type; + + std::vector result_data(1, value_type{0}); + std::span result_span(result_data.data(), result_data.size()); + + std::vector s0_data(1, value_type{3}); + std::span s0_span(s0_data.data(), s0_data.size()); + + std::vector s1_data(1, value_type{5}); + std::span s1_span(s1_data.data(), s1_data.size()); + + TestType result(result_span, shape_type({})); + TestType s0(s0_span, shape_type({})); + TestType s1(s1_span, shape_type({})); + + label_type out(""); + label_type lhs(""); + label_type rhs(""); + the_op(out, lhs, rhs, result, s0, s1); + REQUIRE(result.get_elem({}) == corr_op(s0_data[0], s1_data[0])); +} + +template +void vector_binary_assignment(Fxn1&& the_op, Fxn2&& corr_op) { + using value_type = typename TestType::value_type; + using shape_type = typename TestType::shape_type; + using label_type = typename TestType::label_type; + + std::vector result_data(4, value_type{0}); + std::span result_span(result_data.data(), result_data.size()); + + std::vector s0_data{value_type{1}, value_type{2}, value_type{3}, + value_type{4}}; + std::span s0_span(s0_data.data(), s0_data.size()); + + std::vector s1_data{value_type{5}, value_type{6}, value_type{7}, + value_type{8}}; + std::span s1_span(s1_data.data(), s1_data.size()); + + TestType result(result_span, shape_type({4})); + TestType s0(s0_span, shape_type({4})); + TestType s1(s1_span, shape_type({4})); + + label_type out("i"); + label_type lhs("i"); + label_type rhs("i"); + + the_op(out, lhs, rhs, result, s0, s1); + for(std::size_t i = 0; i < 4; ++i) { + REQUIRE(result.get_elem({i}) == corr_op(s0_data[i], s1_data[i])); + } +} + +template +void matrix_binary_assignment(Fxn1&& the_op, Fxn2&& corr_op) { + using value_type = typename TestType::value_type; + using shape_type = typename TestType::shape_type; + using label_type = typename TestType::label_type; + + std::vector result_data(16, value_type{0}); + std::span result_span(result_data.data(), result_data.size()); + + std::vector s0_data(16); + std::vector s1_data(16); + for(std::size_t i = 0; i < s0_data.size(); ++i) { + s0_data[i] = static_cast(i); + s1_data[i] = static_cast(i * 2); + } + + std::span s0_span(s0_data.data(), s0_data.size()); + std::span s1_span(s1_data.data(), s1_data.size()); + + TestType result(result_span, shape_type({4, 4})); + TestType s0(s0_span, shape_type({4, 4})); + TestType s1(s1_span, shape_type({4, 4})); + + label_type ij("i,j"); + label_type ji("j,i"); + + SECTION("No permutation") { + the_op(ij, ij, ij, result, s0, s1); + for(std::size_t i = 0; i < 4; ++i) { + for(std::size_t j = 0; j < 4; ++j) { + std::size_t idx = i * 4 + j; + auto corr = corr_op(s0_data[idx], s1_data[idx]); + REQUIRE(result.get_elem({i, j}) == corr); + } + } + } + + SECTION("Permute lhs") { + the_op(ij, ji, ij, result, s0, s1); + for(std::size_t i = 0; i < 4; ++i) { + for(std::size_t j = 0; j < 4; ++j) { + std::size_t lhs_idx = j * 4 + i; + std::size_t rhs_idx = i * 4 + j; + auto corr = corr_op(s0_data[lhs_idx], s1_data[rhs_idx]); + REQUIRE(result.get_elem({i, j}) == corr); + } + } + } + + SECTION("Permute rhs") { + the_op(ij, ij, ji, result, s0, s1); + for(std::size_t i = 0; i < 4; ++i) { + for(std::size_t j = 0; j < 4; ++j) { + std::size_t lhs_idx = i * 4 + j; + std::size_t rhs_idx = j * 4 + i; + auto corr = corr_op(s0_data[lhs_idx], s1_data[rhs_idx]); + REQUIRE(result.get_elem({i, j}) == corr); + } + } + } + + SECTION("Permute result") { + the_op(ji, ij, ij, result, s0, s1); + for(std::size_t i = 0; i < 4; ++i) { + for(std::size_t j = 0; j < 4; ++j) { + std::size_t lhs_idx = i * 4 + j; + std::size_t rhs_idx = i * 4 + j; + auto corr = corr_op(s0_data[lhs_idx], s1_data[rhs_idx]); + REQUIRE(result.get_elem({j, i}) == corr); + } + } + } +} + +template +void tensor3_binary_assignment(Fxn1&& the_op, Fxn2&& corr_op) { + using value_type = typename TestType::value_type; + using shape_type = typename TestType::shape_type; + using label_type = typename TestType::label_type; + + const auto n_elements = 8; + std::vector result_data(n_elements, value_type{0}); + std::span result_span(result_data.data(), result_data.size()); + + std::vector t0_data(n_elements); + std::vector t1_data(n_elements); + for(std::size_t i = 0; i < n_elements; ++i) { + t0_data[i] = static_cast(i); + t1_data[i] = static_cast(i * 2); + } + + std::span t0_span(t0_data.data(), t0_data.size()); + std::span t1_span(t1_data.data(), t1_data.size()); + + using rank3_index = std::array; + std::vector tensor3_indices; + for(std::size_t i = 0; i < 2; ++i) { + for(std::size_t j = 0; j < 2; ++j) { + for(std::size_t k = 0; k < 2; ++k) + tensor3_indices.push_back(rank3_index{i, j, k}); + } + } + + TestType result(result_span, shape_type({2, 2, 2})); + TestType t0(t0_span, shape_type({2, 2, 2})); + TestType t1(t1_span, shape_type({2, 2, 2})); + + label_type ijk("i,j,k"); + label_type jik("j,i,k"); + + SECTION("No permutation") { + the_op(ijk, ijk, ijk, result, t0, t1); + for(auto [i, j, k] : tensor3_indices) { + std::size_t lhs_idx = i * 4 + j * 2 + k; + std::size_t rhs_idx = i * 4 + j * 2 + k; + auto corr = corr_op(t0_data[lhs_idx], t1_data[rhs_idx]); + REQUIRE(result.get_elem({i, j, k}) == corr); + } + } + + SECTION("Permute lhs") { + the_op(ijk, jik, ijk, result, t0, t1); + for(auto [i, j, k] : tensor3_indices) { + std::size_t lhs_idx = j * 4 + i * 2 + k; + std::size_t rhs_idx = i * 4 + j * 2 + k; + auto corr = corr_op(t0_data[lhs_idx], t1_data[rhs_idx]); + REQUIRE(result.get_elem({i, j, k}) == corr); + } + } + + SECTION("Permute rhs") { + the_op(ijk, ijk, jik, result, t0, t1); + for(auto [i, j, k] : tensor3_indices) { + std::size_t lhs_idx = i * 4 + j * 2 + k; + std::size_t rhs_idx = j * 4 + i * 2 + k; + auto corr = corr_op(t0_data[lhs_idx], t1_data[rhs_idx]); + REQUIRE(result.get_elem({i, j, k}) == corr); + } + } + + SECTION("Permute result") { + the_op(jik, ijk, ijk, result, t0, t1); + for(auto [i, j, k] : tensor3_indices) { + std::size_t lhs_idx = i * 4 + j * 2 + k; + std::size_t rhs_idx = i * 4 + j * 2 + k; + auto corr = corr_op(t0_data[lhs_idx], t1_data[rhs_idx]); + REQUIRE(result.get_elem({j, i, k}) == corr); + } + } +} + +template +void tensor4_binary_assignment(Fxn1&& the_op, Fxn2&& corr_op) { + using value_type = typename TestType::value_type; + using shape_type = typename TestType::shape_type; + using label_type = typename TestType::label_type; + + const auto n_elements = 16; + std::vector result_data(n_elements, value_type{0}); + std::span result_span(result_data.data(), result_data.size()); + + std::vector t0_data(n_elements); + std::vector t1_data(n_elements); + for(std::size_t i = 0; i < n_elements; ++i) { + t0_data[i] = static_cast(i); + t1_data[i] = static_cast(i * 2); + } + + std::span t0_span(t0_data.data(), t0_data.size()); + std::span t1_span(t1_data.data(), t1_data.size()); + + using rank4_index = std::array; + std::vector tensor4_indices; + for(std::size_t i = 0; i < 2; ++i) { + for(std::size_t j = 0; j < 2; ++j) { + for(std::size_t k = 0; k < 2; ++k) { + for(std::size_t l = 0; l < 2; ++l) { + tensor4_indices.emplace_back(rank4_index{i, j, k, l}); + } + } + } + } + + TestType result(result_span, shape_type({2, 2, 2, 2})); + TestType t0(t0_span, shape_type({2, 2, 2, 2})); + TestType t1(t1_span, shape_type({2, 2, 2, 2})); + + label_type ijkl("i,j,k,l"); + label_type jilk("j,i,l,k"); + + const auto stride0 = 8; + const auto stride1 = 4; + const auto stride2 = 2; + + SECTION("No permutation") { + the_op(ijkl, ijkl, ijkl, result, t0, t1); + for(auto [i, j, k, l] : tensor4_indices) { + std::size_t lhs_idx = i * stride0 + j * stride1 + k * stride2 + l; + std::size_t rhs_idx = i * stride0 + j * stride1 + k * stride2 + l; + auto corr = corr_op(t0_data[lhs_idx], t1_data[rhs_idx]); + REQUIRE(result.get_elem({i, j, k, l}) == corr); + } + } + + SECTION("Permute lhs") { + the_op(ijkl, jilk, ijkl, result, t0, t1); + for(auto [i, j, k, l] : tensor4_indices) { + std::size_t lhs_idx = j * stride0 + i * stride1 + l * stride2 + k; + std::size_t rhs_idx = i * stride0 + j * stride1 + k * stride2 + l; + auto corr = corr_op(t0_data[lhs_idx], t1_data[rhs_idx]); + REQUIRE(result.get_elem({i, j, k, l}) == corr); + } + } + + SECTION("Permute rhs") { + the_op(ijkl, ijkl, jilk, result, t0, t1); + for(auto [i, j, k, l] : tensor4_indices) { + std::size_t lhs_idx = i * stride0 + j * stride1 + k * stride2 + l; + std::size_t rhs_idx = j * stride0 + i * stride1 + l * stride2 + k; + auto corr = corr_op(t0_data[lhs_idx], t1_data[rhs_idx]); + REQUIRE(result.get_elem({i, j, k, l}) == corr); + } + } + + SECTION("Permute result") { + the_op(jilk, ijkl, ijkl, result, t0, t1); + for(auto [i, j, k, l] : tensor4_indices) { + std::size_t lhs_idx = i * stride0 + j * stride1 + k * stride2 + l; + std::size_t rhs_idx = i * stride0 + j * stride1 + k * stride2 + l; + auto corr = corr_op(t0_data[lhs_idx], t1_data[rhs_idx]); + REQUIRE(result.get_elem({j, i, l, k}) == corr); + } + } +} + +} // namespace tensorwrapper::testing diff --git a/tests/cxx/unit_tests/tensorwrapper/backends/testing/hadamard_assignment.hpp b/tests/cxx/unit_tests/tensorwrapper/backends/testing/hadamard_assignment.hpp new file mode 100644 index 00000000..b8e44285 --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/backends/testing/hadamard_assignment.hpp @@ -0,0 +1,72 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include "elementwise_op.hpp" + +namespace tensorwrapper::testing { + +template +void scalar_hadamard_assignment() { + auto the_op = [](auto&& out_idx, auto&& lhs_idx, auto&& rhs_idx, + auto&& result, auto&& t0, auto&& t1) { + result.hadamard_assignment(out_idx, lhs_idx, rhs_idx, t0, t1); + }; + auto corr_op = [](auto a, auto b) { return a * b; }; + scalar_binary_assignment(the_op, corr_op); +} + +template +void vector_hadamard_assignment() { + auto the_op = [](auto&& out_idx, auto&& lhs_idx, auto&& rhs_idx, + auto&& result, auto&& t0, auto&& t1) { + result.hadamard_assignment(out_idx, lhs_idx, rhs_idx, t0, t1); + }; + auto corr_op = [](auto a, auto b) { return a * b; }; + vector_binary_assignment(the_op, corr_op); +} + +template +void matrix_hadamard_assignment() { + auto the_op = [](auto&& out_idx, auto&& lhs_idx, auto&& rhs_idx, + auto&& result, auto&& t0, auto&& t1) { + result.hadamard_assignment(out_idx, lhs_idx, rhs_idx, t0, t1); + }; + auto corr_op = [](auto a, auto b) { return a * b; }; + matrix_binary_assignment(the_op, corr_op); +} + +template +void tensor3_hadamard_assignment() { + auto the_op = [](auto&& out_idx, auto&& lhs_idx, auto&& rhs_idx, + auto&& result, auto&& t0, auto&& t1) { + result.hadamard_assignment(out_idx, lhs_idx, rhs_idx, t0, t1); + }; + auto corr_op = [](auto a, auto b) { return a * b; }; + tensor3_binary_assignment(the_op, corr_op); +} + +template +void tensor4_hadamard_assignment() { + auto the_op = [](auto&& out_idx, auto&& lhs_idx, auto&& rhs_idx, + auto&& result, auto&& t0, auto&& t1) { + result.hadamard_assignment(out_idx, lhs_idx, rhs_idx, t0, t1); + }; + auto corr_op = [](auto a, auto b) { return a * b; }; + tensor4_binary_assignment(the_op, corr_op); +} + +} // namespace tensorwrapper::testing diff --git a/tests/cxx/unit_tests/tensorwrapper/backends/testing/permute_assignment.hpp b/tests/cxx/unit_tests/tensorwrapper/backends/testing/permute_assignment.hpp new file mode 100644 index 00000000..3d2c4bf5 --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/backends/testing/permute_assignment.hpp @@ -0,0 +1,67 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include "unary_op.hpp" + +namespace tensorwrapper::testing { + +template +void scalar_permute_assignment() { + auto the_op = [](auto&& out_idx, auto&& rhs_idx, auto&& result, auto&& t0) { + result.permute_assignment(out_idx, rhs_idx, t0); + }; + auto corr_op = [](auto a) { return a; }; + scalar_unary_assignment(the_op, corr_op); +} + +template +void vector_permute_assignment() { + auto the_op = [](auto&& out_idx, auto&& rhs_idx, auto&& result, auto&& t0) { + result.permute_assignment(out_idx, rhs_idx, t0); + }; + auto corr_op = [](auto a) { return a; }; + vector_unary_assignment(the_op, corr_op); +} + +template +void matrix_permute_assignment() { + auto the_op = [](auto&& out_idx, auto&& rhs_idx, auto&& result, auto&& t0) { + result.permute_assignment(out_idx, rhs_idx, t0); + }; + auto corr_op = [](auto a) { return a; }; + matrix_unary_assignment(the_op, corr_op); +} + +template +void tensor3_permute_assignment() { + auto the_op = [](auto&& out_idx, auto&& rhs_idx, auto&& result, auto&& t0) { + result.permute_assignment(out_idx, rhs_idx, t0); + }; + auto corr_op = [](auto a) { return a; }; + tensor3_unary_assignment(the_op, corr_op); +} + +template +void tensor4_permute_assignment() { + auto the_op = [](auto&& out_idx, auto&& rhs_idx, auto&& result, auto&& t0) { + result.permute_assignment(out_idx, rhs_idx, t0); + }; + auto corr_op = [](auto a) { return a; }; + tensor4_unary_assignment(the_op, corr_op); +} + +} // namespace tensorwrapper::testing diff --git a/tests/cxx/unit_tests/tensorwrapper/backends/testing/scalar_multiplication.hpp b/tests/cxx/unit_tests/tensorwrapper/backends/testing/scalar_multiplication.hpp new file mode 100644 index 00000000..e700fc1c --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/backends/testing/scalar_multiplication.hpp @@ -0,0 +1,83 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include "unary_op.hpp" + +namespace tensorwrapper::testing { + +template +void scalar_scalar_multiplication() { + using value_type = typename TestType::value_type; + value_type scalar{42.0}; + auto the_op = [=](auto&& out_idx, auto&& rhs_idx, auto&& result, + auto&& t0) { + result.scalar_multiplication(out_idx, rhs_idx, scalar, t0); + }; + + auto corr_op = [=](auto a) { return a * scalar; }; + scalar_unary_assignment(the_op, corr_op); +} + +template +void vector_scalar_multiplication() { + using value_type = typename TestType::value_type; + value_type scalar{42.0}; + auto the_op = [=](auto&& out_idx, auto&& rhs_idx, auto&& result, + auto&& t0) { + result.scalar_multiplication(out_idx, rhs_idx, scalar, t0); + }; + auto corr_op = [=](auto a) { return scalar * a; }; + vector_unary_assignment(the_op, corr_op); +} + +template +void matrix_scalar_multiplication() { + using value_type = typename TestType::value_type; + value_type scalar{42.0}; + auto the_op = [=](auto&& out_idx, auto&& rhs_idx, auto&& result, + auto&& t0) { + result.scalar_multiplication(out_idx, rhs_idx, scalar, t0); + }; + auto corr_op = [=](auto a) { return a * scalar; }; + matrix_unary_assignment(the_op, corr_op); +} + +template +void tensor3_scalar_multiplication() { + using value_type = typename TestType::value_type; + value_type scalar{42.0}; + auto the_op = [=](auto&& out_idx, auto&& rhs_idx, auto&& result, + auto&& t0) { + result.scalar_multiplication(out_idx, rhs_idx, scalar, t0); + }; + auto corr_op = [=](auto a) { return scalar * a; }; + tensor3_unary_assignment(the_op, corr_op); +} + +template +void tensor4_scalar_multiplication() { + using value_type = typename TestType::value_type; + value_type scalar{42.0}; + auto the_op = [=](auto&& out_idx, auto&& rhs_idx, auto&& result, + auto&& t0) { + result.scalar_multiplication(out_idx, rhs_idx, scalar, t0); + }; + auto corr_op = [=](auto a) { return a * scalar; }; + tensor4_unary_assignment(the_op, corr_op); +} + +} // namespace tensorwrapper::testing diff --git a/tests/cxx/unit_tests/tensorwrapper/backends/testing/subtraction_assignment.hpp b/tests/cxx/unit_tests/tensorwrapper/backends/testing/subtraction_assignment.hpp new file mode 100644 index 00000000..baf54910 --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/backends/testing/subtraction_assignment.hpp @@ -0,0 +1,72 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include "elementwise_op.hpp" + +namespace tensorwrapper::testing { + +template +void scalar_subtraction_assignment() { + auto the_op = [](auto&& out_idx, auto&& lhs_idx, auto&& rhs_idx, + auto&& result, auto&& t0, auto&& t1) { + result.subtraction_assignment(out_idx, lhs_idx, rhs_idx, t0, t1); + }; + auto corr_op = [](auto a, auto b) { return a - b; }; + scalar_binary_assignment(the_op, corr_op); +} + +template +void vector_subtraction_assignment() { + auto the_op = [](auto&& out_idx, auto&& lhs_idx, auto&& rhs_idx, + auto&& result, auto&& t0, auto&& t1) { + result.subtraction_assignment(out_idx, lhs_idx, rhs_idx, t0, t1); + }; + auto corr_op = [](auto a, auto b) { return a - b; }; + vector_binary_assignment(the_op, corr_op); +} + +template +void matrix_subtraction_assignment() { + auto the_op = [](auto&& out_idx, auto&& lhs_idx, auto&& rhs_idx, + auto&& result, auto&& t0, auto&& t1) { + result.subtraction_assignment(out_idx, lhs_idx, rhs_idx, t0, t1); + }; + auto corr_op = [](auto a, auto b) { return a - b; }; + matrix_binary_assignment(the_op, corr_op); +} + +template +void tensor3_subtraction_assignment() { + auto the_op = [](auto&& out_idx, auto&& lhs_idx, auto&& rhs_idx, + auto&& result, auto&& t0, auto&& t1) { + result.subtraction_assignment(out_idx, lhs_idx, rhs_idx, t0, t1); + }; + auto corr_op = [](auto a, auto b) { return a - b; }; + tensor3_binary_assignment(the_op, corr_op); +} + +template +void tensor4_subtraction_assignment() { + auto the_op = [](auto&& out_idx, auto&& lhs_idx, auto&& rhs_idx, + auto&& result, auto&& t0, auto&& t1) { + result.subtraction_assignment(out_idx, lhs_idx, rhs_idx, t0, t1); + }; + auto corr_op = [](auto a, auto b) { return a - b; }; + tensor4_binary_assignment(the_op, corr_op); +} + +} // namespace tensorwrapper::testing diff --git a/tests/cxx/unit_tests/tensorwrapper/backends/testing/unary_op.hpp b/tests/cxx/unit_tests/tensorwrapper/backends/testing/unary_op.hpp new file mode 100644 index 00000000..6c8f8a05 --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/backends/testing/unary_op.hpp @@ -0,0 +1,214 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +#pragma once +#include +#include + +namespace tensorwrapper::testing { + +template +void scalar_unary_assignment(Fxn1&& the_op, Fxn2&& corr_op) { + using value_type = typename TestType::value_type; + using shape_type = typename TestType::shape_type; + using label_type = typename TestType::label_type; + + std::vector result_data(1, value_type{0}); + std::span result_span(result_data.data(), result_data.size()); + + std::vector s0_data(1, value_type{3}); + std::span s0_span(s0_data.data(), s0_data.size()); + + TestType result(result_span, shape_type({})); + TestType s0(s0_span, shape_type({})); + + label_type out(""); + label_type rhs(""); + the_op(out, rhs, result, s0); + REQUIRE(result.get_elem({}) == corr_op(s0_data[0])); +} + +template +void vector_unary_assignment(Fxn1&& the_op, Fxn2&& corr_op) { + using value_type = typename TestType::value_type; + using shape_type = typename TestType::shape_type; + using label_type = typename TestType::label_type; + + const auto n_elements = 4; + std::vector result_data(n_elements, value_type{0}); + std::span result_span(result_data.data(), result_data.size()); + + std::vector s0_data(n_elements, value_type{3}); + std::span s0_span(s0_data.data(), s0_data.size()); + + TestType result(result_span, shape_type({4})); + TestType s0(s0_span, shape_type({4})); + + label_type out("i"); + label_type rhs("i"); + the_op(out, rhs, result, s0); + for(std::size_t i = 0; i < n_elements; ++i) + REQUIRE(result.get_elem({i}) == corr_op(s0_data[i])); +} + +template +void matrix_unary_assignment(Fxn1&& the_op, Fxn2&& corr_op) { + using value_type = typename TestType::value_type; + using shape_type = typename TestType::shape_type; + using label_type = typename TestType::label_type; + + const auto n_elements = 16; + std::vector result_data(n_elements, value_type{0}); + std::span result_span(result_data.data(), result_data.size()); + + std::vector s0_data(n_elements, value_type{0}); + for(std::size_t i = 0; i < n_elements; ++i) + s0_data[i] = static_cast(i); + + std::span s0_span(s0_data.data(), s0_data.size()); + + TestType result(result_span, shape_type({4, 4})); + TestType s0(s0_span, shape_type({4, 4})); + + label_type ij("i,j"); + label_type ji("j,i"); + + SECTION("No permutation") { + the_op(ij, ij, result, s0); + for(std::size_t i = 0; i < 4; ++i) + for(std::size_t j = 0; j < 4; ++j) { + std::size_t idx = i * 4 + j; + REQUIRE(result.get_elem({i, j}) == corr_op(s0_data[idx])); + } + } + + SECTION("Permute rhs") { + the_op(ij, ji, result, s0); + for(std::size_t i = 0; i < 4; ++i) + for(std::size_t j = 0; j < 4; ++j) { + std::size_t idx = j * 4 + i; + REQUIRE(result.get_elem({i, j}) == corr_op(s0_data[idx])); + } + } + + SECTION("Permute result") { + the_op(ji, ij, result, s0); + for(std::size_t i = 0; i < 4; ++i) + for(std::size_t j = 0; j < 4; ++j) { + std::size_t idx = i * 4 + j; + REQUIRE(result.get_elem({j, i}) == corr_op(s0_data[idx])); + } + } +} + +template +void tensor3_unary_assignment(Fxn1&& the_op, Fxn2&& corr_op) { + using value_type = typename TestType::value_type; + using shape_type = typename TestType::shape_type; + using label_type = typename TestType::label_type; + + const auto n_elements = 8; + std::vector result_data(n_elements, value_type{0}); + std::span result_span(result_data.data(), result_data.size()); + + std::vector s0_data(n_elements, value_type{0}); + for(std::size_t i = 0; i < n_elements; ++i) + s0_data[i] = static_cast(i); + + std::span s0_span(s0_data.data(), s0_data.size()); + + TestType result(result_span, shape_type({2, 2, 2})); + TestType s0(s0_span, shape_type({2, 2, 2})); + + label_type ijk("i,j,k"); + label_type jik("j,i,k"); + + using rank3_index = std::array; + std::vector tensor3_indices; + for(std::size_t i = 0; i < 2; ++i) { + for(std::size_t j = 0; j < 2; ++j) { + for(std::size_t k = 0; k < 2; ++k) + tensor3_indices.push_back(rank3_index{i, j, k}); + } + } + + SECTION("No permutation") { + the_op(ijk, ijk, result, s0); + for(const auto [i, j, k] : tensor3_indices) { + std::size_t idx = i * 4 + j * 2 + k; + REQUIRE(result.get_elem({i, j, k}) == corr_op(s0_data[idx])); + } + } + + SECTION("Permute rhs") { + the_op(ijk, jik, result, s0); + for(const auto [i, j, k] : tensor3_indices) { + std::size_t idx = j * 4 + i * 2 + k; + REQUIRE(result.get_elem({i, j, k}) == corr_op(s0_data[idx])); + } + } +} + +template +void tensor4_unary_assignment(Fxn1&& the_op, Fxn2&& corr_op) { + using value_type = typename TestType::value_type; + using shape_type = typename TestType::shape_type; + using label_type = typename TestType::label_type; + + const auto n_elements = 16; + std::vector result_data(n_elements, value_type{0}); + std::span result_span(result_data.data(), result_data.size()); + + std::vector s0_data(n_elements, value_type{0}); + for(std::size_t i = 0; i < n_elements; ++i) + s0_data[i] = static_cast(i); + + std::span s0_span(s0_data.data(), s0_data.size()); + + TestType result(result_span, shape_type({2, 2, 2, 2})); + TestType s0(s0_span, shape_type({2, 2, 2, 2})); + + label_type ijkl("i,j,k,l"); + label_type jikl("j,i,k,l"); + + using rank4_index = std::array; + std::vector tensor4_indices; + for(std::size_t i = 0; i < 2; ++i) { + for(std::size_t j = 0; j < 2; ++j) { + for(std::size_t k = 0; k < 2; ++k) + for(std::size_t l = 0; l < 2; ++l) + tensor4_indices.push_back(rank4_index{i, j, k, l}); + } + } + + SECTION("No permutation") { + the_op(ijkl, ijkl, result, s0); + for(const auto [i, j, k, l] : tensor4_indices) { + std::size_t idx = i * 8 + j * 4 + k * 2 + l; + REQUIRE(result.get_elem({i, j, k, l}) == corr_op(s0_data[idx])); + } + } + + SECTION("Permute rhs") { + the_op(ijkl, jikl, result, s0); + for(const auto [i, j, k, l] : tensor4_indices) { + std::size_t idx = j * 8 + i * 4 + k * 2 + l; + REQUIRE(result.get_elem({i, j, k, l}) == corr_op(s0_data[idx])); + } + } +} + +} // namespace tensorwrapper::testing diff --git a/tests/cxx/unit_tests/tensorwrapper/buffer/detail_/addition_visitor.cpp b/tests/cxx/unit_tests/tensorwrapper/buffer/detail_/addition_visitor.cpp new file mode 100644 index 00000000..d7b46618 --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/buffer/detail_/addition_visitor.cpp @@ -0,0 +1,38 @@ +/* + * Copyright 2025 NWChemEx-Project + * + * 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 + * + * 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. + */ + +// #include +// #include + +// using namespace tensorwrapper; + +// TEMPLATE_LIST_TEST_CASE("AdditionVisitor", "[buffer][detail_]", +// types::floating_point_types) { +// using VisitorType = buffer::detail_::AdditionVisitor; + +// VisitorType visitor; + +// SECTION("vectors") { +// std::vector lhs{1.0, 2.0, 3.0}; +// std::vector rhs{4.0, 5.0, 6.0}; + +// visitor(std::span(lhs), std::span(rhs)); + +// REQUIRE(lhs[0] == Approx(5.0).epsilon(1e-10)); +// REQUIRE(lhs[1] == Approx(7.0).epsilon(1e-10)); +// REQUIRE(lhs[2] == Approx(9.0).epsilon(1e-10)); +// } +// } diff --git a/tests/cxx/unit_tests/tensorwrapper/symmetry/permutation.cpp b/tests/cxx/unit_tests/tensorwrapper/symmetry/permutation.cpp index a99489e4..e28d4ce2 100644 --- a/tests/cxx/unit_tests/tensorwrapper/symmetry/permutation.cpp +++ b/tests/cxx/unit_tests/tensorwrapper/symmetry/permutation.cpp @@ -58,6 +58,12 @@ TEST_CASE("Permutation") { REQUIRE(p5.size() == mode_index_type(0)); REQUIRE(p5.rank() == mode_index_type(5)); + // One cycle and a fix-point via one-line + Permutation one_line({1, 0, 2}); + REQUIRE(one_line.size() == mode_index_type(1)); + REQUIRE(one_line.rank() == mode_index_type(3)); + REQUIRE(one_line.at(0) == c01); + // Two cycles via one-line Permutation p01_23{1, 0, 3, 2}; REQUIRE(p01_23.size() == mode_index_type(2)); @@ -139,6 +145,22 @@ TEST_CASE("Permutation") { REQUIRE(two_cycles.size() == 2); } + SECTION("apply") { + REQUIRE(defaulted.apply(cycle_type{}) == cycle_type{}); + + // One cycle (0 1) + REQUIRE(one_cycle.apply(cycle_type{0, 1}) == cycle_type{1, 0}); + REQUIRE(one_cycle.apply(cycle_type{1, 2}) == cycle_type{2, 1}); + + // Two cycles (1 3 2)(4 5) + REQUIRE(two_cycles.apply(cycle_type{0, 1, 2, 3, 4, 5}) == + cycle_type{0, 2, 3, 1, 5, 4}); + + Permutation one_line({1, 0, 2}); + REQUIRE(one_line.apply(cycle_type{0, 1, 2}) == cycle_type{1, 0, 2}); + REQUIRE(one_line.apply(cycle_type{2, 3, 4}) == cycle_type{3, 2, 4}); + } + SECTION("swap") { Permutation copy_defaulted(defaulted); Permutation copy_one_cycle(one_cycle);