Merge branch 'dev' into HEAD

src/include/arith.h

@@ -48,6 +48,7 @@ namespace MultiArrayTools
     {
 	static constexpr bool FISSTATIC = true;
 	typedef T value_type;
+	typedef F function;
 
 	template <class... Ops>
 	static auto mk(const Ops&... ops)
@@ -74,6 +75,11 @@ namespace MultiArrayTools
 	{
 	    return a;
 	}
+
+	static inline T selfApply(T& a1, const T& a2)
+	{
+	    return a1 = a2;
+	}
     };
 
     template <typename T, typename U>

src/include/dynamic_operation.h

@@ -17,6 +17,7 @@ namespace MultiArrayTools
 	
 	static constexpr size_t SIZE = 1;
 	static constexpr bool CONT = true;
+	static constexpr bool VABLE = false;
 	
 	DynamicOperationBase() = default;
 	DynamicOperationBase(const DynamicOperationBase& in) = default;
@@ -125,6 +126,7 @@ namespace MultiArrayTools
 	
 	static constexpr size_t SIZE = 1;
 	static constexpr bool CONT = true;
+	static constexpr bool VABLE = false;
 	
 	DynamicO() = default;
 	DynamicO(const DynamicO& in) : mOp(in.mOp ? in.mOp->deepCopy() : nullptr) {}
@@ -144,6 +146,11 @@ namespace MultiArrayTools
 	
         template <class X>
 	inline T get(const DExtTX<X>& pos) const { return mOp->get(pos.reduce()); }
+
+	template <typename V,class X>
+	inline auto vget(const DExtTX<X>& pos) const
+	{ return mOp->template vget<V>(pos.reduce()); }
+	
         template <class X>
 	inline DynamicO& set(const DExtTX<X>& pos) { mOp->set(pos.reduce()); return *this; }
 	inline DExtT rootSteps(std::intptr_t iPtrNum = 0) const { return mOp->rootSteps(iPtrNum); }

src/include/map_range.cc.h

@@ -27,6 +27,18 @@ namespace MultiArrayTools
         assert(mIndPtr != nullptr);
     }
 
+    template <class Op, class Index, class Expr, SpaceType STYPE>
+    std::shared_ptr<ExpressionBase> OpExpr<Op,Index,Expr,STYPE>::deepCopy() const
+    {
+	return std::make_shared<OpExpr<Op,Index,Expr,STYPE>>(*this);
+    }
+    
+    template <class Op, class Index, class Expr, SpaceType STYPE>
+    inline void OpExpr<Op,Index,Expr,STYPE>::operator()(size_t mlast, DExt last)
+    {
+	operator()(mlast, std::dynamic_pointer_cast<ExtT<ExtType>>(last)->ext());
+    }
+
     template <class Op, class Index, class Expr, SpaceType STYPE>
     inline void OpExpr<Op,Index,Expr,STYPE>::operator()(size_t mlast,
                                                             ExtType last)
@@ -63,6 +75,17 @@ namespace MultiArrayTools
 	//return mExpr.rootSteps(iPtrNum).extend( mOp.rootSteps(iPtrNum) );
     }
 
+    template <class Op, class Index, class Expr, SpaceType STYPE>
+    DExt OpExpr<Op,Index,Expr,STYPE>::dRootSteps(std::intptr_t iPtrNum) const
+    {
+	return std::make_shared<ExtT<ExtType>>(rootSteps(iPtrNum));
+    }
+    
+    template <class Op, class Index, class Expr, SpaceType STYPE>
+    DExt OpExpr<Op,Index,Expr,STYPE>::dExtension() const
+    {
+	return std::make_shared<ExtT<ExtType>>(mExt);
+    }
 
     // -> define in range_base.cc
     //std::shared_ptr<RangeFactoryBase> mkMULTI(const char** dp);

src/include/map_range.h

@@ -50,13 +50,14 @@ namespace MultiArrayTools
     
     template <class Op, class Index, class Expr, SpaceType STYPE = SpaceType::ANY>
     //template <class MapF, class IndexPack, class Expr, SpaceType STYPE = SpaceType::ANY>
-    class OpExpr
+    class OpExpr : public ExpressionBase
     {
     public:
         //typedef typename Index::OIType OIType;
 	//typedef SingleIndex<typename Op::value_type,STYPE> OIType;
 	static constexpr size_t LAYER = Expr::LAYER + 1;
 	static constexpr size_t SIZE = Expr::SIZE + Op::SIZE;
+	static constexpr size_t NHLAYER = Expr::NHLAYER + 1;
 	
     private:
 	OpExpr() = default;
@@ -80,11 +81,20 @@ namespace MultiArrayTools
 
         OpExpr(const Op& mapf, const Index* ind, size_t step, Expr ex);
 
+	virtual std::shared_ptr<ExpressionBase> deepCopy() const override final;
+
+	template <size_t VS>
+	inline auto vec() const { return *this; }
+	
+	inline void operator()(size_t mlast, DExt last) override final;
 	inline void operator()(size_t mlast, ExtType last);
-	inline void operator()(size_t mlast = 0);
+	inline void operator()(size_t mlast = 0) override final;
 
 	auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
 
+	virtual DExt dRootSteps(std::intptr_t iPtrNum = 0) const override final;
+	virtual DExt dExtension() const override final;
+	
     };
     
     template <class OIType, class Op, SpaceType XSTYPE, class... Indices>

src/include/multi_array_operation.cc.h

@@ -173,47 +173,45 @@ namespace MultiArrayTools
     };
 
     
-    template <typename T, class Target, class OpClass, OpIndexAff OIA>
+    template <typename T, class IOp, class Target, class OpClass, OpIndexAff OIA>
-    AssignmentExpr2<T,Target,OpClass,OIA>::AssignmentExpr2(T* dataPtr, const Target& tar, const OpClass& sec) :
+    AssignmentExpr<T,IOp,Target,OpClass,OIA>::AssignmentExpr(T* dataPtr, const Target& tar, const OpClass& sec) :
         mTar(tar), mSec(sec), mDataPtr(dataPtr) {}
     
-    template <typename T, class Target, class OpClass, OpIndexAff OIA>
+    template <typename T, class IOp, class Target, class OpClass, OpIndexAff OIA>
-    inline void AssignmentExpr2<T,Target,OpClass,OIA>::operator()(size_t start)
+    inline void AssignmentExpr<T,IOp,Target,OpClass,OIA>::operator()(size_t start)
     {
         ExtType last = rootSteps();
 	last.zero();
-        mDataPtr[OpIndexResolve<OIA>::get(start,last)] = mSec.get(last.next());
+	IOp::f(mDataPtr,OpIndexResolve<OIA>::get(start,last),mSec,last.next());
-        //mDataPtr[last.val()] = mSec.get(last.next());
     }
 
-    template <typename T, class Target, class OpClass, OpIndexAff OIA>
+    template <typename T, class IOp, class Target, class OpClass, OpIndexAff OIA>
-    inline void AssignmentExpr2<T,Target,OpClass,OIA>::operator()(size_t start, ExtType last)
+    inline void AssignmentExpr<T,IOp,Target,OpClass,OIA>::operator()(size_t start, ExtType last)
     {
-        //CHECK;
+	IOp::f(mDataPtr,OpIndexResolve<OIA>::get(start,last),mSec,last.next());
-        mDataPtr[OpIndexResolve<OIA>::get(start,last)] = mSec.get(last.next());
     }
     
-    template <typename T, class Target, class OpClass, OpIndexAff OIA>
+    template <typename T, class IOp, class Target, class OpClass, OpIndexAff OIA>
-    typename AssignmentExpr2<T,Target,OpClass,OIA>::ExtType
+    typename AssignmentExpr<T,IOp,Target,OpClass,OIA>::ExtType
-    AssignmentExpr2<T,Target,OpClass,OIA>::rootSteps(std::intptr_t iPtrNum) const
+    AssignmentExpr<T,IOp,Target,OpClass,OIA>::rootSteps(std::intptr_t iPtrNum) const
     {
 	return mTar.rootSteps(iPtrNum).extend( mSec.rootSteps(iPtrNum) );
     }
 
-    template <typename T, class Target, class OpClass, OpIndexAff OIA>
+    template <typename T, class IOp, class Target, class OpClass, OpIndexAff OIA>
-    inline void AssignmentExpr2<T,Target,OpClass,OIA>::operator()(size_t mlast, DExt last)
+    inline void AssignmentExpr<T,IOp,Target,OpClass,OIA>::operator()(size_t mlast, DExt last)
     {
         (*this)(mlast, std::dynamic_pointer_cast<ExtT<ExtType>>(last)->ext());
     }
 
-    template <typename T, class Target, class OpClass, OpIndexAff OIA>
+    template <typename T, class IOp, class Target, class OpClass, OpIndexAff OIA>
-    inline DExt AssignmentExpr2<T,Target,OpClass,OIA>::dRootSteps(std::intptr_t iPtrNum) const
+    inline DExt AssignmentExpr<T,IOp,Target,OpClass,OIA>::dRootSteps(std::intptr_t iPtrNum) const
     {
         return std::make_shared<ExtT<ExtType>>(rootSteps(iPtrNum));
     }
 
-    template <typename T, class Target, class OpClass, OpIndexAff OIA>
+    template <typename T, class IOp, class Target, class OpClass, OpIndexAff OIA>
-    inline DExt AssignmentExpr2<T,Target,OpClass,OIA>::dExtension() const
+    inline DExt AssignmentExpr<T,IOp,Target,OpClass,OIA>::dExtension() const
     {
         CHECK;
         return nullptr; //???!!!
@@ -304,52 +302,6 @@ namespace MultiArrayTools
         return nullptr; //???!!!
     }
 
-    template <typename T, class Target, class OpClass, OpIndexAff OIA>
-    AddExpr<T,Target,OpClass,OIA>::AddExpr(T* dataPtr, const Target& tar, const OpClass& sec) :
-        mTar(tar), mSec(sec), mDataPtr(dataPtr) {}
-    
-    template <typename T, class Target, class OpClass, OpIndexAff OIA>
-    inline void AddExpr<T,Target,OpClass,OIA>::operator()(size_t start, ExtType last)
-    {
-        mDataPtr[OpIndexResolve<OIA>::get(start,last)] += mSec.get(last.next());
-        //mDataPtr[start] += mSec.get(last.next());
-        //mDataPtr[start] += mSec.template get<ExtType>(last);
-    }
-
-    template <typename T, class Target, class OpClass, OpIndexAff OIA>
-    inline void AddExpr<T,Target,OpClass,OIA>::operator()(size_t start)
-    {
-        ExtType last = rootSteps();
-	last.zero();
-        mDataPtr[OpIndexResolve<OIA>::get(start,last)] += mSec.get(last.next());
-    }
-    
-    template <typename T, class Target, class OpClass, OpIndexAff OIA>
-    typename AddExpr<T,Target,OpClass,OIA>::ExtType AddExpr<T,Target,OpClass,OIA>::rootSteps(std::intptr_t iPtrNum) const
-    {
-	return mTar.rootSteps(iPtrNum).extend( mSec.rootSteps(iPtrNum) );
-	//return mSec.rootSteps(iPtrNum);
-    }
-
-    template <typename T, class Target, class OpClass, OpIndexAff OIA>
-    inline void AddExpr<T,Target,OpClass,OIA>::operator()(size_t mlast, DExt last)
-    {
-        (*this)(mlast, std::dynamic_pointer_cast<ExtT<ExtType>>(last)->ext());
-    }
-
-    template <typename T, class Target, class OpClass, OpIndexAff OIA>
-    inline DExt AddExpr<T,Target,OpClass,OIA>::dRootSteps(std::intptr_t iPtrNum) const
-    {
-        return std::make_shared<ExtT<ExtType>>(rootSteps(iPtrNum));
-    }
-
-    template <typename T, class Target, class OpClass, OpIndexAff OIA>
-    inline DExt AddExpr<T,Target,OpClass,OIA>::dExtension() const
-    {
-        CHECK;
-        return nullptr; //???!!!
-    }
-
     /****************************
      *   ConstOperationRoot     *
      ****************************/
@@ -398,6 +350,14 @@ namespace MultiArrayTools
 	return mDataPtr[pos.val()/*+mOff*/];
     }
 
+    template <typename T, class... Ranges>
+    template <typename V, class ET>
+    inline const V& ConstOperationRoot<T,Ranges...>::vget(ET pos) const
+    {
+	//VCHECK(pos.val());
+	return *(reinterpret_cast<const V*>(mDataPtr+pos.val()));
+    }
+
     template <typename T, class... Ranges>
     template <class ET>
     inline ConstOperationRoot<T,Ranges...>& ConstOperationRoot<T,Ranges...>::set(ET pos)
@@ -445,6 +405,14 @@ namespace MultiArrayTools
 	return static_cast<T>( mOp.get(pos) );
     }
 
+    template <typename T, class Op>
+    template <typename V, class ET>
+    inline V StaticCast<T,Op>::vget(ET pos) const
+    {
+	assert(0); // !!!
+	return V();
+    }
+
     template <typename T, class Op>
     template <class ET>
     inline StaticCast<T,Op>& StaticCast<T,Op>::set(ET pos)
@@ -489,6 +457,15 @@ namespace MultiArrayTools
 	return (mWorkIndex = pos.val()).meta();
     }
 
+    template <class Range>
+    template <typename V, class ET>
+    inline V
+    MetaOperationRoot<Range>::vget(ET pos) const
+    {
+	assert(0); // !!!
+	return V();
+    }
+
     template <class Range>
     template <class ET>
     inline MetaOperationRoot<Range>& MetaOperationRoot<Range>::set(ET pos)
@@ -551,66 +528,116 @@ namespace MultiArrayTools
     }
 
     template <typename T, class... Ranges>
-    template <class OpClass>
+    template <class IOp, class OpClass>
-    auto OperationRoot<T,Ranges...>::assign(const OpClass& in) const
+    auto OperationRoot<T,Ranges...>::asx(const OpClass& in) const
-        -> decltype(mIndex.ifor(1,in.loop(AssignmentExpr2<T,OperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
+        -> decltype(mIndex.ifor(1,in.loop(AssignmentExpr<T,IOp,OperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
-                                          (mOrigDataPtr,*this,in))))
+                                          (mOrigDataPtr,*this,in))).template vec<IOp::VSIZE>())
 
     {
         static_assert( OpClass::SIZE == decltype(in.rootSteps())::SIZE, "Ext Size mismatch" );
-        return mIndex.ifor(1,in.loop(AssignmentExpr2<T,OperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
+        return mIndex.ifor(1,in.loop(AssignmentExpr<T,IOp,OperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
-                                     (mOrigDataPtr,*this,in)));
+                                     (mOrigDataPtr,*this,in))).template vec<IOp::VSIZE>();
+    }
+
+    template <typename T, class... Ranges>
+    template <class IOp, class OpClass>
+    auto OperationRoot<T,Ranges...>::asxExpr(const OpClass& in) const
+        -> decltype(in.loop(AssignmentExpr<T,IOp,OperationRoot<T,Ranges...>,OpClass>
+                            (mOrigDataPtr,*this,in)))
+    {
+        static_assert( OpClass::SIZE == decltype(in.rootSteps())::SIZE, "Ext Size mismatch" );
+        return in.loop(AssignmentExpr<T,IOp,OperationRoot<T,Ranges...>,OpClass>
+                       (mOrigDataPtr,*this,in));
+    }
+    
+    template <typename T, class... Ranges>
+    template <class IOp, class OpClass, class Index>
+    auto OperationRoot<T,Ranges...>::asx(const OpClass& in, const std::shared_ptr<Index>& i) const
+        -> decltype(i->ifor(1,in.loop(AssignmentExpr<T,IOp,OperationRoot<T,Ranges...>,OpClass>
+                                      (mOrigDataPtr,*this,in))).template vec<IOp::VSIZE>())
+    {
+        static_assert( OpClass::SIZE == decltype(in.rootSteps())::SIZE, "Ext Size mismatch" );
+        return i->ifor(1,in.loop(AssignmentExpr<T,IOp,OperationRoot<T,Ranges...>,OpClass>
+                                 (mOrigDataPtr,*this,in))).template vec<IOp::VSIZE>();
+    }
+
+    template <typename T, class... Ranges>
+    template <class OpClass>
+    auto OperationRoot<T,Ranges...>::assign(const OpClass& in) const
+        -> decltype(this->template asx<IAssign<T>>(in))
+    {
+        return this->template asx<IAssign<T>>(in);
     }
 
     template <typename T, class... Ranges>
     template <class OpClass>
     auto OperationRoot<T,Ranges...>::assignExpr(const OpClass& in) const
-        -> decltype(in.loop(AssignmentExpr2<T,OperationRoot<T,Ranges...>,OpClass>
+        -> decltype(this->template asxExpr<IAssign<T>>(in))
-                            (mOrigDataPtr,*this,in)))
     {
-        static_assert( OpClass::SIZE == decltype(in.rootSteps())::SIZE, "Ext Size mismatch" );
+        return this->template asxExpr<IAssign<T>>(in);
-        return in.loop(AssignmentExpr2<T,OperationRoot<T,Ranges...>,OpClass>
-                       (mOrigDataPtr,*this,in));
     }
     
     template <typename T, class... Ranges>
     template <class OpClass, class Index>
     auto OperationRoot<T,Ranges...>::assign(const OpClass& in, const std::shared_ptr<Index>& i) const
-        -> decltype(i->ifor(1,in.loop(AssignmentExpr2<T,OperationRoot<T,Ranges...>,OpClass>
+        -> decltype(this->template asx<IAssign<T>>(in,i))
-                                      (mOrigDataPtr,*this,in))))
     {
-        static_assert( OpClass::SIZE == decltype(in.rootSteps())::SIZE, "Ext Size mismatch" );
+        return this->template asx<IAssign<T>>(in,i);
-        return i->ifor(1,in.loop(AssignmentExpr2<T,OperationRoot<T,Ranges...>,OpClass>
-                                 (mOrigDataPtr,*this,in)));
     }
 
     template <typename T, class... Ranges>
     template <class OpClass>
     auto OperationRoot<T,Ranges...>::plus(const OpClass& in) const
-        -> decltype(mIndex.ifor(1,in.loop(AddExpr<T,OperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
+        -> decltype(this->template asx<IPlus<T>>(in))
-                                          (mOrigDataPtr,*this,in))))
     {
-        static_assert( OpClass::SIZE == decltype(in.rootSteps())::SIZE, "Ext Size mismatch" );
+        return this->template asx<IPlus<T>>(in);
-        return mIndex.ifor(1,in.loop(AddExpr<T,OperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
-                                     (mOrigDataPtr,*this,in)));
     }
 
     template <typename T, class... Ranges>
     template <class OpClass, class Index>
     auto OperationRoot<T,Ranges...>::plus(const OpClass& in, const std::shared_ptr<Index>& i) const
-        -> decltype(i->ifor(1,in.loop(AddExpr<T,OperationRoot<T,Ranges...>,OpClass>
+        -> decltype(this->template asx<IPlus<T>>(in,i))
-                                      (mOrigDataPtr,*this,in))))
     {
-        static_assert( OpClass::SIZE == decltype(in.rootSteps())::SIZE, "Ext Size mismatch" );
+        return this->template asx<IPlus<T>>(in,i);
-        return i->ifor(1,in.loop(AddExpr<T,OperationRoot<T,Ranges...>,OpClass>
-                                 (mOrigDataPtr,*this,in)));
     }
 
+    template <bool VABLE = false>
+    struct VExec
+    {
+	template <template <typename> class F, typename TarOp, class OpClass>
+	static inline void exec(TarOp& th, const OpClass& in)
+	{
+	    typedef typename TarOp::value_type T;
+	    IAccess<T,F<T>> tmp;
+	    th.template asx<decltype(tmp)>(in)();
+	}
+    };
+
+    template <>
+    struct VExec<true>
+    {
+	template <template <typename> class F, typename TarOp, class OpClass>
+	static inline void exec(TarOp& th, const OpClass& in)
+	{
+	    CHECK;
+	    typedef typename TarOp::value_type T;
+	    auto x = th.template asx<IVAccess<T,F<T>>>(in);
+	    if(x.rootSteps(x.vI()) == 1){
+	    //if(0){
+		CHECK;
+		x();
+	    }
+	    else {
+		th.template asx<IAccess<T,F<T>>>(in)();
+	    }
+	}
+    };
+
     template <typename T, class... Ranges>
     template <class OpClass>
     OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::operator=(const OpClass& in)
     {
-        assign(in)();
+	VExec<OpClass::VABLE>::template exec<identity>(*this,in);
         return *this;
     }
 
@@ -618,7 +645,8 @@ namespace MultiArrayTools
     template <class OpClass>
     OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::operator+=(const OpClass& in)
     {
-        plus(in)();
+	VExec<OpClass::VABLE>::template exec<xxxplus>(*this,in);
+        //plus(in)();
         return *this;
     }
 
@@ -641,6 +669,14 @@ namespace MultiArrayTools
 	return mDataPtr[pos.val()];
     }
 
+    template <typename T, class... Ranges>
+    template <typename V, class ET>
+    inline V& OperationRoot<T,Ranges...>::vget(ET pos) const
+    {
+	//VCHECK(pos.val());
+	return *(reinterpret_cast<V*>(mDataPtr+pos.val()));
+    }
+
     template <typename T, class... Ranges>
     template <class ET>
     inline OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::set(ET pos)
@@ -709,53 +745,84 @@ namespace MultiArrayTools
     }
 
     template <typename T, class... Ranges>
-    template <class OpClass>
+    template <class IOp, class OpClass>
-    auto ParallelOperationRoot<T,Ranges...>::assign(const OpClass& in)
+    auto ParallelOperationRoot<T,Ranges...>::asx(const OpClass& in) const
-        -> decltype(mIndex.pifor(1,in.loop(AssignmentExpr2<T,ParallelOperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
+        -> decltype(mIndex.pifor(1,in.loop(AssignmentExpr<T,IOp,ParallelOperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
-                                          (mOrigDataPtr,*this,in))))
+                                          (mOrigDataPtr,*this,in))).template vec<IOp::VSIZE>())
 
     {
-        return mIndex.pifor(1,in.loop(AssignmentExpr2<T,ParallelOperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
+        static_assert( OpClass::SIZE == decltype(in.rootSteps())::SIZE, "Ext Size mismatch" );
-                                      (mOrigDataPtr,*this,in)));
+        return mIndex.pifor(1,in.loop(AssignmentExpr<T,IOp,ParallelOperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
+                                     (mOrigDataPtr,*this,in))).template vec<IOp::VSIZE>();
+    }
+
+    template <typename T, class... Ranges>
+    template <class IOp, class OpClass>
+    auto ParallelOperationRoot<T,Ranges...>::asxExpr(const OpClass& in) const
+        -> decltype(in.loop(AssignmentExpr<T,IOp,ParallelOperationRoot<T,Ranges...>,OpClass>
+                            (mOrigDataPtr,*this,in)))
+    {
+        static_assert( OpClass::SIZE == decltype(in.rootSteps())::SIZE, "Ext Size mismatch" );
+        return in.loop(AssignmentExpr<T,IOp,ParallelOperationRoot<T,Ranges...>,OpClass>
+                       (mOrigDataPtr,*this,in));
+    }
+    
+    template <typename T, class... Ranges>
+    template <class IOp, class OpClass, class Index>
+    auto ParallelOperationRoot<T,Ranges...>::asx(const OpClass& in, const std::shared_ptr<Index>& i) const
+        -> decltype(i->pifor(1,in.loop(AssignmentExpr<T,IOp,ParallelOperationRoot<T,Ranges...>,OpClass>
+                                      (mOrigDataPtr,*this,in))).template vec<IOp::VSIZE>())
+    {
+        static_assert( OpClass::SIZE == decltype(in.rootSteps())::SIZE, "Ext Size mismatch" );
+        return i->pifor(1,in.loop(AssignmentExpr<T,IOp,ParallelOperationRoot<T,Ranges...>,OpClass>
+                                 (mOrigDataPtr,*this,in))).template vec<IOp::VSIZE>();
+    }
+
+    template <typename T, class... Ranges>
+    template <class OpClass>
+    auto ParallelOperationRoot<T,Ranges...>::assign(const OpClass& in) const
+        -> decltype(this->template asx<IAssign<T>>(in))
+    {
+        return this->template asx<IAssign<T>>(in);
+    }
+
+    template <typename T, class... Ranges>
+    template <class OpClass>
+    auto ParallelOperationRoot<T,Ranges...>::assignExpr(const OpClass& in) const
+        -> decltype(this->template asxExpr<IAssign<T>>(in))
+    {
+        return this->template asxExpr<IAssign<T>>(in);
     }
     
     template <typename T, class... Ranges>
     template <class OpClass, class Index>
     auto ParallelOperationRoot<T,Ranges...>::assign(const OpClass& in, const std::shared_ptr<Index>& i) const
-        -> decltype(i->pifor(1,in.loop(AssignmentExpr2<T,ParallelOperationRoot<T,Ranges...>,OpClass>
+        -> decltype(this->template asx<IAssign<T>>(in,i))
-				       (mOrigDataPtr,*this,in))))
     {
-        static_assert( OpClass::SIZE == decltype(in.rootSteps())::SIZE, "Ext Size mismatch" );
+        return this->template asx<IAssign<T>>(in,i);
-        return i->pifor(1,in.loop(AssignmentExpr2<T,ParallelOperationRoot<T,Ranges...>,OpClass>
-				  (mOrigDataPtr,*this,in)));
     }
 
     template <typename T, class... Ranges>
     template <class OpClass>
-    auto ParallelOperationRoot<T,Ranges...>::plus(const OpClass& in)
+    auto ParallelOperationRoot<T,Ranges...>::plus(const OpClass& in) const
-        -> decltype(mIndex.pifor(1,in.loop(AddExpr<T,ParallelOperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
+        -> decltype(this->template asx<IPlus<T>>(in))
-                                           (mOrigDataPtr,*this,in))))
     {
-        return mIndex.pifor(1,in.loop(AddExpr<T,ParallelOperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
+        return this->template asx<IPlus<T>>(in);
-                                      (mOrigDataPtr,*this,in)));
     }
 
     template <typename T, class... Ranges>
     template <class OpClass, class Index>
     auto ParallelOperationRoot<T,Ranges...>::plus(const OpClass& in, const std::shared_ptr<Index>& i) const
-        -> decltype(i->pifor(1,in.loop(AddExpr<T,ParallelOperationRoot<T,Ranges...>,OpClass>
+        -> decltype(this->template asx<IPlus<T>>(in,i))
-				       (mOrigDataPtr,*this,in))))
     {
-        static_assert( OpClass::SIZE == decltype(in.rootSteps())::SIZE, "Ext Size mismatch" );
+        return this->template asx<IPlus<T>>(in,i);
-        return i->pifor(1,in.loop(AddExpr<T,ParallelOperationRoot<T,Ranges...>,OpClass>
-				  (mOrigDataPtr,*this,in)));
     }
 
     template <typename T, class... Ranges>
     template <class OpClass>
     ParallelOperationRoot<T,Ranges...>& ParallelOperationRoot<T,Ranges...>::operator=(const OpClass& in)
     {
-        assign(in)();
+	VExec<OpClass::VABLE>::template exec<identity>(*this,in);
         return *this;
     }
 
@@ -763,7 +830,7 @@ namespace MultiArrayTools
     template <class OpClass>
     ParallelOperationRoot<T,Ranges...>& ParallelOperationRoot<T,Ranges...>::operator+=(const OpClass& in)
     {
-        plus(in)();
+	VExec<OpClass::VABLE>::template exec<xxxplus>(*this,in);
         return *this;
     }
 
@@ -781,6 +848,13 @@ namespace MultiArrayTools
 	return mDataPtr[pos.val()/*+mOff*/];
     }
 
+    template <typename T, class... Ranges>
+    template <typename V, class ET>
+    inline V& ParallelOperationRoot<T,Ranges...>::vget(ET pos) const
+    {
+	return *(reinterpret_cast<V*>(mDataPtr+pos.val()));
+    }
+
     template <typename T, class... Ranges>
     template <class ET>
     inline ParallelOperationRoot<T,Ranges...>& ParallelOperationRoot<T,Ranges...>::set(ET pos)
@@ -826,6 +900,13 @@ namespace MultiArrayTools
 	return mVal;
     }
 
+    template <typename T>
+    template <typename V, class ET>
+    inline V OperationValue<T>::vget(ET pos) const
+    {
+	return static_cast<V>(mVal); // implement???!!!
+    }
+
     template <typename T>
     template <class ET>
     inline OperationValue<T>& OperationValue<T>::set(ET pos)
@@ -910,6 +991,15 @@ namespace MultiArrayTools
 	    template mkOpExpr<SIZE,ET,OpTuple,OpFunction>(mF, pos, mOps);
     }
 
+    template <typename T, class OpFunction, class... Ops>
+    template <typename V, class ET>
+    inline auto Operation<T,OpFunction,Ops...>::vget(ET pos) const
+    {
+	typedef std::tuple<Ops...> OpTuple;
+	return PackNum<sizeof...(Ops)-1>::
+	    template mkVOpExpr<SIZE,V,ET,OpTuple,VFunc<OpFunction>>(mkVFuncPtr(mF), pos, mOps); // implement!!!
+    }
+
     template <typename T, class OpFunction, class... Ops>
     template <class ET>
     inline Operation<T,OpFunction,Ops...>& Operation<T,OpFunction,Ops...>::set(ET pos)
@@ -953,6 +1043,14 @@ namespace MultiArrayTools
 	return mOp.template get<ET>(pos);
     }
 
+    template <typename T, class Op, class IndexType>
+    template <typename V, class ET>
+    inline auto Contraction<T,Op,IndexType>::vget(ET pos) const
+        -> decltype(mOp.template vget<V,ET>(pos))
+    {
+	return mOp.template vget<V,ET>(pos);
+    }
+
     template <typename T, class Op, class IndexType>
     template <class ET>
     inline Contraction<T,Op,IndexType>& Contraction<T,Op,IndexType>::set(ET pos)

src/include/multi_array_operation.h

@@ -8,6 +8,7 @@
 #include <cmath>
 #include <map>
 #include <utility>
+#include <type_traits>
 
 #include "base_def.h"
 #include "mbase_def.h"
@@ -16,8 +17,8 @@
 #include "pack_num.h"
 
 #include "arith.h"
-
 #include "xfor/xfor.h"
+#include "type_operations.h"
 
 namespace MultiArrayTools
 {
@@ -207,11 +208,126 @@ namespace MultiArrayTools
                            TARGET = 1
     };
 
-    template <typename T, class Target, class OpClass, OpIndexAff OIA=OpIndexAff::EXTERN>
+    template <class T>
-    class AssignmentExpr2 : public ExpressionBase
+    struct VType
+    {
+	typedef T type;
+	static constexpr size_t MULT = sizeof(type)/sizeof(T);
+    };
+
+    template <>
+    struct VType<double>
+    {
+	typedef v256 type;
+	static constexpr size_t MULT = sizeof(type)/sizeof(double);
+    };
+
+    template <template <typename...> class F,typename... Ts>
+    inline auto mkVFuncPtr(const std::shared_ptr<F<Ts...>>& f)
+    {
+	return std::shared_ptr<F<typename VType<Ts>::type...>>();
+	// empty, implement corresponding constructors...!!!
+    }
+
+    template <template <typename...> class F,typename... Ts>
+    inline auto mkVFunc(const F<Ts...>& f)
+    {
+	return F<typename VType<Ts>::type...>();
+	// empty, implement corresponding constructors...!!!
+    }
+
+    template <class F>
+    using VFunc = decltype(mkVFunc(std::declval<F>()));
+
+    template <typename T, class F>
+    struct IAccess
+    {
+	typedef T value_type;
+	typedef T in_type;
+	static constexpr size_t VSIZE = sizeof(value_type) / sizeof(in_type);
+
+	template <typename Op, class ExtType>
+	static inline void f(T*& t, size_t pos, const Op& op, ExtType e)
+	{
+	    F::selfApply(t[pos],op.get(e));
+	}
+    };
+
+    template <typename T, class F>
+    struct IVAccess
+    {
+	typedef typename VType<T>::type value_type;
+	typedef T in_type;
+	static constexpr size_t VSIZE = sizeof(value_type) / sizeof(in_type);
+
+	template <typename Op, class ExtType>
+	static inline void f(T*& t, size_t pos, const Op& op, ExtType e)
+	{
+	    //VCHECK(pos);
+	    VFunc<F>::selfApply(*reinterpret_cast<value_type*>(t+pos),op.template vget<value_type>(e));
+	}
+    };
+
+    template <typename T>
+    using xxxplus = plus<T>;
+    
+    template <typename T>
+    using IAssign = IAccess<T,identity<T>>;
+
+    template <typename T>
+    using IPlus = IAccess<T,plus<T>>;
+
+    template <typename T>
+    using IVAssign = IVAccess<T,identity<T>>;
+
+    template <typename T>
+    using IVPlus = IVAccess<T,plus<T>>;
+    /*
+    struct IAssign
+    {
+	template <typename T, typename Op, class ExtType>
+	static inline void f(T*& t, size_t pos, const Op& op, ExtType e)
+	{
+	    t[pos] = op.get(e);
+	}
+    };
+
+    template <typename V>
+    struct IVAssign
+    {
+	template <typename T, typename Op, class ExtType>
+	static inline void f(T*& t, size_t pos, const Op& op, ExtType e)
+	{
+	    reinterpret_cast<V*>(t)[pos] = op.template vget<V>(e);
+	}
+    };
+    
+    struct IPlus
+    {
+	template <typename T, typename Op, class ExtType>
+	static inline void f(T*& t, size_t pos, const Op& op, ExtType e)
+	{
+	    t[pos] += op.get(e);
+	}
+    };
+
+    template <typename V>
+    struct IVPlus
+    {
+	template <typename T, typename Op, class ExtType>
+	static inline void f(T*& t, size_t pos, const Op& op, ExtType e)
+	{
+	    reinterpret_cast<V*>(t)[pos] += op.template vget<V>(e);
+	}
+    };
+    */
+
+    
+    template <typename T, class IOp, class Target, class OpClass, OpIndexAff OIA=OpIndexAff::EXTERN>
+    class AssignmentExpr : public ExpressionBase
     {
     private:
-        AssignmentExpr2() = default;
+        AssignmentExpr() = default;
 
         Target mTar;
         OpClass mSec;
@@ -220,18 +336,19 @@ namespace MultiArrayTools
     public:
 
         static constexpr size_t LAYER = 0;
+        static constexpr size_t NHLAYER = 0;
         static constexpr size_t SIZE = Target::SIZE + OpClass::SIZE;
         typedef decltype(mTar.rootSteps(0).extend( mSec.rootSteps(0) )) ExtType;
 	    
-        AssignmentExpr2(T* dataPtr, const Target& tar, const OpClass& sec);
+        AssignmentExpr(T* dataPtr, const Target& tar, const OpClass& sec);
-        AssignmentExpr2(const AssignmentExpr2& in) = default;
+        AssignmentExpr(const AssignmentExpr& in) = default;
-        AssignmentExpr2(AssignmentExpr2&& in) = default;
+        AssignmentExpr(AssignmentExpr&& in) = default;
-        AssignmentExpr2& operator=(const AssignmentExpr2& in) = default;
+        AssignmentExpr& operator=(const AssignmentExpr& in) = default;
-        AssignmentExpr2& operator=(AssignmentExpr2&& in) = default;
+        AssignmentExpr& operator=(AssignmentExpr&& in) = default;
 
 	virtual std::shared_ptr<ExpressionBase> deepCopy() const override final
 	{
-	    return std::make_shared<AssignmentExpr2<T,Target,OpClass,OIA>>(*this);
+	    return std::make_shared<AssignmentExpr<T,IOp,Target,OpClass,OIA>>(*this);
 	}
 	
         inline void operator()(size_t start = 0); 
@@ -244,6 +361,12 @@ namespace MultiArrayTools
         inline DExt dExtension() const override final;
     };
 
+    template <typename T, class Target, class OpClass, OpIndexAff OIA=OpIndexAff::EXTERN>
+    using AssignmentExpr2 = AssignmentExpr<T,IAssign<T>,Target,OpClass,OIA>;
+
+    template <typename T, class Target, class OpClass, OpIndexAff OIA=OpIndexAff::EXTERN>
+    using AddExpr = AssignmentExpr<T,IPlus<T>,Target,OpClass,OIA>;
+
     template <typename T, class... Ops>
     class MOp
     {
@@ -253,6 +376,7 @@ namespace MultiArrayTools
 
     public:
         static constexpr size_t LAYER = 0;
+	static constexpr size_t NHLAYER = 0;
 	static constexpr size_t SIZE = RootSum<Ops...>::SIZE;
         typedef decltype(RootSumN<sizeof...(Ops)-1>::rootSteps(mOps,0) ) ExtType;
 
@@ -264,6 +388,10 @@ namespace MultiArrayTools
 	MOp& operator=(MOp&& in) = default;
 
         inline size_t get(ExtType last) const;
+
+	template <typename V>
+	inline size_t vget(ExtType last) const { return get(last); }
+	
 	inline MOp& set(ExtType last);
         auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
 
@@ -288,6 +416,7 @@ namespace MultiArrayTools
     public:
 
         static constexpr size_t LAYER = 0;
+	static constexpr size_t NHLAYER = 0;
         static constexpr size_t SIZE = OpClass::SIZE + NextExpr::SIZE;
         typedef decltype(mSec.rootSteps(0).extend( mNExpr.rootSteps(0) ) ) ExtType;
 	    
@@ -304,6 +433,10 @@ namespace MultiArrayTools
 
         inline void operator()(size_t start = 0);
         inline void get(ExtType last);
+
+	template <typename V>
+        inline void vget(ExtType last) { get(last); }
+	
         inline void operator()(size_t start, ExtType last);
         auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
 
@@ -325,48 +458,6 @@ namespace MultiArrayTools
         return MOp<T,Ops...>(exprs...);
     }
 
-    //template <typename T, class OpClass>
-    template <typename T, class Target, class OpClass, OpIndexAff OIA=OpIndexAff::EXTERN>
-    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<ExpressionBase> deepCopy() const override final
-	{
-	    return std::make_shared<AddExpr<T,Target,OpClass,OIA>>(*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 <typename T, class... Ranges>
     class ConstOperationRoot : public OperationTemplate<T,ConstOperationRoot<T,Ranges...> >
     {
@@ -379,6 +470,7 @@ namespace MultiArrayTools
 
 	static constexpr size_t SIZE = 1;
         static constexpr bool CONT = true;
+        static constexpr bool VABLE = true;
 	
 	ConstOperationRoot(const MultiArrayBase<T,Ranges...>& ma,
 			   const std::shared_ptr<typename Ranges::IndexType>&... indices);
@@ -391,6 +483,9 @@ namespace MultiArrayTools
 	template <class ET>
 	inline const T& get(ET pos) const;
 
+	template <typename V, class ET>
+	inline const V& vget(ET pos) const;
+	
 	template <class ET>
 	inline ConstOperationRoot& set(ET pos);
 
@@ -424,12 +519,16 @@ namespace MultiArrayTools
 
 	static constexpr size_t SIZE = Op::SIZE;
         static constexpr bool CONT = false;
+        static constexpr bool VABLE = false;
 	
 	StaticCast(const Op& op);
 
 	template <class ET>
 	inline T get(ET pos) const;
 
+	template <typename V, class ET>
+	inline V vget(ET pos) const;
+
 	template <class ET>
 	inline StaticCast& set(ET pos);
 
@@ -459,12 +558,16 @@ namespace MultiArrayTools
 	
 	static constexpr size_t SIZE = 1;
         static constexpr bool CONT = false;
+        static constexpr bool VABLE = false;
         
 	MetaOperationRoot(const std::shared_ptr<IndexType>& ind);
 
 	template <class ET>
 	inline value_type get(ET pos) const;
 
+	template <typename V, class ET>
+	inline V vget(ET pos) const;
+
 	template <class ET>
 	inline MetaOperationRoot& set(ET pos);
 
@@ -491,6 +594,7 @@ namespace MultiArrayTools
 
 	static constexpr size_t SIZE = 1;
         static constexpr bool CONT = true;
+        static constexpr bool VABLE = true;
 
     private:
 
@@ -507,29 +611,39 @@ namespace MultiArrayTools
 
 	OperationRoot(T* data, const IndexType& ind);
 
+        template <class IOp, class OpClass>
+        auto asx(const OpClass& in) const
+            -> decltype(mIndex.ifor(1,in.loop(AssignmentExpr<T,IOp,OperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
+                                              (mOrigDataPtr,*this,in))).template vec<IOp::VSIZE>());
+
+        template <class IOp, class OpClass>
+        auto asxExpr(const OpClass& in) const
+            -> decltype(in.loop(AssignmentExpr<T,IOp,OperationRoot<T,Ranges...>,OpClass>(mOrigDataPtr,*this,in)));
+            
+        template <class IOp, class OpClass, class Index>
+        auto asx(const OpClass& in, const std::shared_ptr<Index>& i) const
+            -> decltype(i->ifor(1,in.loop(AssignmentExpr<T,IOp,OperationRoot<T,Ranges...>,OpClass>
+                                          (mOrigDataPtr,*this,in))).template vec<IOp::VSIZE>());
+            
         template <class OpClass>
         auto assign(const OpClass& in) const
-            -> decltype(mIndex.ifor(1,in.loop(AssignmentExpr2<T,OperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
+	    -> decltype(this->template asx<IAssign<T>>(in));
-                                              (mOrigDataPtr,*this,in))));
 
 	template <class OpClass>
         auto assignExpr(const OpClass& in) const
-            -> decltype(in.loop(AssignmentExpr2<T,OperationRoot<T,Ranges...>,OpClass>(mOrigDataPtr,*this,in)));
+	    -> decltype(this->template asxExpr<IAssign<T>>(in));
 	    
         template <class OpClass, class Index>
         auto assign(const OpClass& in, const std::shared_ptr<Index>& i) const
-            -> decltype(i->ifor(1,in.loop(AssignmentExpr2<T,OperationRoot<T,Ranges...>,OpClass>
+	    -> decltype(this->template asx<IAssign<T>>(in,i));
-                                          (mOrigDataPtr,*this,in))));
 	
 	template <class OpClass>
         auto plus(const OpClass& in) const
-            -> decltype(mIndex.ifor(1,in.loop(AddExpr<T,OperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
+	    -> decltype(this->template asx<IPlus<T>>(in));
-                                              (mOrigDataPtr,*this,in))));
 
         template <class OpClass, class Index>
         auto plus(const OpClass& in, const std::shared_ptr<Index>& i) const
-            -> decltype(i->ifor(1,in.loop(AddExpr<T,OperationRoot<T,Ranges...>,OpClass>
+	    -> decltype(this->template asx<IPlus<T>>(in,i));
-                                          (mOrigDataPtr,*this,in))));
             
         template <class OpClass>
         OperationRoot& operator=(const OpClass& in);
@@ -544,6 +658,9 @@ namespace MultiArrayTools
 	template <class ET>
 	inline T& get(ET pos) const;
 
+	template <typename V, class ET>
+	inline V& vget(ET pos) const;
+	
 	template <class ET>
 	inline OperationRoot& set(ET pos);
 
@@ -572,6 +689,7 @@ namespace MultiArrayTools
 
 	static constexpr size_t SIZE = 1;
         static constexpr bool CONT = true;
+        static constexpr bool VABLE = true;
 
     private:
 
@@ -585,25 +703,39 @@ namespace MultiArrayTools
 
 	ParallelOperationRoot(T* data, const IndexType& ind);
 
+        template <class IOp, class OpClass>
+        auto asx(const OpClass& in) const
+            -> decltype(mIndex.pifor(1,in.loop(AssignmentExpr<T,IOp,ParallelOperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
+                                              (mOrigDataPtr,*this,in))).template vec<IOp::VSIZE>());
+
+        template <class IOp, class OpClass>
+        auto asxExpr(const OpClass& in) const
+            -> decltype(in.loop(AssignmentExpr<T,IOp,ParallelOperationRoot<T,Ranges...>,OpClass>(mOrigDataPtr,*this,in)));
+            
+        template <class IOp, class OpClass, class Index>
+        auto asx(const OpClass& in, const std::shared_ptr<Index>& i) const
+            -> decltype(i->pifor(1,in.loop(AssignmentExpr<T,IOp,ParallelOperationRoot<T,Ranges...>,OpClass>
+                                          (mOrigDataPtr,*this,in))).template vec<IOp::VSIZE>());
+            
         template <class OpClass>
-        auto assign(const OpClass& in)
+        auto assign(const OpClass& in) const
-            -> decltype(mIndex.pifor(1,in.loop(AssignmentExpr2<T,ParallelOperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
+	    -> decltype(this->template asx<IAssign<T>>(in));
-                                               (mOrigDataPtr,*this,in))));
+
+	template <class OpClass>
+        auto assignExpr(const OpClass& in) const
+	    -> decltype(this->template asxExpr<IAssign<T>>(in));
 	    
         template <class OpClass, class Index>
         auto assign(const OpClass& in, const std::shared_ptr<Index>& i) const
-            -> decltype(i->pifor(1,in.loop(AssignmentExpr2<T,ParallelOperationRoot<T,Ranges...>,OpClass>
+	    -> decltype(this->template asx<IAssign<T>>(in,i));
-                                          (mOrigDataPtr,*this,in))));
 	
 	template <class OpClass>
-        auto plus(const OpClass& in)
+        auto plus(const OpClass& in) const
-            -> decltype(mIndex.pifor(1,in.loop(AddExpr<T,ParallelOperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
+	    -> decltype(this->template asx<IPlus<T>>(in));
-                                               (mOrigDataPtr,*this,in))));
 
         template <class OpClass, class Index>
         auto plus(const OpClass& in, const std::shared_ptr<Index>& i) const
-            -> decltype(i->pifor(1,in.loop(AddExpr<T,ParallelOperationRoot<T,Ranges...>,OpClass>
+	    -> decltype(this->template asx<IPlus<T>>(in,i));
-					   (mOrigDataPtr,*this,in))));
 
         template <class OpClass>
         ParallelOperationRoot& operator=(const OpClass& in);
@@ -616,6 +748,9 @@ namespace MultiArrayTools
 	template <class ET>
 	inline T& get(ET pos) const;
 	 
+	template <typename V, class ET>
+	inline V& vget(ET pos) const;
+
 	template <class ET>
 	inline ParallelOperationRoot& set(ET pos);
 
@@ -639,12 +774,16 @@ namespace MultiArrayTools
 
 	static constexpr size_t SIZE = 0;
         static constexpr bool CONT = true;
+        static constexpr bool VABLE = false;
 
 	OperationValue(const T& val);
 
 	template <class ET>
 	inline const T& get(ET pos) const;
 
+	template <typename V, class ET>
+	inline V vget(ET pos) const;
+
 	template <class ET>
 	inline OperationValue& set(ET pos);
 
@@ -661,7 +800,6 @@ namespace MultiArrayTools
     class OperationPointer : public OperationTemplate<const T*,OperationPointer<T,Op>>
     {
     public:
-
 	typedef T value_type;
 	typedef OperationTemplate<const T*,OperationPointer<T,Op>> OT;
 	
@@ -695,6 +833,18 @@ namespace MultiArrayTools
 	T const** data() const { assert(0); return nullptr; }
     };
 
+    template <typename T, class Op> 
+    inline constexpr bool isVAble()
+    {
+	return Op::VABLE and std::is_same<T,typename Op::value_type>::value;
+    }
+
+    template <typename T, class Op1, class Op2, class... Ops>
+    inline constexpr bool isVAble()
+    {
+	return Op1::VABLE and std::is_same<T,typename Op1::value_type>::value and isVAble<T,Op2,Ops...>();
+    }
+    
     template <typename T, class OpFunction, class... Ops>
     class Operation : public OperationTemplate<T,Operation<T,OpFunction,Ops...> >
     {
@@ -707,6 +857,7 @@ namespace MultiArrayTools
 	static constexpr size_t SIZE = RootSum<Ops...>::SIZE;
 	static constexpr bool FISSTATIC = OpFunction::FISSTATIC;
         static constexpr bool CONT = false;
+        static constexpr bool VABLE = isVAble<T,Ops...>();
 
     private:
 	std::tuple<Ops...> mOps;
@@ -721,6 +872,9 @@ namespace MultiArrayTools
 	template <class ET>
 	inline auto get(ET pos) const;
 
+	template <typename V, class ET>
+	inline auto vget(ET pos) const;
+
 	template <class ET>
 	inline Operation& set(ET pos);
 
@@ -778,6 +932,7 @@ namespace MultiArrayTools
 
 	static constexpr size_t SIZE = Op::SIZE;
         static constexpr bool CONT = Op::CONT;
+        static constexpr bool VABLE = Op::VABLE;
         
     private:
 
@@ -793,6 +948,10 @@ namespace MultiArrayTools
 	inline auto get(ET pos) const
             -> decltype(mOp.template get<ET>(pos));
 
+	template <typename V, class ET>
+	inline auto vget(ET pos) const
+            -> decltype(mOp.template vget<V,ET>(pos));
+
 	template <class ET>
 	inline Contraction& set(ET pos);
 
@@ -818,6 +977,7 @@ namespace MultiArrayTools
 
 	static constexpr size_t SIZE = Op::SIZE;
         static constexpr bool CONT = false;
+	static constexpr bool VABLE = false;
 
     private:
 

src/include/pack_num.h

@@ -81,6 +81,17 @@ namespace MultiArrayHelper
 		( f, pos, ops, std::get<N>(ops).get(Getter<NEXT>::template getX<ETuple>( pos )), args...);
 	}
 
+	template <size_t LAST, typename V, class ETuple, class OpTuple, class OpFunction, typename... Args>
+	static inline auto mkVOpExpr(std::shared_ptr<OpFunction> f, const ETuple& pos, const OpTuple& ops, Args... args)
+	{
+	    typedef typename std::remove_reference<decltype(std::get<N>(ops))>::type NextOpType;
+	    static_assert(LAST >= NextOpType::SIZE, "inconsistent array positions");
+	    static constexpr size_t NEXT = LAST - NextOpType::SIZE;
+	    typedef decltype(std::get<N>(ops).template vget<V>(Getter<NEXT>::template getX<ETuple>( pos ))) ArgT;
+	    return PackNum<N-1>::template mkVOpExpr<NEXT,V,ETuple,OpTuple,OpFunction,ArgT,Args...>
+		( f, pos, ops, std::get<N>(ops).template vget<V>(Getter<NEXT>::template getX<ETuple>( pos )), args...);
+	}
+
 	template <class OpTuple, class Expr>
 	static auto mkLoop( const OpTuple& ot, Expr exp )
 	    -> decltype(std::get<N>(ot).loop( PackNum<N-1>::mkLoop(ot,exp) ))
@@ -168,6 +179,16 @@ namespace MultiArrayHelper
 	    //return OpFunction::apply(std::get<0>(ops).get(Getter<0>::template getX<ETuple>( pos )), args...);
 	}
 
+	template <size_t LAST, typename V, class ETuple, class OpTuple, class OpFunction, typename... Args>
+	static inline auto mkVOpExpr(std::shared_ptr<OpFunction> f, const ETuple& pos, const OpTuple& ops, const Args&... args)
+	{
+	    typedef typename std::remove_reference<decltype(std::get<0>(ops))>::type NextOpType;
+	    static constexpr size_t NEXT = LAST - NextOpType::SIZE;
+	    static_assert(NEXT == 0, "inconsistent array positions");
+	    typedef decltype(std::get<0>(ops).template vget<V>(Getter<0>::template getX<ETuple>( pos ))) ArgT;
+	    return Application<OpFunction::FISSTATIC>::template apply<OpFunction,ArgT,Args...>(f, std::get<0>(ops).template vget<V>(Getter<0>::template getX<ETuple>( pos )), args...);
+	}
+
 	template <class OpTuple, class Expr>
 	static auto mkLoop( const OpTuple& ot, Expr exp )
 	    -> decltype(std::get<0>(ot).loop( exp ))

src/include/type_operations.h

@@ -10,6 +10,7 @@
 #include "mbase_def.h"
 
 #include "pack_num.h"
+#include <cmath>
 
 namespace MultiArrayTools
 {
@@ -471,6 +472,48 @@ namespace MultiArrayTools
 	return o;
     }
 
+    inline double xpow(const double& b, const double& e)
+    {
+	return pow(b,e);
+    }
+
+    inline v256 pow(const v256& b, const v256& e)
+    {
+	v256 out;
+	for(int i = 0; i < 4; i++){
+	    out._x[i] = xpow(b._x[i],e._x[i]);
+	}
+	return out;
+    }
+
+    inline double xexp(const double& a)
+    {
+	return exp(a);
+    }
+
+    inline v256 exp(const v256& a)
+    {
+	v256 out;
+	for(int i = 0; i < 4; i++){
+	    out._x[i] = xexp(a._x[i]);
+	}
+	return out;
+    }
+
+    inline double xsqrt(const double& a)
+    {
+	return sqrt(a);
+    }
+
+    inline v256 sqrt(const v256& a)
+    {
+	v256 out;
+	for(int i = 0; i < 4; i++){
+	    out._x[i] = xsqrt(a._x[i]);
+	}
+	return out;
+    }
+
 } // namespace MultiArrayTools
 
 #endif

src/include/xfor/exttype.h

@@ -86,6 +86,16 @@ namespace MultiArrayHelper
             -> decltype(Getter<N>::getX(*this))
         { return Getter<N>::getX(*this); }
 
+	inline bool operator==(const MExt& in) const
+	{
+	    return mExt == in.mExt and mNext == in.mNext;
+	}
+
+	inline bool operator==(size_t in) const
+	{
+	    return mExt == in and mNext == in;
+	}
+
 	inline MExt operator+(const MExt& in) const;
 	inline MExt operator+(const None& in) const;
 	inline MExt operator*(size_t in) const;
@@ -111,6 +121,16 @@ namespace MultiArrayHelper
         
         static constexpr size_t SIZE = 0;
 
+	inline bool operator==(const None& in) const
+	{
+	    return true;
+	}
+
+	inline bool operator==(size_t in) const
+	{
+	    return true; // CHECK!!!
+	}
+
         inline None operator+(const None& in) const { return None(); }
 	inline None operator*(size_t in) const { return None(); }
 
@@ -158,6 +178,16 @@ namespace MultiArrayHelper
         template <class Y>
         inline MExt(const MExt<Y>& y);
 
+	inline bool operator==(const MExt& in) const
+	{
+	    return mExt == in.mExt;
+	}
+
+	inline bool operator==(size_t in) const
+	{
+	    return mExt == in;
+	}
+
 	inline const size_t& val() const;
 	inline None next() const { return None(); }
 

src/include/xfor/for_type.h

@@ -10,9 +10,12 @@ namespace MultiArrayHelper
 	HIDDEN = 1	
     };
 
-    template <class IndexClass, class Expr, ForType FT = ForType::DEFAULT>
+    template <class IndexClass, class Expr, ForType FT = ForType::DEFAULT, size_t DIV = 1>
     class For;
 
+    template <class IndexClass, class Expr, size_t DIV = 1>
+    class PFor;
+
     
 } // end namespace MultiArrayHelper
 

src/include/xfor/iloop.h

@@ -116,6 +116,7 @@ namespace MultiArrayHelper
             
     public:
         //static constexpr size_t SIZE = NN<LTpSize-1>::lsize(std::declval<LTp>());
+	static constexpr size_t NHLAYER = 10; // some large value
         static constexpr size_t SIZE = NN<LTpSize-1>::template LSIZE<LTp>();
         static constexpr bool CONT = false;
         typedef decltype(NN<LTpSize-1>::rootSteps(mLTp)) ExtType;
@@ -248,6 +249,8 @@ namespace MultiArrayHelper
         static constexpr size_t VarTpSize = LType::VarTpSize;
             
     public:
+
+	static constexpr size_t NHLAYER = 10; // some large value
         static constexpr size_t SIZE = LType::SIZE;
         static constexpr bool CONT = LType::CONT;
         typedef typename LType::ExtType ExtType;

src/include/xfor/xfor.h

@@ -35,6 +35,9 @@ namespace MultiArrayHelper
 	virtual size_t size() const = 0;
 	virtual const size_t& val() const = 0;
 	//virtual size_t rootSteps() const = 0;
+	virtual bool operator==(const ExtBase& in) const = 0;
+	virtual bool operator==(size_t in) 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;
@@ -75,6 +78,12 @@ namespace MultiArrayHelper
 	virtual const size_t& val() const override final { return mExt.val(); }
 	virtual void zero() override final { mExt.zero(); }
 
+	virtual bool operator==(const ExtBase& in) const override final
+	{ return mExt == dynamic_cast<const ExtT<ExtType>&>(in).mExt; }
+	
+	virtual bool operator==(size_t in) const override final
+	{ return mExt == in; }
+
 	virtual DExt operator+(const ExtBase& in) const override final
 	{ return std::make_shared<ExtT<ExtType>>( mExt + dynamic_cast<const ExtT<ExtType>&>(in).mExt ); }
 	virtual DExt operator*(size_t in) const override final
@@ -122,6 +131,11 @@ namespace MultiArrayHelper
 	template <class Y>
 	DExtTX(const Y& y) : mDExt(std::make_shared<ExtT<Y>>(y)) {}
         */
+	bool operator==(const DExtTX& in) const
+	{ return *mDExt == *in.mDExt and mNext == in.mNext; }
+	
+	bool operator==(size_t in) const
+	{ return *mDExt == in and mNext == in; }
 
         template <class Y>
         DExtTX(const DExtTX<Y>& in) : mDExt(in.mDExt), mNext(in.mNext) {}
@@ -173,6 +187,9 @@ namespace MultiArrayHelper
 	ExpressionBase& operator=(const ExpressionBase& in) = default;
 	ExpressionBase& operator=(ExpressionBase&& in) = default;
 
+	//virtual size_t divResid() const { return 1; }
+	virtual std::intptr_t vI() const { return 0; }
+	
 	virtual std::shared_ptr<ExpressionBase> deepCopy() const = 0;
 	
 	virtual void operator()(size_t mlast, DExt last) = 0;
@@ -215,20 +232,20 @@ namespace MultiArrayHelper
     template <size_t ISSTATIC>
     struct ForBound
     {
-	template <size_t BOUND>
+	template <size_t BOUND, size_t DIV = 1>
 	static inline size_t bound(size_t bound)
 	{
-	    return bound;
+	    return bound / DIV;
 	}
     };
 
     template <>
     struct ForBound<1>
     {
-	template <size_t BOUND>
+	template <size_t BOUND, size_t DIV = 1>
 	static constexpr size_t bound(size_t bound)
 	{
-	    return BOUND;
+	    return BOUND / DIV;
 	}
     };
 
@@ -254,6 +271,7 @@ namespace MultiArrayHelper
 	
 	static constexpr size_t LAYER = Expr::LAYER + 1;
 	static constexpr size_t SIZE = Expr::SIZE;
+	static constexpr size_t NHLAYER = Expr::NHLAYER + 1;
 
 	SingleExpression(const SingleExpression& in) = default;
 	SingleExpression& operator=(const SingleExpression& in) = default;
@@ -271,6 +289,9 @@ namespace MultiArrayHelper
 	    return std::make_shared<SingleExpression<IndexClass,Expr>>(*this);
 	}
 
+	template <size_t VS>
+	inline auto vec() const { return *this; }
+
 	inline void operator()(size_t mlast, DExt last) override final;
 	inline void operator()(size_t mlast, ExtType last);
 	inline void operator()(size_t mlast = 0) override final;
@@ -306,6 +327,7 @@ namespace MultiArrayHelper
 	
 	static constexpr size_t LAYER = Expr::LAYER + 1;
 	static constexpr size_t SIZE = Expr::SIZE + 1;
+	static constexpr size_t NHLAYER = Expr::NHLAYER + 1;
 
 	SubExpr(const SubExpr& in) = default;
 	SubExpr& operator=(const SubExpr& in) = default;
@@ -324,6 +346,9 @@ namespace MultiArrayHelper
 	    return std::make_shared<SubExpr<IndexClass,Expr>>(*this);
 	}
 
+	template <size_t VS>
+	inline auto vec() const { return *this; }
+
 	inline void operator()(size_t mlast, DExt last) override final;
 	inline void operator()(size_t mlast, ExtType last) ;
 	inline void operator()(size_t mlast = 0) override final;
@@ -335,15 +360,95 @@ namespace MultiArrayHelper
 	auto extension() const -> ExtType;
     };
 
-    template <class IndexClass, class Expr>
+    template <size_t LAYER, bool ISV>
-    class PFor;
+    struct MkVFor
+    {
+	template <size_t DIV, class IndexClass, class Expr>
+	using ptype = PFor<IndexClass,Expr,1>;
 
-    template <class IndexClass, class Expr, ForType FT>
+	template <size_t DIV, class IndexClass, class Expr, ForType FT>
+	using type = For<IndexClass,Expr,FT,1>;
+    };
+
+    template <>
+    struct MkVFor<1,true>
+    {
+	template <size_t DIV, class IndexClass, class Expr>
+	using ptype = PFor<IndexClass,Expr,DIV>;
+
+	template <size_t DIV, class IndexClass, class Expr, ForType FT>
+	using type = For<IndexClass,Expr,FT,DIV>;
+    };
+
+    template <size_t LAYER>
+    struct MkVExpr
+    {
+	template <size_t VS, class Expr>
+	static auto mk(const Expr& e)
+	{
+	    return e.template vec<VS>();
+	}
+
+	template <class Expr>
+	static inline size_t divResid(const Expr& e)
+	{
+	    return e.divResid();
+	}
+    };
+
+    template <>
+    struct MkVExpr<1>
+    {
+	template <size_t VS, class Expr>
+	static auto mk(const Expr& e)
+	{
+	    return e; // terminate
+	}
+
+	template <class Expr>
+	static inline size_t divResid(const Expr& e)
+	{
+	    return 0;
+	}
+    };
+
+    template <ForType FT, size_t LAYER>
+    struct NHLayer
+    {
+	template <class Expr>
+	static constexpr size_t get()
+	{
+	    return Expr::NHLAYER + 1;
+	}
+    };
+
+    template <size_t LAYER>
+    struct NHLayer<ForType::HIDDEN,LAYER>
+    {
+	template <class Expr>
+	static constexpr size_t get()
+	{
+	    return 0;
+	}
+    };
+
+    template <>
+    struct NHLayer<ForType::DEFAULT,1>
+    {
+	template <class Expr>
+	static constexpr size_t get()
+	{
+	    return Expr::LAYER;
+	}
+    };
+
+    template <class IndexClass, class Expr, ForType FT, size_t DIV>
     class For : public ExpressionBase
     {
     private:
 	For() = default;
 		
+	typedef typename IndexClass::RangeType RangeType;
 	const IndexClass* mIndPtr;
 	size_t mSPos;
 	size_t mMax;
@@ -360,6 +465,8 @@ namespace MultiArrayHelper
 
 	static constexpr size_t LAYER = Expr::LAYER + 1;
 	static constexpr size_t SIZE = Expr::SIZE;
+	static constexpr size_t MAX = RangeType::SIZE / DIV;
+	static constexpr size_t NHLAYER = (FT == ForType::HIDDEN) ? 0 : Expr::NHLAYER + 1;
 	
 	For(const For& in) = default;
 	For& operator=(const For& in) = default;
@@ -374,14 +481,35 @@ namespace MultiArrayHelper
 
 	virtual std::shared_ptr<ExpressionBase> deepCopy() const override final
 	{
-	    return std::make_shared<For<IndexClass,Expr,FT>>(*this);
+	    return std::make_shared<For<IndexClass,Expr,FT,DIV>>(*this);
+	}
+
+	//virtual size_t divResid() const override final { return mMax % DIV + MkVExpr<LAYER>::divResid(mExpr); }
+
+	virtual std::intptr_t vI() const override final
+	{
+	    if(mStep == 1 and NHLAYER == 1 and mMax % DIV == 0){
+	    //if(mStep == 1 and mMax % DIV == 0){
+		//VCHECK(LAYER);
+		//VCHECK(NHLAYER);
+		return reinterpret_cast<std::intptr_t>(mIndPtr);
+	    }
+	    return mExpr.vI();
+	}
+
+	template <size_t VS>
+	auto vec() const
+	{
+	    typedef typename MkVFor<NHLAYER,RangeType::SIZE % DIV == 0 or RangeType::SIZE == static_cast<size_t>(-1)>::
+		template type<VS,IndexClass,decltype(MkVExpr<NHLAYER>::template mk<VS>(mExpr)),FT> oType;
+	    return oType(mIndPtr,mStep,MkVExpr<NHLAYER>::template mk<VS>(mExpr));
 	}
 
 	inline void operator()(size_t mlast, DExt last) override final;
 	inline void operator()(size_t mlast, ExtType last) ;
 	inline void operator()(size_t mlast = 0) override final;
 
-        PFor<IndexClass,Expr> parallel() const;
+        PFor<IndexClass,Expr,DIV> parallel() const;
         
         DExt dRootSteps(std::intptr_t iPtrNum = 0) const override final;
         DExt dExtension() const override final;
@@ -391,12 +519,14 @@ namespace MultiArrayHelper
 
     };
 
-    template <class IndexClass, class Expr>
+   
+    template <class IndexClass, class Expr, size_t DIV>
     class PFor : public ExpressionBase
     {
     private:
 	PFor() = default;
 		
+	typedef typename IndexClass::RangeType RangeType;
 	const IndexClass* mIndPtr;
 	size_t mSPos;
 	size_t mMax;
@@ -413,6 +543,8 @@ namespace MultiArrayHelper
 	
 	static constexpr size_t LAYER = Expr::LAYER + 1;
 	static constexpr size_t SIZE = Expr::SIZE;
+	static constexpr size_t MAX = RangeType::SIZE / DIV;
+	static constexpr size_t NHLAYER = Expr::NHLAYER + 1;
 
 	PFor(const PFor& in) = default;
 	PFor& operator=(const PFor& in) = default;
@@ -425,6 +557,26 @@ namespace MultiArrayHelper
 	PFor(const IndexClass* indPtr,
 	    size_t step, Expr expr);
 
+	//virtual size_t divResid() const override final { return mMax % DIV + MkVExpr<LAYER>::divResid(mExpr); }
+	virtual std::intptr_t vI() const override final
+	{
+	    if(mStep == 1 and NHLAYER == 1 and mMax % DIV == 0){
+	    //if(mStep == 1 and mMax % DIV == 0){
+		//VCHECK(LAYER);
+		//VCHECK(LAYER);
+		return reinterpret_cast<std::intptr_t>(mIndPtr);
+	    }
+	    return mExpr.vI();
+	}
+
+	template <size_t VS>
+	auto vec() const
+	{
+	    typedef typename MkVFor<NHLAYER,RangeType::SIZE % DIV == 0 or RangeType::SIZE == static_cast<size_t>(-1)>::
+		template ptype<VS,IndexClass,decltype(MkVExpr<NHLAYER>::template mk<VS>(mExpr))> oType;
+	    return oType(mIndPtr,mStep,MkVExpr<NHLAYER>::template mk<VS>(mExpr));
+	}
+
 	virtual std::shared_ptr<ExpressionBase> deepCopy() const override final
 	{
 	    return std::make_shared<PFor<IndexClass,Expr>>(*this);
@@ -460,6 +612,7 @@ namespace MultiArrayHelper
 	
 	static constexpr size_t LAYER = 0;
 	static constexpr size_t SIZE = 0;
+	static constexpr size_t NHLAYER = 0;
 
 	DynamicExpression(const DynamicExpression& in) :
 	    mThreadId(omp_get_thread_num()),
@@ -502,6 +655,9 @@ namespace MultiArrayHelper
 	    return std::make_shared<DynamicExpression>(*this);
 	}
 
+	template <size_t VS>
+	inline auto vec() const { return *this; }
+
 	inline void operator()(size_t mlast, DExt last) override final;
         inline void operator()(size_t mlast, DExtT last) { (*this)(mlast,last.get()); }
 	inline void operator()(size_t mlast = 0) override final;
@@ -532,6 +688,7 @@ namespace MultiArrayHelper
 	
 	static constexpr size_t LAYER = Expr::LAYER + 1;
 	static constexpr size_t SIZE = Expr::SIZE;
+	static constexpr size_t NHLAYER = Expr::NHLAYER + 1;
 
 	ExpressionHolder(const ExpressionHolder& in) = default;
 	ExpressionHolder(ExpressionHolder&& in) = default;
@@ -577,50 +734,54 @@ namespace MultiArrayHelper
      *     F o r     *
      *****************/
     
-    template <class IndexClass, class Expr, ForType FT>
+    template <class IndexClass, class Expr, ForType FT, size_t DIV>
-    For<IndexClass,Expr,FT>::For(const std::shared_ptr<IndexClass>& indPtr,
+    For<IndexClass,Expr,FT,DIV>::For(const std::shared_ptr<IndexClass>& indPtr,
 				 size_t step, Expr expr) :
 	mIndPtr(indPtr.get()), mSPos(mIndPtr->pos()), mMax(mIndPtr->max()), mStep(step),
         mExpr(expr), mExt(mExpr.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr )))
     {
+	assert(mMax % DIV == 0);
 	assert(mIndPtr != nullptr);
     }
 
-    template <class IndexClass, class Expr, ForType FT>
+    template <class IndexClass, class Expr, ForType FT, size_t DIV>
-    For<IndexClass,Expr,FT>::For(const IndexClass* indPtr,
+    For<IndexClass,Expr,FT,DIV>::For(const IndexClass* indPtr,
 				 size_t step, Expr expr) :
 	mIndPtr(indPtr), mSPos(mIndPtr->pos()), mMax(mIndPtr->max()), mStep(step),
         mExpr(expr), mExt(mExpr.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr )))
     {
+	//VCHECK(mMax);
+	//VCHECK(DIV);
+	//assert(mMax % DIV == 0);
 	assert(mIndPtr != nullptr);
     }
 
-    template <class IndexClass, class Expr, ForType FT>
+    template <class IndexClass, class Expr, ForType FT, size_t DIV>
-    inline void For<IndexClass,Expr,FT>::operator()(size_t mlast, DExt last)
+    inline void For<IndexClass,Expr,FT,DIV>::operator()(size_t mlast, DExt last)
     {
         operator()(mlast, std::dynamic_pointer_cast<ExtT<ExtType>>(last)->ext());
         //operator()(mlast, *reinterpret_cast<ExtType const*>(last.first));
     }
     
-    template <class IndexClass, class Expr, ForType FT>
+    template <class IndexClass, class Expr, ForType FT, size_t DIV>
-    inline void For<IndexClass,Expr,FT>::operator()(size_t mlast,
+    inline void For<IndexClass,Expr,FT,DIV>::operator()(size_t mlast,
 						    ExtType last)
     {
 	typedef typename IndexClass::RangeType RangeType;
-	for(size_t pos = 0u; pos != ForBound<RangeType::ISSTATIC>::template bound<RangeType::SIZE>(mMax); ++pos){
+	for(size_t pos = 0u; pos != ForBound<RangeType::ISSTATIC>::template bound<RangeType::SIZE,DIV>(mMax); ++pos){
 	    const size_t mnpos = PosForward<FT>::valuex(mlast, mStep, pos);
 	    const ExtType npos = last + mExt*pos;
 	    mExpr(mnpos, npos);
 	}
     }
 
-    template <class IndexClass, class Expr, ForType FT>
+    template <class IndexClass, class Expr, ForType FT, size_t DIV>
-    inline void For<IndexClass,Expr,FT>::operator()(size_t mlast)
+    inline void For<IndexClass,Expr,FT,DIV>::operator()(size_t mlast)
     {
 	typedef typename IndexClass::RangeType RangeType;
 	ExtType last = rootSteps();
 	last.zero();
-	for(size_t pos = 0u; pos != ForBound<RangeType::ISSTATIC>::template bound<RangeType::SIZE>(mMax); ++pos){
+	for(size_t pos = 0u; pos != ForBound<RangeType::ISSTATIC>::template bound<RangeType::SIZE,DIV>(mMax); ++pos){
 	    const size_t mnpos = PosForward<FT>::valuex(mlast, mStep, pos);
 	    const ExtType npos = last + mExt*pos;
 	    mExpr(mnpos, npos);
@@ -628,22 +789,22 @@ namespace MultiArrayHelper
     }
     
     
-    template <class IndexClass, class Expr, ForType FT>
+    template <class IndexClass, class Expr, ForType FT, size_t DIV>
-    auto For<IndexClass,Expr,FT>::rootSteps(std::intptr_t iPtrNum) const
+    auto For<IndexClass,Expr,FT,DIV>::rootSteps(std::intptr_t iPtrNum) const
 	-> ExtType
     {
 	return mExpr.rootSteps(iPtrNum);
     }
 
-    template <class IndexClass, class Expr, ForType FT>
+    template <class IndexClass, class Expr, ForType FT, size_t DIV>
-    auto For<IndexClass,Expr,FT>::extension() const
+    auto For<IndexClass,Expr,FT,DIV>::extension() const
 	-> ExtType
     {
 	return mExt;
     }
 
-    template <class IndexClass, class Expr, ForType FT>
+    template <class IndexClass, class Expr, ForType FT, size_t DIV>
-    DExt For<IndexClass,Expr,FT>::dRootSteps(std::intptr_t iPtrNum) const
+    DExt For<IndexClass,Expr,FT,DIV>::dRootSteps(std::intptr_t iPtrNum) const
     {
         return std::make_shared<ExtT<ExtType>>(rootSteps(iPtrNum));
         //mRootSteps = rootSteps(iPtrNum);
@@ -651,52 +812,54 @@ namespace MultiArrayHelper
 	//					    sizeof(ExtType)/sizeof(size_t));
     }
 
-    template <class IndexClass, class Expr, ForType FT>
+    template <class IndexClass, class Expr, ForType FT, size_t DIV>
-    DExt For<IndexClass,Expr,FT>::dExtension() const
+    DExt For<IndexClass,Expr,FT,DIV>::dExtension() const
     {
         return std::make_shared<ExtT<ExtType>>(mExt);
         //return std::make_pair<size_t const*,size_t>(reinterpret_cast<size_t const*>(&mExt),
 	//					    sizeof(ExtType)/sizeof(size_t));
     }
 
-    template <class IndexClass, class Expr, ForType FT>
+    template <class IndexClass, class Expr, ForType FT, size_t DIV>
-    PFor<IndexClass,Expr> For<IndexClass,Expr,FT>::parallel() const
+    PFor<IndexClass,Expr,DIV> For<IndexClass,Expr,FT,DIV>::parallel() const
     {
         static_assert(FT == ForType::DEFAULT, "hidden for not parallelizable");
-        return PFor<IndexClass,Expr>(mIndPtr, mStep, mExpr);
+        return PFor<IndexClass,Expr,DIV>(mIndPtr, mStep, mExpr);
     }
     
     /******************
      *    P F o r     *
      ******************/
     
-    template <class IndexClass, class Expr>
+    template <class IndexClass, class Expr, size_t DIV>
-    PFor<IndexClass,Expr>::PFor(const std::shared_ptr<IndexClass>& indPtr,
+    PFor<IndexClass,Expr,DIV>::PFor(const std::shared_ptr<IndexClass>& indPtr,
 				 size_t step, Expr expr) :
 	mIndPtr(indPtr.get()), mSPos(mIndPtr->pos()), mMax(mIndPtr->max()), mStep(step),
         mExpr(expr), mExt(mExpr.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr )))
     {
+	//assert(mMax % DIV == 0);
 	assert(mIndPtr != nullptr);
     }
 
-    template <class IndexClass, class Expr>
+    template <class IndexClass, class Expr, size_t DIV>
-    PFor<IndexClass,Expr>::PFor(const IndexClass* indPtr,
+    PFor<IndexClass,Expr,DIV>::PFor(const IndexClass* indPtr,
 				 size_t step, Expr expr) :
 	mIndPtr(indPtr), mSPos(mIndPtr->pos()), mMax(mIndPtr->max()), mStep(step),
         mExpr(expr), mExt(mExpr.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr )))
     {
+	assert(mMax % DIV == 0);
 	assert(mIndPtr != nullptr);
     }
 
-    template <class IndexClass, class Expr>
+    template <class IndexClass, class Expr, size_t DIV>
-    inline void PFor<IndexClass,Expr>::operator()(size_t mlast, DExt last)
+    inline void PFor<IndexClass,Expr,DIV>::operator()(size_t mlast, DExt last)
     {
         operator()(mlast, std::dynamic_pointer_cast<ExtT<ExtType>>(last)->ext());
         //operator()(mlast, *reinterpret_cast<ExtType const*>(last.first));
     }
     
-    template <class IndexClass, class Expr>
+    template <class IndexClass, class Expr, size_t DIV>
-    inline void PFor<IndexClass,Expr>::operator()(size_t mlast,
+    inline void PFor<IndexClass,Expr,DIV>::operator()(size_t mlast,
 						    ExtType last)
     {
         CHECK;
@@ -708,7 +871,7 @@ namespace MultiArrayHelper
         {
             auto expr = mExpr;
 #pragma omp for nowait
-            for(pos = 0; pos < static_cast<int>(ForBound<RangeType::ISSTATIC>::template bound<RangeType::SIZE>(mMax)); pos++){
+            for(pos = 0; pos < static_cast<int>(ForBound<RangeType::ISSTATIC>::template bound<RangeType::SIZE,DIV>(mMax)); pos++){
                 mnpos = PosForward<ForType::DEFAULT>::valuex(mlast, mStep, pos);
                 npos = last + mExt*static_cast<size_t>(pos);
                 expr(mnpos, npos);
@@ -716,11 +879,10 @@ namespace MultiArrayHelper
         }
     }
 
-    template <class IndexClass, class Expr>
+    template <class IndexClass, class Expr, size_t DIV>
-    inline void PFor<IndexClass,Expr>::operator()(size_t mlast)
+    inline void PFor<IndexClass,Expr,DIV>::operator()(size_t mlast)
     {
         CHECK;
-	typedef typename IndexClass::RangeType RangeType;
 	ExtType last = rootSteps();
 	last.zero();
         int pos = 0;
@@ -731,7 +893,7 @@ namespace MultiArrayHelper
         {
             auto expr = mExpr;
 #pragma omp for nowait
-            for(pos = 0; pos < static_cast<int>(ForBound<RangeType::ISSTATIC>::template bound<RangeType::SIZE>(mMax)); pos++){
+            for(pos = 0; pos < static_cast<int>(ForBound<RangeType::ISSTATIC>::template bound<RangeType::SIZE,DIV>(mMax)); pos++){
                 mnpos = PosForward<ForType::DEFAULT>::valuex(mlast, mStep, pos);
                 npos = last + mExt*static_cast<size_t>(pos);
                 expr(mnpos, npos);
@@ -740,22 +902,22 @@ namespace MultiArrayHelper
     }
     
     
-    template <class IndexClass, class Expr>
+    template <class IndexClass, class Expr, size_t DIV>
-    auto PFor<IndexClass,Expr>::rootSteps(std::intptr_t iPtrNum) const
+    auto PFor<IndexClass,Expr,DIV>::rootSteps(std::intptr_t iPtrNum) const
 	-> ExtType
     {
 	return mExpr.rootSteps(iPtrNum);
     }
 
-    template <class IndexClass, class Expr>
+    template <class IndexClass, class Expr, size_t DIV>
-    auto PFor<IndexClass,Expr>::extension() const
+    auto PFor<IndexClass,Expr,DIV>::extension() const
 	-> ExtType
     {
 	return mExt;
     }
 
-    template <class IndexClass, class Expr>
+    template <class IndexClass, class Expr, size_t DIV>
-    DExt PFor<IndexClass,Expr>::dRootSteps(std::intptr_t iPtrNum) const
+    DExt PFor<IndexClass,Expr,DIV>::dRootSteps(std::intptr_t iPtrNum) const
     {
         return std::make_shared<ExtT<ExtType>>(rootSteps(iPtrNum));
         //mRootSteps = rootSteps(iPtrNum);
@@ -763,8 +925,8 @@ namespace MultiArrayHelper
 	//					    sizeof(ExtType)/sizeof(size_t));
     }
 
-    template <class IndexClass, class Expr>
+    template <class IndexClass, class Expr, size_t DIV>
-    DExt PFor<IndexClass,Expr>::dExtension() const
+    DExt PFor<IndexClass,Expr,DIV>::dExtension() const
     {
         return std::make_shared<ExtT<ExtType>>(mExt);
         //return std::make_pair<size_t const*,size_t>(reinterpret_cast<size_t const*>(&mExt),

src/tests/op3_unit_test.cc

@@ -45,6 +45,7 @@ namespace
 
 	std::clock_t begin = std::clock();
 	res1(delta, deltap) += ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
+	//res1(delta, deltap) += ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap);
 	std::clock_t end = std::clock();
 	std::cout << "MultiArray time: " << static_cast<double>( end - begin ) / CLOCKS_PER_SEC
 		  << std::endl;

src/tests/op_perf_test.cc

@@ -174,8 +174,8 @@ namespace {
 	auto beta = MAT::getIndex<SR>();
 	auto gamma = MAT::getIndex<SR>();
 	auto delta = MAT::getIndex<SR>();
-	//auto deltap = MAT::getIndex<SR>();
+	auto deltap = MAT::getIndex<SR>();
-	auto deltap = MAT::getIndex<GenSingleRange<size_t,SpaceType::NONE,1>>();
+	//auto deltap = MAT::getIndex<GenSingleRange<size_t,SpaceType::NONE,1>>();
 	
 	auto mix = MAT::mkMIndex( jj, alpha, beta, gamma );
 
@@ -217,8 +217,8 @@ namespace {
         
         auto begin = std::chrono::system_clock::now();
 	//for(size_t i = 0; i != os; ++i){
-        //res1(ii ,delta, deltap).par() += ma(ii, delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
+        res1(ii ,delta, deltap).par() += ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
-        tcast<v256>(res1)(ii ,delta, deltap).par() += tcast<v256>(ma)(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
+        //tcast<v256>(res1)(ii ,delta, deltap).par() += tcast<v256>(ma)(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
 	//}
         auto end = std::chrono::system_clock::now();
 	std::cout << "MultiArray time: " << std::chrono::duration<double>(end-begin).count()