// -*- 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 "ranges/rheader.h" #include "pack_num.h" #include "arith.h" #include "xfor/xfor.h" namespace MultiArrayTools { namespace { using namespace MultiArrayHelper; } template class OperationBase { public: OperationClass& THIS() { return static_cast(*this); } const OperationClass& THIS() const { return static_cast(*this); } template auto operator+(const OperationBase& in) const; template auto operator-(const OperationBase& in) const; template auto operator*(const OperationBase& in) const; template auto operator/(const OperationBase& in) const; template auto c(const std::shared_ptr& ind) const -> Contraction; template auto sl(const std::shared_ptr&... inds) const -> ConstSlice; template auto slc(const std::shared_ptr&... inds) const -> SliceContraction; template auto p(const std::shared_ptr&... inds) const -> ConstOperationRoot; template auto to(const std::shared_ptr&... inds) const -> MultiArray; template auto addto(const std::shared_ptr&... inds) const -> MultiArray; template auto pto(const std::shared_ptr&... inds) const -> MultiArray; template auto paddto(const std::shared_ptr&... inds) const -> MultiArray; template // Args = Operation Classes auto a(const std::shared_ptr>& ll, const Args&... args) const -> Operation,OperationClass, Args...>; private: friend OperationClass; friend OperationTemplate; OperationBase() = default; }; template class OperationTemplate : public OperationBase { /* empty per default; specialize if needed */ private: OperationTemplate() = default; friend OperationClass; }; template size_t sumRootNum() { return typename Op::rootNum(); } template size_t sumRootNum() { return typename Op1::rootNum() + sumRootNum(); } template struct RootSumN { template struct rs { static constexpr size_t SIZE = Op1::SIZE + RootSumN::template rs::SIZE; }; template static inline auto rootSteps(const std::tuple& etp, std::intptr_t i) { return RootSumN::rootSteps(etp,i).extend( std::get(etp).rootSteps(i) ); } template static inline void exec( size_t start, ExtType last, std::tuple& etp) { std::get(etp)(start,last); RootSumN::exec(start,last.next(),etp); } template static inline size_t get( ExtType last, const std::tuple& etp) { std::get(etp).get(last); return RootSumN::get(last.next(),etp); } template static inline void set( ExtType last, std::tuple& etp) { std::get(etp).set(last); RootSumN::set(last.next(),etp); } }; template <> struct RootSumN<0> { template struct rs { static constexpr size_t SIZE = Op1::SIZE; }; template static inline auto rootSteps(const std::tuple& etp, std::intptr_t i) { return std::get<0>(etp).rootSteps(i); } template static inline void exec( size_t start, ExtType last, std::tuple& etp) { std::get(etp)(start,last); } template static inline size_t get( ExtType last, const std::tuple& etp) { std::get(etp).get(last); return 0; } template static inline void set( ExtType last, std::tuple& etp) { std::get(etp).set(last); } }; template struct RootSum { static constexpr size_t SIZE = RootSumN::template rs::SIZE; }; template struct SelfIdentity { static inline T& sapply(T& a, T b) { return a = b; } }; enum class OpIndexAff { EXTERN = 0, TARGET = 1 }; struct IAssign { template static inline void f(T*& t, size_t pos, const Op& op, ExtType e) { //a = b; t[pos] = op.get(e); } }; template struct IVAssign { template static inline void f(T1& a, const T2& b) { *reinterpret_cast(&a) = b; } }; struct IPlus { template static inline void f(T*& t, size_t pos, const Op& op, ExtType e) { t[pos] += op.get(e); } }; template struct IVPlus { template static inline void f(T1& a, const T2& b) { *reinterpret_cast(&a) += b; } }; template class AssignmentExpr : public ExpressionBase { private: AssignmentExpr() = default; Target mTar; OpClass mSec; T* mDataPtr; public: static constexpr size_t LAYER = 0; static constexpr size_t SIZE = Target::SIZE + OpClass::SIZE; typedef decltype(mTar.rootSteps(0).extend( mSec.rootSteps(0) )) ExtType; AssignmentExpr(T* dataPtr, const Target& tar, const OpClass& sec); AssignmentExpr(const AssignmentExpr& in) = default; AssignmentExpr(AssignmentExpr&& in) = default; AssignmentExpr& operator=(const AssignmentExpr& in) = default; AssignmentExpr& operator=(AssignmentExpr&& in) = default; virtual std::shared_ptr deepCopy() const override final { return std::make_shared>(*this); } inline void operator()(size_t start = 0); inline void operator()(size_t start, ExtType last); auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType; inline void operator()(size_t mlast, DExt last) override final; inline DExt dRootSteps(std::intptr_t iPtrNum = 0) const override final; inline DExt dExtension() const override final; }; template using AssignmentExpr2 = AssignmentExpr; template using AddExpr = AssignmentExpr; template class MOp { private: MOp() = default; std::tuple mOps; public: static constexpr size_t LAYER = 0; static constexpr size_t SIZE = RootSum::SIZE; typedef decltype(RootSumN::rootSteps(mOps,0) ) ExtType; MOp(const Ops&... exprs); MOp(const MOp& in) = default; MOp(MOp&& in) = default; MOp& operator=(const MOp& in) = default; MOp& operator=(MOp&& in) = default; inline size_t get(ExtType last) const; template inline size_t vget(ExtType last) const { return get(last); } inline MOp& set(ExtType last); auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType; template auto loop(Expr exp) const -> decltype(PackNum::mkLoop( mOps, exp)); T* data() const { assert(0); return nullptr; } }; template class GetExpr : public ExpressionBase { private: GetExpr() = default; OpClass mSec; NextExpr mNExpr; public: static constexpr size_t LAYER = 0; static constexpr size_t SIZE = OpClass::SIZE + NextExpr::SIZE; typedef decltype(mSec.rootSteps(0).extend( mNExpr.rootSteps(0) ) ) ExtType; GetExpr(const OpClass& sec, const NextExpr& nexpr); GetExpr(const GetExpr& in) = default; GetExpr(GetExpr&& in) = default; GetExpr& operator=(const GetExpr& in) = default; GetExpr& operator=(GetExpr&& in) = default; virtual std::shared_ptr deepCopy() const override final { return std::make_shared>(*this); } inline void operator()(size_t start = 0); inline void get(ExtType last); template inline void vget(ExtType last) { get(last); } inline void operator()(size_t start, ExtType last); auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType; inline void operator()(size_t mlast, DExt last) override final; inline DExt dRootSteps(std::intptr_t iPtrNum = 0) const override final; inline DExt dExtension() const override final; }; template auto mkGetExpr(const OpClass& op, const NextExpr& nexpr) { return GetExpr(op, nexpr); } template auto mkMOp(const Ops&... exprs) { return MOp(exprs...); } //template /* template class AddExpr : public ExpressionBase { private: AddExpr() = default; Target mTar; OpClass mSec; T* mDataPtr; public: static constexpr size_t LAYER = 0; static constexpr size_t SIZE = Target::SIZE + OpClass::SIZE; typedef decltype(mTar.rootSteps(0).extend( mSec.rootSteps(0) )) ExtType; // static constexpr size_t LAYER = 0; //static constexpr size_t SIZE = OpClass::SIZE; //typedef decltype(mSec.rootSteps()) ExtType; //AddExpr(T* dataPtr, const OpClass& sec); AddExpr(T* dataPtr, const Target& tar, const OpClass& sec); AddExpr(const AddExpr& in) = default; AddExpr(AddExpr&& in) = default; AddExpr& operator=(const AddExpr& in) = default; AddExpr& operator=(AddExpr&& in) = default; virtual std::shared_ptr deepCopy() const override final { return std::make_shared>(*this); } inline void operator()(size_t start = 0); inline void operator()(size_t start, ExtType last); auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType; inline void operator()(size_t mlast, DExt last) override final; inline DExt dRootSteps(std::intptr_t iPtrNum = 0) const override final; inline DExt dExtension() const override final; }; */ template class ConstOperationRoot : public OperationTemplate > { public: typedef T value_type; typedef OperationBase > OT; typedef ContainerRange CRange; typedef ContainerIndex IndexType; static constexpr size_t SIZE = 1; static constexpr bool CONT = true; ConstOperationRoot(const MultiArrayBase& ma, const std::shared_ptr&... indices); ConstOperationRoot(std::shared_ptr > maptr, const std::shared_ptr&... indices); ConstOperationRoot(const T* data, const IndexType& ind); template inline const T& get(ET pos) const; template inline const V& vget(ET pos) const; template inline ConstOperationRoot& set(ET pos); MExt rootSteps(std::intptr_t iPtrNum = 0) const; // nullptr for simple usage with decltype template Expr loop(Expr exp) const; const T* data() const; private: const T* mDataPtr; const T* mOrigDataPtr; IndexType mIndex; std::shared_ptr > mMaPtr; // never remove this ptr, otherwise we lose temporary container instances! }; template class StaticCast : public OperationTemplate > { private: Op mOp; public: typedef T value_type; typedef OperationBase > OT; typedef typename Op::CRange CRange; typedef typename Op::IndexType IndexType; static constexpr size_t SIZE = Op::SIZE; static constexpr bool CONT = false; StaticCast(const Op& op); template inline T get(ET pos) const; template inline V vget(ET pos) const; template inline StaticCast& set(ET pos); auto rootSteps(std::intptr_t iPtrNum = 0) const -> decltype(mOp.rootSteps(iPtrNum)); template Expr loop(Expr exp) const; }; template StaticCast staticcast(const Op& op) { return StaticCast(op); } template class MetaOperationRoot : public OperationTemplate > { public: typedef typename Range::IndexType IndexType; typedef typename IndexType::MetaType value_type; typedef OperationBase > OT; static constexpr size_t SIZE = 1; static constexpr bool CONT = false; MetaOperationRoot(const std::shared_ptr& ind); template inline value_type get(ET pos) const; template inline V vget(ET pos) const; template inline MetaOperationRoot& set(ET pos); MExt rootSteps(std::intptr_t iPtrNum = 0) const; // nullptr for simple usage with decltype template Expr loop(Expr exp) const; private: mutable IndexType mWorkIndex; std::shared_ptr mIndex; }; template class OperationRoot : public OperationTemplate > { public: typedef T value_type; typedef OperationBase > OT; typedef ContainerRange CRange; typedef ContainerIndex IndexType; static constexpr size_t SIZE = 1; static constexpr bool CONT = true; private: T* mDataPtr; T* mOrigDataPtr; IndexType mIndex; public: OperationRoot(MutableMultiArrayBase& ma, const std::shared_ptr&... indices); OperationRoot(MutableMultiArrayBase& ma, const std::tuple...>& indices); OperationRoot(T* data, const IndexType& ind); template auto assign(const OpClass& in) const -> decltype(mIndex.ifor(1,in.loop(AssignmentExpr2,OpClass,OpIndexAff::TARGET> (mOrigDataPtr,*this,in)))); template auto assignExpr(const OpClass& in) const -> decltype(in.loop(AssignmentExpr2,OpClass>(mOrigDataPtr,*this,in))); template auto assign(const OpClass& in, const std::shared_ptr& i) const -> decltype(i->ifor(1,in.loop(AssignmentExpr2,OpClass> (mOrigDataPtr,*this,in)))); template auto plus(const OpClass& in) const -> decltype(mIndex.ifor(1,in.loop(AddExpr,OpClass,OpIndexAff::TARGET> (mOrigDataPtr,*this,in)))); template auto plus(const OpClass& in, const std::shared_ptr& i) const -> decltype(i->ifor(1,in.loop(AddExpr,OpClass> (mOrigDataPtr,*this,in)))); template OperationRoot& operator=(const OpClass& in); template OperationRoot& operator+=(const OpClass& in); OperationRoot& operator=(const OperationRoot& in); ParallelOperationRoot par(); template inline T& get(ET pos) const; template inline V& vget(ET pos) const; template inline OperationRoot& set(ET pos); MExt rootSteps(std::intptr_t iPtrNum = 0) const; // nullptr for simple usage with decltype template Expr loop(Expr exp) const; T* data() const; template auto sl(const std::shared_ptr&... inds) -> Slice; }; template class ParallelOperationRoot : public OperationTemplate > { public: typedef T value_type; typedef OperationBase > OT; typedef ContainerRange CRange; typedef ContainerIndex IndexType; static constexpr size_t SIZE = 1; static constexpr bool CONT = true; private: T* mDataPtr; T* mOrigDataPtr; IndexType mIndex; public: ParallelOperationRoot(MutableMultiArrayBase& ma, const std::shared_ptr&... indices); ParallelOperationRoot(T* data, const IndexType& ind); template auto assign(const OpClass& in) -> decltype(mIndex.pifor(1,in.loop(AssignmentExpr2,OpClass,OpIndexAff::TARGET> (mOrigDataPtr,*this,in)))); template auto assign(const OpClass& in, const std::shared_ptr& i) const -> decltype(i->pifor(1,in.loop(AssignmentExpr2,OpClass> (mOrigDataPtr,*this,in)))); template auto plus(const OpClass& in) -> decltype(mIndex.pifor(1,in.loop(AddExpr,OpClass,OpIndexAff::TARGET> (mOrigDataPtr,*this,in)))); template auto plus(const OpClass& in, const std::shared_ptr& i) const -> decltype(i->pifor(1,in.loop(AddExpr,OpClass> (mOrigDataPtr,*this,in)))); template ParallelOperationRoot& operator=(const OpClass& in); template ParallelOperationRoot& operator+=(const OpClass& in); ParallelOperationRoot& operator=(const ParallelOperationRoot& in); template inline T& get(ET pos) const; template inline V& vget(ET pos) const; template inline ParallelOperationRoot& set(ET pos); MExt rootSteps(std::intptr_t iPtrNum = 0) const; // nullptr for simple usage with decltype template Expr loop(Expr exp) const; T* data() const; }; template class OperationValue : public OperationTemplate > { public: typedef T value_type; typedef OperationBase > OT; typedef ContainerRange CRange; typedef ContainerIndex IndexType; static constexpr size_t SIZE = 0; static constexpr bool CONT = true; OperationValue(const T& val); template inline const T& get(ET pos) const; template inline V vget(ET pos) const; template inline OperationValue& set(ET pos); None rootSteps(std::intptr_t iPtrNum = 0) const; // nullptr for simple usage with decltype template Expr loop(Expr exp) const; private: T mVal; }; template class Operation : public OperationTemplate > { public: typedef T value_type; typedef OperationBase > OT; typedef OpFunction F; static constexpr size_t SIZE = RootSum::SIZE; static constexpr bool FISSTATIC = OpFunction::FISSTATIC; static constexpr bool CONT = false; private: std::tuple mOps; std::shared_ptr mF; // only if non-static public: typedef decltype(PackNum::template mkSteps(0, mOps)) ETuple; Operation(const Ops&... ops); Operation(std::shared_ptr ff, const Ops&... ops); template inline auto get(ET pos) const; template inline auto vget(ET pos) const; template inline Operation& set(ET pos); auto rootSteps(std::intptr_t iPtrNum = 0) const // nullptr for simple usage with decltype -> decltype(PackNum::mkSteps(iPtrNum, mOps)); template auto loop(Expr exp) const -> decltype(PackNum::mkLoop( mOps, exp)); T* data() const { assert(0); return nullptr; } }; namespace { template struct OpMaker { template static inline auto mkOperation(const std::shared_ptr& f, const Ops&... ops) -> Operation { return Operation(f,ops...); } }; template <> struct OpMaker { template static inline auto mkOperation(const std::shared_ptr& f, const Ops&... ops) -> Operation { return Operation(ops...); } }; } template auto mkOperation(const std::shared_ptr& f, const Ops&... ops) -> Operation { return OpMaker::mkOperation(f, ops...); } template class Contraction : public OperationTemplate > { public: typedef T value_type; typedef OperationBase > OT; static constexpr size_t SIZE = Op::SIZE; static constexpr bool CONT = Op::CONT; private: Op mOp; std::shared_ptr mInd; public: typedef decltype(mOp.rootSteps(0)) ETuple; Contraction(const Op& op, std::shared_ptr ind); template inline auto get(ET pos) const -> decltype(mOp.template get(pos)); template inline auto vget(ET pos) const -> decltype(mOp.template vget(pos)); template inline Contraction& set(ET pos); T* data() const { assert(0); return nullptr; } auto rootSteps(std::intptr_t iPtrNum = 0) const // nullptr for simple usage with decltype -> decltype(mOp.rootSteps(iPtrNum)); template auto loop(Expr exp) const -> decltype(mInd->iforh(1,mOp.loop(exp))); }; template // class SliceContraction : public OperationTemplate //, //SliceContraction,Op,Indices...> > class SliceContraction : public OperationTemplate > { public: typedef MultiArray value_type; typedef OperationTemplate > OT; static constexpr size_t SIZE = Op::SIZE; static constexpr bool CONT = false; private: mutable Op mOp; mutable std::shared_ptr > mCont; mutable OperationRoot mTarOp; public: typedef decltype(mOp.rootSteps(0)) ETuple; SliceContraction(const Op& op, std::shared_ptr... ind); template inline const value_type& get(ET pos) const; template inline SliceContraction& set(ET pos); auto rootSteps(std::intptr_t iPtrNum = 0) const // nullptr for simple usage with decltype -> decltype(mOp.rootSteps(iPtrNum)); template auto loop(Expr exp) const -> decltype(mOp.loop(exp)); // no loop }; } #include "type_operations.h" #endif