compiles, basically works correctly (except for 'for' extensions setting, which is not finished)

src/multi_array_operation.h

@@ -166,8 +166,8 @@ namespace MultiArrayTools
 	MBlock<T>& get();
 	const Block<T>& get() const;
 
-	T& get(size_t pos);
+	inline T& get(size_t pos);
-	const T& get(size_t pos) const;
+	inline const T& get(size_t pos) const;
 	
 	std::vector<BTSS> block(const IndexInfo* blockIndex, bool init = false) const;
 	const OperationMaster& block() const;
@@ -178,7 +178,7 @@ namespace MultiArrayTools
 	void performAssignment(std::intptr_t blockIndexNum);
 	OpClass const& mSecond;
 	MutableMultiArrayBase<T,Ranges...>& mArrayRef;
-	T* mData;
+	std::vector<T>& mData;
 	std::shared_ptr<IndexType> mIndex;
 	IndexInfo mIInfo;
 	mutable bType mBlock;
@@ -224,7 +224,7 @@ namespace MultiArrayTools
 		const std::shared_ptr<typename Ranges::IndexType>&... indices);
 
 	MultiArrayBase<T,Ranges...> const& mArrayRef;
-	const T* mData;
+	const std::vector<T>& mData;
 	std::shared_ptr<IndexType> mIndex;
 	IndexInfo mIInfo;
  	mutable bType mBlock;
@@ -276,7 +276,7 @@ namespace MultiArrayTools
 		const std::shared_ptr<typename Ranges::IndexType>&... indices);
 
 	MutableMultiArrayBase<T,Ranges...>& mArrayRef;
-	T* mData;
+	std::vector<T>& mData;
 	std::shared_ptr<IndexType> mIndex;
 	IndexInfo mIInfo;
 	mutable bType mBlock;
@@ -399,7 +399,7 @@ namespace MultiArrayTools
 	    -> decltype(mOp.rootSteps(iPtrNum));
 
 	template <class Expr>
-	auto loop(Expr&& exp) const -> decltype(mInd->ifor(exp))&&;
+	auto loop(Expr&& exp) const -> decltype(mInd->iforh(exp))&&;
     };
     
 }
@@ -547,7 +547,7 @@ namespace MultiArrayTools
     OperationMaster<T,OpClass,Ranges...>::
     OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
 		    std::shared_ptr<typename CRange::IndexType>& index) :
-	mSecond(second), mArrayRef(ma), mData(mArrayRef.data()),
+	mSecond(second), mArrayRef(ma), mData(mArrayRef.datav()),
 	mIndex(mkIndex(index)), mIInfo(*mIndex)
     {
 	auto blockIndex = seekBlockIndex( &mIInfo, second);
@@ -564,7 +564,7 @@ namespace MultiArrayTools
     OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
 		    std::shared_ptr<typename CRange::IndexType>& index,
 		    const IndexInfo* blockIndex) :
-	mSecond(second), mArrayRef(ma), mData(mArrayRef.data()),
+	mSecond(second), mArrayRef(ma), mData(mArrayRef.datav()),
 	mIndex(mkIndex(index)), mIInfo(*mIndex)
     {
 	std::intptr_t blockIndexNum = blockIndex->getPtrNum();
@@ -620,16 +620,18 @@ namespace MultiArrayTools
     }
 
     template <typename T, class OpClass, class... Ranges>
-    T& OperationMaster<T,OpClass,Ranges...>::get(size_t pos)
+    inline T& OperationMaster<T,OpClass,Ranges...>::get(size_t pos)
     {
-	VCHECK(pos);
+	//assert(pos < mIndex->max());
+	//if(pos >= mIndex->max()) { VCHECK(pos); VCHECK(mIndex->max()); assert(0); }
+	//VCHECK(pos);
 	return mData[pos];
     }
 
     template <typename T, class OpClass, class... Ranges>
-    const T& OperationMaster<T,OpClass,Ranges...>::get(size_t pos) const
+    inline const T& OperationMaster<T,OpClass,Ranges...>::get(size_t pos) const
     {
-	VCHECK(pos);
+	//VCHECK(pos);
 	return mData[pos];
     }
     
@@ -661,7 +663,7 @@ namespace MultiArrayTools
     ConstOperationRoot(const MultiArrayBase<T,Ranges...>& ma,
 		       const std::shared_ptr<typename Ranges::IndexType>&... indices) :
 	//OperationTemplate<T,ConstOperationRoot<T,Ranges...> >(this),
-	mArrayRef(ma), mData(mArrayRef.data()),
+	mArrayRef(ma), mData(mArrayRef.datav()),
 	mIndex( mkIndex(ma,indices...) ), mIInfo(*mIndex)
     {}
 
@@ -730,7 +732,7 @@ namespace MultiArrayTools
     OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
 		  const std::shared_ptr<typename Ranges::IndexType>&... indices) :
 	//OperationTemplate<T,OperationRoot<T,Ranges...> >(this),
-	mArrayRef(ma), mData(mArrayRef.data()),
+	mArrayRef(ma), mData(mArrayRef.datav()),
 	mIndex( mkIndex( ma, indices... ) ), mIInfo(*mIndex),
 	mBlockII(nullptr)
     {}
@@ -939,10 +941,10 @@ namespace MultiArrayTools
 
     template <typename T, class Op, class IndexType>
     template <class Expr>
-    auto Contraction<T,Op,IndexType>::loop(Expr&& exp) const -> decltype(mInd->ifor(exp))&&
+    auto Contraction<T,Op,IndexType>::loop(Expr&& exp) const -> decltype(mInd->iforh(exp))&&
     {
-	typedef decltype(mInd->ifor(exp)) LType;
+	typedef decltype(mInd->iforh(exp)) LType;
-	LType&& loop = mInd->ifor(exp);
+	LType&& loop = mInd->iforh(exp);
 	return std::forward<LType>( loop );
     }
 

src/ranges/container_range.h

@@ -92,7 +92,11 @@ namespace MultiArrayTools
 	template <class Exprs>
 	auto ifor(Exprs&& exs) const
 	    -> decltype(RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, exs));
-	
+
+	template <class Exprs>
+	auto iforh(Exprs&& exs) const
+	    -> decltype(RPackNum<sizeof...(Indices)-1>::mkForh(mIPack, exs));
+
     };
 
     
@@ -379,6 +383,14 @@ namespace MultiArrayTools
 	return RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, exs);
     }
 
+    template <class... Indices>
+    template <class Exprs>
+    auto ContainerIndex<Indices...>::iforh(Exprs&& exs) const
+	-> decltype(RPackNum<sizeof...(Indices)-1>::mkForh(mIPack, exs))
+    {
+	return RPackNum<sizeof...(Indices)-1>::mkForh(mIPack, exs);
+    }
+
     
     /*****************************
      *   ContainerRangeFactory   *

src/ranges/index_base.h

@@ -77,7 +77,11 @@ namespace MultiArrayTools
 	template <class Expr>
 	auto ifor(const Expr&& ex) const -> decltype(THIS().template ifor<Expr>(ex))
 	{ return THIS().template ifor<Expr>(ex); }
-	
+
+	template <class Expr>
+	auto iforh(const Expr&& ex) const -> decltype(THIS().template iforh<Expr>(ex))
+	{ return THIS().template iforh<Expr>(ex); }
+
     private:
 
 	friend I;

src/ranges/multi_range.h

@@ -108,6 +108,11 @@ namespace MultiArrayTools
 	template <class Exprs>
 	auto ifor(Exprs&& exs) const
 	    -> decltype(RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, exs));
+
+	template <class Exprs>
+	auto iforh(Exprs&& exs) const
+	    -> decltype(RPackNum<sizeof...(Indices)-1>::mkForh(mIPack, exs));
+
     };
 
     /*************************
@@ -422,7 +427,15 @@ namespace MultiArrayTools
     {
 	return RPackNum<sizeof...(Indices)-1>::mkFor(mIPack, exs);
     }
-    
+
+    template <class... Indices>
+    template <class Exprs>
+    auto MultiIndex<Indices...>::iforh(Exprs&& exs) const
+	-> decltype(RPackNum<sizeof...(Indices)-1>::mkForh(mIPack, exs))
+    {
+	return RPackNum<sizeof...(Indices)-1>::mkForh(mIPack, exs);
+    }
+
     /*************************
      *   MultiRangeFactory   *
      *************************/

src/ranges/rpack_num.h

@@ -235,13 +235,22 @@ namespace MultiArrayHelper
 
 	template <class IndexPack, class Exprs>
 	static auto mkFor(const IndexPack& ipack, Exprs&& exs)
-	    -> decltype(std::get<std::tuple_size<IndexPack>::value-N>(ipack)
+	    -> decltype(std::get<std::tuple_size<IndexPack>::value-N-1>(ipack)
 			->ifor( RPackNum<N-1>::mkFor(ipack, exs) ) )
 	{
-	    return std::get<std::tuple_size<IndexPack>::value-N>(ipack)
+	    return std::get<std::tuple_size<IndexPack>::value-N-1>(ipack)
 		->ifor( RPackNum<N-1>::mkFor(ipack, exs) );
 	}
 
+	template <class IndexPack, class Exprs>
+	static auto mkForh(const IndexPack& ipack, Exprs&& exs)
+	    -> decltype(std::get<std::tuple_size<IndexPack>::value-N-1>(ipack)
+			->iforh( RPackNum<N-1>::mkForh(ipack, exs) ) )
+	{
+	    return std::get<std::tuple_size<IndexPack>::value-N-1>(ipack)
+		->iforh( RPackNum<N-1>::mkForh(ipack, exs) );
+	}
+
     };
 
     
@@ -412,6 +421,15 @@ namespace MultiArrayHelper
 	{
 	    return std::get<std::tuple_size<IndexPack>::value-1>(ipack)->ifor(exs);
 	}
+
+	template <class IndexPack, class Exprs>
+	static auto mkForh(const IndexPack& ipack, Exprs&& exs)
+	    -> decltype(std::get<std::tuple_size<IndexPack>::value-1>(ipack)
+			->iforh(exs) )
+	{
+	    return std::get<std::tuple_size<IndexPack>::value-1>(ipack)->iforh(exs);
+	}
+
     };
 	
 

src/ranges/single_range.h

@@ -68,6 +68,11 @@ namespace MultiArrayTools
 	template <class Expr>
 	auto ifor(Expr&& ex) const
 	    -> For<SingleIndex<U,TYPE>,Expr>;
+
+	template <class Expr>
+	auto iforh(Expr&& ex) const
+	    -> For<SingleIndex<U,TYPE>,Expr,ForType::HIDDEN>;
+
     };
 
     template <typename U, SpaceType TYPE>
@@ -257,6 +262,15 @@ namespace MultiArrayTools
 	return For<SingleIndex<U,TYPE>,Expr>(this, std::forward<Expr>( ex ));
     }
 
+    template <typename U, SpaceType TYPE>
+    template <class Expr>
+    auto SingleIndex<U,TYPE>::iforh(Expr&& ex) const
+	-> For<SingleIndex<U,TYPE>,Expr,ForType::HIDDEN>
+    {
+	//static const size_t LAYER = typename Expr::LAYER; 
+	return For<SingleIndex<U,TYPE>,Expr,ForType::HIDDEN>(this, std::forward<Expr>( ex ));
+    }
+
     
     /********************
      *   SingleRange    *

src/tests/op_perf_test.cc

@@ -183,7 +183,7 @@ namespace {
 		for(size_t a = 0; a != 4; ++a){
 		    for(size_t b = 0; b != 4; ++b){
 			for(size_t c = 0; c != 4; ++c){
-			    const size_t sidx = d*4*4*4*4*4*4*4 + a*5*4*4*4*4 + b*5*4*4*4 + c*5*4 + p;
+			    const size_t sidx = d*4*4*4*4*4*4*4 + a*5*4*4*4*4*4 + b*5*4*4*4 + + c*5*4  + p;
 			    vres[tidx] += data[sidx];
 			}
 		    }

src/xfor/for_type.h

@@ -0,0 +1,18 @@
+
+#ifndef __for_type_h__
+#define __for_type_h__
+
+namespace MultiArrayHelper
+{
+
+    
+    enum class ForType {
+	DEFAULT = 0,
+	HIDDEN = 1	
+    };
+
+    
+} // end namespace MultiArrayHelper
+
+
+#endif

src/xfor/xfor.h

@@ -6,14 +6,33 @@
 #include <memory>
 #include <tuple>
 #include "xfor/for_utils.h"
+#include "xfor/for_type.h"
 
 namespace MultiArrayHelper
 {
 
     // 'HIDDEN FOR' CLASS for nested for loops in contractions a.s.o.
     // (NO COUNTING OF MASTER POSITION !!!!!)
+
+    template <ForType FT = ForType::DEFAULT>
+    struct PosForward
+    {
+	static inline size_t value(size_t last, size_t max, size_t pos)
+	{
+	    return last * max + pos;
+	}
+    };
+
+    template <>
+    struct PosForward<ForType::HIDDEN>
+    {
+	static inline size_t value(size_t last, size_t max, size_t pos)
+	{
+	    return last;
+	}
+    };
     
-    template <class IndexClass, class Expr>
+    template <class IndexClass, class Expr, ForType FT = ForType::DEFAULT>
     class For
     {
     private:
@@ -56,7 +75,7 @@ namespace MultiArrayHelper
     template <>
     size_t exceptMax<1>(size_t max) { return 1; }
 
-    
+   
 } // namespace MultiArrayHelper
 
 /* ========================= *
@@ -68,48 +87,60 @@ namespace MultiArrayHelper
 namespace MultiArrayHelper
 {
 
-    template <class IndexClass, class Expr>
+    template <class IndexClass, class Expr, ForType FT>
-    For<IndexClass,Expr>::For(const std::shared_ptr<IndexClass>& indPtr,
+    For<IndexClass,Expr,FT>::For(const std::shared_ptr<IndexClass>& indPtr,
 			      Expr&& expr) :
 	mIndPtr(indPtr.get()), mExpr(expr),
-	mExt(expr.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr.get() ))) { assert(mIndPtr != nullptr); }
+	mExt(expr.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr.get() )))
+    {
+	assert(mIndPtr != nullptr);
+	//VCHECK(mIndPtr->id());
+	//VCHECK(mIndPtr->max());
+    }
 
-    template <class IndexClass, class Expr>
+    template <class IndexClass, class Expr, ForType FT>
-    For<IndexClass,Expr>::For(const IndexClass* indPtr,
+    For<IndexClass,Expr,FT>::For(const IndexClass* indPtr,
 			      Expr&& expr) :
 	mIndPtr(indPtr), mExpr(std::forward<Expr>( expr )),
-	mExt(expr.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr ) )) { assert(mIndPtr != nullptr); }
+	mExt(expr.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr ) ))
+    {
+	assert(mIndPtr != nullptr);
+	//VCHECK(mIndPtr->id());
+	//VCHECK(mIndPtr->max());
+    }
     
-    template <class IndexClass, class Expr>
+    template <class IndexClass, class Expr, ForType FT>
-    inline void For<IndexClass,Expr>::operator()(size_t mlast,
+    inline void For<IndexClass,Expr,FT>::operator()(size_t mlast,
-						 const ETuple& last) const
+						    const ETuple& last) const
     {
 	auto& ind = *mIndPtr;
-	std::cout << mIndPtr << std::endl;
+	//std::cout << mIndPtr << std::endl;
 	const size_t max = ind.max(); // blocking
 	for(size_t pos = ind.pos(); pos != max; ++pos){
-	    const size_t mnpos = mlast * max + pos;
+	    //const size_t mnpos = mlast * max + pos;
+	    const size_t mnpos = PosForward<FT>::value(mlast, max, pos);
 	    const ETuple npos = std::move( XFPackNum<SIZE-1>::mkPos(pos, mExt, last) );
 	    mExpr(mnpos, npos);
 	}
     }
 
-    template <class IndexClass, class Expr>
+    template <class IndexClass, class Expr, ForType FT>
-    inline void For<IndexClass,Expr>::operator()(size_t mlast) const
+    inline void For<IndexClass,Expr,FT>::operator()(size_t mlast) const
     {
 	const ETuple last;
 	auto& ind = *mIndPtr;
-	std::cout << mIndPtr << std::endl;
+	//std::cout << mIndPtr << std::endl;
 	const size_t max = ind.max(); // blocking
 	for(size_t pos = ind.pos(); pos != max; ++pos){
-	    const size_t mnpos = mlast * max + pos;
+	    //const size_t mnpos = mlast * max + pos;
+	    const size_t mnpos = PosForward<FT>::value(mlast, max, pos);
 	    const ETuple npos = std::move( XFPackNum<SIZE-1>::mkPos(pos, mExt, last) );
 	    mExpr(mnpos, npos);
 	}
     }
     
-    template <class IndexClass, class Expr>
+    template <class IndexClass, class Expr, ForType FT>
-    typename For<IndexClass,Expr>::ETuple For<IndexClass,Expr>::rootSteps(std::intptr_t iPtrNum) const
+    typename For<IndexClass,Expr,FT>::ETuple For<IndexClass,Expr,FT>::rootSteps(std::intptr_t iPtrNum) const
     {
 	return mExpr.rootSteps(iPtrNum);
     }