// -*- C++ -*- #ifndef __multi_array_operation_h__ #define __multi_array_operation_h__ #include #include #include #include #include #include "base_def.h" #include "mbase_def.h" #include "block/block.h" #include "operation_utils.h" #include "ranges/rheader.h" #include "pack_num.h" #include "ranges/index_info.h" namespace MultiArrayTools { namespace { using namespace MultiArrayHelper; } template Block makeBlock(const T* vec, size_t stepSize, size_t blockSize); template MBlock makeBlock(T* vec, size_t stepSize, size_t blockSize); // dont use this for now !! template std::shared_ptr seekBlockIndex(std::shared_ptr ownIdx, const OpClass& second); template const IndexInfo* seekBlockIndex(const IndexInfo* ownII, const OpClass& second); template class OperationTemplate { public: OperationClass& THIS() { return static_cast(*this); } const OperationClass& THIS() const { return static_cast(*this); } template auto operator+(const Second& in) const -> Operation,OperationClass,Second>; template auto operator-(const Second& in) const -> Operation,OperationClass,Second>; template auto operator*(const Second& in) const -> Operation,OperationClass,Second>; template auto operator/(const Second& in) const -> Operation,OperationClass,Second>; template auto c(std::shared_ptr& ind) const -> Contraction; private: friend OperationClass; OperationTemplate() = default; }; template class OperationMaster { public: typedef T value_type; typedef OperationBase OB; typedef ContainerRange CRange; typedef typename MultiRange::IndexType IndexType; typedef MBlock bType; OperationMaster(MutableMultiArrayBase& ma, const OpClass& second, std::shared_ptr& index); OperationMaster(MutableMultiArrayBase& ma, const OpClass& second, std::shared_ptr& index, const IndexInfo* blockIndex); MBlock& get(); const Block& get() const; std::vector block(const IndexInfo* blockIndex, bool init = false) const; const OperationMaster& block() const; protected: std::shared_ptr mkIndex(std::shared_ptr& index); void performAssignment(std::intptr_t blockIndexNum); OpClass const& mSecond; MutableMultiArrayBase& mArrayRef; std::shared_ptr mIndex; IndexInfo mIInfo; mutable bType mBlock; }; template class ConstOperationRoot : /*public OperationBase,*/ public OperationTemplate > { public: typedef T value_type; typedef OperationBase OB; typedef OperationTemplate > OT; typedef ContainerRange CRange; typedef typename CRange::IndexType IndexType; typedef Block bType; ConstOperationRoot(const MultiArrayBase& ma, const std::shared_ptr&... indices); const Block& get() const; std::vector block(const IndexInfo* blockIndex, bool init = false) const; const ConstOperationRoot& block() const; protected: std::shared_ptr mkIndex(const MultiArrayBase& ma, const std::shared_ptr&... indices); MultiArrayBase const& mArrayRef; std::shared_ptr mIndex; IndexInfo mIInfo; mutable bType mBlock; }; template class OperationRoot : public OperationTemplate > { public: typedef T value_type; typedef OperationBase OB; typedef OperationTemplate > OT; typedef ContainerRange CRange; typedef typename CRange::IndexType IndexType; typedef MBlock bType; OperationRoot(MutableMultiArrayBase& ma, const std::shared_ptr&... indices); template OperationMaster operator=(const OpClass& in); const MBlock& get() const; MBlock& get(); OperationRoot& set(const IndexInfo* blockIndex); std::vector block(const IndexInfo* blockIndex, bool init = false) const; const OperationRoot& block() const; protected: std::shared_ptr mkIndex(const MultiArrayBase& ma, const std::shared_ptr&... indices); MutableMultiArrayBase& mArrayRef; std::shared_ptr mIndex; IndexInfo mIInfo; mutable bType mBlock; const IndexInfo* mBlockII; // predefine to save time }; template class Operation : public OperationTemplate > { public: typedef T value_type; typedef OperationBase OB; typedef OperationTemplate > OT; typedef OpFunction F; typedef BlockResult bType; Operation(const Ops&... ops); const BlockResult& get() const; std::vector block(const IndexInfo* blockIndex, bool init = false) const; const Operation& block() const; protected: std::tuple mOps; mutable bType mRes; }; template class Contraction : public OperationTemplate > { public: typedef T value_type; typedef OperationTemplate > OT; typedef BlockResult bType; Contraction(const Op& op, std::shared_ptr ind); const BlockResult& get() const; std::vector block(const IndexInfo* blockIndex, bool init = false) const; const Contraction& block() const; protected: const Op& mOp; std::shared_ptr mInd; mutable bType mRes; }; } /* ========================= * * --- TEMPLATE CODE --- * * ========================= */ namespace MultiArrayTools { namespace { using namespace MultiArrayHelper; } template Block makeBlock(const T* vec, size_t stepSize, size_t blockSize) { return Block(vec, 0, blockSize, stepSize); } template MBlock makeBlock(T* vec, size_t stepSize, size_t blockSize) { return MBlock(vec, 0, blockSize, stepSize); } // dont use this for now !! template std::shared_ptr seekBlockIndex(std::shared_ptr ownIdx, const OpClass& second) { assert(0); // dont use this for now !! std::vector > ivec; seekIndexInst(ownIdx, ivec); std::map, std::vector > mp; for(auto& xx: ivec){ mp[xx] = second.block(xx); } // seek minimal number of VALUEs => guarantees absence of conflicting blocks minimizeAppearanceOfType(mp, BlockType::VALUE); // seek mininmal number of SPLITs => maximal vectorization possible minimizeAppearanceOfType(mp, BlockType::SPLIT); return mp.begin()->first; } template const IndexInfo* seekBlockIndex(const IndexInfo* ownII, const OpClass& second) { const IndexInfo* ii = ownII; while(ii->type() == IndexType::CONT or ii->type() == IndexType::MULTI){ ii = ii->getPtr(ii->dim()-1); } return ii; } /*************************** * OperationTemplate * ***************************/ template template auto OperationTemplate::operator+(const Second& in) const -> Operation,OperationClass,Second> { return Operation,OperationClass,Second>(THIS(), in); } template template auto OperationTemplate::operator-(const Second& in) const -> Operation,OperationClass,Second> { return Operation,OperationClass,Second>(THIS(), in); } template template auto OperationTemplate::operator*(const Second& in) const -> Operation,OperationClass,Second> { return Operation,OperationClass,Second>(THIS(), in); } template template auto OperationTemplate::operator/(const Second& in) const -> Operation,OperationClass,Second> { return Operation,OperationClass,Second>(THIS(), in); } template template auto OperationTemplate::c(std::shared_ptr& ind) const -> Contraction { return Contraction(THIS(), ind); } /************************* * OperationMaster * *************************/ template OperationMaster:: OperationMaster(MutableMultiArrayBase& ma, const OpClass& second, std::shared_ptr& index) : mSecond(second), mArrayRef(ma), mIndex(mkIndex(index)), mIInfo(*mIndex) { auto blockIndex = seekBlockIndex( &mIInfo, second); std::intptr_t blockIndexNum = blockIndex->getPtrNum(); block(blockIndex, true); second.block(blockIndex, true); performAssignment(blockIndexNum); } template OperationMaster:: OperationMaster(MutableMultiArrayBase& ma, const OpClass& second, std::shared_ptr& index, const IndexInfo* blockIndex) : mSecond(second), mArrayRef(ma), mIndex(mkIndex(index)), mIInfo(*mIndex) { std::intptr_t blockIndexNum = blockIndex->getPtrNum(); second.block(blockIndex, true); performAssignment(blockIndexNum); } template std::shared_ptr::IndexType> OperationMaster:: mkIndex(std::shared_ptr& index) { MultiRangeFactory mrf( index->range() ); std::shared_ptr > mr = std::dynamic_pointer_cast >( mrf.create() ); auto i = std::make_shared( mr->begin() ); (*i) = *index; return i; } template void OperationMaster::performAssignment(std::intptr_t blockIndexNum) { //static auto loop = mkLoop(mIndex, *this, mSecond); //loop(); for(*mIndex = 0; mIndex->pos() != mIndex->max(); mIndex->pp(blockIndexNum) ){ block(); get() = mSecond.get(); } } template MBlock& OperationMaster::get() { return mBlock; } template const Block& OperationMaster::get() const { return mBlock; } template std::vector OperationMaster::block(const IndexInfo* blockIndex, bool init) const { std::vector btv(1, getBlockType( &mIInfo, blockIndex, true) ); if(init){ mBlock = makeBlock(mArrayRef.data(), btv[0].second, blockIndex->max()); } return btv; } template const OperationMaster& OperationMaster::block() const { mBlock.set( mIndex->pos() ); return *this; } /**************************** * ConstOperationRoot * ****************************/ template ConstOperationRoot:: ConstOperationRoot(const MultiArrayBase& ma, const std::shared_ptr&... indices) : //OperationTemplate >(this), mArrayRef(ma), mIndex( mkIndex(ma,indices...) ), mIInfo(*mIndex) {} template std::shared_ptr::IndexType> ConstOperationRoot:: mkIndex(const MultiArrayBase& ma, const std::shared_ptr&... indices) { auto i = std::make_shared( ma.range() ); (*mIndex)(indices...); return i; } template const Block& ConstOperationRoot::get() const { block(); return mBlock; } template std::vector ConstOperationRoot::block(const IndexInfo* blockIndex, bool init) const { std::vector btv(1, getBlockType( &mIInfo, blockIndex, true) ); if(init){ mBlock = makeBlock(mArrayRef.data(), btv[0].second, blockIndex->max()); } return btv; } template const ConstOperationRoot& ConstOperationRoot::block() const { mBlock.set( (*mIndex)().pos() ); return *this; } /*********************** * OperationRoot * ***********************/ template OperationRoot:: OperationRoot(MutableMultiArrayBase& ma, const std::shared_ptr&... indices) : //OperationTemplate >(this), mArrayRef(ma), mIndex( mkIndex( ma, indices... ) ), mIInfo(*mIndex), mBlockII(nullptr) {} template std::shared_ptr::IndexType> OperationRoot:: mkIndex(const MultiArrayBase& ma, const std::shared_ptr&... indices) { auto i = std::make_shared( ma.range() ); (*mIndex)(indices...); return i; } template template OperationMaster OperationRoot::operator=(const OpClass& in) { if(mBlockII != nullptr){ return OperationMaster(mArrayRef, in, mIndex, mBlockII); } else { return OperationMaster(mArrayRef, in, mIndex); } } template const MBlock& OperationRoot::get() const { block(); return mBlock; } template MBlock& OperationRoot::get() { block(); return mBlock; } template OperationRoot& OperationRoot::set(const IndexInfo* blockIndex) { mBlockII = blockIndex; return *this; } template std::vector OperationRoot::block(const IndexInfo* blockIndex, bool init) const { std::vector btv(1, getBlockType( &mIInfo, blockIndex, true) ); if(init){ mBlock = makeBlock(mArrayRef.data(), btv[0].second, blockIndex->max()); } return btv; } template const OperationRoot& OperationRoot::block() const { mBlock.set( (*mIndex)().pos() ); return *this; } /******************* * Operation * *******************/ template Operation::Operation(const Ops&... ops) : //OperationTemplate >(this), mOps(ops...) {} template const BlockResult& Operation::get() const { PackNum::template unpackArgs(mRes, mOps); return mRes; } template std::vector Operation::block(const IndexInfo* blockIndex, bool init) const { std::vector btv; PackNum::makeBlockTypeVec(btv, mOps, blockIndex, init); if(init){ mRes.init(blockIndex->max()); } return btv; } template const Operation& Operation::block() const { //mBlock.set( mIndex->pos() ); return *this; } /********************* * Contraction * *********************/ template Contraction::Contraction(const Op& op, std::shared_ptr ind) : //OperationTemplate >(this), mOp(op), mInd(ind) {} template const BlockResult& Contraction::get() const { BlockBinaryOpSelf,decltype(mOp.get())> f(mRes); for(*mInd = 0; mInd->pos() != mInd->max(); ++(*mInd)){ f(mOp.get()); } return mRes; } template std::vector Contraction::block(const IndexInfo* blockIndex, bool init) const { if(init){ mRes.init(blockIndex->max()); } return mOp.block(blockIndex, init); } template const Contraction& Contraction::block() const { return *this; } } #endif