WIP: see last commit (there is still a non-initialized dynamic ext object)

src/include/dynamic_operation.cc.h

@@ -6,9 +6,6 @@ namespace MultiArrayTools
     template <typename T, class Operation>
     const T& DynamicOperation<T,Operation>::get(const DExtT& pos) const
     {
-	VCHECK(pos.size());
-	VCHECK(sizeof(ET)/sizeof(size_t));
-	assert(0);
 	return mOp.get(pos.expl<ET>());
     }
     

src/include/multi_array_operation.cc.h

@@ -172,7 +172,8 @@ namespace MultiArrayTools
     template <typename T, class Target, class OpClass, OpIndexAff OIA>
     inline void AssignmentExpr2<T,Target,OpClass,OIA>::operator()(size_t start)
     {
-        ExtType last;
+        ExtType last = rootSteps();
+	last.zero();
         mDataPtr[OpIndexResolve<OIA>::get(start,last)] = mSec.get(last.next());
         //mDataPtr[last.val()] = mSec.get(last.next());
     }
@@ -224,7 +225,8 @@ namespace MultiArrayTools
     template <typename T, class Target, class OpClass, OpIndexAff OIA>
     inline void AddExpr<T,Target,OpClass,OIA>::operator()(size_t start)
     {
-        ExtType last;
+        ExtType last = rootSteps();
+	last.zero();
         mDataPtr[OpIndexResolve<OIA>::get(start,last)] += mSec.get(last.next());
     }
     

src/include/xfor/exttype.h

@@ -75,6 +75,8 @@ namespace MultiArrayHelper
         template <class Y>
         inline MExt(const MExt<Y>& y);
 
+	inline void zero();
+	
 	inline const size_t& val() const;
 	inline const X& next() const;
 
@@ -108,6 +110,8 @@ namespace MultiArrayHelper
         inline None operator+(const None& in) const { return None(); }
 	inline None operator*(size_t in) const { return None(); }
 
+	inline void zero() {}
+	
         template <class Y>
         Y extend(const Y& y) const
         { return y; }
@@ -150,6 +154,8 @@ namespace MultiArrayHelper
 	inline const size_t& val() const;
 	inline None next() const { return None(); }
 
+	inline void zero();
+	
         template <size_t N>
         inline auto nn() const
             -> decltype(Getter<N>::getX(*this))
@@ -210,6 +216,13 @@ namespace MultiArrayHelper
 	return mNext;
     }
 
+    template <class X>
+    inline void MExt<X>::zero()
+    {
+	mExt = 0u;
+	mNext.zero();
+    }
+
     template <class X>
     inline MExt<X> MExt<X>::operator+(const MExt<X>& in) const
     {
@@ -248,6 +261,11 @@ namespace MultiArrayHelper
         mExt(y.val()) {}
 
     //template <>
+    inline void MExt<None>::zero()
+    {
+	mExt = 0u;
+    }
+
     inline const size_t& MExt<None>::val() const
     {
 	return mExt;

src/include/xfor/xfor.h

@@ -37,6 +37,9 @@ namespace MultiArrayHelper
 	//virtual size_t rootSteps() const = 0;
 	virtual std::shared_ptr<ExtBase> operator+(const ExtBase& in) const = 0;
 	virtual std::shared_ptr<ExtBase> operator*(size_t in) const = 0;
+	virtual void zero() = 0;
+
+	virtual std::shared_ptr<ExtBase> deepCopy() const = 0;
 	
 	template <class ExtType>
 	const ExtType& expl() const;
@@ -62,11 +65,14 @@ namespace MultiArrayHelper
 
 	ExtT(const ExtType& in) : mExt(in) {}
 
+	virtual std::shared_ptr<ExtBase> deepCopy() const override final { return std::make_shared<ExtT<ExtType>>(mExt); }
+
 	virtual size_t size() const override final { return sizeof(ExtType)/sizeof(size_t); }
 	//virtual size_t size() const override final { return ExtType::MExtSize(); }
 	//virtual size_t rootSteps() const override final;
 	const ExtType& ext() const { return mExt; }
 	virtual const size_t& val() const override final { return mExt.val(); }
+	virtual void zero() override final { mExt.zero(); }
 
 	virtual DExt operator+(const ExtBase& in) const override final
 	{ return std::make_shared<ExtT>( mExt + dynamic_cast<const ExtT&>(in).mExt ); }
@@ -79,17 +85,17 @@ namespace MultiArrayHelper
     class DExtTX
     {
     private:
-	DExt mDExt = nullptr;
+	mutable DExt mDExt = nullptr;
 	X mNext;
     public:
 	static constexpr size_t NUM = X::SIZE;
 	static constexpr size_t SIZE = NUM + 1;
 	
         DExtTX() = default;
-        DExtTX(const DExtTX& in) = default;
+        DExtTX(const DExtTX& in) : mDExt(in.mDExt->deepCopy()), mNext(in.mNext) {}
-        DExtTX(DExtTX&& in) = default;
+        DExtTX(DExtTX&& in) : mDExt(in.mDExt->deepCopy()), mNext(in.mNext) {}
-        DExtTX& operator=(const DExtTX& in) = default;
+        DExtTX& operator=(const DExtTX& in) { mNext = in.mNext; mDExt = in.mDExt->deepCopy(); }
-        DExtTX& operator=(DExtTX&& in) = default;
+        DExtTX& operator=(DExtTX&& in) { mNext = in.mNext; mDExt = in.mDExt->deepCopy(); }
         DExtTX(const DExt& in) : mDExt(in) {}
 
 	template <class Y>
@@ -110,7 +116,8 @@ namespace MultiArrayHelper
         { if (not mDExt) return *this; else return DExtTX( (*mDExt) * in ) ; }
 
 	template <class ExtType>
-	inline const ExtType& expl() const { return mDExt->expl<ExtType>(); }
+	inline const ExtType& expl() const
+	{ if(mDExt == nullptr) mDExt = std::make_shared<ExtT<ExtType>>(); assert(mDExt != nullptr); return mDExt->expl<ExtType>(); }
 
 	template <class Y>
 	inline auto extend(const Y& y) const -> DExtTX<decltype(mNext.extend(y))>
@@ -119,6 +126,8 @@ namespace MultiArrayHelper
 	inline const size_t& val() const { return mDExt->val(); }
 	inline const X& next() const { return mNext; }
 
+	inline void zero() { mDExt->zero(); }
+	
 	template <size_t N>
 	inline auto nn() const -> decltype(Getter<N>::getX(*this))
 	    { return Getter<N>::getX(*this); }
@@ -533,7 +542,8 @@ namespace MultiArrayHelper
     inline void For<IndexClass,Expr,FT>::operator()(size_t mlast)
     {
 	typedef typename IndexClass::RangeType RangeType;
-	const ExtType last;
+	ExtType last = rootSteps();
+	last.zero();
 	for(size_t pos = 0u; pos != ForBound<RangeType::ISSTATIC>::template bound<RangeType::SIZE>(mMax); ++pos){
 	    const size_t mnpos = PosForward<FT>::valuex(mlast, mStep, pos);
 	    const ExtType npos = last + mExt*pos;
@@ -635,10 +645,12 @@ namespace MultiArrayHelper
     {
         CHECK;
 	typedef typename IndexClass::RangeType RangeType;
-	const ExtType last;
+	ExtType last = rootSteps();
+	last.zero();
         int pos = 0;
         size_t mnpos = 0;
-        ExtType npos;
+        ExtType npos = rootSteps();
+	npos.zero();
 #pragma omp parallel shared(mExpr) private(pos,mnpos,npos)
         {
             auto expr = mExpr;
@@ -731,7 +743,8 @@ namespace MultiArrayHelper
     inline void SingleExpression<IndexClass,Expr>::operator()(size_t mlast)
     {
 	//typedef typename IndexClass::RangeType RangeType;
-	const ExtType last;
+	ExtType last = rootSteps();
+	last.zero();
 	const size_t pos = mIndPtr->pos();
 	const size_t mnpos = PosForward<ForType::DEFAULT>::value(mlast, mMax, pos);
 	const ExtType npos = last + mExt*pos;
@@ -816,7 +829,8 @@ namespace MultiArrayHelper
     template <class IndexClass, class Expr>
     inline void SubExpr<IndexClass,Expr>::operator()(size_t mlast)
     {
-	const ExtType last;
+	ExtType last = rootSteps();
+	last.zero();
         const size_t pos = (*mSubSet)[last.val()];
         const size_t mnpos = mlast;
         const ExtType npos = last + mExt*pos;