#include "high_level_operation.h" namespace MultiArrayTools { template DynamicO mkDynOp1(const Op& op) { return DynamicO(op); } template template template void HighLevelOpBase::RetT::appendOuterM(const Op& op, const Ops&... ops) { // does not check anything regarding input !!! if(outer.init()){ outer = mkDynOp1(mkMOp(outer,op,ops...)); } else { outer = mkDynOp1(mkMOp(op,ops...)); } } template template void HighLevelOpBase::RetT::appendOuterM() {} template template void HighLevelOpBase::RetT::appendOuter(const DynamicO& in) { if(in.init()){ if(outer.init()){ outer = mkDynOp1(mkMOp(outer,in)); } else { outer = in; } } } template template void HighLevelOpBase::RetT::appendOuter(const RetT& in) { appendOuter(in.outer); } template HighLevelOpRoot::HighLevelOpRoot(const ROP& op) : mOp(op) {} template bool HighLevelOpRoot::root() const { return true; } template template auto HighLevelOpRoot::xcreate(const std::shared_ptr&... inds) -> typename B::template RetT { assert(0); return typename B::template RetT(); } template ROP* HighLevelOpRoot::get() { return &mOp; } namespace { template struct Create { template struct cx { template struct ccx { template static inline void cccx(typename HighLevelOpBase::template RetT& res, const std::array>,M>& in, const std::shared_ptr&... inds, const OPs&... ops, const DOPs&... dops) { static_assert(N > 0, "N > 0 failed"); auto& inn = std::get(in); if(not inn->root()){ auto dop = inn->create(inds...); auto op = *dop.op.data()->mOp; typedef decltype(op) OP; res.appendOuter(dop); assert(dop.op.init()); Create::template cx::template ccx::template cccx (res, in, inds..., op, ops..., dop, dops...); } else { auto& op = *inn->get(); typedef typename std::remove_reference::type OP; Create::template cx::template ccx::template cccx (res, in, inds..., op, ops..., dops...); } } }; }; }; template <> struct Create<0> { template struct cx { template struct ccx { template static inline void cccx(typename HighLevelOpBase::template RetT& res, const std::array>,M>& in, const std::shared_ptr&... inds, const OPs&... ops, const DOPs&... dops) { auto& inn = std::get<0>(in); if(not inn->root()){ auto dop = inn->create(inds...); auto op = *dop.op.data()->mOp; res.appendOuter(dop); res.op = mkDynOutOp(mkFOp(op,ops...), inds...); assert(dop.op.init()); res.appendOuterM(dop.op,dops.op...); } else { auto& op = *inn->get(); res.op = mkDynOutOp(mkFOp(op,ops...), inds...); res.appendOuterM(dops.op...); } } }; }; }; } template HighLevelOp::HighLevelOp(std::array>,N> in) : mIn(in) {} template bool HighLevelOp::root() const { return false; } template ROP* HighLevelOp::get() { assert(0); return nullptr; } template template auto HighLevelOp::xcreate(const std::shared_ptr&... inds) -> typename B::template RetT { typename B::template RetT res; Create::template cx::template ccx::template cccx (res,mIn,inds...); return res; } template HighLevelOpHolder::HighLevelOpHolder(const std::shared_ptr>& op) : mOp(op) {} template bool HighLevelOpHolder::root() const { return mOp->root(); } template template auto HighLevelOpHolder::create(const std::shared_ptr&... inds) const -> decltype(mOp->create(inds...)) { return mOp->create(inds...); } template auto HighLevelOpHolder::get() -> decltype(mOp->get()) { return mOp->get(); } template std::shared_ptr> HighLevelOpHolder::op() const { return mOp; } template HighLevelOpHolder HighLevelOpHolder::operator*(const HighLevelOpHolder& in) const { return HighLevelOpHolder ( std::make_shared,2>> ( std::array>,2>({mOp, in.mOp}) ) ); } template HighLevelOpHolder HighLevelOpHolder::operator+(const HighLevelOpHolder& in) const { return HighLevelOpHolder ( std::make_shared,2>> ( std::array>,2>({mOp, in.mOp}) ) ); } template HighLevelOpHolder HighLevelOpHolder::operator-(const HighLevelOpHolder& in) const { return HighLevelOpHolder ( std::make_shared,2>> ( std::array>,2>({mOp, in.mOp}) ) ); } template HighLevelOpHolder HighLevelOpHolder::operator/(const HighLevelOpHolder& in) const { return HighLevelOpHolder ( std::make_shared,2>> ( std::array>,2>({mOp, in.mOp}) ) ); } template HighLevelOpHolder mkSFunc(const HighLevelOpHolder& a, const HighLevelOpHolder&... as) { constexpr size_t N = sizeof...(ROPs)+1; return HighLevelOpHolder ( std::make_shared> ( std::array>,N>({a.op(), as.op()...}) ) ); } template template HighLevelOpHolder& HighLevelOpHolder::xassign(const HighLevelOpHolder& in, const std::shared_ptr& di, const std::shared_ptr&... is) { const size_t dim = di->dim(); if(dim > 2){ auto ci1 = di->getP(dim-2)->reduced(); auto ci2 = di->getP(dim-1)->reduced(); assert(ci1 != nullptr); assert(ci2 != nullptr); auto odi = mkSubSpaceX(di, dim-2); auto mi = mkMIndex(is..., odi); this->assign(in, mi, ci1, ci2); } else { assert(dim == 2 or dim == 1); auto ci1 = di->getP(dim-1)->reduced(); assert(ci1 != nullptr); auto odi = mkSubSpaceX(di, dim-1); auto mi = mkMIndex(is..., odi); this->assign(in, mi, ci1); } return *this; } template template HighLevelOpHolder& HighLevelOpHolder::xplus(const HighLevelOpHolder& in, const std::shared_ptr& di, const std::shared_ptr&... is) { const size_t dim = di->dim(); if(dim > 2){ auto ci1 = di->getP(dim-2)->reduced(); auto ci2 = di->getP(dim-1)->reduced(); assert(ci1 != nullptr); assert(ci2 != nullptr); auto odi = mkSubSpaceX(di, dim-2); auto mi = mkMIndex(is..., odi); this->plus(in, mi, ci1, ci2); } else { assert(dim == 2 or dim == 1); auto ci1 = di->getP(dim-1)->reduced(); assert(ci1 != nullptr); auto odi = mkSubSpaceX(di, dim-1); auto mi = mkMIndex(is..., odi); this->plus(in, mi, ci1); } return *this; } template std::string printInd(const std::shared_ptr& ind1, const std::shared_ptr& ind2, const std::shared_ptr&... inds) { return std::to_string(reinterpret_cast(ind1.get())) + "(" + std::to_string(ind1->max()) + "), " + printInd(ind2, inds...); } template std::string printInd(const std::shared_ptr& ind1) { return std::to_string(reinterpret_cast(ind1.get())) + "(" + std::to_string(ind1->max()) + ")"; } template template HighLevelOpHolder& HighLevelOpHolder::assign(const HighLevelOpHolder& in, const std::shared_ptr& mi, const std::shared_ptr&... inds) { auto xx = mkArrayPtr(nullr()); ROP& opr = *mOp->get(); if(in.root()){ auto inx = in; opr.par().assign( *inx.get(), mkMIndex(mi,inds...) )(); return *this; } auto loop = mkPILoop ( [&opr,&in,&xx,&inds...,this](){ auto inx = in; auto dop = inx.create(inds...); DynamicO gexp; if(dop.outer.init()){ gexp = mkDynOp1(mkMOp(dop.outer,dop.op)); } else { gexp = mkDynOp1(mkMOp(dop.op)); } auto xloop = mkILoop(std::make_tuple(*dop.op.data()->mOp), std::make_tuple(inds...), std::make_tuple(xx), std::make_tuple(opr.assign( *dop.op.data()->mOp, mkMIndex(inds...) )), std::array({1}), std::array({0})); return mkGetExpr(gexp, xloop); }); mi->pifor(1,loop)(); return *this; } template template HighLevelOpHolder& HighLevelOpHolder::plus(const HighLevelOpHolder& in, const std::shared_ptr& mi, const std::shared_ptr&... inds) { auto xx = mkArrayPtr(nullr()); ROP& opr = *mOp->get(); if(in.root()){ auto inx = in; opr.par().plus( *inx.get(), mkMIndex(mi,inds...) )(); return *this; } auto loop = mkPILoop ( [&opr,&in,&xx,&inds...,this](){ auto inx = in; auto dop = inx.create(inds...); DynamicO gexp; if(dop.outer.init()){ gexp = mkDynOp1(mkMOp(dop.outer,dop.op)); } else { gexp = mkDynOp1(mkMOp(dop.op)); } auto xloop = mkILoop(std::make_tuple(*dop.op.data()->mOp), std::make_tuple(inds...), std::make_tuple(xx), std::make_tuple(opr.plus( *dop.op.data()->mOp, mkMIndex(inds...) )), std::array({1}), std::array({0})); return mkGetExpr(gexp, xloop); }); mi->pifor(1,loop)(); return *this; } template HighLevelOpHolder mkHLO(const ROP& op) { return HighLevelOpHolder(std::make_shared>( op ) ); } #define SP " " #define regFunc1(fff) template \ HighLevelOpHolder hl_##fff (const HighLevelOpHolder& in) \ { return HighLevelOpHolder( std::make_shared,1>> \ ( std::array>,1>( {in.op()} ) ) ); } \ #include "extensions/math.h" #undef regFunc1 #undef SP /* template template inline void SetLInds::mkLIT(const ITuple& itp, const std::shared_ptr& di) { constexpr size_t NN = std::tuple_size::value-N-1; const size_t nn = di->dim()-N-1; typedef typename std::remove_reference(itp))>::type T; std::get(itp) = std::dynamic_pointer_cast(di->get(nn))->getIndex(); SetLInds::mkLIT(itp, di); } template template template inline void SetLInds::xx:: assign(Tar& tar, const Args&... args, const ITp& itp, const std::shared_ptr&... is) { SetLInds::template xx::assign(tar, args..., itp, std::get(itp), is...); } template template template inline void SetLInds::xx:: plus(Tar& tar, const Args&... args, const ITp& itp, const std::shared_ptr&... is) { SetLInds::template xx::plus(tar, args..., itp, std::get(itp), is...); } //template <> template inline void SetLInds<0>::mkLIT(const ITuple& itp, const std::shared_ptr& di) { constexpr size_t NN = std::tuple_size::value-1; const size_t nn = di->dim()-1; typedef typename std::remove_reference(itp))>::type T; std::get(itp) = std::dynamic_pointer_cast(di->get(nn))->getIndex(); } //template <> template template inline void SetLInds<0>::xx:: assign(Tar& tar, const Args&... args, const ITp& itp, const std::shared_ptr&... is) { tar.assign(args..., std::get<0>(itp), is...); } //template <> template template inline void SetLInds<0>::xx:: plus(Tar& tar, const Args&... args, const ITp& itp, const std::shared_ptr&... is) { tar.plus(args..., std::get<0>(itp), is...); } template size_t INDS::CallHLOpBase::depth() const { return mDepth; } template template void INDS::CallHLOp:: assign(HighLevelOpHolder& target, const HighLevelOpHolder& source, const std::shared_ptr&... is, const std::shared_ptr& di) const { auto ip = di->get(di->dim() - this->depth()); auto iregn = ip->regN(); if(iregn.type >= 0 and iregn.depth > sizeof...(LIndices)){ sNext[iregn.type]->assign(target, source, is..., di); } else { ITuple itp; SetLInds::mkLIT(itp,di); auto mi = mkIndex(is...,mkSubSpaceX(di, di->dim() - this->depth())); SetLInds:: template xx,ITuple,HighLevelOpHolder,decltype(mi)>:: assign(target, source, mi, itp); } } template template void INDS::CallHLOp:: plus(HighLevelOpHolder& target, const HighLevelOpHolder& source, const std::shared_ptr&... is, const std::shared_ptr& di) const { auto ip = di->get(di->dim() - this->depth()); auto iregn = ip->regN(); if(iregn.type >= 0 and iregn.depth > sizeof...(LIndices)){ sNext[iregn.type]->plus(target, source, is..., di); } else { ITuple itp; SetLInds::mkLIT(itp,di); auto mi = mkIndex(is...,mkSubSpaceX(di, di->dim() - this->depth())); SetLInds:: template xx,ITuple,HighLevelOpHolder,decltype(mi)>:: plus(target, source, mi, itp); } } */ }