diff --git a/src/include/ranges/rbase_def.h b/src/include/ranges/rbase_def.h
index fba3ff9..3cffe08 100644
--- a/src/include/ranges/rbase_def.h
+++ b/src/include/ranges/rbase_def.h
@@ -54,6 +54,18 @@ namespace MultiArrayTools
     // single_range.h
     template <typename U, SpaceType TYPE>
     class SingleIndex;
+
+    // subrange.h
+    template <class Index>
+    class SubIndex;
+
+    // subrange.h
+    template <class Range>
+    class SubRangeFactory;
+
+    // subrange.h
+    template <class Range>
+    class SubRange;
     
     // multi_range.h
     template <class... Ranges>
diff --git a/src/include/ranges/rpheader.h b/src/include/ranges/rpheader.h
index c3c4142..4dfd831 100644
--- a/src/include/ranges/rpheader.h
+++ b/src/include/ranges/rpheader.h
@@ -5,6 +5,7 @@
 #include "single_range.h"
 #include "multi_range.h"
 #include "container_range.h"
+#include "subrange.h"
 //#include "anonymous_range.h"
 
 //#endif
diff --git a/src/include/ranges/subrange.h b/src/include/ranges/subrange.h
index 8ded014..24cbb4c 100644
--- a/src/include/ranges/subrange.h
+++ b/src/include/ranges/subrange.h
@@ -24,12 +24,14 @@ namespace MultiArrayTools
     }
 
     template <class Index>
-    class SubIndex : public IndexInterface<SubIndex<Index>>
+    class SubIndex : public IndexInterface<SubIndex<Index>,typename Index::MetaType>
     {
-        typedef IndexInterface<SubIndex<Index>> IB;
-        typedef typename SubIndex<Index>::MetaType MetaType;
-        typedef SubRange<typename SubIndex<Index>::RangeType> RangeType;
-        typedef SubIndex IndexType;
+    public:
+        
+        typedef IndexInterface<SubIndex<Index>,typename Index::MetaType> IB;
+        typedef typename Index::MetaType MetaType;
+        typedef SubRange<typename Index::RangeType> RangeType;
+        typedef SubIndex IType;
 
         SubIndex(const std::shared_ptr<RangeType>& range);
 
@@ -37,7 +39,7 @@ namespace MultiArrayTools
         static constexpr size_t totalDim() { return typename Index::totalDim(); }
         static constexpr size_t sDim() { return typename Index::sDim(); }
 
-        static constexpr SpaceType STYPE = typename Index::STYPE;
+        static constexpr SpaceType STYPE = Index::STYPE;
 
         IndexType type() const;
         
@@ -45,16 +47,18 @@ namespace MultiArrayTools
 	SubIndex& operator++();
 	SubIndex& operator--();
 
+        SubIndex& operator()(const std::shared_ptr<Index>& ind); // set full index
+        
 	int pp(std::intptr_t idxPtrNum);
 	int mm(std::intptr_t idxPtrNum);
 
 	std::string stringMeta() const;
-	U meta() const;
-	const U* metaPtr() const;
-	SubIndex& at(const U& metaPos);
-	size_t posAt(const U& metaPos) const;
+	MetaType meta() const;
+	const MetaType* metaPtr() const;
+	SubIndex& at(const MetaType& metaPos);
+	size_t posAt(const MetaType& metaPos) const;
 
-	bool isMeta(const U& metaPos) const;
+	bool isMeta(const MetaType& metaPos) const;
 	
 	size_t dim(); // = 1
 	bool last();
@@ -95,7 +99,7 @@ namespace MultiArrayTools
         SubRangeFactory(const std::shared_ptr<Range>& fullRange,
                         const std::vector<size_t>& subset);
         std::shared_ptr<RangeBase> create();
-    }
+    };
     
     template <class Range>
     class SubRange : public RangeInterface<SubIndex<typename Range::IndexType>>
@@ -116,16 +120,17 @@ namespace MultiArrayTools
 	virtual std::string stringMeta(size_t pos) const final;
 	virtual std::vector<char> data() const final;
 	
-	bool isMeta(const U& metaPos) const;
+	bool isMeta(const MetaType& metaPos) const;
 	
-	const U& get(size_t pos) const;
-	size_t getMeta(const U& metaPos) const;
+	const MetaType& get(size_t pos) const;
+	size_t getMeta(const MetaType& metaPos) const;
 	
 	virtual IndexType begin() const final;
 	virtual IndexType end() const final;
 
         std::shared_ptr<Range> fullRange() const;
         const std::vector<size_t>& subset() const;
+        std::shared_ptr<SingleRange<MetaType,SpaceType::ANY>> outRange() const;
         
         friend SubRangeFactory<Range>;
 
@@ -156,8 +161,8 @@ namespace MultiArrayTools
 
     template <class Index>
     SubIndex<Index>::SubIndex(const std::shared_ptr<RangeType>& range) :
-        IndexInterface<SubIndex<Index>>(range, 0),
-        mExplicitRangePtr(std::dynamic_pointer_cast<RangeType>(IB::mRangePtr)),
+        IndexInterface<SubIndex<Index>,typename Index::MetaType>(range, 0),
+        mExplicitRangePtr(std::dynamic_pointer_cast<RangeType>(IB::mRangePtr))
     {
         mFullIndex = std::make_shared<Index>(mExplicitRangePtr->fullRange());
     }
@@ -170,7 +175,7 @@ namespace MultiArrayTools
     }
 
     template <class Index>
-    SubIndex& SubIndex<Index>::operator=(size_t pos)
+    SubIndex<Index>& SubIndex<Index>::operator=(size_t pos)
     {
         IB::mPos = pos;
         (*mFullIndex) = mExplicitRangePtr->subset()[IB::mPos];
@@ -178,7 +183,7 @@ namespace MultiArrayTools
     }
 
     template <class Index>
-    SubIndex& SubIndex<Index>::operator++()
+    SubIndex<Index>& SubIndex<Index>::operator++()
     {
         ++IB::mPos;
         (*mFullIndex) = mExplicitRangePtr->subset()[IB::mPos];
@@ -186,13 +191,21 @@ namespace MultiArrayTools
     }
 
     template <class Index>
-    SubIndex& SubIndex<Index>::operator--()
+    SubIndex<Index>& SubIndex<Index>::operator--()
     {
         --IB::mPos;
         (*mFullIndex) = mExplicitRangePtr->subset()[IB::mPos];
         return *this;
     }
 
+    template <class Index>
+    SubIndex<Index>& SubIndex<Index>::operator()(const std::shared_ptr<Index>& ind)
+    {
+        assert(mFullIndex->range() == ind->range());
+        mFullIndex = ind;
+        return *this;
+    }
+    
     template <class Index>
     int SubIndex<Index>::pp(std::intptr_t idxPtrNum)
     {
@@ -210,38 +223,39 @@ namespace MultiArrayTools
     template <class Index>
     std::string SubIndex<Index>::stringMeta() const
     {
-        return std::dynamic_pointer_cast<SingleRange<U,TYPE> const>( IB::mRangePtr )->stringMeta(IB::mPos);
+        return std::dynamic_pointer_cast<SingleRange<MetaType,STYPE> const>( IB::mRangePtr )->stringMeta(IB::mPos);
     }
 
     template <class Index>
-    U SubIndex<Index>::meta() const
+    typename SubIndex<Index>::MetaType SubIndex<Index>::meta() const
     {
+        MetaType* x = nullptr;
         return MetaPtrHandle<SubIndex<Index>::RangeType::HASMETACONT>::getMeta
-	    ( mMetaPtr, IB::mPos, mExplicitRangePtr );
+	    ( x, IB::mPos, mExplicitRangePtr );
     }
 
     template <class Index>
-    const U* SubIndex<Index>::metaPtr() const
+    const typename  SubIndex<Index>::MetaType* SubIndex<Index>::metaPtr() const
     {
         assert(0); // not sure where it is used
         return mFullIndex->metaPtr();
     }
 
     template <class Index>
-    SubIndex& SubIndex<Index>::at(const U& metaPos)
+    SubIndex<Index>& SubIndex<Index>::at(const MetaType& metaPos)
     {
         (*this) = mExplicitRangePtr->getMeta( metaPos );
         return *this;
     }
 
     template <class Index>
-    size_t SubIndex<Index>::posAt(const U& metaPos) const
+    size_t SubIndex<Index>::posAt(const MetaType& metaPos) const
     {
         return mExplicitRangePtr->getMeta( metaPos );
     }
 
     template <class Index>
-    bool SubIndex<Index>::isMeta(const U& metaPos) const
+    bool SubIndex<Index>::isMeta(const MetaType& metaPos) const
     {
         return mExplicitRangePtr->isMeta( metaPos );
     }
@@ -265,7 +279,7 @@ namespace MultiArrayTools
     }
 
     template <class Index>
-    std::shared_ptr<SubIndex<Index>::RangeType> SubIndex<Index>::range()
+    std::shared_ptr<typename SubIndex<Index>::RangeType> SubIndex<Index>::range()
     {
         return mExplicitRangePtr;
     }
@@ -303,7 +317,7 @@ namespace MultiArrayTools
         -> For<SubIndex<Index>,SubExpr<Index,Expr>>
     {
         return For<SubIndex<Index>,SubExpr<Index,Expr>>
-            (this, /**/, SubExpr( mFullIndex, /**/ ) );
+            (this, 1, SubExpr<Index,Expr>( mFullIndex, &mExplicitRangePtr->subset(), step, ex ) );
     }
 
     /************************
@@ -331,7 +345,7 @@ namespace MultiArrayTools
     template <class Range>
     SubRange<Range>::SubRange(const std::shared_ptr<Range>& fullRange,
                               const std::vector<size_t>& subset) :
-        RangeInterface<SubRange<Range>>(),
+        RangeInterface<SubIndex<typename Range::IndexType>>(),
         mFullRange(fullRange), mSubSet(subset) {}
 
     template <class Range>
@@ -362,7 +376,8 @@ namespace MultiArrayTools
         h.multiple = 0;
         return h;
     }
-        
+
+    template <class Range>
     std::string SubRange<Range>::stringMeta(size_t pos) const
     {
         return xToString(get(pos));
@@ -386,7 +401,7 @@ namespace MultiArrayTools
     }
 
     template <class Range>
-    bool SubRange<Range>::isMeta(const U& metaPos) const
+    bool SubRange<Range>::isMeta(const MetaType& metaPos) const
     {
         for(size_t i = 0; i != size(); ++i){
             if(get(i) == metaPos){
@@ -397,13 +412,13 @@ namespace MultiArrayTools
     }
 
     template <class Range>
-    const SubRange<Range>::MetaType& SubRange<Range>::get(size_t pos) const
+    const typename SubRange<Range>::MetaType& SubRange<Range>::get(size_t pos) const
     {
         return mFullRange->get( mSubSet[pos] );
     }
 
     template <class Range>
-    size_t SubRange<Range>::getMeta(const U& metaPos) const
+    size_t SubRange<Range>::getMeta(const MetaType& metaPos) const
     {
         for(size_t i = 0; i != size(); ++i){
             if(get(i) == metaPos){
@@ -414,7 +429,7 @@ namespace MultiArrayTools
     }
 
     template <class Range>
-    SubRange<Range>::IndexType SubRange<Range>::begin() const
+    typename SubRange<Range>::IndexType SubRange<Range>::begin() const
     {
         SubRange<Range>::IndexType i( std::dynamic_pointer_cast<SubRange<Range>>
                                       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
@@ -423,7 +438,7 @@ namespace MultiArrayTools
     }
 
     template <class Range>
-    SubRange<Range>::IndexType SubRange<Range>::end() const
+    typename SubRange<Range>::IndexType SubRange<Range>::end() const
     {
         SubRange<Range>::IndexType i( std::dynamic_pointer_cast<SubRange<Range>>
                                       ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
@@ -443,6 +458,17 @@ namespace MultiArrayTools
         return mSubSet;
     }
 
+    template <class Range>
+    std::shared_ptr<SingleRange<typename SubRange<Range>::MetaType,SpaceType::ANY>> SubRange<Range>::outRange() const
+    {
+        std::vector<MetaType> ometa(mSubSet.size());
+        size_t i = 0;
+        for(auto& x: mSubSet){
+            ometa[i++] = mFullRange->get(x);
+        }
+        SingleRangeFactory<MetaType,SpaceType::ANY> srf(ometa);
+        return std::dynamic_pointer_cast<SingleRange<MetaType,SpaceType::ANY>>( srf.create() );
+    }
 }
 
 #endif
diff --git a/src/include/xfor/xfor.h b/src/include/xfor/xfor.h
index 4b11d76..c8a8b4f 100644
--- a/src/include/xfor/xfor.h
+++ b/src/include/xfor/xfor.h
@@ -132,7 +132,57 @@ namespace MultiArrayHelper
         auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
 	auto extension() const -> ExtType;
     };
-    
+
+    template <class IndexClass, class Expr>
+    class SubExpr : public ExpressionBase
+    {
+    private:
+	SubExpr() = default;
+
+	const IndexClass* mIndPtr;
+	size_t mSPos;
+	size_t mMax;
+
+        Expr mExpr;
+	typedef decltype(mExpr.rootSteps()) ExtType;
+	ExtType mExt;
+
+        const std::vector<size_t>* mSubSet;
+        size_t mStep;
+        
+        mutable ExtType mRootSteps;
+
+    public:
+	typedef ExpressionBase EB;
+	
+	static constexpr size_t LAYER = Expr::LAYER + 1;
+	static constexpr size_t SIZE = Expr::SIZE;
+
+	SubExpr(const SubExpr& in) = default;
+	SubExpr& operator=(const SubExpr& in) = default;
+	SubExpr(SubExpr&& in) = default;
+	SubExpr& operator=(SubExpr&& in) = default;
+	
+	SubExpr(const std::shared_ptr<IndexClass>& indPtr,
+                const std::vector<size_t>* subset,
+                size_t step, Expr expr);
+
+	SubExpr(const IndexClass* indPtr,
+                const std::vector<size_t>* subset,
+                size_t step, Expr expr);
+
+	
+	inline void operator()(size_t mlast, DExt last) const override final;
+	inline void operator()(size_t mlast, ExtType last) const;
+	inline void operator()(size_t mlast = 0) const override final;
+
+        DExt dRootSteps(std::intptr_t iPtrNum = 0) const override final;
+        DExt dExtension() const override final;
+
+        auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
+	auto extension() const -> ExtType;
+    };
+
     template <class IndexClass, class Expr, ForType FT>
     class For : public ExpressionBase
     {
@@ -270,7 +320,7 @@ namespace MultiArrayHelper
     template <class IndexClass, class Expr, ForType FT>
     For<IndexClass,Expr,FT>::For(const std::shared_ptr<IndexClass>& indPtr,
 				 size_t step, Expr expr) :
-	mIndPtr(indPtr), mSPos(mIndPtr->pos()), mMax(mIndPtr->max()), mStep(step),
+	mIndPtr(indPtr.get()), mSPos(mIndPtr->pos()), mMax(mIndPtr->max()), mStep(step),
         mExpr(expr), mExt(mExpr.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr )))
     {
 	assert(mIndPtr != nullptr);
@@ -432,6 +482,96 @@ namespace MultiArrayHelper
 						    sizeof(ExtType)/sizeof(size_t));
     }
     
+    /****************
+     *   SubExpr    *
+     ****************/
+    
+    template <class IndexClass, class Expr>
+    SubExpr<IndexClass,Expr>::SubExpr(const std::shared_ptr<IndexClass>& indPtr,
+                                      const std::vector<size_t>* subset,
+                                      size_t step,
+                                      Expr expr) :
+	mIndPtr(indPtr.get()), mSPos(mIndPtr->pos()), mMax(mIndPtr->max()),
+        mExpr(expr), mExt(mExpr.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr ))),
+        mSubSet(subset), mStep(step)
+    {
+	assert(mIndPtr != nullptr);
+    }
+
+    template <class IndexClass, class Expr>
+    SubExpr<IndexClass,Expr>::SubExpr(const IndexClass* indPtr,
+                                      const std::vector<size_t>* subset,
+                                      size_t step,
+                                      Expr expr) :
+	mIndPtr(indPtr), mSPos(mIndPtr->pos()), mMax(mIndPtr->max()),
+        mExpr(expr), mExt(mExpr.rootSteps( reinterpret_cast<std::intptr_t>( mIndPtr ))),
+        mSubSet(subset), mStep(step)
+    {
+	assert(mIndPtr != nullptr);
+    }
+
+    template <class IndexClass, class Expr>
+    inline void SubExpr<IndexClass,Expr>::operator()(size_t mlast, DExt last) const
+    {
+        operator()(mlast, *reinterpret_cast<ExtType const*>(last.first));
+    }
+    
+    template <class IndexClass, class Expr>
+    inline void SubExpr<IndexClass,Expr>::operator()(size_t mlast,
+                                                     ExtType last) const
+    {
+        // INCLUDE FOR LOOP HERE AGIAN !!!!
+        const size_t pos = (*mSubSet)[mlast];
+        const size_t mnpos = mlast * mStep;
+        VCHECK(mlast);
+        VCHECK(pos);
+        VCHECK(mnpos);
+        const ExtType npos = last + mExt*pos;
+        VCHECK(npos.val());
+        VCHECK(npos.next().val());
+        mExpr(mnpos, npos);
+    }
+
+    template <class IndexClass, class Expr>
+    inline void SubExpr<IndexClass,Expr>::operator()(size_t mlast) const
+    {
+	const ExtType last;
+        const size_t pos = (*mSubSet)[mlast];
+        const size_t mnpos = mlast * mStep;
+        const ExtType npos = last + mExt*pos;
+        mExpr(mnpos, npos);
+    }
+    
+    
+    template <class IndexClass, class Expr>
+    auto SubExpr<IndexClass,Expr>::rootSteps(std::intptr_t iPtrNum) const
+	-> ExtType
+    {
+	return mExpr.rootSteps(iPtrNum);
+    }
+
+    template <class IndexClass, class Expr>
+    auto SubExpr<IndexClass,Expr>::extension() const
+	-> ExtType
+    {
+	return mExt;
+    }
+
+    template <class IndexClass, class Expr>
+    DExt SubExpr<IndexClass,Expr>::dRootSteps(std::intptr_t iPtrNum) const
+    {
+        mRootSteps = rootSteps(iPtrNum);
+        return std::make_pair<size_t const*,size_t>(reinterpret_cast<size_t const*>(&mRootSteps),
+						    sizeof(ExtType)/sizeof(size_t));
+    }
+
+    template <class IndexClass, class Expr>
+    DExt SubExpr<IndexClass,Expr>::dExtension() const
+    {
+        return std::make_pair<size_t const*,size_t>(reinterpret_cast<size_t const*>(&mExt),
+						    sizeof(ExtType)/sizeof(size_t));
+    }
+
     /***************************
      *   DynamicExpression   *
      ***************************/
diff --git a/src/tests/op_unit_test.cc b/src/tests/op_unit_test.cc
index 5917d02..a2b002d 100644
--- a/src/tests/op_unit_test.cc
+++ b/src/tests/op_unit_test.cc
@@ -159,6 +159,38 @@ namespace {
 				   -14.364, -1.868, -25.703, 13.836, 23.563, 41.339 };
     };
 
+    class OpTest_Sub : public ::testing::Test
+    {
+    protected:
+
+	typedef SingleRangeFactory<char,SpaceType::ANY> SRF;
+	typedef SRF::oType SRange;
+
+        OpTest_Sub()
+	{
+	    swapFactory<SRF>(rfbptr, {'x', 'l'} );
+	    sr1ptr = std::dynamic_pointer_cast<SRange>( rfbptr->create() );
+
+	    swapFactory<SRF>(rfbptr, {'1', '2', '3'} );
+	    sr2ptr = std::dynamic_pointer_cast<SRange>( rfbptr->create() );
+
+	    swapFactory<SRF>(rfbptr, {'a', 'b'} );
+	    sr3ptr = std::dynamic_pointer_cast<SRange>( rfbptr->create() );
+
+        }
+
+	std::shared_ptr<RangeFactoryBase> rfbptr;
+        std::shared_ptr<SRange> sr1ptr;
+	std::shared_ptr<SRange> sr2ptr;
+	std::shared_ptr<SRange> sr3ptr;
+
+        std::vector<double> v1 = { 2.917, 9.436, 0.373, 0.353, 4.005, 1.070,
+                                   -14.364, -1.868, -25.703, 13.836, 23.563, 41.339 };
+
+        std::vector<double> v2 = { 0.353, 4.005, 1.070, 2.310, 9.243, 2.911 };
+    };
+
+    
     class MapTest : public ::testing::Test
     {
     protected:
@@ -506,6 +538,59 @@ namespace {
 	EXPECT_EQ( fabs( res.at('g') - (7.192+5.063) ) < 0.0001, true );
     }
 
+    TEST_F(OpTest_Sub, Exec)
+    {
+        MultiArray<double,SRange,SRange,SRange> ma1(sr1ptr, sr2ptr, sr3ptr, v1);
+        MultiArray<double,SRange,SRange> ma2(sr3ptr, sr2ptr, v2);
+
+        SubRangeFactory<SRange> subf(sr2ptr, std::vector<size_t>({0,2}));
+        auto subptr = MAT::createExplicit(subf);
+        
+        MultiArray<double,SubRange<SRange>,SRange,SRange> res(subptr,sr3ptr,sr1ptr,0.);
+
+        auto i1 = MAT::getIndex( sr1ptr );
+        auto i2 = MAT::getIndex( sr2ptr );
+        auto i3 = MAT::getIndex( sr3ptr );
+        auto si = MAT::getIndex( subptr );
+        (*si)(i2);
+
+        res(si,i3,i1) = ma2(i3,i2) - ma1(i1,i2,i3);
+
+        EXPECT_EQ( res.size(), 8 );
+        EXPECT_EQ( res.vdata().size(), 8 );
+        EXPECT_EQ( MAT::rptr<0>( res )->size(), 2 );
+
+        EXPECT_EQ( MAT::rptr<0>( res )->isMeta('1'), true );
+        EXPECT_EQ( MAT::rptr<0>( res )->isMeta('3'), true );
+        
+        EXPECT_EQ( xround( res.at(mkt('1','a','x')) ), xround(ma2.at(mkt('a','1')) - ma1.at(mkt('x','1','a'))) );
+        EXPECT_EQ( xround( res.at(mkt('1','a','l')) ), xround(ma2.at(mkt('a','1')) - ma1.at(mkt('l','1','a'))) );
+
+        EXPECT_EQ( xround( res.at(mkt('1','b','x')) ), xround(ma2.at(mkt('b','1')) - ma1.at(mkt('x','1','b'))) );
+        EXPECT_EQ( xround( res.at(mkt('1','b','l')) ), xround(ma2.at(mkt('b','1')) - ma1.at(mkt('l','1','b'))) );
+
+        EXPECT_EQ( xround( res.at(mkt('3','a','x')) ), xround(ma2.at(mkt('a','3')) - ma1.at(mkt('x','3','a'))) );
+        EXPECT_EQ( xround( res.at(mkt('3','a','l')) ), xround(ma2.at(mkt('a','3')) - ma1.at(mkt('l','3','a'))) );
+
+        EXPECT_EQ( xround( res.at(mkt('3','b','x')) ), xround(ma2.at(mkt('b','3')) - ma1.at(mkt('x','3','b'))) );
+        EXPECT_EQ( xround( res.at(mkt('3','b','l')) ), xround(ma2.at(mkt('b','3')) - ma1.at(mkt('l','3','b'))) );
+
+        auto res2 = res.format( subptr->outRange(), sr3ptr, sr1ptr );
+
+        EXPECT_EQ( xround( res2.at(mkt('1','a','x')) ), xround(ma2.at(mkt('a','1')) - ma1.at(mkt('x','1','a'))) );
+        EXPECT_EQ( xround( res2.at(mkt('1','a','l')) ), xround(ma2.at(mkt('a','1')) - ma1.at(mkt('l','1','a'))) );
+
+        EXPECT_EQ( xround( res2.at(mkt('1','b','x')) ), xround(ma2.at(mkt('b','1')) - ma1.at(mkt('x','1','b'))) );
+        EXPECT_EQ( xround( res2.at(mkt('1','b','l')) ), xround(ma2.at(mkt('b','1')) - ma1.at(mkt('l','1','b'))) );
+
+        EXPECT_EQ( xround( res2.at(mkt('3','a','x')) ), xround(ma2.at(mkt('a','3')) - ma1.at(mkt('x','3','a'))) );
+        EXPECT_EQ( xround( res2.at(mkt('3','a','l')) ), xround(ma2.at(mkt('a','3')) - ma1.at(mkt('l','3','a'))) );
+
+        EXPECT_EQ( xround( res2.at(mkt('3','b','x')) ), xround(ma2.at(mkt('b','3')) - ma1.at(mkt('x','3','b'))) );
+        EXPECT_EQ( xround( res2.at(mkt('3','b','l')) ), xround(ma2.at(mkt('b','3')) - ma1.at(mkt('l','3','b'))) );
+        
+    }
+    
     TEST_F(OpTest_MDim, ExecOp1)
     {
 	MultiArray<double,SRange,SRange> res(sr2ptr,sr4ptr);