diff --git a/cmake/check_user_env.cmake b/cmake/check_user_env.cmake index 52c987aec..24fe5b9bd 100644 --- a/cmake/check_user_env.cmake +++ b/cmake/check_user_env.cmake @@ -45,25 +45,15 @@ if(DEFINED GA_RUNTIME) endif() endif() -check_compiler_version(C Clang 9) -check_compiler_version(CXX Clang 9) +check_compiler_version(C Clang 18) +check_compiler_version(CXX Clang 18) -check_compiler_version(C AppleClang 15) -check_compiler_version(CXX AppleClang 15) +check_compiler_version(C AppleClang 19) +check_compiler_version(CXX AppleClang 19) -check_compiler_version(C GNU 9.1) -check_compiler_version(CXX GNU 9.1) -check_compiler_version(Fortran GNU 9.1) - -#TODO:Check for GCC>=9 compatibility -# check_compiler_version(C Intel 19) -# check_compiler_version(CXX Intel 19) -# check_compiler_version(Fortran Intel 19) - -#TODO:Check for GCC>=9 compatibility -check_compiler_version(C PGI 20) -check_compiler_version(CXX PGI 20) -check_compiler_version(Fortran PGI 20) +check_compiler_version(C GNU 14.1) +check_compiler_version(CXX GNU 14.1) +check_compiler_version(Fortran GNU 14.1) find_package(MPI REQUIRED) @@ -97,7 +87,7 @@ if(${PROJECT_NAME}_ENABLE_CUDA) message(FATAL_ERROR "CUDA Toolkit not found.") endif() - set(_CUDA_MIN "11.7") + set(_CUDA_MIN "12.8") if(CMAKE_CUDA_COMPILER_VERSION VERSION_LESS ${_CUDA_MIN}) message(FATAL_ERROR "CUDA version provided \ (${CMAKE_CUDA_COMPILER_VERSION}) \ diff --git a/src/tamm/CMakeLists.txt b/src/tamm/CMakeLists.txt index 497b00e8d..391cb10bd 100644 --- a/src/tamm/CMakeLists.txt +++ b/src/tamm/CMakeLists.txt @@ -67,6 +67,7 @@ set(TAMM_INCLUDES eigen_includes.hpp runtime_engine.hpp block_buffer.hpp + block_scratch.hpp lru_cache.hpp kernels/assign.hpp kernels/multiply.hpp @@ -92,7 +93,6 @@ set(TAMM_INCLUDES tensor_variant.hpp op_executor.hpp symbol.hpp - op_cost.hpp block_operations.hpp setop.hpp scanop.hpp diff --git a/src/tamm/addop.hpp b/src/tamm/addop.hpp index 82214bba3..8fa9fd33e 100644 --- a/src/tamm/addop.hpp +++ b/src/tamm/addop.hpp @@ -10,6 +10,7 @@ #include "tamm/block_assign_plan.hpp" #include "tamm/block_operations.hpp" +#include "tamm/block_scratch.hpp" #include "tamm/boundvec.hpp" #include "tamm/errors.hpp" #include "tamm/kernels/assign.hpp" @@ -44,134 +45,47 @@ namespace tamm::internal { template struct AddOpPlanBase { using AddOpT = AddOp; + + // writes()/accumulates()/reads() are identical for every plan: an assign + // writes the LHS, an update accumulates into it, and the RHS is always read. + // (The former local/global split was a no-op for scheduling — the variants + // were concatenated and only membership matters for dependency tests — and no + // code queried the local/global variants individually.) TensorBase* writes(const AddOpT& addop) const { - auto ret1 = local_writes(addop); - auto ret2 = global_writes(addop); - ret1.insert(ret1.end(), ret2.begin(), ret2.end()); - return !ret1.empty() ? ret1[0] : nullptr; + return addop.is_assign() ? addop.lhs().base_ptr() : nullptr; } - TensorBase* accumulates(const AddOpT& addop) const { - auto ret1 = local_accumulates(addop); - auto ret2 = global_accumulates(addop); - ret1.insert(ret1.end(), ret2.begin(), ret2.end()); - return !ret1.empty() ? ret1[0] : nullptr; + return addop.is_assign() ? nullptr : addop.lhs().base_ptr(); } - std::vector reads(const AddOpT& addop) const { - auto ret1 = local_reads(addop); - auto ret2 = global_reads(addop); - ret1.insert(ret1.end(), ret2.begin(), ret2.end()); - return !ret1.empty() ? ret1 : std::vector{}; + return {addop.rhs().base_ptr()}; } - virtual std::vector global_writes(const AddOpT& addop) const = 0; - virtual std::vector global_accumulates(const AddOpT& addop) const = 0; - virtual std::vector global_reads(const AddOpT& addop) const = 0; - virtual std::vector local_writes(const AddOpT& addop) const = 0; - virtual std::vector local_accumulates(const AddOpT& addop) const = 0; - virtual std::vector local_reads(const AddOpT& addop) const = 0; - virtual void apply(const AddOpT& addop, ExecutionContext& ec, ExecutionHW hw) = 0; + virtual ~AddOpPlanBase() = default; }; // AddOpPlanBase template struct FlatAddPlan: public AddOpPlanBase { using AddOpT = AddOp; - std::vector global_writes(const AddOpT& addop) const override { return {}; } - std::vector global_accumulates(const AddOpT& addop) const override { return {}; } - - std::vector global_reads(const AddOpT& addop) const override { return {}; } - - std::vector local_writes(const AddOpT& addop) const override { - if(addop.is_assign()) { return {addop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_accumulates(const AddOpT& addop) const override { - if(!addop.is_assign()) { return {addop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_reads(const AddOpT& addop) const override { - return {addop.rhs().base_ptr()}; - } void apply(const AddOpT& addop, ExecutionContext& ec, ExecutionHW hw) override; }; // FlatAddPlan template struct LHSAddPlan: public AddOpPlanBase { using AddOpT = AddOp; - std::vector global_writes(const AddOpT& addop) const override { return {}; } - std::vector global_accumulates(const AddOpT& addop) const override { return {}; } - std::vector global_reads(const AddOpT& addop) const override { return {}; } - - std::vector local_writes(const AddOpT& addop) const override { - if(addop.is_assign()) { return {addop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_accumulates(const AddOpT& addop) const override { - if(!addop.is_assign()) { return {addop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_reads(const AddOpT& addop) const override { - return {addop.rhs().base_ptr()}; - } void apply(const AddOpT& addop, ExecutionContext& ec, ExecutionHW hw) override; }; // LHSAddPlan template struct GeneralFlatAddPlan: public AddOpPlanBase { using AddOpT = AddOp; - std::vector global_writes(const AddOpT& addop) const override { - if(addop.is_assign()) { return {addop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector global_accumulates(const AddOpT& addop) const override { - if(!addop.is_assign()) { return {addop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector global_reads(const AddOpT& addop) const override { - return {addop.rhs().base_ptr()}; - } - - std::vector local_writes(const AddOpT& addop) const override { - if(addop.is_assign()) { return {addop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_accumulates(const AddOpT& addop) const override { - if(!addop.is_assign()) { return {addop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_reads(const AddOpT& addop) const override { - return {addop.rhs().base_ptr()}; - } void apply(const AddOpT& addop, ExecutionContext& ec, ExecutionHW hw) override; }; // GeneralFlatAddPlan template struct GeneralLHSAddPlan: public AddOpPlanBase { using AddOpT = AddOp; - std::vector global_writes(const AddOpT& addop) const override { - if(addop.is_assign()) { return {addop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector global_accumulates(const AddOpT& addop) const override { - if(!addop.is_assign()) { return {addop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector global_reads(const AddOpT& addop) const override { - return {addop.rhs().base_ptr()}; - } - std::vector local_writes(const AddOpT& addop) const override { - if(addop.is_assign()) { return {addop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_accumulates(const AddOpT& addop) const override { - if(!addop.is_assign()) { return {addop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_reads(const AddOpT& addop) const override { - return {addop.rhs().base_ptr()}; - } void apply(const AddOpT& addop, ExecutionContext& ec, ExecutionHW hw) override; }; // GeneralLHSAddPlan @@ -198,8 +112,8 @@ class AddOp: public Op { rhs_lbls = IndexLabelVec(labels.begin() + lhs.labels().size(), labels.begin() + lhs.labels().size() + rhs.labels().size()); - lhs_.set_labels(lhs_lbls); - rhs_.set_labels(rhs_lbls); + lhs_.set_labels(std::move(lhs_lbls)); + rhs_.set_labels(std::move(rhs_lbls)); } if(lhs.has_str_lbl()) { fillin_labels(); } @@ -208,7 +122,7 @@ class AddOp: public Op { validate(); } - AddOp(const AddOp&) = default; + // Copy/move are implicitly generated (Rule of Zero); clone() copies. T alpha() const { return alpha_; } @@ -236,7 +150,7 @@ class AddOp: public Op { } std::shared_ptr clone() const override { - return std::shared_ptr(new AddOp{*this}); + return std::make_shared(*this); } void execute(ExecutionContext& ec, ExecutionHW hw = ExecutionHW::CPU) override { @@ -361,31 +275,8 @@ class AddOp: public Op { tamm_terminate(os.str()); } - IndexLabelVec ilv{lhs_.labels()}; - ilv.insert(ilv.end(), rhs_.labels().begin(), rhs_.labels().end()); - - for(size_t i = 0; i < ilv.size(); i++) { - for(const auto& dl: ilv[i].secondary_labels()) { - size_t j; - for(j = 0; j < ilv.size(); j++) { - if(dl.tiled_index_space() == ilv[j].tiled_index_space() && dl.label() == ilv[j].label()) { - break; - } - } - EXPECTS(j < ilv.size()); - } - } - - for(size_t i = 0; i < ilv.size(); i++) { - const auto& ilbl = ilv[i]; - for(size_t j = i + 1; j < ilv.size(); j++) { - const auto& jlbl = ilv[j]; - if(ilbl.tiled_index_space() == jlbl.tiled_index_space() && ilbl.label() == jlbl.label() && - ilbl.label_str() == jlbl.label_str()) { - EXPECTS(ilbl == jlbl); - } - } - } + const auto ilv = internal::merge_vector(lhs_.labels(), rhs_.labels()); + internal::validate_index_labels(ilv); } void fillin_int_labels() { @@ -540,8 +431,8 @@ void GeneralFlatAddPlan::apply(const AddOpT // EXPECTS(pg_lhs.size() == pg_rhs.size()); // EXPECTS(pg_lhs.size() == pg_ec.size()); - BlockAssignPlan::OpType optype = is_assign ? optype = BlockAssignPlan::OpType::set - : BlockAssignPlan::OpType::update; + BlockAssignPlan::OpType optype = + is_assign ? BlockAssignPlan::OpType::set : BlockAssignPlan::OpType::update; BlockAssignPlan plan{lhs_lt.labels(), rhs_lt.labels(), optype}; std::vector pg_lhs_in_ec = pg_lhs.rank_translate(pg_ec); @@ -550,6 +441,10 @@ void GeneralFlatAddPlan::apply(const AddOpT Proc round_robin_counter = 0; Proc ec_pg_size = Proc{ec.pg().size()}; + // Reused (grow-only) scratch for the remote-copy case; RAII-owned. + internal::BlockScratch lhs_scratch; + internal::BlockScratch rhs_scratch; + for(size_t i = 0; i < pg_lhs_in_ec.size(); i++) { Proc assigned_proc; if(pg_lhs_in_ec[i] >= Proc{0}) { @@ -567,34 +462,24 @@ void GeneralFlatAddPlan::apply(const AddOpT } if(proc_me_in_ec == assigned_proc) { - bool alloced_lhs_buf{false}; - bool alloced_rhs_buf{false}; - T1* lhs_buf{nullptr}; - T2* rhs_buf{nullptr}; + T1* lhs_buf{nullptr}; + T2* rhs_buf{nullptr}; /// get total buffer size for a given Proc size_t lhs_size = lhs_tensor.total_buf_size(i); size_t rhs_size = rhs_tensor.total_buf_size(i); EXPECTS(lhs_size == rhs_size); if(lhs_size <= 0) continue; - if(proc_me_in_ec == pg_lhs_in_ec[i]) { - lhs_buf = lhs_tensor.access_local_buf(); - alloced_lhs_buf = false; - } + if(proc_me_in_ec == pg_lhs_in_ec[i]) { lhs_buf = lhs_scratch.view(lhs_tensor.access_local_buf()); } else { - lhs_buf = new T1[lhs_size]; - alloced_lhs_buf = true; + lhs_buf = lhs_scratch.owned(lhs_size); auto* lhs_mem_region = lhs_tensor.memory_region(); /// get all of lhs's buf at i-th proc to lhs_buf lhs_mem_region->get(Proc{i}, Offset{0}, Size{lhs_size}, lhs_buf); } - if(proc_me_in_ec == pg_rhs_in_ec[i]) { - rhs_buf = rhs_tensor.access_local_buf(); - alloced_rhs_buf = false; - } + if(proc_me_in_ec == pg_rhs_in_ec[i]) { rhs_buf = rhs_scratch.view(rhs_tensor.access_local_buf()); } else { - rhs_buf = new T2[rhs_size]; - alloced_rhs_buf = true; + rhs_buf = rhs_scratch.owned(rhs_size); EXPECTS(rhs_buf != nullptr); auto* rhs_mem_region = rhs_tensor.memory_region(); /// get all of rhs's buf at i-th proc to rhs_buf @@ -612,8 +497,6 @@ void GeneralFlatAddPlan::apply(const AddOpT auto* lhs_mem_region = lhs_tensor.memory_region(); lhs_mem_region->put(Proc{i}, Offset{0}, Size{lhs_size}, lhs_buf); } - if(alloced_lhs_buf) { delete[] lhs_buf; } - if(alloced_rhs_buf) { delete[] rhs_buf; } } } } @@ -656,21 +539,22 @@ void GeneralLHSAddPlan::apply(const AddOpT& LabelLoopNest loop_nest{merged_use_labels}; + // Reused (grow-only) scratch for the non-local / view / dense cases; RAII-owned. + internal::BlockScratch lhs_scratch; + internal::BlockScratch rhs_scratch; + auto lambda = [&](const IndexVector& l_blockid, const IndexVector& r_blockid) { auto [lhs_proc, lhs_offset] = ldist.locate(l_blockid); auto lhs_blocksize = lhs_tensor.block_size(l_blockid); auto lhs_blockdims = lhs_tensor.block_dims(l_blockid); T1* lhs_buf{nullptr}; - bool lhs_alloced{false}; if(proc_lhs_to_ec[lhs_proc.value()] == proc_me_in_ec && lhs_tensor.kind() != TensorBase::TensorKind::view && lhs_tensor.kind() != TensorBase::TensorKind::dense) { - lhs_buf = lhs_tensor.access_local_buf() + lhs_offset.value(); - lhs_alloced = false; + lhs_buf = lhs_scratch.view(lhs_tensor.access_local_buf() + lhs_offset.value()); } else { - lhs_buf = new T1[lhs_blocksize]; - lhs_alloced = true; + lhs_buf = lhs_scratch.owned(lhs_blocksize); span lhs_span{lhs_buf, lhs_blocksize}; lhs_tensor.get(l_blockid, lhs_span); } @@ -679,17 +563,14 @@ void GeneralLHSAddPlan::apply(const AddOpT& auto rhs_blocksize = rhs_tensor.block_size(r_blockid); auto rhs_blockdims = rhs_tensor.block_dims(r_blockid); T2* rhs_buf{nullptr}; - bool rhs_alloced{false}; if(proc_rhs_to_ec[rhs_proc.value()] == proc_me_in_ec && rhs_tensor.kind() != TensorBase::TensorKind::view && rhs_tensor.kind() != TensorBase::TensorKind::lambda && rhs_tensor.kind() != TensorBase::TensorKind::dense) { - rhs_buf = rhs_tensor.access_local_buf() + rhs_offset.value(); - rhs_alloced = false; + rhs_buf = rhs_scratch.view(rhs_tensor.access_local_buf() + rhs_offset.value()); } else { - rhs_buf = new T2[rhs_blocksize]; - rhs_alloced = true; + rhs_buf = rhs_scratch.owned(rhs_blocksize); span rhs_span{rhs_buf, rhs_blocksize}; rhs_tensor.get(r_blockid, rhs_span); } @@ -705,8 +586,6 @@ void GeneralLHSAddPlan::apply(const AddOpT& span lhs_span{lhs_buf, lhs_blocksize}; lhs_tensor.put(l_blockid, lhs_span); } - if(lhs_alloced) { delete[] lhs_buf; } - if(rhs_alloced) { delete[] rhs_buf; } }; Proc round_robin_counter = 0; diff --git a/src/tamm/allocop.hpp b/src/tamm/allocop.hpp index 179f11832..3dfcccd51 100644 --- a/src/tamm/allocop.hpp +++ b/src/tamm/allocop.hpp @@ -18,15 +18,15 @@ class AllocOp: public Op { public: AllocOp(TensorType tensor, ExecutionContext& ec): tensor_{tensor}, ec_{ec} {} - AllocOp(const AllocOp&) = default; - TensorType tensor() const { return tensor_; } OpList canonicalize() const override { return OpList{(*this)}; } OpType op_type() const override { return OpType::alloc; } - std::shared_ptr clone() const override { return std::shared_ptr(new AllocOp{*this}); } + std::shared_ptr clone() const override { + return std::make_shared(*this); + } void execute(ExecutionContext& ec, ExecutionHW hw = ExecutionHW::CPU) override { tensor_.allocate(&ec_); diff --git a/src/tamm/attribute.hpp b/src/tamm/attribute.hpp index e4be515d7..544c2fa8f 100644 --- a/src/tamm/attribute.hpp +++ b/src/tamm/attribute.hpp @@ -102,7 +102,7 @@ class Attribute { * * @returns empty */ - bool empty() const { return attr_map_.find(T{0}) != attr_map_.end(); } + bool empty() const { return attr_map_.contains(T{0}); } protected: AttributeToRangeMap attr_map_; // @@ -116,8 +116,8 @@ using SpatialAttribute = Attribute; namespace std { template<> struct hash { - typedef tamm::SpinAttribute argument_type; - typedef std::size_t result_type; + using argument_type = tamm::SpinAttribute; + using result_type = std::size_t; result_type operator()(argument_type const& attribute) const noexcept { using tamm::internal::hash_combine; @@ -136,8 +136,8 @@ struct hash { template<> struct hash { - typedef tamm::SpatialAttribute argument_type; - typedef std::size_t result_type; + using argument_type = tamm::SpatialAttribute; + using result_type = std::size_t; result_type operator()(argument_type const& attribute) const noexcept { using tamm::internal::hash_combine; diff --git a/src/tamm/block_buffer.hpp b/src/tamm/block_buffer.hpp index bd8cdd1f7..3b127bbda 100644 --- a/src/tamm/block_buffer.hpp +++ b/src/tamm/block_buffer.hpp @@ -3,94 +3,152 @@ #include "tamm/tensor.hpp" #include "tamm/types.hpp" #include +#include +#include namespace tamm { class RuntimeEngine; /** - * @brief The class used to pass block buffers to user functions + * @brief Non-owning view + optional owner over a contiguous block of T. * - * @tparam T + * C++20 rewrite: + * - Ownership managed by std::vector storage_ instead of raw new[]/delete[]. + * - Exposed as std::span buf_span_ (non-owning view into storage_). + * - The 'allocated' bool flag is eliminated entirely. + * - Copy / move / dtor are Rule-of-Zero (compiler-generated, always correct). + * + * @tparam T Element type of the block */ template class BlockBuffer { public: + // ------------------------------------------------------------------- + // Constructors + // ------------------------------------------------------------------- + BlockBuffer() = default; - BlockBuffer(span buf_span, IndexedTensor indexedTensor, RuntimeEngine* re, - bool allocated = false): - buf_span{buf_span}, allocated{allocated}, indexedTensor{indexedTensor}, re{re} {} - BlockBuffer(const BlockBuffer& block_buffer): - indexedTensor{block_buffer.indexedTensor}, re{block_buffer.re} { - if(allocated) delete[] buf_span.data(); - allocated = true; - const auto size = block_buffer.buf_span.size(); - T* buffer = new T[size]; // that will need to be more complicated once we get device buffers - std::copy(block_buffer.buf_span.begin(), block_buffer.buf_span.end(), buffer); - buf_span = span{buffer, size}; + + /// Wrap an externally-owned span (no allocation, not an owner). + BlockBuffer(std::span buf_span, IndexedTensor indexedTensor, + RuntimeEngine* re) + : buf_span_{buf_span}, indexedTensor_{indexedTensor}, re_{re} {} + + /// Allocate a fresh buffer, fetch the block from the tensor. + BlockBuffer(Tensor tensor, IndexVector blockid) + : storage_(tensor.block_size(blockid)), + indexedTensor_{tensor, blockid} { + buf_span_ = std::span{storage_.data(), storage_.size()}; + tensor.get(blockid, buf_span_); } - BlockBuffer(BlockBuffer&& block_buffer) { - buf_span = std::move(block_buffer.buf_span); - allocated = false; - indexedTensor = std::move(block_buffer.indexedTensor); - re = block_buffer.re; - block_buffer.re = nullptr; + + // Rule-of-Zero: copy / move / dtor are all compiler-generated. + // - Copy: deep-copies storage_, buf_span_ reconstructed from it below + // via the copy ctor body (storage_ copy + span retarget). + // - Move: storage_ and buf_span_ are both moved correctly. + // - Dtor: storage_ RAII-cleans up automatically. + // + // We need a custom copy constructor only to retarget buf_span_ after + // the vector copy (span still points at the source's allocation). + BlockBuffer(const BlockBuffer& other) + : storage_{other.storage_}, + indexedTensor_{other.indexedTensor_}, + re_{other.re_} { + // If other was owning (storage_ has data), point our span at our copy. + // If other was non-owning (storage_ empty, span points externally), + // copy the span view as-is. + if (!storage_.empty()) + buf_span_ = std::span{storage_.data(), storage_.size()}; + else + buf_span_ = other.buf_span_; } - BlockBuffer& operator=(const BlockBuffer& block_buffer) { - indexedTensor = block_buffer.indexedTensor; - re = block_buffer.indexedTensor; - if(allocated) delete[] buf_span.data(); - allocated = true; - const auto size = block_buffer.buf_span.size(); - T* buffer = new T[size]; // that will need to be more complicated once we get device buffers - std::copy(block_buffer.buf_span.begin(), block_buffer.buf_span.end(), buffer); - buf_span = span{buffer, size}; + + BlockBuffer& operator=(const BlockBuffer& other) { + if (this == &other) return *this; + storage_ = other.storage_; + indexedTensor_ = other.indexedTensor_; + re_ = other.re_; + if (!storage_.empty()) + buf_span_ = std::span{storage_.data(), storage_.size()}; + else + buf_span_ = other.buf_span_; return *this; } - BlockBuffer(Tensor tensor, IndexVector blockid): - allocated{true}, indexedTensor{tensor, blockid} { - const size_t size = tensor.block_size(blockid); - T* buffer = new T[size]; - buf_span = span{buffer, size}; - tensor.get(blockid, buf_span); + + // Move ctor: after std::vector move, storage_ is valid in *this and + // empty in other; retarget span. + BlockBuffer(BlockBuffer&& other) noexcept + : storage_{std::move(other.storage_)}, + buf_span_{other.buf_span_}, + indexedTensor_{std::move(other.indexedTensor_)}, + re_{other.re_} { + if (!storage_.empty()) + buf_span_ = std::span{storage_.data(), storage_.size()}; + other.buf_span_ = {}; + other.re_ = nullptr; } - ~BlockBuffer() { - if(allocated) delete[] buf_span.data(); + + BlockBuffer& operator=(BlockBuffer&& other) noexcept { + if (this == &other) return *this; + storage_ = std::move(other.storage_); + buf_span_ = other.buf_span_; + indexedTensor_ = std::move(other.indexedTensor_); + re_ = other.re_; + if (!storage_.empty()) + buf_span_ = std::span{storage_.data(), storage_.size()}; + other.buf_span_ = {}; + other.re_ = nullptr; + return *this; } - // Whatever else is necessary to make the type regular - - auto begin() { return buf_span.begin(); } - auto begin() const { return buf_span.begin(); } - auto end() { return buf_span.end(); } - auto end() const { return buf_span.end(); } - auto get_span() { return buf_span; } - const auto get_span() const { return buf_span; } - auto data() { return buf_span.data(); } - const auto data() const { return buf_span.data(); } - void release_put() { - indexedTensor.put(buf_span); - release(); + ~BlockBuffer() = default; + + // ------------------------------------------------------------------- + // Iterators / data access + // ------------------------------------------------------------------- + auto begin() { return buf_span_.begin(); } + auto begin() const { return buf_span_.begin(); } + auto end() { return buf_span_.end(); } + auto end() const { return buf_span_.end(); } + + [[nodiscard]] std::span get_span() { return buf_span_; } + [[nodiscard]] std::span get_span() const { return buf_span_; } + [[nodiscard]] T* data() { return buf_span_.data(); } + [[nodiscard]] const T* data() const { return buf_span_.data(); } + + // ------------------------------------------------------------------- + // Release helpers (write-back and free) + // ------------------------------------------------------------------- + void release_put() { + indexedTensor_.put(buf_span_); + storage_.clear(); + buf_span_ = {}; } void release_put(Tensor tensor, IndexVector blockid) { - tensor.put(blockid, buf_span); - release(); + tensor.put(blockid, buf_span_); + storage_.clear(); + buf_span_ = {}; } void release_add() { - indexedTensor.add(buf_span); - release(); + indexedTensor_.add(buf_span_); + storage_.clear(); + buf_span_ = {}; } void release_add(Tensor tensor, IndexVector blockid) { - tensor.add(blockid, buf_span); - release(); + tensor.add(blockid, buf_span_); + storage_.clear(); + buf_span_ = {}; } void release() { - if(allocated) { - delete[] buf_span.data(); - allocated = false; - } + storage_.clear(); + buf_span_ = {}; } - std::vector block_dims() { return indexedTensor.first.block_dims(indexedTensor.second); } + + [[nodiscard]] std::vector block_dims() { + return indexedTensor_.first.block_dims(indexedTensor_.second); + } + template BlockBuffer& operator=(const V val) { std::fill(begin(), end(), val); @@ -98,25 +156,22 @@ class BlockBuffer { } private: - span buf_span; - bool allocated = false; - IndexedTensor indexedTensor; - // re is a pointer to allow it to be uninitialized. - RuntimeEngine* re; + std::vector storage_; ///< owns memory when non-empty + std::span buf_span_; ///< non-owning view (into storage_ or external) + IndexedTensor indexedTensor_; + RuntimeEngine* re_{nullptr}; }; template -bool operator==(const BlockBuffer lhs, const BlockBuffer rhs) { - return lhs.size = rhs.size && - std::equal(lhs.get_data(), lhs.get_data() + lhs.get_size(), rhs.get_data()) && - lhs.get_tensor() == rhs.get_tensor() && - lhs.get_block_id() == rhs.get_block_id(); +bool operator==(const BlockBuffer& lhs, const BlockBuffer& rhs) { + return lhs.get_span().size() == rhs.get_span().size() && + std::equal(lhs.get_span().begin(), lhs.get_span().end(), + rhs.get_span().begin()); } template inline auto& operator<<(Stream& os, BlockBuffer bf) { - // std::copy(bf.begin(), bf.end(), std::ostream_iterator(os, " ")); - for(auto it = bf.begin(); it != bf.end(); ++it) { os << *it << " "; } + for (auto it = bf.begin(); it != bf.end(); ++it) { os << *it << " "; } return os; } diff --git a/src/tamm/block_lambda_plan.hpp b/src/tamm/block_lambda_plan.hpp index 5e7ff339e..c0e3f4aee 100644 --- a/src/tamm/block_lambda_plan.hpp +++ b/src/tamm/block_lambda_plan.hpp @@ -2,6 +2,8 @@ #include "tamm/blockops_cpu.hpp" #include "tamm/types.hpp" +#include +#include namespace tamm { /////////////////////////////////////////////////////////////////////////////// @@ -46,23 +48,17 @@ class BlockLambdaPlan { void prep_flat_plan() { plan_ = Plan::flat; std::array labels_list{labels_}; - for(size_t i = 1; i < labels_list.size(); i++) { - if(labels_list[0].size() != labels_list[1].size()) { - plan_ = Plan::invalid; - return; - } - } - if(labels_list.size() > 0 && - internal::unique_entries(labels_list[0]).size() != labels_list[0].size()) { - plan_ = Plan::invalid; - return; - } - for(size_t i = 1; i < labels_list.size(); i++) { - if(!std::equal(labels_list[0].begin(), labels_list[0].end(), labels_list[i].begin())) { - plan_ = Plan::invalid; - return; - } - } + if(labels_list.empty()) return; + + // The flat plan applies only when every operand has the exact same label + // list (same size and same labels), and those labels are unique. + // Bug fix: the size check previously compared against labels_list[1] for + // every i (a fixed index) instead of labels_list[i]. + const auto& ref = labels_list[0]; + const bool all_same = + internal::unique_entries(ref).size() == ref.size() && + std::ranges::all_of(labels_list, [&](const auto& l) { return std::ranges::equal(l, ref); }); + if(!all_same) { plan_ = Plan::invalid; } } void prep_ipgen_loop_plan() { diff --git a/src/tamm/block_mult_plan.hpp b/src/tamm/block_mult_plan.hpp index e1a0e36a4..3688d6b69 100644 --- a/src/tamm/block_mult_plan.hpp +++ b/src/tamm/block_mult_plan.hpp @@ -4,40 +4,26 @@ #include "tamm/block_span.hpp" #include "tamm/blockops_blas.hpp" #include "tamm/errors.hpp" +#include "tamm/scalar.hpp" #include "tamm/tiled_index_space.hpp" #include "tamm/types.hpp" #include +#include // std::byte +#include /** - * @brief Block multiply plan selection logic. - * - * - Terms: - * - Reduction index: an index in RHS but not in LHS. e.g., j in A(i) += - * B(i,j) - * - Hadamard index: an index in LHS and both RHS tensors. e.g., l in A(l,i,j) - * += B(l,i,k) * C(l,k,j) - * - * - Choose FLAT plan if: - * - rhs are scalars or 1d - * - lhs is 1d - * - lhs and non-scalar rhs have the same label - * - * - else, choose LOOP GEMM plan if: - * - No repeated labels in any labeled tensor - * - No reduction labels - * - Hadamard labels (if any) are outermost in all tensors - * - * - else, choose LOOP TTGT plan if: - * - No repeated labels in any labeled tensors - * - No reduction labels - * - * - else, choose general plan if: - * - No repeated labels in any labeled tensors - * - * - else: - * - NOT_IMPLEMENTED() for now + * @brief Block multiply plan selection and dispatch. * + * C++20 / perf design: + * - BlockMultPlan stores a cached std::variant of the selected sub-plan, + * built once in the constructor. apply_impl() is a zero-copy std::visit + * dispatch — no per-call plan reconstruction or label-vector copies. + * - GeneralMultPlan caches permuted block dims (recomputed only on shape + * change) and intermediate byte-buffers (resized only on growth), + * eliminating all hot-path allocations in the CCSD contraction loop. + * - GemmPlan stores raw data pointers (void*) directly, not + * pointer-to-pointer, avoiding the dangling-address UB of &lhs.buf(). */ namespace tamm::internal { @@ -60,13 +46,7 @@ class FlatBlockMultPlan { switch(op_type_) { case FlatOpType::scalar_scalar: scalar_mult_assign(lhs, rhs1, rhs2); break; case FlatOpType::scalar_vector: scalar_vec_mult_assign(lhs, rhs1, rhs2); break; - case FlatOpType::vector_vector: - vec_vec_mult_assign(lhs, rhs1, rhs2); - break; - /// @bug: clang doen't like to have defaults when enum type is - /// used in switch cases - // default: - // break; + case FlatOpType::vector_vector: vec_vec_mult_assign(lhs, rhs1, rhs2); break; } } @@ -76,13 +56,7 @@ class FlatBlockMultPlan { switch(op_type_) { case FlatOpType::scalar_scalar: scalar_mult_assign(lhs, alpha, rhs1, rhs2); break; case FlatOpType::scalar_vector: scalar_vec_mult_assign(lhs, alpha, rhs1, rhs2); break; - case FlatOpType::vector_vector: - vec_vec_mult_assign(lhs, alpha, rhs1, rhs2); - break; - /// @bug: clang doen't like to have defaults when enum type is - /// used in switch cases - // default: - // break; + case FlatOpType::vector_vector: vec_vec_mult_assign(lhs, alpha, rhs1, rhs2); break; } } @@ -91,13 +65,7 @@ class FlatBlockMultPlan { switch(op_type_) { case FlatOpType::scalar_scalar: scalar_mult_update(lhs, rhs1, rhs2); break; case FlatOpType::scalar_vector: scalar_vec_mult_update(lhs, rhs1, rhs2); break; - case FlatOpType::vector_vector: - vec_vec_mult_update(lhs, rhs1, rhs2); - break; - /// @bug: clang doen't like to have defaults when enum type is - /// used in switch cases - // default: - // break; + case FlatOpType::vector_vector: vec_vec_mult_update(lhs, rhs1, rhs2); break; } } @@ -107,13 +75,7 @@ class FlatBlockMultPlan { switch(op_type_) { case FlatOpType::scalar_scalar: scalar_mult_update(lhs, alpha, rhs1, rhs2); break; case FlatOpType::scalar_vector: scalar_vec_mult_update(lhs, alpha, rhs1, rhs2); break; - case FlatOpType::vector_vector: - vec_vec_mult_update(lhs, alpha, rhs1, rhs2); - break; - /// @bug: clang doen't like to have defaults when enum type is - /// used in switch cases - // default: - // break; + case FlatOpType::vector_vector: vec_vec_mult_update(lhs, alpha, rhs1, rhs2); break; } } @@ -123,22 +85,12 @@ class FlatBlockMultPlan { switch(op_type_) { case FlatOpType::scalar_scalar: scalar_mult_update(beta, lhs, alpha, rhs1, rhs2); break; case FlatOpType::scalar_vector: scalar_vec_mult_update(beta, lhs, alpha, rhs1, rhs2); break; - case FlatOpType::vector_vector: - vec_vec_mult_update(beta, lhs, alpha, rhs1, rhs2); - break; - /// @bug: clang doen't like to have defaults when enum type is - /// used in switch cases - // default: - // break; + case FlatOpType::vector_vector: vec_vec_mult_update(beta, lhs, alpha, rhs1, rhs2); break; } } private: - enum class FlatOpType { - scalar_scalar, - scalar_vector, - vector_vector, - }; + enum class FlatOpType { scalar_scalar, scalar_vector, vector_vector }; void prep() { EXPECTS(lhs_labels_.size() < 2); @@ -146,114 +98,111 @@ class FlatBlockMultPlan { EXPECTS(rhs2_labels_.size() < 2); if(rhs1_labels_.size() == 0 && rhs2_labels_.size() == 0 && lhs_labels_.size() == 0) { - op_type_ = FlatOpType::scalar_scalar; - valid_ = true; - return; + op_type_ = FlatOpType::scalar_scalar; valid_ = true; return; } else if(rhs1_labels_.size() == 1 && rhs2_labels_.size() == 1 && lhs_labels_.size() == 1) { - op_type_ = FlatOpType::vector_vector; - valid_ = true; - return; + op_type_ = FlatOpType::vector_vector; valid_ = true; return; } else if((lhs_labels_.size() == 1 && rhs1_labels_.size() == 1) || rhs2_labels_.size() == 1) { - op_type_ = FlatOpType::scalar_vector; - valid_ = true; + op_type_ = FlatOpType::scalar_vector; valid_ = true; } else { valid_ = false; } } template - void scalar_mult_update(T beta, BlockSpan& lhs, T alpha, const BlockSpan& rhs1, - const BlockSpan& rhs2) { + void scalar_mult_update(T beta, BlockSpan& lhs, T alpha, + const BlockSpan& rhs1, const BlockSpan& rhs2) { lhs[0] = (beta * lhs[0]) + (alpha * rhs1[0] * rhs2[0]); } - template - void scalar_mult_update(BlockSpan& lhs, T alpha, const BlockSpan& rhs1, - const BlockSpan& rhs2) { + void scalar_mult_update(BlockSpan& lhs, T alpha, + const BlockSpan& rhs1, const BlockSpan& rhs2) { lhs[0] += alpha * rhs1[0] * rhs2[0]; } - template - void scalar_mult_update(BlockSpan& lhs, const BlockSpan& rhs1, const BlockSpan& rhs2) { + void scalar_mult_update(BlockSpan& lhs, + const BlockSpan& rhs1, const BlockSpan& rhs2) { lhs[0] += rhs1[0] * rhs2[0]; } - template - void scalar_mult_assign(BlockSpan& lhs, T alpha, const BlockSpan& rhs1, - const BlockSpan& rhs2) { + void scalar_mult_assign(BlockSpan& lhs, T alpha, + const BlockSpan& rhs1, const BlockSpan& rhs2) { lhs[0] = alpha * rhs1[0] * rhs2[0]; } - template - void scalar_mult_assign(BlockSpan& lhs, const BlockSpan& rhs1, const BlockSpan& rhs2) { + void scalar_mult_assign(BlockSpan& lhs, + const BlockSpan& rhs1, const BlockSpan& rhs2) { lhs[0] = rhs1[0] * rhs2[0]; } - template void scalar_vec_mult_update(T beta, BlockSpan& lhs_vec, T alpha, - const BlockSpan& rhs1_scalar, const BlockSpan& rhs2_vec) { + const BlockSpan& rhs1_scalar, + const BlockSpan& rhs2_vec) { EXPECTS(lhs_vec.num_elements() == rhs2_vec.num_elements()); - for(size_t i = 0; i < lhs_vec.num_elements(); i++) { - auto new_alpha = alpha * rhs1_scalar[0]; - blockops::cpu::flat_update(beta, lhs_vec, new_alpha, rhs2_vec); - } + blockops::cpu::flat_update(beta, lhs_vec, alpha * rhs1_scalar[0], rhs2_vec); } - template - void scalar_vec_mult_update(BlockSpan& lhs_vec, T alpha, const BlockSpan& rhs1_scalar, + void scalar_vec_mult_update(BlockSpan& lhs_vec, T alpha, + const BlockSpan& rhs1_scalar, const BlockSpan& rhs2_vec) { EXPECTS(lhs_vec.num_elements() == rhs2_vec.num_elements()); - for(size_t i = 0; i < lhs_vec.num_elements(); i++) { - auto new_alpha = alpha * rhs1_scalar[0]; - blockops::cpu::flat_update(lhs_vec, new_alpha, rhs2_vec); - } + blockops::cpu::flat_update(lhs_vec, alpha * rhs1_scalar[0], rhs2_vec); } - template - void scalar_vec_mult_assign(BlockSpan& lhs_vec, T alpha, const BlockSpan& rhs1_scalar, + void scalar_vec_mult_assign(BlockSpan& lhs_vec, T alpha, + const BlockSpan& rhs1_scalar, const BlockSpan& rhs2_vec) { EXPECTS(lhs_vec.num_elements() == rhs2_vec.num_elements()); - for(size_t i = 0; i < lhs_vec.num_elements(); i++) { - auto new_alpha = alpha * rhs1_scalar[0]; - blockops::cpu::flat_assign(lhs_vec, new_alpha, rhs2_vec); - } + blockops::cpu::flat_assign(lhs_vec, alpha * rhs1_scalar[0], rhs2_vec); } - template - void vec_vec_mult_update(T beta, BlockSpan& lhs_vec, T alpha, const BlockSpan& rhs1_vec, + void scalar_vec_mult_update(BlockSpan& lhs_vec, + const BlockSpan& rhs1_scalar, + const BlockSpan& rhs2_vec) { + EXPECTS(lhs_vec.num_elements() == rhs2_vec.num_elements()); + blockops::cpu::flat_update(lhs_vec, rhs1_scalar[0], rhs2_vec); + } + template + void scalar_vec_mult_assign(BlockSpan& lhs_vec, + const BlockSpan& rhs1_scalar, + const BlockSpan& rhs2_vec) { + EXPECTS(lhs_vec.num_elements() == rhs2_vec.num_elements()); + blockops::cpu::flat_assign(lhs_vec, rhs1_scalar[0], rhs2_vec); + } + template + void vec_vec_mult_update(T beta, BlockSpan& lhs_vec, T alpha, + const BlockSpan& rhs1_vec, const BlockSpan& rhs2_vec) { - for(size_t i = 0; i < lhs_vec.num_elements(); i++) { + for(size_t i = 0; i < lhs_vec.num_elements(); i++) lhs_vec[i] = (beta * lhs_vec[i]) + (alpha * rhs1_vec[i] * rhs2_vec[i]); - } } - template - void vec_vec_mult_update(BlockSpan& lhs_vec, T alpha, const BlockSpan& rhs1_vec, + void vec_vec_mult_update(BlockSpan& lhs_vec, T alpha, + const BlockSpan& rhs1_vec, const BlockSpan& rhs2_vec) { - for(size_t i = 0; i < lhs_vec.num_elements(); i++) { + for(size_t i = 0; i < lhs_vec.num_elements(); i++) lhs_vec[i] += alpha * rhs1_vec[i] * rhs2_vec[i]; - } } - template - void vec_vec_mult_update(BlockSpan& lhs_vec, const BlockSpan& rhs1_vec, + void vec_vec_mult_update(BlockSpan& lhs_vec, + const BlockSpan& rhs1_vec, const BlockSpan& rhs2_vec) { - for(size_t i = 0; i < lhs_vec.num_elements(); i++) { lhs_vec[i] += rhs1_vec[i] * rhs2_vec[i]; } + for(size_t i = 0; i < lhs_vec.num_elements(); i++) + lhs_vec[i] += rhs1_vec[i] * rhs2_vec[i]; } - template - void vec_vec_mult_assign(BlockSpan& lhs_vec, T alpha, const BlockSpan& rhs1_vec, + void vec_vec_mult_assign(BlockSpan& lhs_vec, T alpha, + const BlockSpan& rhs1_vec, const BlockSpan& rhs2_vec) { - for(size_t i = 0; i < lhs_vec.num_elements(); i++) { + for(size_t i = 0; i < lhs_vec.num_elements(); i++) lhs_vec[i] = alpha * rhs1_vec[i] * rhs2_vec[i]; - } } - template - void vec_vec_mult_assign(BlockSpan& lhs_vec, const BlockSpan& rhs1_vec, + void vec_vec_mult_assign(BlockSpan& lhs_vec, + const BlockSpan& rhs1_vec, const BlockSpan& rhs2_vec) { - for(size_t i = 0; i < lhs_vec.num_elements(); i++) { lhs_vec[i] = rhs1_vec[i] * rhs2_vec[i]; } + for(size_t i = 0; i < lhs_vec.num_elements(); i++) + lhs_vec[i] = rhs1_vec[i] * rhs2_vec[i]; } IndexLabelVec lhs_labels_; @@ -261,16 +210,10 @@ class FlatBlockMultPlan { IndexLabelVec rhs2_labels_; FlatOpType op_type_; bool valid_; -}; // class FlatBlockMultPlan -/** - * @brief basic ttgt plan - * - * @todo - * - Loop ttgt - * - Avoid transpose when not needed - * - Separate prep() from apply() - * - */ +}; + +// --------------------------------------------------------------------------- + class TTGTPlan { public: TTGTPlan(): valid_{false} {} @@ -284,8 +227,8 @@ class TTGTPlan { } template - void apply(T beta, BlockSpan& lhs, T alpha, const BlockSpan& rhs1, - const BlockSpan& rhs2) { + void apply(T beta, BlockSpan& lhs, T alpha, + const BlockSpan& rhs1, const BlockSpan& rhs2) { NOT_IMPLEMENTED(); } @@ -302,26 +245,24 @@ class TTGTPlan { IndexLabelVec rhs1_labels_; IndexLabelVec rhs2_labels_; bool valid_; -}; // class TTGTPlan +}; + +// --------------------------------------------------------------------------- class GemmPlan { public: GemmPlan(): valid_{false} {} - GemmPlan(const GemmPlan&) = default; GemmPlan& operator=(const GemmPlan&) = default; GemmPlan(const IndexLabelVec& lhs_labels, const IndexLabelVec& rhs1_labels, const IndexLabelVec& rhs2_labels) { if(has_repeated_indices(lhs_labels) || has_repeated_indices(rhs1_labels) || - has_repeated_indices(rhs2_labels)) { - return; - } + has_repeated_indices(rhs2_labels)) { return; } auto find = [](const auto& collection, const auto& element) { return std::find(collection.begin(), collection.end(), element); }; - auto has = [find](const auto& collection, const auto& element) { return find(collection, element) != collection.end(); }; @@ -337,9 +278,7 @@ class GemmPlan { if(num_mindices < 0 || num_nindices < 0 || num_kindices < 0) { return; } if((lhs_labels.size() != num_mindices + num_nindices) || (arg1_labels.size() != num_mindices + num_kindices) || - (arg2_labels.size() != num_kindices + num_nindices)) { - return; - } + (arg2_labels.size() != num_kindices + num_nindices)) { return; } bool transpose_arg1 = (num_mindices > 0 && lhs_labels[0] != arg1_labels[0]); bool transpose_arg2 = (num_nindices > 0 && lhs_labels[num_mindices] == arg2_labels[0]); @@ -349,18 +288,15 @@ class GemmPlan { int arg1_kpos = transpose_arg1 ? 0 : num_mindices; int arg2_kpos = transpose_arg2 ? num_nindices : 0; int arg2_npos = transpose_arg2 ? 0 : num_kindices; - if(!std::equal(lhs_labels.begin() + lhs_mpos, lhs_labels.begin() + lhs_mpos + num_mindices, - arg1_labels.begin() + arg1_mpos)) { - return; - } - if(!std::equal(lhs_labels.begin() + lhs_npos, lhs_labels.begin() + lhs_npos + num_nindices, - arg2_labels.begin() + arg2_npos)) { - return; - } - if(!std::equal(arg1_labels.begin() + arg1_kpos, arg1_labels.begin() + arg1_kpos + num_kindices, - arg2_labels.begin() + arg2_kpos)) { - return; - } + if(!std::equal(lhs_labels.begin() + lhs_mpos, + lhs_labels.begin() + lhs_mpos + num_mindices, + arg1_labels.begin() + arg1_mpos)) { return; } + if(!std::equal(lhs_labels.begin() + lhs_npos, + lhs_labels.begin() + lhs_npos + num_nindices, + arg2_labels.begin() + arg2_npos)) { return; } + if(!std::equal(arg1_labels.begin() + arg1_kpos, + arg1_labels.begin() + arg1_kpos + num_kindices, + arg2_labels.begin() + arg2_kpos)) { return; } num_mindices_ = num_mindices; num_nindices_ = num_nindices; @@ -373,27 +309,27 @@ class GemmPlan { bool is_valid() const { return valid_; } - template - void apply(T0 beta, BlockSpan& lhs, T1 alpha, BlockSpan& rhs1, BlockSpan& rhs2) { + template + void apply(T0 beta, BlockSpan& lhs, T1 alpha, + BlockSpan& rhs1, BlockSpan& rhs2) { EXPECTS(valid_); update_plan(lhs, rhs1, rhs2); size_t loff = 0, r1off = 0, r2off = 0; - for(size_t i = 0; i < num_batches_; i++) { + for(size_t i = 0; i < static_cast(num_batches_); i++) { auto TransA = transpose_arg1_ ? blas::Op::Trans : blas::Op::NoTrans; auto TransB = transpose_arg2_ ? blas::Op::Trans : blas::Op::NoTrans; int lda = !transpose_arg1_ ? K_ : M_; int ldb = !transpose_arg2_ ? N_ : K_; int ldc = N_; - - auto* bufa = rhs1_is_arg1_ ? reinterpret_cast(&bufa_) : reinterpret_cast(&bufa_); - auto* bufb = rhs1_is_arg1_ ? reinterpret_cast(&bufb_) : reinterpret_cast(&bufb_); - auto* bufc = reinterpret_cast(&bufc_); - - blas::gemm(blas::Layout::RowMajor, TransA, TransB, M_, N_, K_, alpha, &bufa[r1off], lda, - &bufb[r2off], ldb, beta, &bufc[loff], ldc); - loff += M_ * N_; - r1off += M_ * K_; - r2off += K_ * N_; + // Fix: store and use raw data pointers directly (not pointer-to-pointer). + auto* bufa = static_cast(bufa_); + auto* bufb = static_cast(bufb_); + auto* bufc = static_cast(bufc_); + blas::gemm(blas::Layout::RowMajor, TransA, TransB, M_, N_, K_, + alpha, bufa + r1off, lda, + bufb + r2off, ldb, + beta, bufc + loff, ldc); + loff += M_ * N_; r1off += M_ * K_; r2off += K_ * N_; } } @@ -403,54 +339,47 @@ class GemmPlan { } template - void update_plan(const BlockSpan& lhs, const BlockSpan& rhs1, const BlockSpan& rhs2) { - M_ = 1; - N_ = 1; - K_ = 1; - - const std::vector& lhs_bdims = lhs.block_dims(); - for(size_t i = 0; i < num_mindices_; i++) { M_ *= lhs_bdims[num_batch_indices_ + i]; } - for(size_t i = num_mindices_; i < lhs_bdims.size(); i++) { + void update_plan(const BlockSpan& lhs, + const BlockSpan& rhs1, + const BlockSpan& rhs2) { + M_ = 1; N_ = 1; K_ = 1; + const auto& lhs_bdims = lhs.block_dims(); // const& — no copy + for(size_t i = 0; i < static_cast(num_mindices_); i++) + M_ *= lhs_bdims[num_batch_indices_ + i]; + for(size_t i = static_cast(num_mindices_); i < lhs_bdims.size(); i++) N_ *= lhs_bdims[num_batch_indices_ + i]; - } - size_t kstart_idx = (transpose_arg1_ ? 0 : num_mindices_); + size_t kstart_idx = (transpose_arg1_ ? 0 : static_cast(num_mindices_)); const auto& arg1_bdims = (rhs1_is_arg1_ ? rhs1.block_dims() : rhs2.block_dims()); - for(int k = 0; k < num_kindices_; k++) { + for(int k = 0; k < num_kindices_; k++) K_ *= arg1_bdims[num_batch_indices_ + kstart_idx + k]; - } - bufc_ = &lhs.buf(); - bufa_ = rhs1_is_arg1_ ? &rhs1.buf() : &rhs2.buf(); - bufb_ = rhs1_is_arg1_ ? &rhs2.buf() : &rhs1.buf(); - - int num_batches_ = 1; - for(int i = 0; i < num_batch_indices_; i++) { num_batches_ *= lhs_bdims[i]; } - } + // Fix Bug 2: store the data pointer VALUE (void*), not the address of + // the T* temporary returned by buf(). Previously '&lhs.buf()' took + // the address of a prvalue, which is a dangling pointer / UB. + bufc_ = static_cast(lhs.buf()); + bufa_ = static_cast(rhs1_is_arg1_ ? rhs1.buf() : rhs2.buf()); + bufb_ = static_cast(rhs1_is_arg1_ ? rhs2.buf() : rhs1.buf()); + num_batches_ = 1; + for(int i = 0; i < num_batch_indices_; i++) num_batches_ *= lhs_bdims[i]; + } + + bool valid_{false}; + int num_batch_indices_{0}; // in-class init: was uninitialized (UB) before commit 3ac53c2 + int num_mindices_{}; + int num_nindices_{}; + int num_kindices_{}; + bool rhs1_is_arg1_{}; + bool transpose_arg1_{}; + bool transpose_arg2_{}; + int M_{}, N_{}, K_{}; + int num_batches_{}; + // Fix Bug 2: plain void* storing the data pointer value directly. + void* bufa_{}; + void* bufb_{}; + void* bufc_{}; +}; + +// --------------------------------------------------------------------------- - bool valid_; - int num_batch_indices_; - int num_mindices_; - int num_nindices_; - int num_kindices_; // number of summation indices - bool rhs1_is_arg1_; // true if LHS = RHS1()*RHS2(). False if LHS = - // RHS2()*RHS1() - bool transpose_arg1_; - bool transpose_arg2_; - - int M_, N_, K_; - int num_batches_; - void *bufa_, *bufb_, *bufc_; -}; // class GemmPlan - -/** - * @brief General mult plan with repeated indices. This is not designed for - * sparse blocks (labels with dependent indices). - * - * @todo - * - Avoid block assign plans when not needed - * - There could bt multiple intermediate, one here and one in TTGT. Incorporate - * TTGT functionality here to minimize copies. - * - */ class GeneralMultPlan { public: GeneralMultPlan(): valid_{false} {} @@ -459,50 +388,84 @@ class GeneralMultPlan { bool is_valid() const { return valid_; } - GeneralMultPlan(const IndexLabelVec& lhs_labels, const IndexLabelVec& rhs1_labels, + GeneralMultPlan(const IndexLabelVec& lhs_labels, + const IndexLabelVec& rhs1_labels, const IndexLabelVec& rhs2_labels): - valid_{true}, lhs_labels_{lhs_labels}, rhs1_labels_{rhs1_labels}, rhs2_labels_{rhs2_labels} { + valid_{true}, + lhs_labels_{lhs_labels}, + rhs1_labels_{rhs1_labels}, + rhs2_labels_{rhs2_labels} { for(const auto& lbl: lhs_labels) { - if(std::find(lhs_inter_labels_.begin(), lhs_inter_labels_.end(), lbl) == - lhs_inter_labels_.end()) { + if(std::find(lhs_inter_labels_.begin(), lhs_inter_labels_.end(), lbl) + == lhs_inter_labels_.end()) lhs_inter_labels_.push_back(lbl); - } } for(const auto& lbl: rhs1_labels) { - if(std::find(rhs1_inter_labels_.begin(), rhs1_inter_labels_.end(), lbl) == - rhs1_inter_labels_.end()) { + if(std::find(rhs1_inter_labels_.begin(), rhs1_inter_labels_.end(), lbl) + == rhs1_inter_labels_.end()) rhs1_inter_labels_.push_back(lbl); - } } for(const auto& lbl: rhs2_labels) { - if(std::find(rhs2_inter_labels_.begin(), rhs2_inter_labels_.end(), lbl) == - rhs2_inter_labels_.end()) { + if(std::find(rhs2_inter_labels_.begin(), rhs2_inter_labels_.end(), lbl) + == rhs2_inter_labels_.end()) rhs2_inter_labels_.push_back(lbl); - } } - lhs_ba_plan_ = BlockAssignPlan{lhs_labels_, lhs_inter_labels_, BlockAssignPlan::OpType::set}; - rhs1_ba_plan_ = BlockAssignPlan{rhs1_inter_labels_, rhs1_labels_, BlockAssignPlan::OpType::set}; - rhs2_ba_plan_ = BlockAssignPlan{rhs2_inter_labels_, rhs2_labels_, BlockAssignPlan::OpType::set}; - linter_to_l_perm_ = perm_map_compute(lhs_inter_labels_, lhs_labels_); + lhs_ba_plan_ = BlockAssignPlan{lhs_labels_, lhs_inter_labels_, + BlockAssignPlan::OpType::set}; + rhs1_ba_plan_ = BlockAssignPlan{rhs1_inter_labels_, rhs1_labels_, + BlockAssignPlan::OpType::set}; + rhs2_ba_plan_ = BlockAssignPlan{rhs2_inter_labels_, rhs2_labels_, + BlockAssignPlan::OpType::set}; + linter_to_l_perm_ = perm_map_compute(lhs_inter_labels_, lhs_labels_); r1_to_r1inter_perm_ = perm_map_compute(rhs1_labels_, rhs1_inter_labels_); r2_to_r2inter_perm_ = perm_map_compute(rhs2_labels_, rhs2_inter_labels_); - ttgt_plan_ = TTGTPlan{lhs_inter_labels_, rhs1_inter_labels_, rhs2_inter_labels_}; + ttgt_plan_ = TTGTPlan{lhs_inter_labels_, rhs1_inter_labels_, rhs2_inter_labels_}; } - template - void apply(T0 beta, BlockSpan& lhs, T1 alpha, BlockSpan& rhs1, BlockSpan& rhs2) { + template + void apply(T0 beta, BlockSpan& lhs, T1 alpha, + BlockSpan& rhs1, BlockSpan& rhs2) { EXPECTS(valid_); - std::vector linter_buf(lhs.num_elements()); - std::vector r1inter_buf(rhs1.num_elements()); - std::vector r2inter_buf(rhs2.num_elements()); - - const auto& linter_dims = perm_map_apply(lhs.block_dims(), linter_to_l_perm_); - const auto& r1inter_dims = perm_map_apply(rhs1.block_dims(), r1_to_r1inter_perm_); - const auto& r2inter_dims = perm_map_apply(rhs2.block_dims(), r2_to_r2inter_perm_); - BlockSpan lhs_inter{linter_buf.data(), linter_dims}; - BlockSpan rhs1_inter{r1inter_buf.data(), r1inter_dims}; - BlockSpan rhs2_inter{r2inter_buf.data(), r2inter_dims}; + const size_t lhs_nelems = lhs.num_elements(); + const size_t rhs1_nelems = rhs1.num_elements(); + const size_t rhs2_nelems = rhs2.num_elements(); + + // Fix Bug 1: intermediate buffers are std::vector sized in + // bytes, not std::vector. This handles float, double, and + // complex types uniformly without type-punning or memory corruption. + if(linter_buf_.size() < lhs_nelems * sizeof(T1)) linter_buf_.resize(lhs_nelems * sizeof(T1)); + if(r1inter_buf_.size() < rhs1_nelems * sizeof(T2)) r1inter_buf_.resize(rhs1_nelems * sizeof(T2)); + if(r2inter_buf_.size() < rhs2_nelems * sizeof(T3)) r2inter_buf_.resize(rhs2_nelems * sizeof(T3)); + + // Cache permuted dims — recompute only when block shape changes. + const auto& lhs_dims = lhs.block_dims(); // const& — no copy + const auto& rhs1_dims = rhs1.block_dims(); + const auto& rhs2_dims = rhs2.block_dims(); + + if(linter_dims_.empty() || + !std::equal(lhs_dims.begin(), lhs_dims.end(), linter_dims_src_.begin())) { + linter_dims_src_.assign(lhs_dims.begin(), lhs_dims.end()); + linter_dims_ = perm_map_apply(linter_dims_src_, linter_to_l_perm_); + } + if(r1inter_dims_.empty() || + !std::equal(rhs1_dims.begin(), rhs1_dims.end(), r1inter_dims_src_.begin())) { + r1inter_dims_src_.assign(rhs1_dims.begin(), rhs1_dims.end()); + r1inter_dims_ = perm_map_apply(r1inter_dims_src_, r1_to_r1inter_perm_); + } + if(r2inter_dims_.empty() || + !std::equal(rhs2_dims.begin(), rhs2_dims.end(), r2inter_dims_src_.begin())) { + r2inter_dims_src_.assign(rhs2_dims.begin(), rhs2_dims.end()); + r2inter_dims_ = perm_map_apply(r2inter_dims_src_, r2_to_r2inter_perm_); + } + + // Reinterpret byte buffers as the correct element type. + BlockSpan lhs_inter { + reinterpret_cast(linter_buf_.data()), linter_dims_}; + BlockSpan rhs1_inter{ + reinterpret_cast(r1inter_buf_.data()), r1inter_dims_}; + BlockSpan rhs2_inter{ + reinterpret_cast(r2inter_buf_.data()), r2inter_dims_}; rhs1_ba_plan_.apply(rhs1_inter, rhs1); rhs2_ba_plan_.apply(rhs2_inter, rhs2); @@ -525,16 +488,39 @@ class GeneralMultPlan { BlockAssignPlan rhs2_ba_plan_; BlockAssignPlan lhs_ba_plan_; TTGTPlan ttgt_plan_; -}; // class GeneralMultPlan + // Fix Bug 1: byte buffers — correct for any element type T1/T2/T3. + std::vector linter_buf_; + std::vector r1inter_buf_; + std::vector r2inter_buf_; + // Cached permuted dims (recomputed only when block shape changes). + std::vector linter_dims_src_, linter_dims_; + std::vector r1inter_dims_src_, r1inter_dims_; + std::vector r2inter_dims_src_, r2inter_dims_; +}; } // namespace tamm::internal +// --------------------------------------------------------------------------- + namespace tamm { + +/** + * @brief BlockMultPlan — select and cache the best contraction sub-plan. + * + * The sub-plan is selected once in the constructor; apply_impl() dispatches + * via std::visit with zero label-vector copies at call time. + * + * The public apply() overload (Scalar lscale/rscale) now correctly forwards + * to apply_impl() instead of silently doing nothing. + */ class BlockMultPlan { public: enum class OpType { set, update }; - BlockMultPlan(const IndexLabelVec& lhs_labels, const IndexLabelVec& rhs1_labels, - const IndexLabelVec& rhs2_labels, OpType optype): + + BlockMultPlan(const IndexLabelVec& lhs_labels, + const IndexLabelVec& rhs1_labels, + const IndexLabelVec& rhs2_labels, + OpType optype): lhs_labels_{lhs_labels}, rhs1_labels_{rhs1_labels}, rhs2_labels_{rhs2_labels}, @@ -544,69 +530,96 @@ class BlockMultPlan { has_reduction_index_{has_reduction_index()}, has_hadamard_index_{has_hadamard_index()} { prep_flat_plan(); - if(plan_ == Plan::invalid) { prep_loop_gemm_plan(); } - if(plan_ == Plan::invalid) { prep_loop_ttgt_plan(); } - if(plan_ == Plan::invalid) { prep_general_plan(); } + if(plan_ == Plan::invalid) prep_loop_gemm_plan(); + if(plan_ == Plan::invalid) prep_loop_ttgt_plan(); + if(plan_ == Plan::invalid) prep_general_plan(); if(plan_ == Plan::invalid) { NOT_IMPLEMENTED(); } EXPECTS(plan_ != Plan::invalid); + build_cached_plan(); + } + + // Primary hot-path dispatch: typed scalars, all types known at compile time. + template + void apply_impl(T1 lscale, BlockSpan& lhs, T1 rscale, + BlockSpan& rhs1, BlockSpan& rhs2) { + std::visit([&](auto& p) { + using PT = std::decay_t; + if constexpr (std::is_same_v) { + if(optype_ == OpType::set) + p.apply_assign(lhs, rscale, rhs1, rhs2); + else + p.apply_update(lscale, lhs, rscale, rhs1, rhs2); + } else if constexpr (std::is_same_v) { + p.apply(lscale, lhs, rscale, rhs1, rhs2); + } else if constexpr (std::is_same_v) { + p.apply(lscale, lhs, rscale, rhs1, rhs2); + } else if constexpr (std::is_same_v) { + p.apply(lscale, lhs, rscale, rhs1, rhs2); + } else { + NOT_ALLOWED(); + } + }, cached_plan_); } + // Public Scalar overload forwards to the typed dispatch. + // + // This plan is the single-element-type contraction path: lhs, rhs1 and rhs2 + // must all share the same element type. Mixed-type contractions (e.g. + // real x complex) are handled by kernels::block_multiply in multop.hpp, not + // by BlockMultPlan. We constrain the overload to T1==T2==T3 so the previous + // reinterpret_cast&>(rhs) UB (valid only when the types were + // already identical) can no longer be reached with incompatible types. template - void apply_impl(T1 lscale, BlockSpan& lhs, T1 rscale, BlockSpan& rhs1, - BlockSpan& rhs2) { + requires (std::is_same_v && std::is_same_v) + void apply(Scalar lscale, BlockSpan& lhs, Scalar rscale, + BlockSpan& rhs1, BlockSpan& rhs2) { + apply_impl(lscale.template get(), lhs, + rscale.template get(), rhs1, rhs2); + } + +private: + using CachedPlan = std::variant< + internal::FlatBlockMultPlan, + internal::GemmPlan, + internal::TTGTPlan, + internal::GeneralMultPlan + >; + + void build_cached_plan() { switch(plan_) { - case Plan::invalid: { - NOT_ALLOWED(); + case Plan::flat_assign: + case Plan::flat_update: + cached_plan_ = internal::FlatBlockMultPlan{lhs_labels_, rhs1_labels_, rhs2_labels_}; break; - } - case Plan::flat_assign: { - internal::FlatBlockMultPlan flat_plan{lhs_labels_, rhs1_labels_, rhs2_labels_}; - flat_plan.apply_assign(lscale, lhs, rscale, rhs1, rhs2, optype_); + case Plan::loop_gemm: + cached_plan_ = internal::GemmPlan{lhs_labels_, rhs1_labels_, rhs2_labels_}; break; - } - case Plan::loop_gemm: { - internal::GemmPlan gemm_plan{lhs_labels_, rhs1_labels_, rhs2_labels_}; - gemm_plan.apply(lscale, lhs, rscale, rhs1, rhs2); - break; - } - case Plan::loop_ttgt: { - internal::TTGTPlan ttgt_plan{lhs_labels_, rhs1_labels_, rhs2_labels_}; - ttgt_plan.apply(lscale, lhs, rscale, rhs1, rhs2); + case Plan::loop_ttgt: + cached_plan_ = internal::TTGTPlan{lhs_labels_, rhs1_labels_, rhs2_labels_}; break; - } - case Plan::general: { - internal::GeneralMultPlan general_plan{lhs_labels_, rhs1_labels_, rhs2_labels_}; - general_plan.apply(lscale, lhs, rscale, rhs1, rhs2); + case Plan::general: + cached_plan_ = internal::GeneralMultPlan{lhs_labels_, rhs1_labels_, rhs2_labels_}; break; - } default: UNREACHABLE(); } } - template - void apply(Scalar lscale, BlockSpan& lhs, Scalar rscale, BlockSpan& rhs1, - BlockSpan& rhs2) {} - -private: bool has_reduction_index() { - std::set lhs_labels(lhs_labels_.begin(), lhs_labels_.end()); - std::set rhs_labels(rhs1_labels_.begin(), rhs1_labels_.end()); - for(auto lbl: rhs2_labels_) { rhs_labels.insert(lbl); } - + std::set lhs_set(lhs_labels_.begin(), lhs_labels_.end()); + std::set rhs_set(rhs1_labels_.begin(), rhs1_labels_.end()); + for(auto lbl: rhs2_labels_) rhs_set.insert(lbl); IndexLabelVec reduction_lbls; - - std::set_difference(rhs_labels.begin(), rhs_labels.end(), lhs_labels.begin(), lhs_labels.end(), + std::set_difference(rhs_set.begin(), rhs_set.end(), + lhs_set.begin(), lhs_set.end(), std::back_inserter(reduction_lbls)); - - return reduction_lbls.size() != 0; + return !reduction_lbls.empty(); } bool has_hadamard_index() { - std::set rhs1_labels(rhs1_labels_.begin(), rhs1_labels_.end()); - std::set rhs2_labels(rhs2_labels_.begin(), rhs2_labels_.end()); - for(const auto& lbl: lhs_labels_) { - if(rhs1_labels.count(lbl) != 0 && rhs2_labels.count(lbl) != 0) { return true; } - } + std::set rhs1_set(rhs1_labels_.begin(), rhs1_labels_.end()); + std::set rhs2_set(rhs2_labels_.begin(), rhs2_labels_.end()); + for(const auto& lbl: lhs_labels_) + if(rhs1_set.count(lbl) && rhs2_set.count(lbl)) return true; return false; } @@ -617,34 +630,22 @@ class BlockMultPlan { } IndexLabelVec get_hadamard_labels() { - if(!has_hadamard_index_) { return {}; } + if(!has_hadamard_index_) return {}; IndexLabelVec result; - - std::set rhs1_labels(rhs1_labels_.begin(), rhs1_labels_.end()); - std::set rhs2_labels(rhs2_labels_.begin(), rhs2_labels_.end()); - - for(const auto& lbl: lhs_labels_) { - if(rhs1_labels.count(lbl) != 0 && rhs2_labels.count(lbl) != 0) { result.push_back(lbl); } - } - + std::set rhs1_set(rhs1_labels_.begin(), rhs1_labels_.end()); + std::set rhs2_set(rhs2_labels_.begin(), rhs2_labels_.end()); + for(const auto& lbl: lhs_labels_) + if(rhs1_set.count(lbl) && rhs2_set.count(lbl)) result.push_back(lbl); return result; } - /** - * @brief Choose FLAT plan if: - * - rhs are scalars or 1d - * - lhs is 1d - * - lhs and non-scalar rhs have the same label - */ void prep_flat_plan() { if(rhs1_labels_.size() == 0 && rhs2_labels_.size() == 0 && lhs_labels_.size() == 0) { - plan_ = optype_ == OpType::set ? Plan::flat_assign : Plan::flat_update; - return; + plan_ = optype_ == OpType::set ? Plan::flat_assign : Plan::flat_update; return; } else if(rhs1_labels_.size() == 1 && rhs2_labels_.size() == 1 && lhs_labels_.size() == 1) { - if(lhs_labels_ == rhs1_labels_ && lhs_labels_ == rhs1_labels_) { - plan_ = optype_ == OpType::set ? Plan::flat_assign : Plan::flat_update; - return; + if(lhs_labels_ == rhs1_labels_ && lhs_labels_ == rhs2_labels_) { + plan_ = optype_ == OpType::set ? Plan::flat_assign : Plan::flat_update; return; } } else if(rhs1_labels_.size() == 1 || rhs1_labels_.size() == 0) { @@ -652,73 +653,53 @@ class BlockMultPlan { size_t rhs_size = std::max(rhs1_labels_.size(), rhs2_labels_.size()); if(lhs_labels_.size() == rhs_size) { if((rhs1_labels_.size() == 1 && lhs_labels_ == rhs1_labels_) || - (rhs2_labels_.size() == 1 && lhs_labels_ == rhs2_labels_)) { + (rhs2_labels_.size() == 1 && lhs_labels_ == rhs2_labels_)) plan_ = optype_ == OpType::set ? Plan::flat_assign : Plan::flat_update; - } } } } } - /** - * @brief choose LOOP GEMM plan if: - * - No repeated labels in any labeled tensor - * - No reduction labels - * - Hadamard labels (if any) are outermost in all tensors - */ void prep_loop_gemm_plan() { if(!has_repeated_index_ && !has_reduction_index_) { if(has_hadamard_index_) { auto hadamard_labels = get_hadamard_labels(); - const ptrdiff_t hlabels = hadamard_labels.size(); + const ptrdiff_t hlabels = static_cast(hadamard_labels.size()); for(size_t i = 0; i < hadamard_labels.size(); i++) { - auto lbl = lhs_labels_[i]; - auto rhs1_pos = - std::find(rhs1_labels_.begin(), rhs2_labels_.end(), lbl) - rhs1_labels_.begin(); - auto rhs2_pos = - std::find(rhs1_labels_.begin(), rhs2_labels_.end(), lbl) - rhs1_labels_.begin(); - if(rhs1_pos >= hlabels || rhs2_pos >= hlabels) { return; } + auto lbl = lhs_labels_[i]; + auto rhs1_pos = std::find(rhs1_labels_.begin(), rhs1_labels_.end(), lbl) + - rhs1_labels_.begin(); + auto rhs2_pos = std::find(rhs2_labels_.begin(), rhs2_labels_.end(), lbl) + - rhs2_labels_.begin(); + if(rhs1_pos >= hlabels || rhs2_pos >= hlabels) return; } plan_ = Plan::loop_gemm; + } else { + plan_ = Plan::loop_gemm; } - else { plan_ = Plan::loop_gemm; } } } - /** - * @brief choose LOOP TTGT plan if: - * - No repeated labels in any labeled tensors - * - No reduction labels - */ void prep_loop_ttgt_plan() { - if(!has_repeated_index_ && !has_reduction_index_) { plan_ = Plan::loop_ttgt; } + if(!has_repeated_index_ && !has_reduction_index_) plan_ = Plan::loop_ttgt; } - /** - * @brief choose general plan if: - * - No repeated labels in any labeled tensors - */ void prep_general_plan() { - if(!has_repeated_index_) { plan_ = Plan::general; } + if(!has_repeated_index_) plan_ = Plan::general; } - enum class Plan { - flat_assign, - flat_update, - loop_gemm, - loop_ttgt, - general, - invalid, - }; + enum class Plan { flat_assign, flat_update, loop_gemm, loop_ttgt, general, invalid }; IndexLabelVec lhs_labels_; IndexLabelVec rhs1_labels_; IndexLabelVec rhs2_labels_; OpType optype_; Plan plan_; + CachedPlan cached_plan_; ///< built once in ctor; zero-copy std::visit dispatch bool has_repeated_index_; bool has_reduction_index_; bool has_hadamard_index_; -}; // class BlockMultPlan +}; + } // namespace tamm diff --git a/src/tamm/block_scratch.hpp b/src/tamm/block_scratch.hpp new file mode 100644 index 000000000..9b0f274d0 --- /dev/null +++ b/src/tamm/block_scratch.hpp @@ -0,0 +1,48 @@ +#pragma once + +#include +#include + +namespace tamm::internal { + +/** + * @brief RAII + pooled scratch for the per-block "own-or-view" buffer pattern. + * + * Many op-apply loops (set/add/mult) need a working pointer to a block that is + * EITHER a view into the tensor's local buffer (no ownership, no copy) OR a + * freshly allocated copy fetched from a remote rank. The legacy code expressed + * this with a raw `new T[n]` + a `bool alloced` flag + a matching `delete[]`, + * which leaks if anything between the `new` and the `delete[]` throws. + * + * BlockScratch owns a grow-only `std::vector` that is reused across loop + * iterations (pooling: no per-block malloc/free churn once it reaches the max + * block size) and hands back a raw `T*` for the kernel. Declare ONE instance + * outside the block loop; call view()/owned() per block. + * + * Not thread-safe: use one instance per worker thread. + */ +template +class BlockScratch { +public: + BlockScratch() = default; + BlockScratch(const BlockScratch&) = delete; + BlockScratch& operator=(const BlockScratch&) = delete; + BlockScratch(BlockScratch&&) = default; + BlockScratch& operator=(BlockScratch&&) = default; + ~BlockScratch() = default; + + /// Use an external (tensor-local) buffer directly; no ownership, no copy. + [[nodiscard]] T* view(T* external) noexcept { return external; } + + /// Provide an owned, reusable buffer of at least @p n elements. + /// The storage is grow-only and reused on subsequent calls (pooling). + [[nodiscard]] T* owned(std::size_t n) { + if(storage_.size() < n) { storage_.resize(n); } + return storage_.data(); + } + +private: + std::vector storage_; +}; + +} // namespace tamm::internal diff --git a/src/tamm/block_span.hpp b/src/tamm/block_span.hpp index c78033b79..ca889978b 100644 --- a/src/tamm/block_span.hpp +++ b/src/tamm/block_span.hpp @@ -1,51 +1,173 @@ #pragma once -#include "tamm/scalar.hpp" +// Non-owning view over a contiguous tensor block with runtime extents. +// +// NOTE on std::mdspan: TAMM blocks have a *runtime* rank. std::mdspan (and +// the kokkos reference implementation) fix the rank at compile time — there is +// no standard "dynamic-rank" mdspan (std::dextents is +// ill-formed: it attempts to build an extents object of rank 2^64-1). All +// TAMM block kernels treat a block as a flat, contiguous buffer plus an +// extents list (see block_assign_plan.hpp, blockops_*.hpp), so we model +// BlockSpan directly as pointer + extents and expose a std::span for flat +// element access. Kernel code that knows the rank at compile time can build a +// fixed-rank std::mdspan directly from buf() + block_dims() where beneficial. + +#include +#include +#include +#include +#include + +#include "tamm/errors.hpp" #include "tamm/types.hpp" namespace tamm { +/** + * @brief Non-owning view over a contiguous block of T with runtime-determined + * multi-dimensional extents. + * + * The block is stored contiguously in row-major (layout_right) order. Element + * access is available flat (operator[](i)) or as a std::span (flat_span()). + * + * C++20 features used: + * - std::span for the flat view (flat_span()) and the extents span accessor. + * - std::ranges::contiguous_range constraint on the factory helper. + * - [[nodiscard]] on all pure query accessors. + * + * @tparam T Element type of the block + */ template class BlockSpan { public: - enum class BufKind { cpu, invalid }; - // BlockSpan(const Tensor& tensor, const IndexVector& blockid, T* buf) - // : buf_{buf} { - // EXPECTS(buf != nullptr); - // block_dims_ = tensor.block_dims(blockid); - // num_elements_ = 1; - // for (const auto& bd : block_dims) { - // num_elements_ *= bd; - // } - // } - - BlockSpan(T* buf, const std::vector& block_dims): - buf_kind_{BufKind::cpu}, buf_{buf}, block_dims_{block_dims} { + enum class BufKind : uint8_t { cpu, invalid }; + + // ------------------------------------------------------------------ + // Constructors + // ------------------------------------------------------------------ + + /// Default constructor: produces an invalid/null span. + BlockSpan() noexcept : buf_kind_{BufKind::invalid}, buf_{nullptr}, num_elements_{0} {} + + /** + * @brief Primary constructor: pointer + runtime dimension list. + * + * Extents are cached once (single allocation) so block_dims() can return a + * const reference without re-allocating on every call. + * + * @param[in] buf Pointer to the data buffer (must not be null) + * @param[in] dims Extents of each dimension as a contiguous span + * + * @pre buf != nullptr + */ + BlockSpan(T* buf, std::span dims) + : buf_kind_{BufKind::cpu}, buf_{buf}, extents_(dims.begin(), dims.end()) { EXPECTS(buf != nullptr); num_elements_ = 1; - for(const auto& bd: block_dims) { num_elements_ *= bd; } + for(size_t d: extents_) num_elements_ *= d; } - BlockSpan(): buf_kind_{BufKind::invalid}, buf_{nullptr}, num_elements_{0} {} + /// Convenience constructor: accept an std::vector directly. + BlockSpan(T* buf, const std::vector& block_dims) + : BlockSpan(buf, std::span{block_dims}) {} + + BlockSpan(const BlockSpan&) = default; + BlockSpan(BlockSpan&&) = default; + ~BlockSpan() = default; + BlockSpan& operator=(const BlockSpan&) = default; + BlockSpan& operator=(BlockSpan&&) = default; - BlockSpan(const BlockSpan&) = default; - BlockSpan(BlockSpan&&) = default; - ~BlockSpan() = default; - BlockSpan& operator=(const BlockSpan&) = default; - BlockSpan& operator=(BlockSpan&&) = default; + // ------------------------------------------------------------------ + // Accessors + // ------------------------------------------------------------------ - const std::vector& block_dims() const { return block_dims_; } + /// Raw data pointer (mutable). + [[nodiscard]] T* buf() noexcept { return buf_; } + /// Raw data pointer (const). + [[nodiscard]] const T* buf() const noexcept { return buf_; } + + /// Alias for buf() — matches std::span / std::vector naming. + [[nodiscard]] T* data() noexcept { return buf_; } + [[nodiscard]] const T* data() const noexcept { return buf_; } + + /// Flat element count (product of all extents). + [[nodiscard]] size_t num_elements() const noexcept { return num_elements_; } + + /// Number of dimensions. + [[nodiscard]] size_t rank() const noexcept { return extents_.size(); } + + /// Extent along dimension d. + [[nodiscard]] size_t extent(size_t d) const { + EXPECTS(d < extents_.size()); + return extents_[d]; + } - T* buf() { return buf_; } + /** + * @brief Flat (linear) element access into the contiguous block. + * + * The block is stored contiguously (row-major), so flat indexing is + * well-defined for any rank. This is the access pattern used by the + * element-wise block-multiply kernels (scalar/vector Hadamard paths), + * e.g. lhs[0], lhs_vec[i]. + * + * @param i Flat element index in [0, num_elements()). + */ + [[nodiscard]] T& operator[](size_t i) noexcept { return buf_[i]; } + [[nodiscard]] const T& operator[](size_t i) const noexcept { return buf_[i]; } - const T* buf() const { return buf_; } + /// Flat contiguous view over the whole block. + [[nodiscard]] std::span flat_span() noexcept { return {buf_, num_elements_}; } + [[nodiscard]] std::span flat_span() const noexcept { return {buf_, num_elements_}; } - size_t num_elements() const { return num_elements_; } + /** + * @brief Extents of the block, one per dimension. + * + * Returns a const reference to the cached extents vector so it can be passed + * directly to the block kernels (index_permute_assign, ipgen_assign, ...) + * which take `const std::vector&`. No allocation per call. + * + * @return const std::vector& of length rank() + */ + [[nodiscard]] const std::vector& block_dims() const noexcept { return extents_; } + + /// Zero-allocation span view over the extents (for range algorithms). + [[nodiscard]] std::span block_dims_span() const noexcept { + return std::span{extents_.data(), extents_.size()}; + } + + /// Backward-compatible copy of the extents. + [[nodiscard]] std::vector block_dims_vec() const { return extents_; } + + /// True when the span holds a valid (non-null) buffer. + [[nodiscard]] bool is_valid() const noexcept { + return buf_kind_ != BufKind::invalid && buf_ != nullptr; + } private: BufKind buf_kind_; T* buf_; - std::vector block_dims_; - size_t num_elements_; -}; // class BlockSpan + size_t num_elements_{0}; + std::vector extents_; ///< cached extents (allocated once at construction) +}; + +// --------------------------------------------------------------------------- +// Factory helpers +// --------------------------------------------------------------------------- + +/// Build a BlockSpan from a raw pointer + initializer-list of dimensions. +template +[[nodiscard]] inline BlockSpan +make_block_span(T* buf, std::initializer_list dims) { + return BlockSpan{buf, std::vector{dims}}; +} + +/// Build a BlockSpan from a raw pointer + contiguous range of dimensions. +template + requires std::is_convertible_v, size_t> +[[nodiscard]] inline BlockSpan +make_block_span(T* buf, const R& dims) { + std::vector dv(std::ranges::begin(dims), std::ranges::end(dims)); + return BlockSpan{buf, std::span{dv}}; +} + } // namespace tamm diff --git a/src/tamm/blockops_blas.hpp b/src/tamm/blockops_blas.hpp index f60cc00ab..7b81068b1 100644 --- a/src/tamm/blockops_blas.hpp +++ b/src/tamm/blockops_blas.hpp @@ -18,10 +18,8 @@ prep_buffers(T1 lscale, BlockSpan& lhs, T1 rscale, BlockSpan& rhs1, Bloc auto adims = rhs1.block_dims(); auto bdims = rhs2.block_dims(); auto cdims = lhs.block_dims(); - const size_t asize = - std::accumulate(adims.begin(), adims.end(), (size_t) 1, std::multiplies()); - const size_t bsize = - std::accumulate(bdims.begin(), bdims.end(), (size_t) 1, std::multiplies()); + const size_t asize = std::reduce(adims.begin(), adims.end(), size_t{1}, std::multiplies<>{}); + const size_t bsize = std::reduce(bdims.begin(), bdims.end(), size_t{1}, std::multiplies<>{}); // const size_t csize = // std::accumulate(cdims.begin(), cdims.end(), (size_t) 1, std::multiplies()); diff --git a/src/tamm/blockops_cpu.hpp b/src/tamm/blockops_cpu.hpp index 89928665d..0415902cc 100644 --- a/src/tamm/blockops_cpu.hpp +++ b/src/tamm/blockops_cpu.hpp @@ -1,13 +1,16 @@ #pragma once +#include #include +#include // std::ssize #include +#include #include #include "tamm/block_span.hpp" #include "tamm/iteration.hpp" #include "tamm/perm.hpp" -//#include "tamm/scalar.hpp" +#include "tamm/scalar.hpp" #include "tamm/types.hpp" namespace tamm::blockops::cpu { @@ -20,9 +23,8 @@ namespace tamm::blockops::cpu { template void flat_set(BlockSpan& lhs, const T2& value_) { - auto buf = lhs.buf(); - size_t num_elements = lhs.num_elements(); - for(size_t i = 0; i < num_elements; ++i) { buf[i] = value_; } + auto* buf = lhs.buf(); + std::fill(buf, buf + lhs.num_elements(), static_cast(value_)); } template @@ -176,7 +178,7 @@ void flat_assign(BlockSpan& lhs, const BlockSpan& rhs) { TL* lbuf = lhs.buf(); const TR* rbuf = rhs.buf(); const size_t num_elements = lhs.num_elements(); - for(int i = 0; i < num_elements; i++) { *lbuf++ = *rbuf++; } + std::copy(rbuf, rbuf + num_elements, lbuf); } template @@ -252,8 +254,8 @@ void flat_update(const Scalar& lscale, BlockSpan& lhs, const Scalar& rscale, template void flat_lambda(Func&& func, BlockSpan& block, BlockSpans&&... rest) { EXPECTS(((block.buf() != nullptr) && ... && (rest.buf() != nullptr))); - int num_elements = block.tensor().block_size(block.blockid()); - for(int i = 0; i < num_elements; i++) { + const size_t num_elements = block.tensor().block_size(block.blockid()); + for(size_t i = 0; i < num_elements; i++) { std::forward(func)(block.buf()[i], (std::forward(rest).buf()[i])...); } } @@ -264,9 +266,11 @@ void flat_lambda(Func&& func, BlockSpan& block, BlockSpans&&... rest) { // /////////////////////////////////////////////////////////////////////////////// -inline size_t idx(int n, const size_t* id, const std::vector& ldims, const PermVector& p) { - size_t idx = 0; - for(int i = 0; i < n - 1; i++) { idx = (idx + id[p[i]]) * ldims[p[i + 1]]; } +inline size_t idx(std::span id, const std::vector& ldims, + const PermVector& p) { + const size_t n = id.size(); + size_t idx = 0; + for(size_t i = 0; i + 1 < n; i++) { idx = (idx + id[p[i]]) * ldims[p[i + 1]]; } if(n > 0) { idx += id[p[n - 1]]; } return idx; } @@ -285,30 +289,33 @@ void index_permute_assign(TL* lbuf, const TR* rbuf, const PermVector& perm_to_de for(size_t i = 0; i < ldims[0]; i++) { lbuf[i] = rbuf[i]; } } else if(ndim == 2) { - size_t i[2], c; + std::array i{}; + size_t c; for(c = 0, i[0] = 0; i[0] < ldims[0]; i[0]++) { for(i[1] = 0; i[1] < ldims[1]; i[1]++, c++) { - lbuf[c] = rbuf[idx(2, i, ldims, perm_to_dest)]; + lbuf[c] = rbuf[idx(i, ldims, perm_to_dest)]; } } } else if(ndim == 3) { - size_t i[3], c; + std::array i{}; + size_t c; for(c = 0, i[0] = 0; i[0] < ldims[0]; i[0]++) { for(i[1] = 0; i[1] < ldims[1]; i[1]++) { for(i[2] = 0; i[2] < ldims[2]; i[2]++, c++) { - lbuf[c] = rbuf[idx(3, i, ldims, perm_to_dest)]; + lbuf[c] = rbuf[idx(i, ldims, perm_to_dest)]; } } } } else if(ndim == 4) { - size_t i[4], c; + std::array i{}; + size_t c; for(c = 0, i[0] = 0; i[0] < ldims[0]; i[0]++) { for(i[1] = 0; i[1] < ldims[1]; i[1]++) { for(i[2] = 0; i[2] < ldims[2]; i[2]++) { for(i[3] = 0; i[3] < ldims[3]; i[3]++, c++) { - lbuf[c] = rbuf[idx(4, i, ldims, perm_to_dest)]; + lbuf[c] = rbuf[idx(i, ldims, perm_to_dest)]; } } } @@ -331,30 +338,33 @@ void index_permute_assign(TL* lbuf, TL rscale, const TR* rbuf, const PermVector& for(size_t i = 0; i < ldims[0]; i++) { lbuf[i] = rscale * rbuf[i]; } } else if(ndim == 2) { - size_t i[2], c; + std::array i{}; + size_t c; for(c = 0, i[0] = 0; i[0] < ldims[0]; i[0]++) { for(i[1] = 0; i[1] < ldims[1]; i[1]++, c++) { - lbuf[c] = rscale * rbuf[idx(2, i, ldims, perm_to_dest)]; + lbuf[c] = rscale * rbuf[idx(i, ldims, perm_to_dest)]; } } } else if(ndim == 3) { - size_t i[3], c; + std::array i{}; + size_t c; for(c = 0, i[0] = 0; i[0] < ldims[0]; i[0]++) { for(i[1] = 0; i[1] < ldims[1]; i[1]++) { for(i[2] = 0; i[2] < ldims[2]; i[2]++, c++) { - lbuf[c] = rscale * rbuf[idx(3, i, ldims, perm_to_dest)]; + lbuf[c] = rscale * rbuf[idx(i, ldims, perm_to_dest)]; } } } } else if(ndim == 4) { - size_t i[4], c; + std::array i{}; + size_t c; for(c = 0, i[0] = 0; i[0] < ldims[0]; i[0]++) { for(i[1] = 0; i[1] < ldims[1]; i[1]++) { for(i[2] = 0; i[2] < ldims[2]; i[2]++) { for(i[3] = 0; i[3] < ldims[3]; i[3]++, c++) { - lbuf[c] = rscale * rbuf[idx(4, i, ldims, perm_to_dest)]; + lbuf[c] = rscale * rbuf[idx(i, ldims, perm_to_dest)]; } } } @@ -377,30 +387,33 @@ void index_permute_update(TL* lbuf, const TR* rbuf, const PermVector& perm_to_de for(size_t i = 0; i < ldims[0]; i++) { lbuf[i] += rbuf[i]; } } else if(ndim == 2) { - size_t i[2], c; + std::array i{}; + size_t c; for(c = 0, i[0] = 0; i[0] < ldims[0]; i[0]++) { for(i[1] = 0; i[1] < ldims[1]; i[1]++, c++) { - lbuf[c] += rbuf[idx(2, i, ldims, perm_to_dest)]; + lbuf[c] += rbuf[idx(i, ldims, perm_to_dest)]; } } } else if(ndim == 3) { - size_t i[3], c; + std::array i{}; + size_t c; for(c = 0, i[0] = 0; i[0] < ldims[0]; i[0]++) { for(i[1] = 0; i[1] < ldims[1]; i[1]++) { for(i[2] = 0; i[2] < ldims[2]; i[2]++, c++) { - lbuf[c] += rbuf[idx(3, i, ldims, perm_to_dest)]; + lbuf[c] += rbuf[idx(i, ldims, perm_to_dest)]; } } } } else if(ndim == 4) { - size_t i[4], c; + std::array i{}; + size_t c; for(c = 0, i[0] = 0; i[0] < ldims[0]; i[0]++) { for(i[1] = 0; i[1] < ldims[1]; i[1]++) { for(i[2] = 0; i[2] < ldims[2]; i[2]++) { for(i[3] = 0; i[3] < ldims[3]; i[3]++, c++) { - lbuf[c] += rbuf[idx(4, i, ldims, perm_to_dest)]; + lbuf[c] += rbuf[idx(i, ldims, perm_to_dest)]; } } } @@ -423,30 +436,33 @@ void index_permute_update(TL* lbuf, TL rscale, const TR* rbuf, const PermVector& for(size_t i = 0; i < ldims[0]; i++) { lbuf[i] += rscale * rbuf[i]; } } else if(ndim == 2) { - size_t i[2], c; + std::array i{}; + size_t c; for(c = 0, i[0] = 0; i[0] < ldims[0]; i[0]++) { for(i[1] = 0; i[1] < ldims[1]; i[1]++, c++) { - lbuf[c] += rscale * rbuf[idx(2, i, ldims, perm_to_dest)]; + lbuf[c] += rscale * rbuf[idx(i, ldims, perm_to_dest)]; } } } else if(ndim == 3) { - size_t i[3], c; + std::array i{}; + size_t c; for(c = 0, i[0] = 0; i[0] < ldims[0]; i[0]++) { for(i[1] = 0; i[1] < ldims[1]; i[1]++) { for(i[2] = 0; i[2] < ldims[2]; i[2]++, c++) { - lbuf[c] += rscale * rbuf[idx(3, i, ldims, perm_to_dest)]; + lbuf[c] += rscale * rbuf[idx(i, ldims, perm_to_dest)]; } } } } else if(ndim == 4) { - size_t i[4], c; + std::array i{}; + size_t c; for(c = 0, i[0] = 0; i[0] < ldims[0]; i[0]++) { for(i[1] = 0; i[1] < ldims[1]; i[1]++) { for(i[2] = 0; i[2] < ldims[2]; i[2]++) { for(i[3] = 0; i[3] < ldims[3]; i[3]++, c++) { - lbuf[c] += rscale * rbuf[idx(4, i, ldims, perm_to_dest)]; + lbuf[c] += rscale * rbuf[idx(i, ldims, perm_to_dest)]; } } } @@ -469,30 +485,33 @@ void index_permute_update(TL lscale, TL* lbuf, TL rscale, const TR* rbuf, for(size_t i = 0; i < ldims[0]; i++) { lbuf[i] = lscale * lbuf[i] + rscale * rbuf[i]; } } else if(ndim == 2) { - size_t i[2], c; + std::array i{}; + size_t c; for(c = 0, i[0] = 0; i[0] < ldims[0]; i[0]++) { for(i[1] = 0; i[1] < ldims[1]; i[1]++, c++) { - lbuf[c] = lscale * lbuf[i] + rscale * rbuf[idx(2, i, ldims, perm_to_dest)]; + lbuf[c] = lscale * lbuf[c] + rscale * rbuf[idx(i, ldims, perm_to_dest)]; } } } else if(ndim == 3) { - size_t i[3], c; + std::array i{}; + size_t c; for(c = 0, i[0] = 0; i[0] < ldims[0]; i[0]++) { for(i[1] = 0; i[1] < ldims[1]; i[1]++) { for(i[2] = 0; i[2] < ldims[2]; i[2]++, c++) { - lbuf[c] = lscale * lbuf[i] + rscale * rbuf[idx(3, i, ldims, perm_to_dest)]; + lbuf[c] = lscale * lbuf[c] + rscale * rbuf[idx(i, ldims, perm_to_dest)]; } } } } else if(ndim == 4) { - size_t i[4], c; + std::array i{}; + size_t c; for(c = 0, i[0] = 0; i[0] < ldims[0]; i[0]++) { for(i[1] = 0; i[1] < ldims[1]; i[1]++) { for(i[2] = 0; i[2] < ldims[2]; i[2]++) { for(i[3] = 0; i[3] < ldims[3]; i[3]++, c++) { - lbuf[c] = lscale * lbuf[i] + rscale * rbuf[idx(4, i, ldims, perm_to_dest)]; + lbuf[c] = lscale * lbuf[c] + rscale * rbuf[idx(i, ldims, perm_to_dest)]; } } } @@ -510,7 +529,10 @@ void index_permute_update(TL lscale, TL* lbuf, TL rscale, const TR* rbuf, inline size_t ipgen_idx(const std::vector& index_vec, const std::vector& dims_vec) { size_t ret = 0, ld = 1; EXPECTS(index_vec.size() == dims_vec.size()); - for(int i = index_vec.size(); i >= 0; i--) { + // Iterate the valid range [size-1 .. 0]. The previous start value of + // index_vec.size() read index_vec[size()] / dims_vec[size()] out of bounds + // on the first iteration. + for(std::ptrdiff_t i = std::ssize(index_vec) - 1; i >= 0; i--) { ret += ld * index_vec[i]; ld *= dims_vec[i]; } diff --git a/src/tamm/boundvec.hpp b/src/tamm/boundvec.hpp index b51e2311e..567d41869 100644 --- a/src/tamm/boundvec.hpp +++ b/src/tamm/boundvec.hpp @@ -1,10 +1,10 @@ -// Copyright 2016 Pacific Northwest National Laboratory - #pragma once #include -#include -#include +#include +#include +#include +#include #include "tamm/errors.hpp" @@ -13,295 +13,304 @@ namespace tamm { /** * @brief Vector of bounded length. * - * This class provides a vector-like interface while having bounded size like an - * array. Bounds checks can be added in debug mode. + * A vector-like container whose maximum number of elements is fixed at + * compile time as @p maxsize. Storage is entirely on the stack (a + * @c std::array), so no heap allocation ever occurs. * - * @tparam T Type of each element contained in the vector - * @tparam maxsize Maximum size of the vector + * @tparam T Element type. + * @tparam maxsize Compile-time upper bound on the number of elements. * - * @todo Add bounds checking when BOUNDVEC_DEBUG option is enabled + * @todo Consider replacing with a fully-standard `std::inplace_vector` + * once C++26 support is available. */ -template -class BoundVec: public std::array { +template +class BoundVec { public: - using size_type = typename std::array::size_type; - using value_type = typename std::array::value_type; + // ----------------------------------------------------------------------- + // Type aliases (mirror std::vector interface) + // ----------------------------------------------------------------------- + using value_type = T; + using size_type = std::size_t; + using reference = T&; + using const_reference = const T&; + using pointer = T*; + using const_pointer = const T*; using iterator = typename std::array::iterator; using const_iterator = typename std::array::const_iterator; using reverse_iterator = typename std::array::reverse_iterator; using const_reverse_iterator = typename std::array::const_reverse_iterator; - using reference = typename std::array::reference; - using const_reference = typename std::array::const_reference; - using std::array::begin; - using std::array::rend; + // ----------------------------------------------------------------------- + // Constructors + // ----------------------------------------------------------------------- + + /** @brief Default constructor — produces an empty vector. */ + BoundVec() = default; /** - * @brief Constructor a zero-sized vector + * @brief Construct @p count copies of @p value. + * @param count Number of elements to create. + * @param value Value to copy into each element (defaults to T{}). + * @pre count <= maxsize */ - BoundVec() noexcept: size_{0} {} + explicit BoundVec(size_type count, const T& value = T{}) { + EXPECTS(count <= maxsize); + for(size_type i = 0; i < count; ++i) { data_[size_++] = value; } + } /** - * @brief Construct a vector of specified size, with an optional initial value - * - * @param[in] count size of constructor vector - * @param[in] value initial value of all elements in the constructed vector + * @brief Construct from an initializer list. + * @param list Elements to copy into the vector. + * @pre list.size() <= maxsize */ - explicit BoundVec(size_type count, const T& value = T()): size_{0} { - for(size_type i = 0; i < count; i++) { push_back(value); } + BoundVec(std::initializer_list list) { + EXPECTS(list.size() <= maxsize); + for(const auto& v: list) { data_[size_++] = v; } } - BoundVec(const BoundVec&) = default; - BoundVec(BoundVec&&) = default; - BoundVec& operator=(const BoundVec&) = default; - BoundVec& operator=(BoundVec&&) = default; - ~BoundVec() { - for(auto& v: *this) { v.~value_type(); } + /** + * @brief Construct from a pair of iterators. + * @tparam Iter Input iterator type. + * @param first Beginning of the source range. + * @param last One-past-end of the source range. + * @pre std::distance(first, last) <= maxsize + * + * Constrained to std::input_iterator so it does not compete with the + * (count, value) constructor for calls like BoundVec(n, val). + */ + template + BoundVec(Iter first, Iter last) { + for(; first != last; ++first) { push_back(*first); } } - template - BoundVec(Itr first, Itr last): size_{0} { - for(auto itr = first; itr != last; ++itr) { push_back(*itr); } - } + // ----------------------------------------------------------------------- + // Capacity + // ----------------------------------------------------------------------- - BoundVec(std::initializer_list init): size_{0} { - for(auto v: init) { push_back(v); } - } + /** @brief Return the number of live elements. */ + [[nodiscard]] size_type size() const noexcept { return size_; } - /** - * @brief Size of this vector - * - * @return Size of the vector - */ - constexpr size_type size() const noexcept { return size_; } + /** @brief Return the compile-time maximum number of elements. */ + [[nodiscard]] static constexpr size_type max_size() noexcept { return maxsize; } - /** - * @brief Maximum number of elements this vector can hold - * - * @return Maximum size of this vector - */ - constexpr size_type max_size() const noexcept { return maxsize; } + /** @brief Return true iff the vector holds no elements. */ + [[nodiscard]] bool empty() const noexcept { return size_ == 0; } + + // ----------------------------------------------------------------------- + // Modifiers + // ----------------------------------------------------------------------- + + /** @brief Remove all elements (does not release storage). */ + void clear() noexcept { size_ = 0; } /** - * @brief Is this vector empty - * - * @return True if this vector is empty, false otherwise. + * @brief Append a copy of @p val. + * @param val Value to append. + * @pre size() < maxsize */ - constexpr bool empty() const noexcept { return size_ == 0; } + void push_back(const T& val) { + EXPECTS(size_ < maxsize); + data_[size_++] = val; + } /** - * @brief Clear the contents of this vector + * @brief Append by moving @p val. + * @param val Value to move-append. + * @pre size() < maxsize */ - void clear() noexcept { - for(auto& v: *this) { v.~value_type(); } - size_ = 0; + void push_back(T&& val) { + EXPECTS(size_ < maxsize); + data_[size_++] = std::move(val); } /** - * @brief Push one element to the back of the vector - * - * @param[in] value The value to be pushed - * - * @pre size() < max_size + * @brief Remove the last element. + * @pre !empty() */ - void push_back(const T& value) { - EXPECTS(size() < maxsize); - this->at(size_++) = value; + void pop_back() { + EXPECTS(size_ > 0); + --size_; } /** - * @brief Push one element to the back of the vector + * @brief Resize the vector to @p sz elements. * - * @param[in,out] value The value to be pushed + * If @p sz > size(), new elements are value-initialised. If @p sz < + * size(), excess elements are discarded. * - * @pre size() < max_size + * @param sz New size. + * @pre sz <= maxsize */ - void push_back(T&& value) { - EXPECTS(size() < maxsize); - this->at(size_++) = std::move(value); + void resize(size_type sz) { + // Fix: sz is unsigned so "sz >= 0" is always true and never fires. + // The meaningful guard is an upper-bound check. + EXPECTS(sz <= maxsize); + if(sz > size_) { + for(size_type i = size_; i < sz; ++i) { data_[i] = T{}; } + } + size_ = sz; } /** - * @brief Remove one element from the back of the vector - * - * @pre size() > 0 + * @brief Insert a copy of @p val at the back (alias for push_back). + * @param val Value to insert. + * @pre size() < maxsize + * @return Iterator to the inserted element. */ - void pop_back() noexcept { - EXPECTS(size() > 0); - back().~value_type(); - size_ -= 1; + iterator insert_back(const T& val) { + push_back(val); + return end() - 1; } /** - * @brief Resize vector to desired size - * - * @pre size >=0 - * - * @param[in] size Size desired for the vector + * @brief Insert by moving @p val at the back (alias for push_back). + * @param val Value to move-insert. + * @pre size() < maxsize + * @return Iterator to the inserted element. */ - void resize(size_type size) { - EXPECTS(size >= 0); - EXPECTS(size < maxsize); - for(size_t i = size; i < size_; i++) { this->at(i).~value_type(); } - for(size_t i = size_; i < size; i++) { this->at(i) = value_type(); } - size_ = size; + iterator insert_back(T&& val) { + push_back(std::move(val)); + return end() - 1; } /** - * @brief Insert a sequence of elements, specified using iterators, at the back of the vector - * - * @param[in] first Starting iterator position for tasks to be inserted - * @param[in] last Ending iterator position for the tasks to be inserted - * - * @pre size() + std::distance(@param first, @param last) <= maxsize + * @brief Append the range [first, last) at the back. + * @tparam InputIt Input iterator type. + * @param first Beginning of the source range. + * @param last One-past-end of the source range. + * @pre size() + distance(first,last) <= maxsize */ - template + template void insert_back(InputIt first, InputIt last) { - EXPECTS(size_ + (last - first) <= maxsize); - for(auto itr = first; itr != last; ++itr) { push_back(*itr); } + for(; first != last; ++first) { push_back(*first); } } /** - * @brief Insert a given value multiple times at the back of the vector - * - * @param[in] count Number of times the given value is to be inserted - * @param[in] value Value to be inserted - * - * @pre size() + @param count <= maxsize + * @brief Append @p count copies of @p value at the back. + * @param count Number of copies to append. + * @param value Value to copy. + * @pre size() + count <= maxsize */ void insert_back(size_type count, const T& value) { - EXPECTS(size_ + count <= maxsize); - for(size_type i = 0; i < count; i++) { push_back(value); } + for(size_type i = 0; i < count; ++i) { push_back(value); } } - // BoundVec& operator = (BoundVec& bvec) { - // size_ = bvec.size_; - // std::copy(bvec.begin(), bvec.end(), begin()); - // return *this; - // } - - // BoundVec& operator=(const BoundVec& bvec) { - // size_ = bvec.size_; - // std::copy(bvec.begin(), bvec.end(), begin()); - // return *this; - // } - - // BoundVec& operator=(BoundVec&& bvec) { - // clear(); - // for(auto&& bv : bvec) { push_back(std::move(bv)); } - // return *this; - // } - - /** - * @brief Obtain end iterator past the last element in the vector - * - * @return The end iterator - */ - iterator end() noexcept { return std::array::begin() + size_; } - - /** - * @brief Obtain end iterator past the last element in the vector - * - * @return Const end iterator - */ - const_iterator end() const noexcept { return std::array::begin() + size_; } - - // reverse_iterator rbegin() const { - // return std::array::begin() + size_; - // } + // ----------------------------------------------------------------------- + // Iterators + // ----------------------------------------------------------------------- + + /** @brief Return an iterator to the first live element. */ + iterator begin() noexcept { return data_.begin(); } + /** @brief Return a const iterator to the first live element. */ + const_iterator begin() const noexcept { return data_.begin(); } + + /** @brief Return an iterator one past the last live element. */ + iterator end() noexcept { return data_.begin() + size_; } + /** @brief Return a const iterator one past the last live element. */ + const_iterator end() const noexcept { return data_.begin() + size_; } + + /** @brief Const iterator to the first live element. */ + const_iterator cbegin() const noexcept { return data_.begin(); } + /** @brief Const iterator one past the last live element. */ + const_iterator cend() const noexcept { return data_.begin() + size_; } + + /** @brief Reverse iterator to the last live element. */ + reverse_iterator rbegin() noexcept { return reverse_iterator{end()}; } + const_reverse_iterator rbegin() const noexcept { return const_reverse_iterator{end()}; } + /** @brief Reverse iterator one before the first live element. */ + reverse_iterator rend() noexcept { return reverse_iterator{begin()}; } + const_reverse_iterator rend() const noexcept { return const_reverse_iterator{begin()}; } + + // ----------------------------------------------------------------------- + // Raw storage access + // ----------------------------------------------------------------------- + + /** @brief Pointer to the underlying contiguous storage (mutable). */ + [[nodiscard]] pointer data() noexcept { return data_.data(); } + /** @brief Pointer to the underlying contiguous storage (const). */ + [[nodiscard]] const_pointer data() const noexcept { return data_.data(); } + + // ----------------------------------------------------------------------- + // Element access + // ----------------------------------------------------------------------- /** - * Obtain reference to first element in the vector - * - * @pre size() > 0 - * - * @return Reference to first element + * @brief Unchecked element access. + * @param i Index (0-based). + * @pre i < size() */ - reference front() noexcept { - EXPECTS(size() > 0); - return this->at(0); + reference operator[](size_type i) { + EXPECTS(i < size_); + return data_[i]; } - - /** - * Obtain a const reference to first element in the vector - * - * @pre size() > 0 - * - * @return Reference to first element - */ - const_reference front() const noexcept { - EXPECTS(size() > 0); - return this->at(0); + /** @brief Unchecked const element access. */ + const_reference operator[](size_type i) const { + EXPECTS(i < size_); + return data_[i]; } /** - * Obtain reference to the last element in the vector - * - * @pre size() > 0 - * - * @return Reference to last element + * @brief Return a reference to the first element. + * @pre !empty() */ - reference back() noexcept { - EXPECTS(size() > 0); - return this->at(size_ - 1); + reference front() { + EXPECTS(size_ > 0); + return data_.front(); } - - /** - * Obtain a const reference to the last element in the vector - * - * @pre size() > 0 - * - * @return Const reference to last element - */ - const_reference back() const noexcept { - EXPECTS(size() > 0); - return this->at(size_ - 1); + /** @brief Return a const reference to the first element. */ + const_reference front() const { + EXPECTS(size_ > 0); + return data_.front(); } -private: /** - * @brief Size of the vector + * @brief Return a reference to the last live element. + * @pre !empty() */ - size_type size_; - - /** - * @brief Equality operator to compare two vectors - * - * @param[in] lhs One vector to be compared - * - * @param[in] rhs The other vector to be compared - */ - friend bool operator==(const BoundVec& lhs, const BoundVec& rhs) { - return lhs.size() == rhs.size() && std::equal(lhs.begin(), lhs.end(), rhs.begin(), rhs.end()); + reference back() { + EXPECTS(size_ > 0); + return data_[size_ - 1]; + } + /** @brief Return a const reference to the last live element. */ + const_reference back() const { + EXPECTS(size_ > 0); + return data_[size_ - 1]; } - /** - * @brief Inequality operator to compare two vectors - * - * @param[in] lhs One vector to be compared - * - * @param[in] rhs The other vector to be compared - */ - friend bool operator!=(const BoundVec& lhs, const BoundVec& rhs) { - return !(lhs == rhs); + // ----------------------------------------------------------------------- + // Comparison + // ----------------------------------------------------------------------- + + /** @brief Equality operator — two BoundVecs are equal iff their live + * elements are pairwise equal. */ + bool operator==(const BoundVec& rhs) const noexcept { + if(size_ != rhs.size_) { return false; } + for(size_type i = 0; i < size_; i++) { + if(data_[i] != rhs.data_[i]) { return false; } + } + return true; } -}; // class BoundVec + +private: + /** @brief Size of the vector */ + size_type size_{0}; + std::array data_{}; +}; /** - * @brief Dump a vector to an output stream - * - * @param[in,out] os Output stream to write to - * - * @param[in] bvec The vector to be output - * - * @return The modified output stream @param os + * @brief Dump a vector to an output stream (diagnostic helper). + * @param os Target output stream. + * @param bv Vector to print. + * @return Reference to @p os. */ -template -inline std::ostream& operator<<(std::ostream& os, const BoundVec& bvec) { - os << std::string{"[ "}; - for(auto el: bvec) { os << el << " "; } - os << std::string{"]"}; - return os; +template +std::ostream& operator<<(std::ostream& os, const BoundVec& bv) { + os << '['; + for(std::size_t i = 0; i < bv.size(); ++i) { + if(i) os << ", "; + os << bv[i]; + } + return os << ']'; } } // namespace tamm diff --git a/src/tamm/deallocop.hpp b/src/tamm/deallocop.hpp index c00bfea28..a04d90310 100644 --- a/src/tamm/deallocop.hpp +++ b/src/tamm/deallocop.hpp @@ -19,7 +19,6 @@ class DeallocOp: public Op { public: DeallocOp(TensorType tensor): tensor_{tensor} {} - DeallocOp(const DeallocOp&) = default; TensorType tensor() const { return tensor_; } @@ -27,7 +26,9 @@ class DeallocOp: public Op { OpType op_type() const override { return OpType::dealloc; } - std::shared_ptr clone() const override { return std::shared_ptr(new DeallocOp{*this}); } + std::shared_ptr clone() const override { + return std::make_shared(*this); + } void execute(ExecutionContext& ec, ExecutionHW hw = ExecutionHW::CPU) override { tensor_.deallocate(); diff --git a/src/tamm/distribution.hpp b/src/tamm/distribution.hpp index 6852b8a95..71aadf30d 100644 --- a/src/tamm/distribution.hpp +++ b/src/tamm/distribution.hpp @@ -127,13 +127,13 @@ class Distribution { void set_ga_handle(int ga_handle) { ga_ = ga_handle; } - void set_proc_grid(std::vector pg) { proc_grid_ = pg; } + void set_proc_grid(std::vector pg) { proc_grid_ = std::move(pg); } void set_proc_buf_size(Size proc_buf_size) { proc_buf_size_ = proc_buf_size; } void set_max_proc_buf_size(Size max_proc_buf_size) { max_proc_buf_size_ = max_proc_buf_size; } - std::vector proc_grid() const { return proc_grid_; } + [[nodiscard]] const std::vector& proc_grid() const { return proc_grid_; } /** * @brief Construct a new Distribution object using a TensorBase object and diff --git a/src/tamm/errors.hpp b/src/tamm/errors.hpp index 9dcc6123f..1735fe719 100644 --- a/src/tamm/errors.hpp +++ b/src/tamm/errors.hpp @@ -1,70 +1,114 @@ #pragma once +// C++20: std::source_location replaces __FILE__/__LINE__ macro patchwork. +// Macro wrappers are retained for ABI / call-site compatibility. + #include #include +#include // C++20 +#include +#include +#include namespace tamm { +// --------------------------------------------------------------------------- +// Core diagnostic helpers (inline, no macros) +// --------------------------------------------------------------------------- + +/// Throw a descriptive runtime_error with full source location. +[[noreturn]] inline void +tamm_error(std::string_view msg, + std::source_location loc = std::source_location::current()) { + std::string full = std::string{loc.file_name()} + ':' + + std::to_string(loc.line()) + + " in '" + loc.function_name() + "': " + + std::string{msg}; + std::cerr << "TAMM ERROR: " << full << '\n'; + throw std::runtime_error(full); +} + +/// Precondition check: throws if cond is false. +inline void +tamm_expects(bool cond, std::string_view msg, + std::source_location loc = std::source_location::current()) { + if (!cond) [[unlikely]] { + std::string full = std::string{loc.file_name()} + ':' + + std::to_string(loc.line()) + + " in '" + loc.function_name() + "': EXPECTS failed [" + + std::string{msg} + ']'; + std::cerr << full << '\n'; + throw std::runtime_error(full); + } +} + +/// Precondition check with custom message. +inline void +tamm_expects_str(bool cond, std::string_view cond_str, std::string_view user_msg, + std::source_location loc = std::source_location::current()) { + if (!cond) [[unlikely]] { + std::string full = std::string{loc.file_name()} + ':' + + std::to_string(loc.line()) + + " in '" + loc.function_name() + "': " + + std::string{cond_str} + " -- " + std::string{user_msg}; + std::cerr << full << '\n'; + throw std::runtime_error(full); + } +} + +// --------------------------------------------------------------------------- +// Macro wrappers (preserved for call-site compatibility) +// source_location is captured automatically at each call site. +// --------------------------------------------------------------------------- + // clang-format off +/** + * @brief Wrapper for assertion checking. + * + * This is meant to identify preconditions and possibly include additional + * operations (e.g., an error message). + */ +// static_cast performs a *contextual* conversion at the call site, so +// pointer-like conditions with an explicit operator bool (std::shared_ptr, +// std::unique_ptr, ...) still work, e.g. EXPECTS(distribution_). +#define EXPECTS(cond) tamm_expects(static_cast(cond), #cond) + +/** + * @brief Wrapper for assertion checking with a custom string message. + */ +#define EXPECTS_STR(cond, str) tamm_expects_str(static_cast(cond), #cond, str) + +#define EXPECTS_NOTHROW(cond) assert(cond) + /** * @brief Mark code options that are not yet implemented. */ #define NOT_IMPLEMENTED() \ - do { \ - std::cerr << "ERROR: (Not implemented)" \ - << "file:" << __FILE__ << "function:" << __func__ \ - << " line:" << __LINE__ << std::endl; \ - throw std::runtime_error(""); \ - } while(0) + tamm_error("Not implemented") /** * @brief Mark code options that are not yet allowed. */ #define NOT_ALLOWED() \ - do { \ - std::cerr << "ERROR: (Not allowed)" \ - << "file:" << __FILE__ << "function:" << __func__ \ - << " line:" << __LINE__ << std::endl; \ - throw std::runtime_error(""); \ - } while(0) - -/** - * @brief Mark code paths that should be unreachable - */ -#define UNREACHABLE() \ - do { \ - std::cerr << "ERROR: (unreachable)" \ - << "file:" << __FILE__ << "function:" << __func__ \ - << " line:" << __LINE__ << std::endl; \ - } while(0) + tamm_error("Not allowed") /** - * @brief Wrapper for assertion checking. + * @brief Mark code paths that should be unreachable. * - * This is meant to identify preconditions and possibly include additional - * operations (e.g., an error message). + * Uses C++23 std::unreachable() when available, otherwise aborts via + * tamm_error + __builtin_unreachable for optimizer hints. */ -#define EXPECTS_NOTHROW(cond) assert(cond) -#define EXPECTS(cond) \ - do { \ - if(!(cond)) { \ - std::cerr << "EXPECTS failed. Condition: " << __FILE__ << ":" \ - << __LINE__ << " " << #cond << "\n"; \ - throw std::string{"EXPECT condition failed: "} + \ - std::string{#cond}; \ - } \ - } while(0) - -#define EXPECTS_STR(cond, str) \ - do { \ - if(!(cond)) { \ - std::cerr << "Assertion failed. Condition: " << __FILE__ << ":" \ - << __LINE__ << " " << #cond << " - " \ - << std::string{#str} <<"\n"; \ - throw std::string{"Error: "} + std::string{#str}; \ - } \ +#if defined(__cpp_lib_unreachable) && __cpp_lib_unreachable >= 202202L +# include +# define UNREACHABLE() (std::unreachable()) +#else +# define UNREACHABLE() \ + do { \ + tamm_error("Reached supposedly unreachable code"); \ + __builtin_unreachable(); \ } while(0) +#endif // clang-format on diff --git a/src/tamm/execution_context.hpp b/src/tamm/execution_context.hpp index 4f29dc491..0af50aeb3 100644 --- a/src/tamm/execution_context.hpp +++ b/src/tamm/execution_context.hpp @@ -356,13 +356,11 @@ class ExecutionContext { std::stringstream& get_profile_data() { return profile_data_; } std::string get_profile_header() { - std::string pheader = - "ID;Level;OpType;OP;total_op_time_min;total_op_time_max;total_op_time_avg;"; - pheader += "get_time_min;get_time_max;get_time_avg;"; - pheader += "block_compute_time_min;block_compute_time_max;block_compute_time_avg;"; - pheader += "copy_time_min;copy_time_max;copy_time_avg;"; - pheader += "acc_time_min;acc_time_max;acc_time_avg"; - return pheader; + return "ID;Level;OpType;OP;total_op_time_min;total_op_time_max;total_op_time_avg;" + "get_time_min;get_time_max;get_time_avg;" + "block_compute_time_min;block_compute_time_max;block_compute_time_avg;" + "copy_time_min;copy_time_max;copy_time_avg;" + "acc_time_min;acc_time_max;acc_time_avg"; } template diff --git a/src/tamm/index_space.cpp b/src/tamm/index_space.cpp index aeaf442d7..e429a87b9 100644 --- a/src/tamm/index_space.cpp +++ b/src/tamm/index_space.cpp @@ -160,4 +160,31 @@ bool operator>=(const IndexSpace& lhs, const IndexSpace& rhs) { return (lhs > rhs) || (lhs == rhs); } +// Defined here (rather than in index_space_interface.hpp) so that the +// std::vector members are only ever destroyed/copied where +// TiledIndexSpace is a complete type. Required since C++20. +namespace internal { +const std::vector& empty_tiled_index_space_vec() { + static const std::vector empty; + return empty; +} +} // namespace internal + +DependentIndexSpaceImpl::DependentIndexSpaceImpl( + const std::vector& indep_spaces, + const std::map& dep_space_relation): + dep_spaces_{indep_spaces}, dep_space_relation_{dep_space_relation}, named_ranges_{} { + max_size_ = 0; + for(const auto& pair: dep_space_relation) { + max_size_ = std::max(max_size_, pair.second.num_indices()); + } +} + +DependentIndexSpaceImpl::DependentIndexSpaceImpl( + const std::vector& indep_spaces, const IndexSpace& ref_space, + const std::map& dep_space_relation): + dep_spaces_{indep_spaces}, dep_space_relation_{dep_space_relation}, named_ranges_{} {} + +size_t DependentIndexSpaceImpl::num_key_tiled_index_spaces() const { return dep_spaces_.size(); } + } // namespace tamm diff --git a/src/tamm/index_space.hpp b/src/tamm/index_space.hpp index b7a231a27..cd934cee7 100644 --- a/src/tamm/index_space.hpp +++ b/src/tamm/index_space.hpp @@ -469,7 +469,8 @@ class IndexSpace { protected: std::shared_ptr impl_; /**< shared pointer to the implementation */ - size_t hash_value_; /**< hash value associated with the IndexSpace */ + size_t hash_value_{0}; /**< hash value associated with the IndexSpace (0 for + a default-constructed / empty IndexSpace) */ }; // class IndexSpace } // namespace tamm diff --git a/src/tamm/index_space_interface.hpp b/src/tamm/index_space_interface.hpp index db4e1243d..2b32c64f1 100644 --- a/src/tamm/index_space_interface.hpp +++ b/src/tamm/index_space_interface.hpp @@ -7,6 +7,19 @@ #include namespace tamm { +namespace internal { +/** + * @brief Returns a reference to a shared, always-empty vector of + * TiledIndexSpace. Defined out-of-line (in index_space.cpp) where + * TiledIndexSpace is a complete type, so that impl classes which need to hand + * out an empty dependency list do not have to store a std::vector of the + * (here) incomplete TiledIndexSpace type. Required since C++20 instantiates + * std::vector's (constexpr) destructor eagerly, which needs a complete element + * type. + */ +const std::vector& empty_tiled_index_space_vec(); +} // namespace internal + /** * @ingroup index_space * @class IndexSpaceInterface @@ -233,10 +246,15 @@ class IndexSpaceInterface { * @tparam ContainerType stl container type with iterator * (RandomAccessIterator) support * @param [in] data_vec input vector - * @returns true returned if there are duplicates + * @returns true returned if there are no duplicate elements + * + * NOTE: name deliberately reads "no_duplicate" — it returns true when the + * input is duplicate-free. (The previous name has_duplicate() was inverted + * relative to its return value; callers use it as a uniqueness precondition + * via EXPECTS(has_no_duplicate(...)).) */ template - static bool has_duplicate(const ContainerType& data_vec) { + static bool has_no_duplicate(const ContainerType& data_vec) { ContainerType temp_vec = data_vec; std::sort(temp_vec.begin(), temp_vec.end()); @@ -268,7 +286,7 @@ class IndexSpaceInterface { } // Check no overlap on the ranges std::sort(att_indices.begin(), att_indices.end()); - EXPECTS(has_duplicate(att_indices)); + EXPECTS(has_no_duplicate(att_indices)); // Check for full coverage of the indices EXPECTS(indices.size() == att_indices.size()); @@ -320,7 +338,7 @@ class RangeIndexSpaceImpl: public IndexSpaceInterface { named_subspaces_{construct_subspaces(named_ranges, spin)}, spin_{construct_spin(spin)}, spatial_{construct_spatial(spatial)} { - EXPECTS(has_duplicate(indices_)); + EXPECTS(has_no_duplicate(indices_)); } /// @todo do we need these copy/move constructor/operators? @@ -353,7 +371,9 @@ class RangeIndexSpaceImpl: public IndexSpaceInterface { // Maximum number of indices in this index space std::size_t max_num_indices() const override { return indices_.size(); } - const std::vector& key_tiled_index_spaces() const override { return empty_vec_; } + const std::vector& key_tiled_index_spaces() const override { + return internal::empty_tiled_index_space_vec(); + } const std::map& map_tiled_index_spaces() const override { return empty_map_; @@ -399,8 +419,7 @@ class RangeIndexSpaceImpl: public IndexSpaceInterface { NameToRangeMap named_ranges_; /**< Map from name to subspace ranges*/ std::map named_subspaces_; /**< Map from names to (sub) IndexSpaces */ SpinAttribute spin_; /**< Spin attribute associated with the IndexSpace */ - SpatialAttribute spatial_; /**< Spatial attribute associated with the IndexSpace */ - std::vector empty_vec_; /**< Empty vector for dependencies */ + SpatialAttribute spatial_; /**< Spatial attribute associated with the IndexSpace */ std::map empty_map_; /**< Empty map for dependency relations */ /** @@ -566,7 +585,9 @@ class SubSpaceImpl: public IndexSpaceInterface { // Maximum number of indices in this index space std::size_t max_num_indices() const override { return indices_.size(); } - const std::vector& key_tiled_index_spaces() const override { return empty_vec_; } + const std::vector& key_tiled_index_spaces() const override { + return internal::empty_tiled_index_space_vec(); + } const std::map& map_tiled_index_spaces() const override { return empty_map_; @@ -614,7 +635,6 @@ class SubSpaceImpl: public IndexSpaceInterface { NameToRangeMap named_ranges_; /**< Map from name to subspace ranges*/ std::map named_subspaces_; /**< Map from names to (sub) IndexSpaces */ IndexSpace root_space_; /**< Root IndexSpace */ - std::vector empty_vec_; /**< Empty vector for dependencies */ std::map empty_map_; /**< Empty map for dependency relations */ /** * @brief Helper method for constructing the new set of @@ -679,7 +699,7 @@ class AggregateSpaceImpl: public IndexSpaceInterface { indices_{construct_indices(spaces)}, named_ranges_{named_ranges}, named_subspaces_{construct_subspaces(named_ranges)} { - // EXPECTS(has_duplicate(indices_)); + // EXPECTS(has_no_duplicate(indices_)); if(names.size() > 0) { add_ref_names(spaces, names); } if(subspace_references.size() > 0) { add_subspace_references(subspace_references); } } @@ -714,7 +734,9 @@ class AggregateSpaceImpl: public IndexSpaceInterface { // Maximum number of indices in this index space std::size_t max_num_indices() const override { return indices_.size(); } - const std::vector& key_tiled_index_spaces() const override { return empty_vec_; } + const std::vector& key_tiled_index_spaces() const override { + return internal::empty_tiled_index_space_vec(); + } const std::map& map_tiled_index_spaces() const override { return empty_map_; @@ -795,7 +817,6 @@ class AggregateSpaceImpl: public IndexSpaceInterface { NameToRangeMap named_ranges_; /**< Map from name to subspace ranges*/ std::map named_subspaces_; /**< Map from names to (sub) IndexSpaces */ std::vector empty_range_; /**< Empty range vector for spin relation */ - std::vector empty_vec_; /**< Empty vector for dependencies */ std::map empty_map_; /**< Empty map for dependency relations */ /** @@ -921,15 +942,7 @@ class DependentIndexSpaceImpl: public IndexSpaceInterface { * IndexSpaces */ DependentIndexSpaceImpl(const std::vector& indep_spaces, - const std::map& dep_space_relation): - dep_spaces_{indep_spaces}, dep_space_relation_{dep_space_relation}, named_ranges_{} { - // std::cerr << __FUNCTION__ << " " << __LINE__ << "\n"; - max_size_ = 0; - for(const auto& pair: dep_space_relation) { - max_size_ = std::max(max_size_, pair.second.num_indices()); - } - // std::cerr << __FUNCTION__ << " " << __LINE__ << "\n"; - } + const std::map& dep_space_relation); /*** * @brief Construct a new Dependent Index Space Impl object @@ -942,10 +955,15 @@ class DependentIndexSpaceImpl: public IndexSpaceInterface { */ DependentIndexSpaceImpl(const std::vector& indep_spaces, const IndexSpace& ref_space, - const std::map& dep_space_relation): - dep_spaces_{indep_spaces}, dep_space_relation_{dep_space_relation}, named_ranges_{} {} + const std::map& dep_space_relation); /// @todo do we need these constructor/operators + // Note: kept as in-class '= default' (not forced out-of-line). These are + // only instantiated on actual use, which happens in index_space.cpp where + // TiledIndexSpace is complete. Forcing them out-of-line would eagerly + // instantiate the copy/move assignment of named_ranges_, whose value type is + // a 'const std::vector' (NameToRangeMap) and is therefore not + // assignable -- ill-formed under libc++. DependentIndexSpaceImpl(DependentIndexSpaceImpl&&) = default; DependentIndexSpaceImpl(const DependentIndexSpaceImpl&) = default; DependentIndexSpaceImpl& operator=(DependentIndexSpaceImpl&&) = default; @@ -1006,7 +1024,8 @@ class DependentIndexSpaceImpl: public IndexSpaceInterface { return dep_spaces_; } - size_t num_key_tiled_index_spaces() const override { return dep_spaces_.size(); } + // Out-of-line (index_space.cpp): .size() needs TiledIndexSpace complete. + size_t num_key_tiled_index_spaces() const override; const std::map& map_tiled_index_spaces() const override { return dep_space_relation_; diff --git a/src/tamm/kernels/assign.hpp b/src/tamm/kernels/assign.hpp index bcb1333d1..cbcc3942e 100644 --- a/src/tamm/kernels/assign.hpp +++ b/src/tamm/kernels/assign.hpp @@ -7,9 +7,11 @@ #include "hptt/hptt.h" #include +#include #include #include #include +#include #include namespace tamm { @@ -168,9 +170,10 @@ void ipacc4(const SizeVec& loop_dims, T* dst, const SizeVec& loop_dld, T scale, } } -inline size_t idx(int n, const size_t* id, const Size* sz, const PermVector& p) { - Size idx = 0; - for(int i = 0; i < n - 1; i++) { idx = (idx + id[p[i]]) * sz[p[i + 1]]; } +inline size_t idx(std::span id, const Size* sz, const PermVector& p) { + const size_t n = id.size(); + Size idx = 0; + for(size_t i = 0; i + 1 < n; i++) { idx = (idx + id[p[i]]) * sz[p[i + 1]]; } if(n > 0) { idx += id[p[n - 1]]; } return idx.value(); } @@ -190,35 +193,35 @@ void index_permute_acc(T* dbuf, const T* sbuf, const PermVector& perm_to_dest, c } else if(ndim == 2) { Size sz[] = {ddims[0], ddims[1]}; - size_t i[2]; + std::array i{}; size_t c; for(c = 0, i[0] = 0; i[0] < sz[0]; i[0]++) { for(i[1] = 0; i[1] < sz[1]; i[1]++, c++) { - dbuf[c] += scale * sbuf[idx(2, i, sz, perm_to_dest)]; + dbuf[c] += scale * sbuf[idx(i, sz, perm_to_dest)]; } } } else if(ndim == 3) { Size sz[] = {ddims[0], ddims[1], ddims[2]}; - size_t i[3]; + std::array i{}; size_t c; for(c = 0, i[0] = 0; i[0] < sz[0]; i[0]++) { for(i[1] = 0; i[1] < sz[1]; i[1]++) { for(i[2] = 0; i[2] < sz[2]; i[2]++, c++) { - dbuf[c] += scale * sbuf[idx(3, i, sz, perm_to_dest)]; + dbuf[c] += scale * sbuf[idx(i, sz, perm_to_dest)]; } } } } else if(ndim == 4) { Size sz[] = {ddims[0], ddims[1], ddims[2], ddims[3]}; - size_t i[4]; + std::array i{}; size_t c; for(c = 0, i[0] = 0; i[0] < sz[0]; i[0]++) { for(i[1] = 0; i[1] < sz[1]; i[1]++) { for(i[2] = 0; i[2] < sz[2]; i[2]++) { for(i[3] = 0; i[3] < sz[3]; i[3]++, c++) { - dbuf[c] += scale * sbuf[idx(4, i, sz, perm_to_dest)]; + dbuf[c] += scale * sbuf[idx(i, sz, perm_to_dest)]; } } } @@ -242,35 +245,35 @@ void index_permute(T* dbuf, const T* sbuf, const PermVector& perm_to_dest, const } else if(ndim == 2) { Size sz[] = {ddims[0], ddims[1]}; - size_t i[2]; + std::array i{}; size_t c; for(c = 0, i[0] = 0; i[0] < sz[0]; i[0]++) { for(i[1] = 0; i[1] < sz[1]; i[1]++, c++) { - dbuf[c] = scale * sbuf[idx(2, i, sz, perm_to_dest)]; + dbuf[c] = scale * sbuf[idx(i, sz, perm_to_dest)]; } } } else if(ndim == 3) { Size sz[] = {ddims[0], ddims[1], ddims[2]}; - size_t i[3]; + std::array i{}; size_t c; for(c = 0, i[0] = 0; i[0] < sz[0]; i[0]++) { for(i[1] = 0; i[1] < sz[1]; i[1]++) { for(i[2] = 0; i[2] < sz[2]; i[2]++, c++) { - dbuf[c] = scale * sbuf[idx(3, i, sz, perm_to_dest)]; + dbuf[c] = scale * sbuf[idx(i, sz, perm_to_dest)]; } } } } else if(ndim == 4) { Size sz[] = {ddims[0], ddims[1], ddims[2], ddims[3]}; - size_t i[4]; + std::array i{}; size_t c; for(c = 0, i[0] = 0; i[0] < sz[0]; i[0]++) { for(i[1] = 0; i[1] < sz[1]; i[1]++) { for(i[2] = 0; i[2] < sz[2]; i[2]++) { for(i[3] = 0; i[3] < sz[3]; i[3]++, c++) { - dbuf[c] = scale * sbuf[idx(4, i, sz, perm_to_dest)]; + dbuf[c] = scale * sbuf[idx(i, sz, perm_to_dest)]; } } } diff --git a/src/tamm/label_translator.hpp b/src/tamm/label_translator.hpp index bfa8c080c..bf1023599 100644 --- a/src/tamm/label_translator.hpp +++ b/src/tamm/label_translator.hpp @@ -143,5 +143,3 @@ class LabelTranslator { }; // class LabelTranslator } // namespace tamm::internal - -#pragma once diff --git a/src/tamm/labeled_tensor.hpp b/src/tamm/labeled_tensor.hpp index 0742e6727..167477b54 100644 --- a/src/tamm/labeled_tensor.hpp +++ b/src/tamm/labeled_tensor.hpp @@ -3,6 +3,7 @@ // #include "tamm/ops.hpp" #include "tamm/tensor.hpp" #include +#include namespace tamm { template @@ -40,9 +41,9 @@ class LabeledTensor { const StringLabelVec& str_labels() const { return slv_; } const std::vector& str_map() const { return str_map_; } - void set_labels(const IndexLabelVec& ilv) { + void set_labels(IndexLabelVec ilv) { EXPECTS(ilv_.size() == ilv.size()); - ilv_ = ilv; + ilv_ = std::move(ilv); slv_.clear(); slv_.resize(ilv_.size()); str_map_ = std::vector(ilv_.size(), false); diff --git a/src/tamm/lru_cache.hpp b/src/tamm/lru_cache.hpp index 37805a456..ac676a76b 100644 --- a/src/tamm/lru_cache.hpp +++ b/src/tamm/lru_cache.hpp @@ -56,14 +56,20 @@ class LRUCache { return {hit, cached_value_[key]}; } + /// Return a reference to the cached value for @p key. + /// @pre key is present in the cache (checked via EXPECTS) Value& access(const Key& key) { - EXPECTS(cached_value_.find(key) != cached_value_.end()); - return *cached_value_.find(key); + auto it = cached_value_.find(key); + EXPECTS(it != cached_value_.end()); + return it->second; } + /// Return a const reference to the cached value for @p key. + /// @pre key is present in the cache (checked via EXPECTS) const Value& access(const Key& key) const { - EXPECTS(cached_value_.find(key) != cached_value_.end()); - return *cached_value_.find(key); + auto it = cached_value_.find(key); + EXPECTS(it != cached_value_.end()); + return it->second; } void gather_stats(std::vector& vec) { @@ -84,7 +90,7 @@ class LRUCache { } }; uint32_t max_size_; - uint32_t cycle_; + uint32_t cycle_{0}; // fix: initialize to avoid UB std::map> cache_; std::map cycle_to_key_; std::map cached_value_; diff --git a/src/tamm/mapop.hpp b/src/tamm/mapop.hpp index dcf943714..914beaa60 100644 --- a/src/tamm/mapop.hpp +++ b/src/tamm/mapop.hpp @@ -43,7 +43,7 @@ class MapOp: public Op { OpType op_type() const override { return OpType::map; } std::shared_ptr clone() const override { - return std::shared_ptr(new MapOp{*this}); + return std::make_shared(*this); } void execute(ExecutionContext& ec, ExecutionHW hw = ExecutionHW::CPU) override { @@ -193,29 +193,7 @@ class MapOp: public Op { for(size_t i = 0; i < N; i++) { ilv.insert(ilv.end(), rhs_[i].labels().begin(), rhs_[i].labels().end()); } - - for(size_t i = 0; i < ilv.size(); i++) { - for(const auto& dl: ilv[i].secondary_labels()) { - size_t j; - for(j = 0; j < ilv.size(); j++) { - if(dl.tiled_index_space() == ilv[j].tiled_index_space() && dl.label() == ilv[j].label()) { - break; - } - } - EXPECTS(j < ilv.size()); - } - } - - for(size_t i = 0; i < ilv.size(); i++) { - const auto& ilbl = ilv[i]; - for(size_t j = i + 1; j < ilv.size(); j++) { - const auto& jlbl = ilv[j]; - if(ilbl.tiled_index_space() == jlbl.tiled_index_space() && ilbl.label() == jlbl.label() && - ilbl.label_str() == jlbl.label_str()) { - EXPECTS(ilbl == jlbl); - } - } - } + internal::validate_index_labels(ilv); } LabeledTensorT lhs_; diff --git a/src/tamm/memory_manager_local.hpp b/src/tamm/memory_manager_local.hpp index 75d8867d4..42e67affa 100644 --- a/src/tamm/memory_manager_local.hpp +++ b/src/tamm/memory_manager_local.hpp @@ -267,7 +267,7 @@ class MemoryManagerLocal: public MemoryManager { for(size_t i = 0; i < mp.local_nelements().value(); i++) { switch(mp.eltype_) { case ElementType::double_precision: - os << i << " " << (reinterpret_cast(mp.buf_))[i] << std::endl; + os << i << " " << (reinterpret_cast(mp.buf_))[i] << '\n'; break; default: NOT_IMPLEMENTED(); } diff --git a/src/tamm/multop.hpp b/src/tamm/multop.hpp index 9edd13ef9..ecf270050 100644 --- a/src/tamm/multop.hpp +++ b/src/tamm/multop.hpp @@ -31,77 +31,35 @@ namespace tamm::internal { template struct MultOpPlanBase { using MultOpT = MultOp; + + // writes()/accumulates()/reads() are identical for every plan: an assign + // writes the LHS, an update accumulates into it, and both RHS operands are + // read. (The former local/global split was a no-op for scheduling — the + // variants were concatenated and only membership matters for dependency + // tests — and no code queried the local/global variants individually.) TensorBase* writes(const MultOpT& multop) const { - auto ret1 = local_writes(multop); - auto ret2 = global_writes(multop); - ret1.insert(ret1.end(), ret2.begin(), ret2.end()); - return !ret1.empty() ? ret1[0] : nullptr; + return multop.is_assign() ? multop.lhs().base_ptr() : nullptr; } - TensorBase* accumulates(const MultOpT& multop) const { - auto ret1 = local_accumulates(multop); - auto ret2 = global_accumulates(multop); - ret1.insert(ret1.end(), ret2.begin(), ret2.end()); - return !ret1.empty() ? ret1[0] : nullptr; + return multop.is_assign() ? nullptr : multop.lhs().base_ptr(); } - std::vector reads(const MultOpT& multop) const { - auto ret1 = local_reads(multop); - auto ret2 = global_reads(multop); - ret1.insert(ret1.end(), ret2.begin(), ret2.end()); - return !ret1.empty() ? ret1 : std::vector{}; + return {multop.rhs1().base_ptr(), multop.rhs2().base_ptr()}; } - virtual std::vector global_writes(const MultOpT& multop) const = 0; - virtual std::vector global_accumulates(const MultOpT& multop) const = 0; - virtual std::vector global_reads(const MultOpT& multop) const = 0; - virtual std::vector local_writes(const MultOpT& multop) const = 0; - virtual std::vector local_accumulates(const MultOpT& multop) const = 0; - virtual std::vector local_reads(const MultOpT& multop) const = 0; - virtual void apply(const MultOpT& multop, ExecutionContext& ec, ExecutionHW hw) = 0; + virtual ~MultOpPlanBase() = default; }; // MultOpPlanBase template struct FlatMultPlan: public MultOpPlanBase { using MultOpT = MultOp; - std::vector global_writes(const MultOpT& multop) const override { return {}; } - std::vector global_accumulates(const MultOpT& multop) const override { return {}; } - - std::vector global_reads(const MultOpT& multop) const override { return {}; } - - std::vector local_writes(const MultOpT& multop) const override { - if(multop.is_assign()) { return {multop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_accumulates(const MultOpT& multop) const override { - if(!multop.is_assign()) { return {multop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_reads(const MultOpT& multop) const override { - return {multop.rhs1().base_ptr(), multop.rhs2().base_ptr()}; - } void apply(const MultOpT& multop, ExecutionContext& ec, ExecutionHW hw) override; }; // FlatMultPlan template struct LHSMultPlan: public MultOpPlanBase { using MultOpT = MultOp; - std::vector global_writes(const MultOpT& multop) const override { return {}; } - std::vector global_accumulates(const MultOpT& multop) const override { return {}; } - std::vector global_reads(const MultOpT& multop) const override { return {}; } - - std::vector local_writes(const MultOpT& multop) const override { - if(multop.is_assign()) { return {multop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_accumulates(const MultOpT& multop) const override { - if(!multop.is_assign()) { return {multop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_reads(const MultOpT& multop) const override { - return {multop.rhs1().base_ptr(), multop.rhs2().base_ptr()}; - } void apply(const MultOpT& multop, ExecutionContext& ec, ExecutionHW hw) override; }; // LHSMultPlan @@ -109,29 +67,6 @@ template { using MultOpT = MultOp; - std::vector global_writes(const MultOpT& multop) const override { - if(multop.is_assign()) { return {multop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector global_accumulates(const MultOpT& multop) const override { - if(!multop.is_assign()) { return {multop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector global_reads(const MultOpT& multop) const override { - return {multop.rhs1().base_ptr(), multop.rhs2().base_ptr()}; - } - - std::vector local_writes(const MultOpT& multop) const override { - if(multop.is_assign()) { return {multop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_accumulates(const MultOpT& multop) const override { - if(!multop.is_assign()) { return {multop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_reads(const MultOpT& multop) const override { - return {multop.rhs1().base_ptr(), multop.rhs2().base_ptr()}; - } void apply(const MultOpT& multop, ExecutionContext& ec, ExecutionHW hw) override; }; // GeneralFlatMultPlan @@ -139,28 +74,6 @@ template { using MultOpT = MultOp; - std::vector global_writes(const MultOpT& multop) const override { - if(multop.is_assign()) { return {multop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector global_accumulates(const MultOpT& multop) const override { - if(!multop.is_assign()) { return {multop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector global_reads(const MultOpT& multop) const override { - return {multop.rhs1().base_ptr(), multop.rhs2().base_ptr()}; - } - std::vector local_writes(const MultOpT& multop) const override { - if(multop.is_assign()) { return {multop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_accumulates(const MultOpT& multop) const override { - if(!multop.is_assign()) { return {multop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_reads(const MultOpT& multop) const override { - return {multop.rhs1().base_ptr(), multop.rhs2().base_ptr()}; - } void apply(const MultOpT& multop, ExecutionContext& ec, ExecutionHW hw) override; }; // GeneralLHSMultPlan @@ -172,7 +85,7 @@ template struct AddBuf { AddBuf(T2* ta, T3* tb, T1* cbuf, const IndexVector& blockid): blockid_{blockid}, cbuf_{cbuf}, ta_{ta}, tb_{tb} {} - ~AddBuf() {} + ~AddBuf() = default; T1* cbuf_; T2* abuf_; @@ -227,9 +140,9 @@ class MultOp: public Op { rhs2_lbls = IndexLabelVec(labels.begin() + lhs.labels().size() + rhs1.labels().size(), labels.begin() + lhs.labels().size() + rhs1.labels().size() + rhs2.labels().size()); - lhs_.set_labels(lhs_lbls); - rhs1_.set_labels(rhs1_lbls); - rhs2_.set_labels(rhs2_lbls); + lhs_.set_labels(std::move(lhs_lbls)); + rhs1_.set_labels(std::move(rhs1_lbls)); + rhs2_.set_labels(std::move(rhs2_lbls)); } if(lhs.has_str_lbl()) { fillin_labels(); } @@ -238,7 +151,7 @@ class MultOp: public Op { validate(); } - MultOp(const MultOp&) = default; + // Copy/move are implicitly generated (Rule of Zero); clone() copies. LabeledTensorT1 lhs() const { return lhs_; } @@ -268,7 +181,9 @@ class MultOp: public Op { return result; } - std::shared_ptr clone() const override { return std::shared_ptr(new MultOp{*this}); } + std::shared_ptr clone() const override { + return std::make_shared(*this); + } using TensorElType1 = typename LabeledTensorT1::element_type; using TensorElType2 = typename LabeledTensorT2::element_type; @@ -285,7 +200,7 @@ class MultOp: public Op { all_labels.insert(all_labels.end(), rhs2_.labels().begin(), rhs2_.labels().end()); LabelLoopNest loop_nest{all_labels}; - std::vector*> add_bufs; + std::vector>> add_bufs; // function to compute one block auto lambda = [=, &oprof, &add_bufs, &loop_nest, &ec](const IndexVector itval) { auto ctensor = lhs_.tensor(); @@ -461,13 +376,13 @@ class MultOp: public Op { th_a = static_cast(memDevicePool.allocate(asize * sizeof(TensorElType2))); th_b = static_cast(memDevicePool.allocate(bsize * sizeof(TensorElType3))); - ab = new AddBuf{th_a, th_b, cbuf, - translated_cblockid}; + add_bufs.push_back(std::make_unique>( + th_a, th_b, cbuf, translated_cblockid)); #else - ab = new AddBuf{ctensor, cbuf, - translated_cblockid}; + add_bufs.push_back(std::make_unique>( + ctensor, cbuf, translated_cblockid)); #endif - add_bufs.push_back(ab); + ab = add_bufs.back().get(); { TimerGuard tg_bc{&oprof.multOpBCTime}; @@ -530,7 +445,6 @@ class MultOp: public Op { // add the computed update to the tensor ctensor.add(translated_cblockid, {ab->cbuf_, csize}); } - delete ab; add_bufs.clear(); #endif } @@ -540,13 +454,8 @@ class MultOp: public Op { //@todo make parallel // do_work(ec, loop_nest, lambda); - bool has_sparse_labels = false; - for(auto& lbl: all_labels) { - if(lbl.is_dependent()) { - has_sparse_labels = true; - break; - } - } + const bool has_sparse_labels = + std::ranges::any_of(all_labels, [](const auto& lbl) { return lbl.is_dependent(); }); if(1 && (lhs_.tensor().is_dense() /* && !lhs_.tensor().has_spin() */) && (rhs1_.tensor().is_dense() /* && !rhs1_.tensor().has_spin() */) && @@ -565,7 +474,6 @@ class MultOp: public Op { for(auto& ab: add_bufs) { (ab->tensor_).nb_add(ab->blockid_, ab->cbuf_, &(ab->nbhdl_)); ab->wait(); - delete ab; } add_bufs.clear(); } @@ -662,16 +570,16 @@ class MultOp: public Op { SizeVec cdims_sz; for(const auto v: cdims) { cdims_sz.push_back(v); } - AddBuf* ab{nullptr}; + std::unique_ptr> ab; #if defined(USE_CUDA) || defined(USE_HIP) || defined(USE_DPCPP) TensorElType2* th_a{nullptr}; TensorElType3* th_b{nullptr}; - ab = new AddBuf{ - th_a, th_b, {}, translated_cblockid}; + ab = std::make_unique>( + th_a, th_b, static_cast(nullptr), translated_cblockid); #else - ab = - new AddBuf{ctensor, {}, translated_cblockid}; + ab = std::make_unique>( + ctensor, nullptr, translated_cblockid); #endif { @@ -846,7 +754,7 @@ class MultOp: public Op { #endif } // multoptime - delete ab; + ab.reset(); memHostPool.deallocate(cbuf, csize * sizeof(TensorElType1)); }; //@todo use a scheduler @@ -992,32 +900,9 @@ class MultOp: public Op { tamm_terminate(os.str()); } - IndexLabelVec ilv{lhs_.labels()}; - ilv.insert(ilv.end(), rhs1_.labels().begin(), rhs1_.labels().end()); - ilv.insert(ilv.end(), rhs2_.labels().begin(), rhs2_.labels().end()); - - for(size_t i = 0; i < ilv.size(); i++) { - for(const auto& dl: ilv[i].secondary_labels()) { - size_t j; - for(j = 0; j < ilv.size(); j++) { - if(dl.tiled_index_space() == ilv[j].tiled_index_space() && dl.label() == ilv[j].label()) { - break; - } - } - EXPECTS(j < ilv.size()); - } - } - - for(size_t i = 0; i < ilv.size(); i++) { - const auto& ilbl = ilv[i]; - for(size_t j = i + 1; j < ilv.size(); j++) { - const auto& jlbl = ilv[j]; - if(ilbl.tiled_index_space() == jlbl.tiled_index_space() && ilbl.label() == jlbl.label() && - ilbl.label_str() == jlbl.label_str()) { - EXPECTS(ilbl == jlbl); - } - } - } + const auto ilv = internal::merge_vector(lhs_.labels(), rhs1_.labels(), + rhs2_.labels()); + internal::validate_index_labels(ilv); } LabeledTensorT1 lhs_; diff --git a/src/tamm/op_visitors.hpp b/src/tamm/op_visitors.hpp index 645e2423b..3148d2088 100644 --- a/src/tamm/op_visitors.hpp +++ b/src/tamm/op_visitors.hpp @@ -1,7 +1,6 @@ #pragma once #include "tamm/interfaces.hpp" -#include "tamm/op_cost.hpp" #include "tamm/op_dag.hpp" #include "tamm/tiled_index_space.hpp" #include "tamm/types.hpp" @@ -1924,7 +1923,7 @@ class CanonicalizeVisitor: public VisitorBase { } if(no_slicing) { - canonicalized_ops_.push_back({std::move(parforop.op().clone()), LabelPair{}}); + canonicalized_ops_.push_back({parforop.op().clone(), LabelPair{}}); } } diff --git a/src/tamm/opmin.hpp b/src/tamm/opmin.hpp index 90d6ec77f..a3be195f5 100644 --- a/src/tamm/opmin.hpp +++ b/src/tamm/opmin.hpp @@ -649,7 +649,7 @@ class OpMin { auto optimized_op = optimizer.optimize(); optimized_ops.push_back(std::move(optimized_op)); } - else { optimized_ops.push_back(std::move(op->clone())); } + else { optimized_ops.push_back(op->clone()); } } std::unique_ptr result_op = (*optimized_ops.at(0)).clone(); @@ -661,7 +661,7 @@ class OpMin { result_op->accept(clear_visitor); result_op->set_attribute(lhs_labels); - return std::move(result_op); + return result_op; } protected: diff --git a/src/tamm/ops.hpp b/src/tamm/ops.hpp index 083071ae7..dab89e755 100644 --- a/src/tamm/ops.hpp +++ b/src/tamm/ops.hpp @@ -16,12 +16,8 @@ namespace tamm::internal { template class LabelMap { public: - LabelMap() = default; - LabelMap(const LabelMap&) = default; - LabelMap(LabelMap&&) = default; - LabelMap& operator=(const LabelMap&) = default; - LabelMap& operator=(LabelMap&&) = default; - ~LabelMap() = default; + // Rule of Zero: only member is a std::map, so all special members are correct + // when compiler-generated. LabelMap& update(const IndexLabelVec& labels, const std::vector& vals) { EXPECTS(labels.size() == vals.size()); diff --git a/src/tamm/perm.hpp b/src/tamm/perm.hpp index ef0374c7b..9cde537dc 100644 --- a/src/tamm/perm.hpp +++ b/src/tamm/perm.hpp @@ -27,7 +27,8 @@ PermVector perm_compute(const std::vector& from, const std::vector& to) { for(auto p: to) { auto itr = std::find(from.begin(), from.end(), p); EXPECTS(itr != from.end()); - layout.push_back(itr - from.begin()); + // Use std::distance instead of raw iterator subtraction for generality + layout.push_back(static_cast(std::distance(from.begin(), itr))); } return layout; } @@ -39,8 +40,9 @@ bool are_permutations(const std::vector& vec1, const std::vector& vec2) { for(size_t i = 0; i < vec1.size(); i++) { auto it = std::find(vec2.begin(), vec2.end(), vec1[i]); if(it == vec2.end()) { return false; } - if(taken[std::distance(vec2.begin(), it)] == true) { return false; } - taken[std::distance(vec2.begin(), it)] = true; + auto idx = static_cast(std::distance(vec2.begin(), it)); + if(taken[idx]) { return false; } + taken[idx] = true; } return true; } @@ -50,9 +52,9 @@ PermVector perm_map_compute(const std::vector& unique_vec, const std::vector< PermVector ret; for(const auto& val: vec_required) { auto it = std::find(unique_vec.begin(), unique_vec.end(), val); - EXPECTS(it >= unique_vec.begin()); EXPECTS(it != unique_vec.end()); - ret.push_back(it - unique_vec.begin()); + // Use std::distance instead of raw iterator subtraction + ret.push_back(static_cast(std::distance(unique_vec.begin(), it))); } return ret; } @@ -62,8 +64,10 @@ std::vector perm_map_apply(const std::vector& input_vec, const std::vector& perm_map) { std::vector ret; for(const auto& pm: perm_map) { - EXPECTS(pm < input_vec.size()); - ret.push_back(input_vec[pm]); + // Fix: cast to size_t before comparison to avoid signed/unsigned mismatch + // (a negative pm would silently wrap past this check otherwise) + EXPECTS(static_cast(pm) < input_vec.size()); + ret.push_back(input_vec[static_cast(pm)]); } return ret; } @@ -73,8 +77,9 @@ void perm_map_apply(std::vector& out_vec, const std::vector& input_vec, const std::vector& perm_map) { out_vec.resize(perm_map.size()); for(size_t i = 0; i < perm_map.size(); i++) { - EXPECTS(perm_map[i] < input_vec.size()); - out_vec[i] = input_vec[perm_map[i]]; + // Fix: cast to size_t before comparison to avoid signed/unsigned mismatch + EXPECTS(static_cast(perm_map[i]) < input_vec.size()); + out_vec[i] = input_vec[static_cast(perm_map[i])]; } } diff --git a/src/tamm/proc_grid.hpp b/src/tamm/proc_grid.hpp index 2d7630ff6..23a7b2e46 100644 --- a/src/tamm/proc_grid.hpp +++ b/src/tamm/proc_grid.hpp @@ -35,7 +35,7 @@ static double dd_ev(const int64_t ndims, const std::vector ardims, t = 1.0; for(k = 0; k < ndims; k++) { q = (ardims[k] / pedims[k]) * pedims[k]; - t = t * (q / (double) ardims[k]); + t = t * (q / static_cast(ardims[k])); } return t; } @@ -276,12 +276,12 @@ static std::vector compute_proc_grid(const int64_t ndims, h = istart; q = (tard[istart] < p0 * pedims[istart]) ? 1.1 - : (tard[istart] % (p0 * pedims[istart])) / (double) tard[istart]; + : (tard[istart] % (p0 * pedims[istart])) / static_cast(tard[istart]); for(j = 1; j < ndims; j++) { ilook = (istart + istep * j) % ndims; w = (tard[ilook] < p0 * pedims[ilook]) ? 1.1 - : (tard[ilook] % (p0 * pedims[ilook])) / (double) tard[ilook]; + : (tard[ilook] % (p0 * pedims[ilook])) / static_cast(tard[ilook]); if(w < q) { q = w; h = ilook; diff --git a/src/tamm/proc_group.hpp b/src/tamm/proc_group.hpp index de70d6db4..72d48831a 100644 --- a/src/tamm/proc_group.hpp +++ b/src/tamm/proc_group.hpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -17,14 +18,13 @@ class ProcGroup { * */ ProcGroup(): pginfo_{std::make_shared()} {} - ProcGroup(const ProcGroup&) = default; - ProcGroup(ProcGroup&& pg): pginfo_{std::move(pg.pginfo_)} {} - ProcGroup& operator=(ProcGroup pg) { - using std::swap; - swap(*this, pg); - return *this; - } - ~ProcGroup() = default; + // Rule of Zero: the only member is a shared_ptr, so compiler-generated copy, + // move, and assignment are all correct. + ProcGroup(const ProcGroup&) = default; + ProcGroup(ProcGroup&&) = default; + ProcGroup& operator=(const ProcGroup&) = default; + ProcGroup& operator=(ProcGroup&&) = default; + ~ProcGroup() = default; /** * @brief Construct a new Proc Group object by wrapping the given MPI @@ -155,7 +155,7 @@ class ProcGroup { // Create subgroup from first nranks of parent group static ProcGroup create_subgroup(const ProcGroup& parent_group, int nranks) { std::vector ranks(nranks); - for(int i = 0; i < nranks; i++) ranks[i] = i; + std::iota(ranks.begin(), ranks.end(), 0); return create_subgroup(parent_group, ranks); } @@ -763,7 +763,7 @@ class ProcGroup { // also works when GA is initialized with an existing MPI communicator MPI_Comm_group(GA_MPI_Comm(), &group_world); - for(int i = 0; i < nranks; i++) { ranks[i] = i; } + std::iota(ranks, ranks + nranks, 0); MPI_Group_translate_ranks(group, nranks, ranks, group_world, ranks_world); int ga_pg_default = GA_Pgroup_get_default(); @@ -790,7 +790,7 @@ class ProcGroup { MPI_Comm_group(parent_group.comm(), &group_world); - for(int i = 0; i < nranks; i++) { ranks[i] = i; } + std::iota(ranks, ranks + nranks, 0); MPI_Group_translate_ranks(group, nranks, ranks, group_world, ranks_world); int ga_pg_default = GA_Pgroup_get_default(); diff --git a/src/tamm/range.hpp b/src/tamm/range.hpp index 97da25b68..0b3f2fb5c 100644 --- a/src/tamm/range.hpp +++ b/src/tamm/range.hpp @@ -147,9 +147,12 @@ class Range { } protected: - Index lo_; /**< Low end of range */ - Index hi_; /**< High end of range */ - Index step_; /**< step size for the range */ + // In-class initializers: a default-constructed Range (Range() = default) must + // not leave these members uninitialized — Range is used as a map value and in + // hashing, which would otherwise read indeterminate values (UB). + Index lo_{0}; /**< Low end of range */ + Index hi_{0}; /**< High end of range */ + Index step_{1}; /**< step size for the range */ private: /** @@ -238,8 +241,8 @@ using AttributeToRangeMap = std::map>; namespace std { template<> struct hash { - typedef tamm::Range argument_type; - typedef std::size_t result_type; + using argument_type = tamm::Range; + using result_type = std::size_t; result_type operator()(argument_type const& range) const noexcept { using tamm::internal::hash_combine; diff --git a/src/tamm/scalar.hpp b/src/tamm/scalar.hpp index df21b82fb..43118b6a7 100644 --- a/src/tamm/scalar.hpp +++ b/src/tamm/scalar.hpp @@ -10,6 +10,14 @@ namespace tamm { +namespace detail { +// Local, dependency-free complex trait (avoids pulling in the heavyweight +// tamm/utils.hpp just for internal::is_complex_v). +template struct is_complex: std::false_type {}; +template struct is_complex>: std::true_type {}; +template inline constexpr bool is_complex_v = is_complex::value; +} // namespace detail + /** * @brief Scalar object for representing scalar values from different data types * @@ -24,11 +32,11 @@ class Scalar { Scalar(ElType eltype) { switch(eltype) { - case ElType::inv: value_ = (double) 1.0; break; - case ElType::i32: value_ = (int) 1; break; - case ElType::i64: value_ = (int64_t) 1; break; - case ElType::fp32: value_ = (float) 1.0; break; - case ElType::fp64: value_ = (double) 1.0; break; + case ElType::inv: value_ = 1.0; break; + case ElType::i32: value_ = 1; break; + case ElType::i64: value_ = int64_t{1}; break; + case ElType::fp32: value_ = 1.0f; break; + case ElType::fp64: value_ = 1.0; break; case ElType::cfp32: value_ = std::complex(1.0, 0.0); break; case ElType::cfp64: value_ = std::complex(1.0, 0.0); break; } @@ -81,6 +89,37 @@ class Scalar { ElementType value() const { return value_; } + /** + * @brief Extract the stored value converted to the requested element type T. + * + * Visits the underlying variant and converts the active alternative to T. + * Real target types drop the imaginary part of a stored complex value; + * complex target types are constructed from real alternatives. This is the + * type-safe replacement for `static_cast(scalar)`, which is ill-formed + * because Scalar holds a std::variant with no implicit numeric conversion. + */ + template + [[nodiscard]] T get() const { + return std::visit( + [](auto v) -> T { + using V = std::decay_t; + if constexpr(detail::is_complex_v) { + if constexpr(detail::is_complex_v) + return T{static_cast(v.real()), + static_cast(v.imag())}; + else + return T{static_cast(v)}; + } + else { + if constexpr(detail::is_complex_v) + return static_cast(v.real()); + else + return static_cast(v); + } + }, + value_); + } + private: ElementType value_; }; diff --git a/src/tamm/scanop.hpp b/src/tamm/scanop.hpp index 214a3411b..2d07b3784 100644 --- a/src/tamm/scanop.hpp +++ b/src/tamm/scanop.hpp @@ -24,7 +24,7 @@ class ScanOp: public Op { OpType op_type() const override { return OpType::scan; } std::shared_ptr clone() const override { - return std::shared_ptr(new ScanOp{*this}); + return std::make_shared(*this); } void execute(ExecutionContext& ec, ExecutionHW hw = ExecutionHW::CPU) override { @@ -91,29 +91,7 @@ class ScanOp: public Op { */ void validate() { IndexLabelVec ilv{lhs_.labels()}; - - for(size_t i = 0; i < ilv.size(); i++) { - for(const auto& dl: ilv[i].secondary_labels()) { - size_t j; - for(j = 0; j < ilv.size(); j++) { - if(dl.tiled_index_space() == ilv[j].tiled_index_space() && dl.label() == ilv[j].label()) { - break; - } - } - EXPECTS(j < ilv.size()); - } - } - - for(size_t i = 0; i < ilv.size(); i++) { - const auto& ilbl = ilv[i]; - for(size_t j = i + 1; j < ilv.size(); j++) { - const auto& jlbl = ilv[j]; - if(ilbl.tiled_index_space() == jlbl.tiled_index_space() && ilbl.label() == jlbl.label() && - ilbl.label_str() == jlbl.label_str()) { - EXPECTS(ilbl == jlbl); - } - } - } + internal::validate_index_labels(ilv); } LabeledTensorT lhs_; diff --git a/src/tamm/scheduler.hpp b/src/tamm/scheduler.hpp index 38f323907..a3cf6720e 100644 --- a/src/tamm/scheduler.hpp +++ b/src/tamm/scheduler.hpp @@ -1,5 +1,6 @@ #pragma once +#include #include #include "ga/ga-mpi.h" @@ -208,8 +209,8 @@ class Scheduler { auto misc_start = std::chrono::high_resolution_clock::now(); auto order = levelize_and_order(ops_, start_idx_, ops_.size()); EXPECTS(order.size() == ops_.size() - start_idx_); - size_t lvl = 0; - AtomicCounter* ac = new AtomicCounterGA(ec().pg(), order.size()); + size_t lvl = 0; + auto ac = std::make_unique(ec().pg(), order.size()); ac->allocate(0); auto misc_end = std::chrono::high_resolution_clock::now(); double misc_time = @@ -250,7 +251,7 @@ class Scheduler { // oprof.multOpAddTime = 0; // t1 = t3; } - ec().set_ac(IndexedAC(ac, i)); + ec().set_ac(IndexedAC(ac.get(), i)); if(ops_[order[i].second]->exhw_ != ExecutionHW::DEFAULT) execute_on = ops_[order[i].second]->exhw_; auto t2 = std::chrono::high_resolution_clock::now(); @@ -297,7 +298,7 @@ class Scheduler { ec().set_ac(IndexedAC(nullptr, 0)); misc_start = t3; ac->deallocate(); - delete ac; + ac.reset(); misc_end = std::chrono::high_resolution_clock::now(); misc_time += std::chrono::duration_cast>((misc_end - misc_start)).count(); @@ -428,17 +429,17 @@ class Scheduler { size_t off = start_idx_; for(size_t g: groups) { EXPECTS(g > 0); - AtomicCounter* ac = new AtomicCounterGA(ec().pg(), g); + auto ac = std::make_unique(ec().pg(), g); ac->allocate(0); for(size_t i = off; i < off + g; i++, start_idx_++) { - ec().set_ac(IndexedAC(ac, i - off)); + ec().set_ac(IndexedAC(ac.get(), i - off)); ops_[i]->execute(ec()); } ec().set_ac(IndexedAC(nullptr, 0)); ac->deallocate(); - delete ac; + ac.reset(); // memory fence. for now GA_Sync() // GA_Sync(); diff --git a/src/tamm/setop.hpp b/src/tamm/setop.hpp index 75fd73b63..475002fcf 100644 --- a/src/tamm/setop.hpp +++ b/src/tamm/setop.hpp @@ -6,6 +6,7 @@ #include #include +#include "tamm/block_scratch.hpp" #include "tamm/boundvec.hpp" #include "tamm/errors.hpp" #include "tamm/label_translator.hpp" @@ -26,102 +27,45 @@ namespace tamm::internal { template struct SetOpPlanBase { using SetOpT = SetOp; + + // writes()/accumulates() are identical for every plan: an assign writes the + // LHS, an update accumulates into it. (The former local/global split was a + // no-op — local+global were always concatenated and only element [0] used — + // and no code queried the local/global variants individually.) TensorBase* writes(const SetOpT& setop) const { - auto ret1 = local_writes(setop); - auto ret2 = global_writes(setop); - ret1.insert(ret1.end(), ret2.begin(), ret2.end()); - return !ret1.empty() ? ret1[0] : nullptr; + return setop.is_assign() ? setop.lhs().base_ptr() : nullptr; } - TensorBase* accumulates(const SetOpT& setop) const { - auto ret1 = local_accumulates(setop); - auto ret2 = global_accumulates(setop); - ret1.insert(ret1.end(), ret2.begin(), ret2.end()); - return !ret1.empty() ? ret1[0] : nullptr; + return setop.is_assign() ? nullptr : setop.lhs().base_ptr(); } - virtual std::vector global_writes(const SetOpT& setop) const = 0; - virtual std::vector global_accumulates(const SetOpT& setop) const = 0; - virtual std::vector local_writes(const SetOpT& setop) const = 0; - virtual std::vector local_accumulates(const SetOpT& setop) const = 0; - virtual void apply(const SetOpT& setop, ExecutionContext& ec, ExecutionHW hw) = 0; + virtual void apply(const SetOpT& setop, ExecutionContext& ec, ExecutionHW hw) = 0; + virtual ~SetOpPlanBase() = default; }; // SetOpPlanBase template struct FlatPlan: public SetOpPlanBase { using SetOpT = SetOp; - std::vector global_writes(const SetOpT& setop) const override { return {}; } - std::vector global_accumulates(const SetOpT& setop) const override { return {}; } - std::vector local_writes(const SetOpT& setop) const override { - if(setop.is_assign()) { return {setop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_accumulates(const SetOpT& setop) const override { - if(!setop.is_assign()) { return {setop.lhs().base_ptr()}; } - else { return {}; } - } void apply(const SetOpT& setop, ExecutionContext& ec, ExecutionHW hw) override; }; template struct LHSPlan: public SetOpPlanBase { using SetOpT = SetOp; - std::vector global_writes(const SetOpT& setop) const override { return {}; } - std::vector global_accumulates(const SetOpT& setop) const override { return {}; } - std::vector local_writes(const SetOpT& setop) const override { - if(setop.is_assign()) { return {setop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_accumulates(const SetOpT& setop) const override { - if(!setop.is_assign()) { return {setop.lhs().base_ptr()}; } - else { return {}; } - } void apply(const SetOpT& setop, ExecutionContext& ec, ExecutionHW hw) override; -}; // namespace tamm::internal +}; template struct GeneralFlatPlan: public SetOpPlanBase { using SetOpT = SetOp; - std::vector global_writes(const SetOpT& setop) const override { - if(setop.is_assign()) { return {setop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector global_accumulates(const SetOpT& setop) const override { - if(!setop.is_assign()) { return {setop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_writes(const SetOpT& setop) const override { - if(setop.is_assign()) { return {setop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_accumulates(const SetOpT& setop) const override { - if(!setop.is_assign()) { return {setop.lhs().base_ptr()}; } - else { return {}; } - } void apply(const SetOpT& setop, ExecutionContext& ec, ExecutionHW hw) override; -}; // GeneralFlatPlan +}; template struct GeneralLHSPlan: public SetOpPlanBase { using SetOpT = SetOp; - std::vector global_writes(const SetOpT& setop) const override { - if(setop.is_assign()) { return {setop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector global_accumulates(const SetOpT& setop) const override { - if(!setop.is_assign()) { return {setop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_writes(const SetOpT& setop) const override { - if(setop.is_assign()) { return {setop.lhs().base_ptr()}; } - else { return {}; } - } - std::vector local_accumulates(const SetOpT& setop) const override { - if(!setop.is_assign()) { return {setop.lhs().base_ptr()}; } - else { return {}; } - } void apply(const SetOpT& setop, ExecutionContext& ec, ExecutionHW hw) override; -}; // GeneralLHSPlan +}; } // namespace tamm::internal @@ -140,7 +84,7 @@ class SetOp: public Op { if(!lhs.has_str_lbl() && !lhs.labels().empty()) { auto lbls = lhs.labels(); internal::update_labels(lbls); - lhs_.set_labels(lbls); + lhs_.set_labels(std::move(lbls)); } if(lhs.has_str_lbl()) { fillin_labels(); } @@ -162,7 +106,7 @@ class SetOp: public Op { EXPECTS(plan_obj_ != nullptr); } - SetOp(const SetOp&) = default; + // Copy/move are implicitly generated (Rule of Zero); clone() copies. T alpha() const { return alpha_; } @@ -173,7 +117,7 @@ class SetOp: public Op { OpList canonicalize() const override { return OpList{(*this)}; } std::shared_ptr clone() const override { - return std::shared_ptr(new SetOp{*this}); + return std::make_shared(*this); } OpType op_type() const override { return OpType::set; } @@ -194,17 +138,17 @@ class SetOp: public Op { TensorBase* writes(const ExecutionContext& ec) const { if(lhs_.tensor().pg() == ec.pg()) { return plan_obj_->writes(*this); } - else { general_plan_obj_->writes(*this); } + else { return general_plan_obj_->writes(*this); } } std::vector reads(const ExecutionContext& ec) const { if(lhs_.tensor().pg() == ec.pg()) { return plan_obj_->reads(*this); } - else { general_plan_obj_->reads(*this); } + else { return general_plan_obj_->reads(*this); } } TensorBase* accumulates(const ExecutionContext& ec) const { if(lhs_.tensor().pg() == ec.pg()) { return plan_obj_->accumulates(*this); } - else { general_plan_obj_->accumulates(*this); } + else { return general_plan_obj_->accumulates(*this); } } bool is_memory_barrier() const { return false; } @@ -234,29 +178,7 @@ class SetOp: public Op { */ void validate() { IndexLabelVec ilv{lhs_.labels()}; - - for(size_t i = 0; i < ilv.size(); i++) { - for(const auto& dl: ilv[i].secondary_labels()) { - size_t j; - for(j = 0; j < ilv.size(); j++) { - if(dl.tiled_index_space() == ilv[j].tiled_index_space() && dl.label() == ilv[j].label()) { - break; - } - } - EXPECTS(j < ilv.size()); - } - } - - for(size_t i = 0; i < ilv.size(); i++) { - const auto& ilbl = ilv[i]; - for(size_t j = i + 1; j < ilv.size(); j++) { - const auto& jlbl = ilv[j]; - if(ilbl.tiled_index_space() == jlbl.tiled_index_space() && ilbl.label() == jlbl.label() && - ilbl.label_str() == jlbl.label_str()) { - EXPECTS(ilbl == jlbl); - } - } - } + internal::validate_index_labels(ilv); } LabeledTensorT lhs_; @@ -355,8 +277,8 @@ void GeneralFlatPlan::apply(const SetOp& s // EXPECTS(pg_lhs.size() == pg_ec.size()); - BlockSetPlan::OpType optype = is_assign ? optype = BlockSetPlan::OpType::set - : BlockSetPlan::OpType::update; + BlockSetPlan::OpType optype = + is_assign ? BlockSetPlan::OpType::set : BlockSetPlan::OpType::update; BlockSetPlan plan{lhs_lt.labels(), optype}; std::vector pg_lhs_in_ec = pg_lhs.rank_translate(pg_ec); @@ -364,6 +286,9 @@ void GeneralFlatPlan::apply(const SetOp& s Proc round_robin_counter = 0; Proc ec_pg_size = Proc{ec.pg().size()}; + // Reused (grow-only) scratch for the remote-copy case; RAII-owned. + internal::BlockScratch lhs_scratch; + for(size_t i = 0; i < pg_lhs_in_ec.size(); i++) { Proc assigned_proc; if(pg_lhs_in_ec[i] >= Proc{0}) { @@ -377,19 +302,16 @@ void GeneralFlatPlan::apply(const SetOp& s } if(proc_me_in_ec == assigned_proc) { - bool alloced_lhs_buf{false}; LHS_ElType* lhs_buf{nullptr}; /// get total buffer size for a given Proc size_t lhs_size = lhs_tensor.total_buf_size(i); if(lhs_size <= 0) continue; if(proc_me_in_ec == pg_lhs_in_ec[i]) { - lhs_buf = lhs_tensor.access_local_buf(); - alloced_lhs_buf = false; + lhs_buf = lhs_scratch.view(lhs_tensor.access_local_buf()); } else { - lhs_buf = new LHS_ElType[lhs_size]; - alloced_lhs_buf = true; + lhs_buf = lhs_scratch.owned(lhs_size); auto* lhs_mem_region = lhs_tensor.memory_region(); /// get all of lhs's buf at i-th proc to lhs_buf lhs_mem_region->get(Proc{i}, Offset{0}, Size{lhs_size}, lhs_buf); @@ -406,7 +328,6 @@ void GeneralFlatPlan::apply(const SetOp& s auto* lhs_mem_region = lhs_tensor.memory_region(); lhs_mem_region->put(Proc{i}, Offset{0}, Size{lhs_size}, lhs_buf); } - if(alloced_lhs_buf) { delete[] lhs_buf; } } } } @@ -437,21 +358,21 @@ void GeneralLHSPlan::apply(const SetOp& se LabelLoopNest loop_nest{lhs_lt.labels()}; + // Reused (grow-only) scratch for the non-local / view case; RAII-owned. + internal::BlockScratch lhs_scratch; + auto lambda = [&](const IndexVector& l_blockid) { auto [lhs_proc, lhs_offset] = ldist.locate(l_blockid); auto lhs_blocksize = lhs_tensor.block_size(l_blockid); auto lhs_blockdims = lhs_tensor.block_dims(l_blockid); LHS_ElType* lhs_buf{nullptr}; - bool lhs_alloced{false}; if(proc_lhs_to_ec[lhs_proc.value()] == proc_me_in_ec && lhs_tensor.kind() != TensorBase::TensorKind::view) { - lhs_buf = lhs_tensor.access_local_buf() + lhs_offset.value(); - lhs_alloced = false; + lhs_buf = lhs_scratch.view(lhs_tensor.access_local_buf() + lhs_offset.value()); } else { - lhs_buf = new LHS_ElType[lhs_blocksize]; - lhs_alloced = true; + lhs_buf = lhs_scratch.owned(lhs_blocksize); span lhs_span{lhs_buf, lhs_blocksize}; lhs_tensor.get(l_blockid, lhs_span); } @@ -465,7 +386,6 @@ void GeneralLHSPlan::apply(const SetOp& se span lhs_span{lhs_buf, lhs_blocksize}; lhs_tensor.put(l_blockid, lhs_span); } - if(lhs_alloced) { delete[] lhs_buf; } }; Proc round_robin_counter = 0; diff --git a/src/tamm/strong_num.hpp b/src/tamm/strong_num.hpp index 9b804ba5b..7d584ab1b 100644 --- a/src/tamm/strong_num.hpp +++ b/src/tamm/strong_num.hpp @@ -1,8 +1,12 @@ // Copyright 2016 Pacific Northwest National Laboratory +// C++20 modernization: concepts, operator<=>, [[nodiscard]], requires-clauses. #pragma once -#include +#include +#include // std::strong_ordering, operator<=> +#include // std::integral, std::floating_point +#include // std::hash #include #include #include @@ -18,7 +22,16 @@ namespace tamm { * @todo Add debug mode checks for overflow and underflow */ -#define DEBUG_STRONGNUM +// --------------------------------------------------------------------------- +// Concept: StrongNumeric +// Constrains what underlying types StrongNum may wrap. +// --------------------------------------------------------------------------- +template +concept StrongNumeric = std::integral || std::floating_point; + +// --------------------------------------------------------------------------- +// checked_cast: narrow cast with optional debug assertion. +// --------------------------------------------------------------------------- /** * @brief Return input value in the desired target type after possibly checking @@ -29,18 +42,14 @@ namespace tamm { * @param s Value being type cast * @return Value after type casting */ -template::value>, - typename = std::enable_if_t::value>, - typename = std::enable_if_t::value>> -constexpr Target checked_cast(Source s) { -#if defined(DEBUG_STRONGNUM) +template + requires (!std::is_same_v) +constexpr Target checked_cast(Source s) noexcept(false) { auto r = static_cast(s); - assert(static_cast(r) == s); - return r; -#else - return static_cast(s); +#if defined(TAMM_DEBUG_STRONGNUM) + assert(static_cast(r) == s && "checked_cast: narrowing lost data"); #endif + return r; } /** @@ -50,14 +59,8 @@ constexpr Target checked_cast(Source s) { * @param s input value * @return same input value */ -template::value>> -constexpr T checked_cast(T s) { - return s; -} - -// template Target strongnum_cast(Source s) { -// return checked_cast(s); -// } +template +constexpr T checked_cast(T s) noexcept { return s; } /** * @brief Strongly typed wrapper for a numeric type. @@ -74,10 +77,16 @@ constexpr T checked_cast(T s) { * Checked casts are to be used to convert between types and possibly check the * conversions in debug mode. * + * C++20 changes vs. original: + * - StrongNumeric concept replaces enable_if arithmetic checks. + * - 12 comparison operators replaced by operator<=> (spaceship). + * - [[nodiscard]] applied to value() and all binary operators. + * - requires-clause replaces enable_if on constructor. + * * @tparam Space Unique type name - * @tparam T Numeric typed being wrapper + * @tparam T Numeric type being wrapped */ -template +template struct StrongNum { /** * @brief Type of wrapper number @@ -91,244 +100,161 @@ struct StrongNum { */ using NumType = StrongNum; - StrongNum() = default; - StrongNum(const StrongNum&) = default; - StrongNum& operator=(const StrongNum&) = default; - ~StrongNum() = default; - - template::value>, - typename = std::enable_if_t::value>> - StrongNum(const T2 v1): v{checked_cast(v1)} {} - - template::value>> - NumType& operator=(T2 t) { + // ---- Lifecycle -------------------------------------------------------- + StrongNum() = default; + StrongNum(const StrongNum&) = default; + StrongNum& operator=(const StrongNum&) = default; + ~StrongNum() = default; + + /// Implicit construction from any compatible arithmetic type. + /// + /// NOTE: this is intentionally NOT explicit. TAMM relies pervasively on + /// implicit conversions such as `Proc p = GA_Nodeid();`, `Offset off = 0;`, + /// and `Size sz = block_size(...)`. Making it explicit breaks hundreds of + /// call sites, so the original (implicit) behaviour is preserved here. + template + requires std::is_convertible_v + constexpr StrongNum(T2 v1) noexcept : v{checked_cast(v1)} {} + + // ---- Assignment from raw arithmetic ----------------------------------- + template + NumType& operator=(T2 t) noexcept { v = checked_cast(t); return *this; } - template::value>> - NumType& operator+=(T2 t) { - v += checked_cast(t); - return *this; - } - - NumType& operator+=(NumType d) { - v += d.v; - return *this; - } - - template::value>> - NumType& operator-=(T2 t) { - v -= checked_cast(t); - return *this; - } - - NumType& operator-=(NumType d) { - v -= d.v; - return *this; - } - - template::value>> - NumType& operator*=(T2 t) { - v *= checked_cast(t); - return *this; - } - - NumType& operator*=(NumType d) { - v *= d.v; - return *this; - } - - template::value>> - NumType& operator/=(T2 t) { - v /= checked_cast(t); - return *this; - } - - NumType& operator/=(NumType d) { - v /= d.v; - return *this; - } - - template::value>> - NumType& operator^=(T2 t) { - v ^= checked_cast(t); - return *this; - } - - NumType& operator^=(NumType d) { - v ^= d.v; - return *this; - } - - NumType& operator++() { - v += 1; - return *this; - } - - NumType operator++(int) { - NumType ret{*this}; - v += 1; - return ret; - } - - NumType& operator--() { - v -= 1; - return *this; - } - - NumType operator--(int) { - NumType ret{*this}; - v -= 1; - return ret; - } - - NumType operator+(NumType d) const { return v + d.v; } - NumType operator-(NumType d) const { return v - d.v; } - NumType operator*(NumType d) const { return v * d.v; } - NumType operator/(NumType d) const { return v / d.v; } - NumType operator%(NumType d) const { return v % d.v; } - - template::value>> - NumType operator+(T2 t) const { - return v + checked_cast(t); - } - - template::value>> - NumType operator-(T2 t) const { - return v - checked_cast(t); - } - - template::value>> - NumType operator*(T2 t) const { - return v * checked_cast(t); - } - - template::value>> - NumType operator/(T2 t) const { - return v / checked_cast(t); - } - - template::value>> - NumType operator%(T2 t) const { - return v % checked_cast(t); - } - - bool operator==(NumType d) const { return v == d.v; } - bool operator!=(NumType d) const { return v != d.v; } - bool operator>=(NumType d) const { return v >= d.v; } - bool operator<=(NumType d) const { return v <= d.v; } - bool operator>(NumType d) const { return v > d.v; } - bool operator<(NumType d) const { return v < d.v; } - - template - bool operator==(T2 t) const { - return v == checked_cast(t); - } - - template - bool operator!=(T2 t) const { - return v != checked_cast(t); - } - - template - bool operator>=(T2 t) const { - return v >= checked_cast(t); - } - - template - bool operator<=(T2 t) const { - return v <= checked_cast(t); - } - - template - bool operator>(T2 t) const { - return v > checked_cast(t); - } - - template - bool operator<(T2 t) const { - return v < checked_cast(t); - } - - T value() const { return v; } - T& value() { return v; } - - // template T1 value() const { return checked_cast(v); } + // ---- Compound assignment (NumType operand) ---------------------------- + NumType& operator+=(NumType d) noexcept { v += d.v; return *this; } + NumType& operator-=(NumType d) noexcept { v -= d.v; return *this; } + NumType& operator*=(NumType d) noexcept { v *= d.v; return *this; } + NumType& operator/=(NumType d) noexcept { v /= d.v; return *this; } + NumType& operator%=(NumType d) noexcept { v %= d.v; return *this; } + NumType& operator^=(NumType d) noexcept { v ^= d.v; return *this; } + + // ---- Compound assignment (raw arithmetic operand) --------------------- + template + NumType& operator+=(T2 t) noexcept { v += checked_cast(t); return *this; } + template + NumType& operator-=(T2 t) noexcept { v -= checked_cast(t); return *this; } + template + NumType& operator*=(T2 t) noexcept { v *= checked_cast(t); return *this; } + template + NumType& operator/=(T2 t) noexcept { v /= checked_cast(t); return *this; } + template + NumType& operator^=(T2 t) noexcept { v ^= checked_cast(t); return *this; } + + // ---- Increment / decrement ------------------------------------------- + NumType& operator++() noexcept { v += 1; return *this; } + NumType operator++(int) noexcept { NumType ret{*this}; v += 1; return ret; } + NumType& operator--() noexcept { v -= 1; return *this; } + NumType operator--(int) noexcept { NumType ret{*this}; v -= 1; return ret; } + + // ---- Binary arithmetic (NumType operands) ---------------------------- + [[nodiscard]] NumType operator+(NumType d) const noexcept { return NumType{v + d.v}; } + [[nodiscard]] NumType operator-(NumType d) const noexcept { return NumType{v - d.v}; } + [[nodiscard]] NumType operator*(NumType d) const noexcept { return NumType{v * d.v}; } + [[nodiscard]] NumType operator/(NumType d) const noexcept { return NumType{v / d.v}; } + [[nodiscard]] NumType operator%(NumType d) const noexcept { return NumType{v % d.v}; } + + // ---- Binary arithmetic (raw arithmetic operands) --------------------- + template + [[nodiscard]] NumType operator+(T2 t) const noexcept { return NumType{v + checked_cast(t)}; } + template + [[nodiscard]] NumType operator-(T2 t) const noexcept { return NumType{v - checked_cast(t)}; } + template + [[nodiscard]] NumType operator*(T2 t) const noexcept { return NumType{v * checked_cast(t)}; } + template + [[nodiscard]] NumType operator/(T2 t) const noexcept { return NumType{v / checked_cast(t)}; } + template + [[nodiscard]] NumType operator%(T2 t) const noexcept { return NumType{v % checked_cast(t)}; } + + // ---- Comparison (same-type) ------------------------------------------- + // Defaulted operator== AND operator<=>. operator<=> alone does NOT + // synthesize ==/!=, so both are needed; == synthesizes != and <=> + // synthesizes < <= > >=. + [[nodiscard]] bool operator==(const NumType& d) const noexcept = default; + [[nodiscard]] auto operator<=>(const NumType& d) const noexcept = default; + + // Heterogeneous comparisons against raw arithmetic (non-defaulted). + template + [[nodiscard]] bool operator==(T2 t) const noexcept { return v == checked_cast(t); } + template + [[nodiscard]] auto operator<=>(T2 t) const noexcept { return v <=> checked_cast(t); } + + // ---- Value accessors -------------------------------------------------- + [[nodiscard]] T value() const noexcept { return v; } + T& value() noexcept { return v; } private: - T v; /**< Value wrapped by this object */ + T v{}; /**< Value wrapped by this object */ }; -template -inline StrongNum strongnum_cast(StrongNum v2) { - return {checked_cast(v2.value())}; -} - -template -inline bool operator==(T val, StrongNum snum) { - return snum == val; -} - -template -inline bool operator!=(T val, StrongNum snum) { - return snum != val; -} - -template -inline bool operator>=(T val, StrongNum snum) { - return snum <= val; +// --------------------------------------------------------------------------- +// strongnum_cast: cross-space conversion +// --------------------------------------------------------------------------- +template +[[nodiscard]] inline StrongNum +strongnum_cast(StrongNum v2) noexcept { + return StrongNum{checked_cast(v2.value())}; } -template -inline bool operator<=(T val, StrongNum snum) { - return snum >= val; +// --------------------------------------------------------------------------- +// Non-member operators: raw-lhs op StrongNum-rhs +// (Needed because the member operators only cover StrongNum-lhs.) +// --------------------------------------------------------------------------- +template +[[nodiscard]] inline bool operator==(T val, StrongNum sn) noexcept { + return sn == val; } - -template -inline bool operator>(T val, StrongNum snum) { - return snum < val; +template +[[nodiscard]] inline auto operator<=>(T val, StrongNum sn) noexcept { + return checked_cast(val) <=> sn.value(); } -template -inline bool operator<(T val, StrongNum snum) { - return snum > val; +template +[[nodiscard]] inline StrongNum operator-(T v, StrongNum sn) noexcept { + return StrongNum{checked_cast(v) - sn.value()}; } - -template -StrongNum operator-(T value, StrongNum snum) { - return checked_cast(value) - snum; +template +[[nodiscard]] inline StrongNum operator+(T v, StrongNum sn) noexcept { + return StrongNum{checked_cast(v) + sn.value()}; } - -template -StrongNum operator*(T value, StrongNum snum) { - return checked_cast(value) * snum; +template +[[nodiscard]] inline StrongNum operator*(T v, StrongNum sn) noexcept { + return StrongNum{checked_cast(v) * sn.value()}; } - -template -StrongNum operator/(T value, StrongNum snum) { - return checked_cast(value) / snum; +template +[[nodiscard]] inline StrongNum operator/(T v, StrongNum sn) noexcept { + return StrongNum{checked_cast(v) / sn.value()}; } - -template -StrongNum operator+(T value, StrongNum snum) { - return checked_cast(value) + snum; +template +[[nodiscard]] inline StrongNum operator%(T v, StrongNum sn) noexcept { + return StrongNum{checked_cast(v) % sn.value()}; } -template -StrongNum operator%(T value, StrongNum snum) { - return checked_cast(value) % snum; -} - -template +// --------------------------------------------------------------------------- +// I/O +// --------------------------------------------------------------------------- +template inline std::ostream& operator<<(std::ostream& os, StrongNum s) { return os << s.value(); } - -template +template inline std::istream& operator>>(std::istream& is, StrongNum& s) { is >> s.value(); return is; } } // namespace tamm + +// --------------------------------------------------------------------------- +// std::hash specialization so StrongNum can be used in unordered containers. +// --------------------------------------------------------------------------- +namespace std { +template +struct hash> { + [[nodiscard]] size_t operator()(tamm::StrongNum s) const noexcept { + return std::hash{}(s.value()); + } +}; +} // namespace std diff --git a/src/tamm/tamm_io.hpp b/src/tamm/tamm_io.hpp index 1e6c50b07..3ea3bd994 100644 --- a/src/tamm/tamm_io.hpp +++ b/src/tamm/tamm_io.hpp @@ -195,10 +195,10 @@ hid_t get_hdf5_dt() { else if constexpr(is_same_v) return H5T_NATIVE_FLOAT; else if constexpr(is_same_v) return H5T_NATIVE_DOUBLE; else if constexpr(is_same_v, T>) { - typedef struct { + struct complex_t { float re; /*real part*/ float im; /*imaginary part*/ - } complex_t; + }; hid_t complex_id = H5Tcreate(H5T_COMPOUND, sizeof(complex_t)); H5Tinsert(complex_id, "real", HOFFSET(complex_t, re), H5T_NATIVE_FLOAT); @@ -206,10 +206,10 @@ hid_t get_hdf5_dt() { return complex_id; } else if constexpr(is_same_v, T>) { - typedef struct { + struct complex_t { double re; /*real part*/ double im; /*imaginary part*/ - } complex_t; + }; hid_t complex_id = H5Tcreate(H5T_COMPOUND, sizeof(complex_t)); H5Tinsert(complex_id, "real", HOFFSET(complex_t, re), H5T_NATIVE_DOUBLE); diff --git a/src/tamm/tamm_utils.hpp b/src/tamm/tamm_utils.hpp index f3d2a39da..335a387b8 100644 --- a/src/tamm/tamm_utils.hpp +++ b/src/tamm/tamm_utils.hpp @@ -184,7 +184,7 @@ template void print_vector(std::vector vec, std::string filename = "") { std::stringstream tstring; for(size_t i = 0; i < vec.size(); i++) - tstring << i + 1 << "\t" << std::fixed << std::setprecision(12) << vec[i] << std::endl; + tstring << i + 1 << "\t" << std::fixed << std::setprecision(12) << vec[i] << '\n'; if(!filename.empty()) { std::ofstream tos(filename, std::ios::out); @@ -610,7 +610,7 @@ std::vector diagonal(LabeledTensor ltensor) { } } - int dsize = (int) dvec.size(); + int dsize = static_cast(dvec.size()); ec.pg().broadcast(&dsize, 0); if(ec.pg().rank() != 0) dvec.resize(dsize); ec.pg().broadcast(dvec.data(), dsize, 0); @@ -2044,14 +2044,14 @@ void print_dense_tensor(const Tensor& tensor, std::function& tensor, std::function& tensor, std::function friend bool operator==(const Tensor& lhs, const Tensor& rhs); - bool is_sparse() const { - for(auto& tis: tiled_index_spaces()) { - if(tis.is_dependent()) return true; - } - return false; + [[nodiscard]] bool is_sparse() const { + return std::ranges::any_of(tiled_index_spaces(), + [](const auto& tis) { return tis.is_dependent(); }); } - bool is_dense() const { return !is_sparse(); } + [[nodiscard]] bool is_dense() const { return !is_sparse(); } T* access_local_buf() { return impl_->access_local_buf(); } diff --git a/src/tamm/tensor_base.hpp b/src/tamm/tensor_base.hpp index 35391c232..d325d3946 100644 --- a/src/tamm/tensor_base.hpp +++ b/src/tamm/tensor_base.hpp @@ -125,44 +125,44 @@ class TensorBase { ExecutionContext* execution_context() const { return ec_; } - auto tindices() const { return block_indices_; } + [[nodiscard]] const std::vector& tindices() const { return block_indices_; } - TAMM_SIZE block_size(const IndexVector& blockid) const { + [[nodiscard]] TAMM_SIZE block_size(const IndexVector& blockid) const { size_t ret = 1; EXPECTS(blockid.size() == num_modes()); size_t rank = block_indices_.size(); for(size_t i = 0; i < rank; i++) { IndexVector dep_idx_vals{}; - if(dep_map_.find(i) != dep_map_.end()) { - for(const auto& pos: dep_map_.at(i)) { dep_idx_vals.push_back(blockid[pos]); } + if(auto it = dep_map_.find(i); it != dep_map_.end()) { + for(const auto& pos: it->second) { dep_idx_vals.push_back(blockid[pos]); } } ret *= block_indices_[i](dep_idx_vals).tile_size(blockid[i]); } return ret; } - std::vector block_dims(const IndexVector& blockid) const { + [[nodiscard]] std::vector block_dims(const IndexVector& blockid) const { std::vector ret; EXPECTS(blockid.size() == num_modes()); size_t rank = block_indices_.size(); for(size_t i = 0; i < rank; i++) { IndexVector dep_idx_vals{}; - if(dep_map_.find(i) != dep_map_.end()) { - for(const auto& pos: dep_map_.at(i)) { dep_idx_vals.push_back(blockid[pos]); } + if(auto it = dep_map_.find(i); it != dep_map_.end()) { + for(const auto& pos: it->second) { dep_idx_vals.push_back(blockid[pos]); } } ret.push_back(block_indices_[i](dep_idx_vals).tile_size(blockid[i])); } return ret; } - std::vector block_offsets(const IndexVector& blockid) const { + [[nodiscard]] std::vector block_offsets(const IndexVector& blockid) const { std::vector ret; EXPECTS(blockid.size() == num_modes()); size_t rank = num_modes(); for(size_t i = 0; i < rank; i++) { IndexVector dep_idx_vals{}; - if(dep_map_.find(i) != dep_map_.end()) { - for(const auto& pos: dep_map_.at(i)) { dep_idx_vals.push_back(blockid[pos]); } + if(auto it = dep_map_.find(i); it != dep_map_.end()) { + for(const auto& pos: it->second) { dep_idx_vals.push_back(blockid[pos]); } } ret.push_back(block_indices_[i](dep_idx_vals).tile_offset(blockid[i])); } @@ -239,12 +239,9 @@ class TensorBase { Spin spin_total() const { return spin_total_; } - bool is_dense() const { - bool result = true; - for(const auto& tis: block_indices_) { - if(tis.is_dependent()) { return false; } - } - return result; + [[nodiscard]] bool is_dense() const { + return std::ranges::none_of(block_indices_, + [](const auto& tis) { return tis.is_dependent(); }); } bool is_non_zero(const IndexVector& blockid) const { diff --git a/src/tamm/tensor_impl.hpp b/src/tamm/tensor_impl.hpp index ffd7cb004..0b81901db 100644 --- a/src/tamm/tensor_impl.hpp +++ b/src/tamm/tensor_impl.hpp @@ -9,9 +9,13 @@ #include "tamm/mem_profiler.hpp" #include "tamm/memory_manager_local.hpp" #include "tamm/tensor_base.hpp" +#include #include #include +#include +#include #include +#include #if defined(USE_UPCXX) #include #endif @@ -119,7 +123,7 @@ class TensorImpl: public TensorBase { // for(const auto& tis : t_spaces) { EXPECTS(tis.has_spin()); } - spin_mask_ = spin_mask; + spin_mask_ = std::move(spin_mask); has_spin_symmetry_ = true; has_spatial_symmetry_ = false; // spin_total_ = calculate_spin(); @@ -138,7 +142,7 @@ class TensorImpl: public TensorBase { // for(const auto& tlbl : t_labels) { // EXPECTS(tlbl.tiled_index_space().has_spin()); // } - spin_mask_ = spin_mask; + spin_mask_ = std::move(spin_mask); has_spin_symmetry_ = true; has_spatial_symmetry_ = false; // spin_total_ = calculate_spin(); @@ -162,15 +166,14 @@ class TensorImpl: public TensorBase { SpinMask spin_mask; size_t upper = spin_sizes[0]; size_t lower = spin_sizes.size() > 1 ? spin_sizes[1] : t_spaces.size() - upper; - size_t ignore = spin_sizes.size() > 2 ? spin_sizes[1] : t_spaces.size() - (upper + lower); + size_t ignore = spin_sizes.size() > 2 ? spin_sizes[2] : t_spaces.size() - (upper + lower); - for(size_t i = 0; i < upper; i++) { spin_mask.push_back(SpinPosition::upper); } + spin_mask.reserve(upper + lower + ignore); + spin_mask.insert(spin_mask.end(), upper, SpinPosition::upper); + spin_mask.insert(spin_mask.end(), lower, SpinPosition::lower); + spin_mask.insert(spin_mask.end(), ignore, SpinPosition::ignore); - for(size_t i = 0; i < lower; i++) { spin_mask.push_back(SpinPosition::lower); } - - for(size_t i = 0; i < ignore; i++) { spin_mask.push_back(SpinPosition::ignore); } - - spin_mask_ = spin_mask; + spin_mask_ = std::move(spin_mask); has_spin_symmetry_ = true; has_spatial_symmetry_ = false; // spin_total_ = calculate_spin(); @@ -207,15 +210,14 @@ class TensorImpl: public TensorBase { SpinMask spin_mask; size_t upper = spin_sizes[0]; size_t lower = spin_sizes.size() > 1 ? spin_sizes[1] : t_labels.size() - upper; - size_t ignore = spin_sizes.size() > 2 ? spin_sizes[1] : t_labels.size() - (upper + lower); - - for(size_t i = 0; i < upper; i++) { spin_mask.push_back(SpinPosition::upper); } + size_t ignore = spin_sizes.size() > 2 ? spin_sizes[2] : t_labels.size() - (upper + lower); - for(size_t i = 0; i < lower; i++) { spin_mask.push_back(SpinPosition::lower); } + spin_mask.reserve(upper + lower + ignore); + spin_mask.insert(spin_mask.end(), upper, SpinPosition::upper); + spin_mask.insert(spin_mask.end(), lower, SpinPosition::lower); + spin_mask.insert(spin_mask.end(), ignore, SpinPosition::ignore); - for(size_t i = 0; i < ignore; i++) { spin_mask.push_back(SpinPosition::ignore); } - - spin_mask_ = spin_mask; + spin_mask_ = std::move(spin_mask); has_spin_symmetry_ = true; has_spatial_symmetry_ = false; // spin_total_ = calculate_spin(); @@ -266,11 +268,11 @@ class TensorImpl: public TensorBase { // get memory profiler instance auto& memprof = MemProfiler::instance(); - auto defd = ec->get_default_distribution(); - Distribution* distribution = - ec->distribution(defd->get_tensor_base(), defd->get_dist_proc()); // defd->kind()); - // Distribution* distribution = - // ec->distribution(defd.tensor_base(), nproc ); + // Owning local factory pointers; unique_ptr guarantees release even if + // clone() below throws (RAII, replaces manual delete). + std::unique_ptr defd{ec->get_default_distribution()}; + std::unique_ptr distribution{ + ec->distribution(defd->get_tensor_base(), defd->get_dist_proc())}; // defd->kind()); MemoryManager* memory_manager = ec->memory_manager(); EXPECTS(distribution != nullptr); EXPECTS(memory_manager != nullptr); @@ -282,10 +284,6 @@ class TensorImpl: public TensorBase { distribution_ = std::shared_ptr(distribution->clone(this, memory_manager->pg().size())); } - - // Delete unused pointers - delete defd; - delete distribution; #if 0 auto rank = memory_manager->pg().rank(); auto buf_size = distribution_->buf_size(rank); @@ -335,7 +333,7 @@ class TensorImpl: public TensorBase { if(!is_non_zero(idx_vec)) { Size size = block_size(idx_vec); EXPECTS(size <= buff_span.size()); - for(size_t i = 0; i < size; i++) { buff_span[i] = (T) 0; } + std::fill_n(buff_span.begin(), size.value(), T{0}); return; } @@ -363,7 +361,7 @@ class TensorImpl: public TensorBase { if(!is_non_zero(idx_vec)) { Size size = block_size(idx_vec); EXPECTS(size <= buff_span.size()); - for(size_t i = 0; i < size; i++) { buff_span[i] = (T) 0; } + std::fill_n(buff_span.begin(), size.value(), T{0}); return; } @@ -598,12 +596,12 @@ class LambdaTensorImpl: public TensorImpl { * @param [in] lambda a function for constructing the Tensor */ LambdaTensorImpl(const TiledIndexSpaceVec& tis_vec, Func lambda): - TensorImpl(tis_vec), lambda_{lambda} { + TensorImpl(tis_vec), lambda_{std::move(lambda)} { setKind(TensorBase::TensorKind::lambda); } LambdaTensorImpl(const IndexLabelVec& til_vec, Func lambda): - TensorImpl(til_vec), lambda_{lambda} { + TensorImpl(til_vec), lambda_{std::move(lambda)} { setKind(TensorBase::TensorKind::lambda); } @@ -808,8 +806,9 @@ class DenseTensorImpl: public TensorImpl { ec_ = ec; const int ndims = num_modes(); - auto defd = ec->get_default_distribution(); - Distribution* distribution = ec->distribution(defd->get_tensor_base(), defd->get_dist_proc()); + std::unique_ptr defd{ec->get_default_distribution()}; + std::unique_ptr distribution{ + ec->distribution(defd->get_tensor_base(), defd->get_dist_proc())}; EXPECTS(distribution != nullptr); @@ -818,9 +817,6 @@ class DenseTensorImpl: public TensorImpl { distribution_ = std::shared_ptr(distribution->clone(this, ec->pg().size())); proc_grid_ = distribution_->proc_grid(); - delete defd; - delete distribution; - auto tis_dims = tindices(); std::vector> new_tiles(ndims); @@ -949,7 +945,9 @@ class DenseTensorImpl: public TensorImpl { tile_index++; } - if(local_nelems_ = tile_offsets[my_rank]) + // Assign then test for non-zero (parenthesised to make the intent + // explicit and silence -Wparentheses). + if((local_nelems_ = tile_offsets[my_rank]) != 0) for(int i = 4 - ndims; i < 4; ++i) local_buf_dims_.push_back(local_tiles_.back().lo[i] + local_tiles_.back().dim[i] - local_tiles_.front().lo[i]); @@ -996,7 +994,7 @@ class DenseTensorImpl: public TensorImpl { } if(is_bgd) { - nblocks = std::accumulate(nblock, nblock + ndims, (int) 1, std::multiplies()); + nblocks = std::reduce(nblock, nblock + ndims, 1, std::multiplies<>{}); int proclist_c[nblocks]; std::iota(proclist_c, proclist_c + nblocks, 0); GA_Set_restricted(ga_, proclist_c, nblocks); @@ -1416,7 +1414,7 @@ class DenseTensorImpl: public TensorImpl { int64_t len; NGA_Access_block_segment64(ga_, GA_Pgroup_nodeid(GA_Get_pgroup(ga_)), reinterpret_cast(&ptr), &len); - res = (size_t) len; + res = static_cast(len); #endif return res; } @@ -1496,7 +1494,7 @@ class ViewTensorImpl: public TensorImpl { // Ctors ViewTensorImpl() = default; ViewTensorImpl(Tensor ref_tensor, const TiledIndexSpaceVec& tis_vec, Func ref_map_func): - TensorImpl(tis_vec), ref_tensor_{ref_tensor}, ref_map_func_{ref_map_func} { + TensorImpl(tis_vec), ref_tensor_{ref_tensor}, ref_map_func_{std::move(ref_map_func)} { setKind(TensorBase::TensorKind::view); if(ref_tensor_.is_allocated()) { distribution_ = @@ -1508,7 +1506,7 @@ class ViewTensorImpl: public TensorImpl { } ViewTensorImpl(Tensor ref_tensor, const IndexLabelVec& labels, Func ref_map_func): - TensorImpl(labels), ref_tensor_{ref_tensor}, ref_map_func_{ref_map_func} { + TensorImpl(labels), ref_tensor_{ref_tensor}, ref_map_func_{std::move(ref_map_func)} { setKind(TensorBase::TensorKind::view); if(ref_tensor_.is_allocated()) { @@ -1524,9 +1522,9 @@ class ViewTensorImpl: public TensorImpl { CopyFunc get_func, CopyFunc put_func): TensorImpl(labels), ref_tensor_{ref_tensor}, - ref_map_func_{ref_map_func}, - get_func_{get_func}, - put_func_{put_func} { + ref_map_func_{std::move(ref_map_func)}, + get_func_{std::move(get_func)}, + put_func_{std::move(put_func)} { setKind(TensorBase::TensorKind::view); if(ref_tensor_.is_allocated()) { @@ -1793,8 +1791,9 @@ class TensorUnitTiled: public TensorImpl { EXPECTS(tensor_opt_.is_allocated()); if(!is_allocated()) { - auto defd = ec->get_default_distribution(); - Distribution* distribution = ec->distribution(defd->get_tensor_base(), defd->get_dist_proc()); + std::unique_ptr defd{ec->get_default_distribution()}; + std::unique_ptr distribution{ + ec->distribution(defd->get_tensor_base(), defd->get_dist_proc())}; MemoryManager* memory_manager = ec->memory_manager(); EXPECTS(distribution != nullptr); EXPECTS(memory_manager != nullptr); @@ -1805,8 +1804,8 @@ class TensorUnitTiled: public TensorImpl { EXPECTS(distribution_ != nullptr); - delete defd; - delete distribution; + defd.reset(); + distribution.reset(); delete memory_manager; auto eltype = tensor_element_type(); diff --git a/src/tamm/tensor_variant.hpp b/src/tamm/tensor_variant.hpp index fa9ec1e51..b0ac6ac72 100644 --- a/src/tamm/tensor_variant.hpp +++ b/src/tamm/tensor_variant.hpp @@ -4,6 +4,9 @@ #include "tamm/scheduler.hpp" #include "tamm/tensor.hpp" +#include +#include + namespace tamm { class TensorVariant { public: @@ -110,11 +113,10 @@ class TensorVariant { } size_t mem_size() const { - size_t result = el_size(); - auto tis_list = std::visit( - overloaded{[&](const auto& tensor) { return tensor.tiled_index_spaces(); }}, value_); - for(const auto& tis: tis_list) { result *= tis.max_num_indices(); } - return result; + auto tis_list = std::visit( + overloaded{[&](const auto& tensor) { return tensor.tiled_index_spaces(); }}, value_); + return std::transform_reduce(tis_list.begin(), tis_list.end(), el_size(), std::multiplies<>{}, + [](const auto& tis) { return tis.max_num_indices(); }); } bool is_allocated() const { diff --git a/src/tamm/tiled_index_space.hpp b/src/tamm/tiled_index_space.hpp index a5b869d41..e33040613 100644 --- a/src/tamm/tiled_index_space.hpp +++ b/src/tamm/tiled_index_space.hpp @@ -188,16 +188,14 @@ class TiledIndexSpace { */ TiledIndexSpace operator()(std::string id) const { if(id == "all") { return (*this); } - if(tiled_info_->tiled_named_subspaces_.find(id) == tiled_info_->tiled_named_subspaces_.end()) { - std::cerr << "Named sub-space " + id + " doesn't exist!" << std::endl; - } - if(tiled_info_->tiled_named_subspaces_.find(id) == tiled_info_->tiled_named_subspaces_.end()) { - std::ostringstream os; - os << "[TAMM ERROR] Named sub-space doesn't exist!\n" << __FILE__ << ":L" << __LINE__; - tamm_terminate(os.str()); - } - - return tiled_info_->tiled_named_subspaces_.at(id); + // Single lookup (C++20 if-init) instead of find()+find()+at(). + const auto& named = tiled_info_->tiled_named_subspaces_; + if(auto it = named.find(id); it != named.end()) { return it->second; } + std::cerr << "Named sub-space " + id + " doesn't exist!" << std::endl; + std::ostringstream os; + os << "[TAMM ERROR] Named sub-space doesn't exist!\n" << __FILE__ << ":L" << __LINE__; + tamm_terminate(os.str()); + return {}; // unreachable (tamm_terminate does not return) } /** @@ -210,12 +208,10 @@ class TiledIndexSpace { TiledIndexSpace operator()(const IndexVector& dep_idx_vec = {}) const { if(dep_idx_vec.empty()) { return (*this); } const auto& t_dep_map = tiled_info_->tiled_dep_map_; - EXPECTS(t_dep_map.find(dep_idx_vec) != t_dep_map.end()); - // if(t_dep_map.find(dep_idx_vec) == t_dep_map.end()){ - // return TiledIndexSpace{IndexSpace{{}}}; - // } - - return t_dep_map.at(dep_idx_vec); + // Single lookup instead of find() (in EXPECTS) + at(). + auto it = t_dep_map.find(dep_idx_vec); + EXPECTS(it != t_dep_map.end()); + return it->second; } /** @@ -230,14 +226,11 @@ class TiledIndexSpace { lookup_dependent_space(const IndexVector& dep_idx_vec = {}) const { if(dep_idx_vec.empty()) { return {(*this), true}; } const auto& t_dep_map = tiled_info_->tiled_dep_map_; - if(t_dep_map.find(dep_idx_vec) == t_dep_map.end()) { - return {TiledIndexSpace{IndexSpace{IndexVector{}}}, false}; - } - else if(t_dep_map.at(dep_idx_vec) == TiledIndexSpace{IndexSpace{IndexVector{}}}) { - return {t_dep_map.at(dep_idx_vec), false}; - } - - return {t_dep_map.at(dep_idx_vec), true}; + // Single lookup instead of find() + up to three at() calls. + auto it = t_dep_map.find(dep_idx_vec); + if(it == t_dep_map.end()) { return {TiledIndexSpace{IndexSpace{IndexVector{}}}, false}; } + const bool is_empty = (it->second == TiledIndexSpace{IndexSpace{IndexVector{}}}); + return {it->second, !is_empty}; } /** @@ -1395,7 +1388,8 @@ class TiledIndexSpace { TiledIndexSpaceInfo object*/ std::shared_ptr parent_tis_; /**< Shared pointer to the parent TiledIndexSpace object*/ - size_t hash_value_; + size_t hash_value_{0}; // 0 for a default-constructed TiledIndexSpace; avoids + // reading an indeterminate value in hash()/is_identical. bool is_dense_subspace_ = false; /** diff --git a/src/tamm/types.hpp b/src/tamm/types.hpp index b4506f59a..900d4680b 100644 --- a/src/tamm/types.hpp +++ b/src/tamm/types.hpp @@ -1,4 +1,9 @@ // Copyright 2016 Pacific Northwest National Laboratory +// C++20 modernization: removed pre-C++17 apply shim (C++20 guarantees +// std::apply), inline constexpr maxrank, [[nodiscard]] on queries. +// IndexVector remains std::vector (see note below) because it doubles +// as the storage for full index/tile lists of an IndexSpace, which exceed +// maxrank. Rank-bounded, allocation-sensitive vectors use TensorVec. #pragma once @@ -7,52 +12,71 @@ #include "tamm/boundvec.hpp" #include "tamm/errors.hpp" #include "tamm/strong_num.hpp" +#include // std::atomic (thread-safe make_label) +#include // std::bit_cast (C++20) #include +#include #include #include +#include // C++20 #if defined(USE_UPCXX) #include #endif -//#include - namespace tamm { -// Free functions -#if __cplusplus < 201703L +// --------------------------------------------------------------------------- +// Internal utilities +// --------------------------------------------------------------------------- namespace internal { -template -constexpr decltype(auto) apply_impl(F&& f, Tuple&& t, std::index_sequence) { - return f(std::get(std::forward(t))...); -} -} // namespace internal -template -constexpr decltype(auto) apply(F&& f, Tuple&& t) { - return internal::apply_impl( - std::forward(f), std::forward(t), - std::make_index_sequence>::value>{}); +/// Mix hash of v into seed (boost-style hash_combine). +template +void hash_combine(size_t& seed, const T& v) { + constexpr size_t magic = 0x9e3779b9ULL; + seed ^= std::hash{}(v) + magic + (seed << 6) + (seed >> 2); } -#endif -// End Free functions -namespace internal { -template -void hash_combine(size_t& seed, T const& v) { - seed ^= std::hash{}(v) + 0x9e3779b9 + (seed << 6) + (seed >> 2); +/// Fold a parameter pack into a vector (C++17 fold-expression, kept here). +template +void unfold_vec(std::vector& vec, Args&&... args) { + static_assert((std::is_constructible_v && ...)); + (vec.push_back(std::forward(args)), ...); } + } // namespace internal -using TAMM_SIZE = uint64_t; -// IndexSpace related type definitions +// --------------------------------------------------------------------------- +// Fundamental scalar / index aliases +// --------------------------------------------------------------------------- +using TAMM_SIZE = uint64_t; using Index = uint32_t; -using IndexVector = std::vector; -using IndexIterator = std::vector::const_iterator; using Tile = uint32_t; -// DAG related Hash -using HashData = uint64_t; +using HashData = uint64_t; using StringLabelVec = std::vector; +// TensorRank and the compile-time maximum tensor rank (used by the BoundVec +// aliases TensorVec/BlockDimVec/PermVec below). +using TensorRank = size_t; +inline constexpr TensorRank maxrank{8}; + +// --------------------------------------------------------------------------- +// IndexVector: dynamic block-id / index-list vector. +// +// NOTE: This is std::vector, NOT a fixed-capacity BoundVec. Although +// a *block id* never exceeds maxrank entries, IndexVector is reused throughout +// TAMM to store full index lists and tile-offset arrays of an IndexSpace / +// TiledIndexSpace (see index_space*.hpp, tiled_index_space.hpp, range.hpp), +// which routinely hold thousands of elements. It also relies on std::vector +// API (data(), erase(), range insert(), (count,value) ctor) at many call +// sites. A fixed-capacity container here overflows / fails to compile. +// +// Rank-bounded, allocation-sensitive helpers use TensorVec/BlockDimVec +// (BoundVec) instead — see below. +// --------------------------------------------------------------------------- +using IndexVector = std::vector; +using IndexIterator = IndexVector::const_iterator; + class TiledIndexSpace; using TiledIndexSpaceVec = std::vector; class TiledIndexLabel; @@ -61,38 +85,39 @@ using IndexLabelVec = std::vector; using Perm = int32_t; using PermVector = std::vector; -////////////////////////////////// -struct IrrepSpace; -using Irrep = StrongNum; -struct SpinSpace; -using Spin = StrongNum; -struct SpatialSpace; -using Spatial = StrongNum; -using TensorRank = size_t; -struct OffsetSpace; -using Offset = StrongNum; -struct BlockIndexSpace; -using BlockIndex = StrongNum; -struct ProcSpace; -using Proc = StrongNum; -struct SignSpace; -using Sign = StrongNum; +// --------------------------------------------------------------------------- +// StrongNum aliases +// --------------------------------------------------------------------------- +struct IrrepSpace; using Irrep = StrongNum; +struct SpinSpace; using Spin = StrongNum; +struct SpatialSpace; using Spatial = StrongNum; +struct OffsetSpace; using Offset = StrongNum; +struct BlockIndexSpace; using BlockIndex = StrongNum; +struct ProcSpace; using Proc = StrongNum; +struct SignSpace; using Sign = StrongNum; // these are typedefs for usability using Size = Offset; using BlockCount = BlockIndex; -using Label = int; // needs to support negative values - -using IntLabel = int32_t; // a simple integer label for indices +using Label = int; // needs to support negative values +using IntLabel = int32_t; // a simple integer label for indices using IntLabelVec = std::vector; +using SizeVec = std::vector; +using ProcGrid = std::vector; +using ProcList = std::vector; -using SizeVec = std::vector; -using ProcGrid = std::vector; -using ProcList = std::vector; +// BoundVec-based tensor dimension/permutation helpers +template +using TensorVec = BoundVec; +using BlockDimVec = TensorVec; +using PermVec = TensorVec; -enum class AllocationStatus { invalid, created, attached, deallocated, orphaned }; +// --------------------------------------------------------------------------- +// Enumerations (all scoped enum class) +// --------------------------------------------------------------------------- +enum class AllocationStatus : uint8_t { invalid, created, attached, deallocated, orphaned }; -enum class ElementType { +enum class ElementType : uint8_t { invalid, single_precision, double_precision, @@ -100,74 +125,48 @@ enum class ElementType { double_complex }; -enum class DistributionKind { invalid, nw, dense, simple_round_robin, view, unit_tile }; +enum class DistributionKind : uint8_t { + invalid, nw, dense, simple_round_robin, view, unit_tile +}; -enum class MemoryManagerKind { invalid, ga, local }; +enum class MemoryManagerKind : uint8_t { invalid, ga, local }; +// --------------------------------------------------------------------------- +// Element type queries +// --------------------------------------------------------------------------- template -constexpr ElementType tensor_element_type() { - return ElementType::invalid; -} - -template<> -constexpr ElementType tensor_element_type() { - return ElementType::double_precision; -} - -template<> -constexpr ElementType tensor_element_type() { - return ElementType::single_precision; -} - -template<> -constexpr ElementType tensor_element_type>() { - return ElementType::single_complex; -} - -template<> -constexpr ElementType tensor_element_type>() { - return ElementType::double_complex; -} - -static inline constexpr size_t element_size(ElementType eltype) { - size_t ret = 0; - switch(eltype) { - case ElementType::single_precision: ret = sizeof(float); break; - case ElementType::double_precision: ret = sizeof(double); break; - case ElementType::single_complex: ret = sizeof(std::complex); break; - case ElementType::double_complex: ret = sizeof(std::complex); break; - default: UNREACHABLE(); +constexpr ElementType tensor_element_type() noexcept { return ElementType::invalid; } +template<> constexpr ElementType tensor_element_type() noexcept { return ElementType::double_precision; } +template<> constexpr ElementType tensor_element_type() noexcept { return ElementType::single_precision; } +template<> constexpr ElementType tensor_element_type>() noexcept { return ElementType::single_complex; } +template<> constexpr ElementType tensor_element_type>() noexcept { return ElementType::double_complex; } + +[[nodiscard]] static inline constexpr size_t element_size(ElementType eltype) noexcept { + switch (eltype) { + case ElementType::single_precision: return sizeof(float); + case ElementType::double_precision: return sizeof(double); + case ElementType::single_complex: return sizeof(std::complex); + case ElementType::double_complex: return sizeof(std::complex); + default: return 0; } - return ret; } -const TensorRank maxrank{8}; - -template -using TensorVec = BoundVec; - -using BlockDimVec = TensorVec; -// using DimTypeVec = TensorVec; -using PermVec = TensorVec; - -/////////// - -/////////////////// - +// --------------------------------------------------------------------------- +// Remaining enumerations +// --------------------------------------------------------------------------- using RangeValue = int64_t; -// enum class IndexSpaceType { mo, mso, ao, aso, aux }; -enum class SpinPosition { ignore, upper, lower }; -enum class IndexPosition { upper, lower, neither }; - -enum class SpinType { ao_spin, mo_spin }; - -enum class ExecutionHW { CPU, GPU, DEFAULT }; - -enum class ReduceOp { min, max, sum, maxloc, minloc }; +enum class SpinPosition : uint8_t { ignore, upper, lower }; +enum class IndexPosition : uint8_t { upper, lower, neither }; +enum class SpinType : uint8_t { ao_spin, mo_spin }; +enum class ExecutionHW : uint8_t { CPU, GPU, DEFAULT }; +enum class ReduceOp : uint8_t { min, max, sum, maxloc, minloc }; using SpinMask = std::vector; +// --------------------------------------------------------------------------- +// Runtime data handles (UPC++ vs GA) +// --------------------------------------------------------------------------- #if defined(USE_UPCXX) using rtDataHandlePtr = upcxx::future<>*; using rtDataHandle = upcxx::future<>; @@ -182,7 +181,7 @@ class DataCommunicationHandle { ~DataCommunicationHandle() = default; void waitForCompletion() { - if(!getCompletionStatus()) { + if (!getCompletionStatus()) { #if defined(USE_UPCXX) data_handle_.wait(); #else @@ -191,134 +190,105 @@ class DataCommunicationHandle { setCompletionStatus(); } } - void setCompletionStatus() { status_ = true; } - void resetCompletionStatus() { status_ = false; } - - bool getCompletionStatus() { - /* - if(status_ == false) - status_ = NGA_NbTest(&data_handle_); - */ - - return status_; - } - rtDataHandlePtr getDataHandlePtr() { return &data_handle_; } + void setCompletionStatus() noexcept { status_ = true; } + void resetCompletionStatus() noexcept { status_ = false; } + [[nodiscard]] bool getCompletionStatus() const noexcept { return status_; } + rtDataHandlePtr getDataHandlePtr() noexcept { return &data_handle_; } rtDataHandle data_handle_; - private: bool status_{true}; }; using DataCommunicationHandlePtr = DataCommunicationHandle*; -////////////////// - -// namespace SpinType { -// const Spin alpha{1}; -// const Spin beta{2}; -// }; // namespace SpinType - +// --------------------------------------------------------------------------- +// MPI type helpers +// --------------------------------------------------------------------------- #if !defined(USE_UPCXX) template -static inline MPI_Datatype mpi_type() { +[[nodiscard]] static inline MPI_Datatype mpi_type() { using std::is_same_v; - - if constexpr(is_same_v) return MPI_INT; - else if constexpr(is_same_v) return MPI_C_BOOL; - else if constexpr(is_same_v) return MPI_CHAR; - else if constexpr(is_same_v) return MPI_INT64_T; - else if constexpr(is_same_v) return MPI_UNSIGNED; - else if constexpr(is_same_v) return MPI_UNSIGNED_LONG; - else if constexpr(is_same_v) return MPI_FLOAT; - else if constexpr(is_same_v) return MPI_DOUBLE; - else if constexpr(is_same_v, T>) return MPI_COMPLEX; - else if constexpr(is_same_v, T>) return MPI_DOUBLE_COMPLEX; - else NOT_IMPLEMENTED(); // unhandled type + if constexpr (is_same_v) return MPI_INT; + else if constexpr (is_same_v) return MPI_C_BOOL; + else if constexpr (is_same_v) return MPI_CHAR; + else if constexpr (is_same_v) return MPI_INT64_T; + else if constexpr (is_same_v) return MPI_UNSIGNED; + else if constexpr (is_same_v) return MPI_UNSIGNED_LONG; + else if constexpr (is_same_v) return MPI_FLOAT; + else if constexpr (is_same_v) return MPI_DOUBLE; + else if constexpr (is_same_v,T>) return MPI_COMPLEX; + else if constexpr (is_same_v,T>) return MPI_DOUBLE_COMPLEX; + else NOT_IMPLEMENTED(); } template -static inline MPI_Datatype mpi_type_loc() { +[[nodiscard]] static inline MPI_Datatype mpi_type_loc() { using std::is_same_v; - - if constexpr(is_same_v) return MPI_2INT; - else if constexpr(is_same_v) return MPI_2REAL; - else if constexpr(is_same_v) return MPI_2DOUBLE_PRECISION; - else NOT_IMPLEMENTED(); // unhandled type + if constexpr (is_same_v) return MPI_2INT; + else if constexpr (is_same_v) return MPI_2REAL; + else if constexpr (is_same_v) return MPI_2DOUBLE_PRECISION; + else NOT_IMPLEMENTED(); } -static inline MPI_Op mpi_op(ReduceOp rop) { - if(rop == ReduceOp::min) return MPI_MIN; - else if(rop == ReduceOp::max) return MPI_MAX; - else if(rop == ReduceOp::sum) return MPI_SUM; - else if(rop == ReduceOp::minloc) return MPI_MINLOC; - else if(rop == ReduceOp::maxloc) return MPI_MAXLOC; - else NOT_IMPLEMENTED(); // unhandled op +[[nodiscard]] static inline MPI_Op mpi_op(ReduceOp rop) { + switch (rop) { + case ReduceOp::min: return MPI_MIN; + case ReduceOp::max: return MPI_MAX; + case ReduceOp::sum: return MPI_SUM; + case ReduceOp::minloc: return MPI_MINLOC; + case ReduceOp::maxloc: return MPI_MAXLOC; + default: NOT_IMPLEMENTED(); + } } #endif -namespace internal { -template -void unfold_vec(std::vector& v, Args&&... args) { - static_assert((std::is_constructible_v && ...)); - (v.push_back(std::forward(args)), ...); -} -} // namespace internal - +// --------------------------------------------------------------------------- +// make_label: thread-safe monotonic label generator +// --------------------------------------------------------------------------- inline Label make_label() { - static Label lbl = 0; - return lbl++; + static std::atomic