diff --git a/src/base_def.h b/src/base_def.h
index 321b288..6409a6e 100644
--- a/src/base_def.h
+++ b/src/base_def.h
@@ -113,6 +113,10 @@ namespace MultiArrayTools
     template <typename T, class Operation, class... MAOps>
     class MultiArrayOperation;
 
+    // multi_array_operation.h
+    template <typename T, class Operation, class Range, class MAOp>
+    class MultiArrayContraction;
+    
     // slice.h
     template <typename T, class Range, class MARange>
     class Slice;
diff --git a/src/multi_array_operation.cc b/src/multi_array_operation.cc
index 4f296db..113605a 100644
--- a/src/multi_array_operation.cc
+++ b/src/multi_array_operation.cc
@@ -34,7 +34,7 @@ namespace MultiArrayTools
     template <typename T, class Range>
     void MultiArrayOperationRoot<T,Range>::performAssignment(const MultiArrayOperationBase<T>& in)
     {
-#error "WRITE MAOR INTRINSIC CONTRACT FUNCTION"
+	//#error "WRITE MAOR INTRINSIC CONTRACT FUNCTION"
 	//CHECK;
 	in.linkIndicesTo(MAOB::mIibPtr);
 	//CHECK;
@@ -42,6 +42,7 @@ namespace MultiArrayTools
 	//CHECK;
 	const size_t endPos = mArrayRef.endIndex().pos();
 	std::cout << "assignment: " << endPos << " elements" << std::endl;
+	// assignment loop
 	for(iref = mArrayRef.beginIndex().pos(); iref != mArrayRef.endIndex(); ++iref){
 	    std::cout << iref.pos() << '\r' << std::flush;
 	    get() = in.get();
@@ -50,6 +51,7 @@ namespace MultiArrayTools
 	MAOB::mIibPtr->freeLinked();
     }
 
+    /*
     template <typename T, class Range>
     template <class RangeX>
     MultiArrayOperationRoot<T,Range>&
@@ -70,7 +72,7 @@ namespace MultiArrayTools
 	sl.setConst(in.mArrayRef, name(), dynamic_cast<const typename RangeX::IndexType&>( in.index() ), in.name());
 	return *this;
     }
-    
+    */
     
     // CONST SLICE !!!!!
     
@@ -102,10 +104,9 @@ namespace MultiArrayTools
     MultiArrayOperationRoot<T,Range>::operator=(MultiArrayOperationRoot<T,Range>& in)
     {
 	//CHECK;
-	maketurnSlice(in);
-	if(mArrayRef.isSlice() and not mArrayRef.isInit()){
-	    return makeSlice(in);
-	}
+	//	if(mArrayRef.isSlice() and not mArrayRef.isInit()){
+	//    return makeSlice(in);
+	//}
 	performAssignment(in);
 	freeIndex();
 	return *this;
@@ -117,9 +118,9 @@ namespace MultiArrayTools
     MultiArrayOperationRoot<T,Range>::operator=(MultiArrayOperationRoot<T,Range2>& in)
     {
 	//CHECK;
-	if(mArrayRef.isSlice() and not mArrayRef.isInit()){
-	    return makeSlice(in);
-	}
+	//if(mArrayRef.isSlice() and not mArrayRef.isInit()){
+	//    return makeSlice(in);
+	//}
 	performAssignment(in);
 	freeIndex();
 	return *this;
@@ -132,10 +133,10 @@ namespace MultiArrayTools
     MultiArrayOperationRoot<T,Range>::operator=(const MultiArrayOperationRoot<T,Range2>& in)
     {
 	//CHECK;
-	if(mArrayRef.isSlice() and not mArrayRef.isInit()){
+	//if(mArrayRef.isSlice() and not mArrayRef.isInit()){
 	    //CHECK;
-	    return makeConstSlice(in);
-	}
+	//    return makeConstSlice(in);
+	//}
 	performAssignment(in);
 	freeIndex();
 	return *this;
@@ -147,11 +148,11 @@ namespace MultiArrayTools
     MultiArrayOperationRoot<T,Range>::operator=(const MultiArrayOperation<T,Operation,MAOps...>& in)
     {
 	//CHECK;
-	if(mArrayRef.isSlice() and not mArrayRef.isInit()){
+	//if(mArrayRef.isSlice() and not mArrayRef.isInit()){
 	    // NO SLICE CREATION !!! (total array not initialized!!)
 	    // throw !
-	    assert(0);
-	}
+	//    assert(0);
+	//}
    	performAssignment(in);
 	freeIndex();
 	return *this;
@@ -166,6 +167,31 @@ namespace MultiArrayTools
 	return MultiArrayOperation<T,Operation,MultiArrayOperationRoot<T,Range>, MAOps...>(op, *this, secs...);
     }
 
+    template <typename T, class Range>
+    template <class Range2, class ContractOperation>
+    MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperationRoot<T,Range> >
+    MultiArrayOperationRoot<T,Range>::contract(const ContractOperation& cop,
+					       const std::string& indexName) const
+    {
+	typename Range2::IndexType* ind = dynamic_cast<typename Range2::IndexType*>( mIndex.getLinked(indexName) );
+	return MultiArrayContraction<T,ContractOperation,Range2,
+				     MultiArrayOperationRoot<T,Range> >(cop, *this, *ind);
+    }
+
+    template <typename T, class Range>
+    template <class Range2, class ContractOperation>
+    MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperationRoot<T,Range> >
+    MultiArrayOperationRoot<T,Range>::contract(const ContractOperation& cop,
+					       const std::string& indexName,
+					       const typename Range2::IndexType& begin,
+					       const typename Range2::IndexType& end) const
+    {
+	typename Range2::IndexType* ind = dynamic_cast<typename Range2::IndexType*>( mIndex.getLinked(indexName) );
+	return MultiArrayContraction<T,ContractOperation,Range2,
+				     MultiArrayOperationRoot<T,Range> >(cop, *this, *ind, begin, end);
+    }
+
+    
     template <typename T, class Range>
     template <class MAOp>
     auto MultiArrayOperationRoot<T,Range>::operator+(const MAOp& sec)
@@ -411,6 +437,31 @@ namespace MultiArrayTools
 	return MultiArrayOperation<T,Operation,ConstMultiArrayOperationRoot<T,Range>, MAOps...>(op, *this, secs...);
     }
 
+    template <typename T, class Range>
+    template <class Range2, class ContractOperation>
+    MultiArrayContraction<T,ContractOperation,Range2,ConstMultiArrayOperationRoot<T,Range> >
+    ConstMultiArrayOperationRoot<T,Range>::contract(const ContractOperation& cop,
+						    const std::string& indexName) const
+    {
+	typename Range2::IndexType* ind = dynamic_cast<typename Range2::IndexType*>( mIndex.getLinked(indexName) );
+	return MultiArrayContraction<T,ContractOperation,Range2,
+				     ConstMultiArrayOperationRoot<T,Range> >(cop, *this, *ind);
+    }
+
+    template <typename T, class Range>
+    template <class Range2, class ContractOperation>
+    MultiArrayContraction<T,ContractOperation,Range2,ConstMultiArrayOperationRoot<T,Range> >
+    ConstMultiArrayOperationRoot<T,Range>::contract(const ContractOperation& cop,
+						    const std::string& indexName,
+						    const typename Range2::IndexType& begin,
+						    const typename Range2::IndexType& end) const
+    {
+	typename Range2::IndexType* ind = dynamic_cast<typename Range2::IndexType*>( mIndex.getLinked(indexName) );
+	return MultiArrayContraction<T,ContractOperation,Range2,
+				     ConstMultiArrayOperationRoot<T,Range> >(cop, *this, *ind, begin, end);
+    }
+
+
     template <typename T, class Range>
     template <class MAOp>
     auto ConstMultiArrayOperationRoot<T,Range>::operator+(const MAOp& sec) const
@@ -639,20 +690,55 @@ namespace MultiArrayTools
 	TupleIndicesLinker<sizeof...(MAOps)-1>::linkTupleIndicesTo(mArgs, target);
     }
 
-    /*
-    template <typename T, class Operation, class... MAOps>    
-    T& MultiArrayOperation<T,Operation,MAOps...>::get()
-    {
-	mVal = OperationCall<sizeof...(MAOps)-1>::
-	    template callOperation(mOp, mArgs);
-	return mVal;
-	}*/
-
-    template <typename T, class Operation, class... MAOps>    
+    
+    template <typename T, class Operation, class... MAOps>
     const T& MultiArrayOperation<T,Operation,MAOps...>::get() const
     {
 	mVal = OperationCall<sizeof...(MAOps)-1>::
 	    template callOperation(mOp, mArgs);
 	return mVal;
     }
+
+    /*******************************
+     *   MultiArrayContraction     *
+     *******************************/
+
+    template <typename T, class ContractOperation, class Range, class MAOp>
+    MultiArrayContraction<T,ContractOperation,Range,MAOp>::
+    MultiArrayContraction(const ContractOperation& cop, const MAOp& mao,
+			  const typename Range::IndexType& runIndex) :
+	MultiArrayOperation<T,ContractOperation,MAOp>(cop, mao),
+	mBeginIndex(runIndex), mEndIndex(runIndex),
+	mRunIndex(runIndex)
+    {
+	mBeginIndex.setPos(0);
+	mEndIndex.setPos(mRunIndex.max());
+	MAO::linkIndicesTo(&mRunIndex);
+    }
+
+    template <typename T, class ContractOperation, class Range, class MAOp>
+    MultiArrayContraction<T,ContractOperation,Range,MAOp>::
+    MultiArrayContraction(const ContractOperation& cop, const MAOp& mao,
+			  const typename Range::IndexType& runIndex,
+			  const typename Range::IndexType& beginIndex,
+			  const typename Range::IndexType& endIndex) :
+	MultiArrayOperation<T,ContractOperation,MAOp>(cop, mao),
+	mBeginIndex(beginIndex), mEndIndex(endIndex),
+	mRunIndex(runIndex)
+    {
+	MAO::linkIndicesTo(&mRunIndex);
+    }
+
+    
+    // for the moment simplest case only:
+    template <typename T, class ContractOperation, class Range, class MAOp>
+    const T& MultiArrayContraction<T,ContractOperation,Range,MAOp>::get() const
+    {
+	MAO::mOp.reset();
+	for(mRunIndex.copyPos( mBeginIndex ); mRunIndex.pos() != mEndIndex.pos(); ++mRunIndex){
+	    MAO::mOp(std::get<0>(MAO::mArgs).get() );
+	}
+	MAO::mOp.endOp(MAO::mVal);
+	return MAO::mOp();
+    }
 }
diff --git a/src/multi_array_operation.h b/src/multi_array_operation.h
index 64c7ecc..def9563 100644
--- a/src/multi_array_operation.h
+++ b/src/multi_array_operation.h
@@ -73,6 +73,17 @@ namespace MultiArrayTools
 	template <class Operation, class... MAOps>
 	MultiArrayOperation<T,Operation,MultiArrayOperationRoot<T,Range>, MAOps...>
 	operator()(const Operation& op, const MAOps&... secs) const;
+
+	template < class Range2, class ContractOperation>
+	MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperationRoot<T,Range> >
+	contract(const ContractOperation& cop, const std::string& indexName) const;
+	
+	template <class Range2, class ContractOperation>
+	MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperationRoot<T,Range> >
+	contract(const ContractOperation& cop, const std::string& indexName,
+		 const typename Range2::IndexType& begin,
+		 const typename Range2::IndexType& end) const;
+
 	
 	template <class MAOp>
 	auto operator+(const MAOp& sec) -> decltype(operator()(std::plus<T>(), sec));
@@ -131,13 +142,13 @@ namespace MultiArrayTools
 
 	void performAssignment(const MultiArrayOperationBase<T>& in);
 
+	/*
 	template <class RangeX>
 	MultiArrayOperationRoot& makeSlice(MultiArrayOperationRoot<T,RangeX>& in);
-
 	
 	template <class RangeX>
 	const MultiArrayOperationRoot& makeConstSlice(const MultiArrayOperationRoot<T,RangeX>& in);
-	
+	*/
 	MutableMultiArrayBase<T,Range>& mArrayRef;
 	mutable IndexType mIndex;
 	Name mNm;
@@ -155,23 +166,20 @@ namespace MultiArrayTools
 	ConstMultiArrayOperationRoot(const MultiArrayBase<T,Range>& ma, const Name& nm);
 	ConstMultiArrayOperationRoot(const MultiArrayOperationRoot<T,Range>& in);
 
-	/*
-	const ConstMultiArrayOperationRoot& operator=(const ConstMultiArrayOperationRoot& in);
-
-	template <class Range2>
-	const ConstMultiArrayOperationRoot& operator=(const ConstMultiArrayOperationRoot<T,Range2>& in);
-
-	template <class Range2>
-	const ConstMultiArrayOperationRoot& operator=(const MultiArrayOperationRoot<T,Range2>& in);
-	*/
-
-	//template <class Operation, class... MAOps>
-	//MultiArrayOperation<T,Operation,MultiArrayOperationRoot<T,Range>, MAOps...>
-	//operator()(Operation& op, const MAOps&... secs) const;
-
 	template <class Operation, class... MAOps>
 	MultiArrayOperation<T,Operation,ConstMultiArrayOperationRoot<T,Range>, MAOps...>
 	operator()(const Operation& op, const MAOps&... secs) const;
+
+	template <class Range2, class ContractOperation>
+	MultiArrayContraction<T,ContractOperation,Range2,ConstMultiArrayOperationRoot<T,Range> >
+	contract(const ContractOperation& cop, const std::string& indexName) const;
+
+	
+	template <class Range2, class ContractOperation>
+	MultiArrayContraction<T,ContractOperation,Range2,ConstMultiArrayOperationRoot<T,Range> >
+	contract(const ContractOperation& cop, const std::string& indexName,
+		 const typename Range2::IndexType& begin,
+		 const typename Range2::IndexType& end) const;
 	
 	template <class MAOp>
 	auto operator+(const MAOp& sec) const -> decltype(operator()(std::plus<T>(), sec));
@@ -262,7 +270,6 @@ namespace MultiArrayTools
 	
 	virtual void linkIndicesTo(IndefinitIndexBase* target) const override;
 
-	//virtual T& get() override;
 	virtual const T& get() const override;
 
     protected:
@@ -272,8 +279,30 @@ namespace MultiArrayTools
 	OBT mArgs; // include first arg also here !!!
     };
 
+    template <typename T, class ContractOperation, class Range, class MAOp>
+    class MultiArrayContraction : public MultiArrayOperation<T,ContractOperation,MAOp>
+    {
+    public:
+	typedef MultiArrayOperationBase<T> MAOB;
+	typedef MultiArrayOperation<T,ContractOperation,MAOp> MAO;
 
+	MultiArrayContraction(const ContractOperation& cop, const MAOp& mao,
+			      const typename Range::IndexType& runIndex);
 
+	
+	MultiArrayContraction(const ContractOperation& cop, const MAOp& mao,
+			      const typename Range::IndexType& runIndex,
+			      const typename Range::IndexType& beginIndex,
+			      const typename Range::IndexType& endIndex);
+	
+	virtual const T& get() const override;
+
+    protected:
+	typename Range::IndexType mBeginIndex;
+	typename Range::IndexType mEndIndex;
+        mutable typename Range::IndexType mRunIndex;
+    };
+    
 }
 
 #include "multi_array_operation.cc"
diff --git a/src/unit_test.cc b/src/unit_test.cc
index 769e71a..f1c7298 100644
--- a/src/unit_test.cc
+++ b/src/unit_test.cc
@@ -9,6 +9,33 @@
 namespace MAT = MultiArrayTools;
 
 namespace {
+
+    template <typename T>
+    struct sum
+    {
+    public:
+	sum() = default;
+
+	T& operator()() const
+	{
+	    return res;
+	}
+	
+	T& operator()(const T& a) const
+	{
+	    return res += a;
+	}
+	
+	void endOp(T& res) const {}
+
+	void reset() const
+	{
+	    res = static_cast<T>(0);
+	}
+	
+    private:
+	mutable T res = static_cast<T>(0);
+    };
     
     class OneDimTest : public ::testing::Test
     {
@@ -135,6 +162,32 @@ namespace {
 	MultiArray3dAny ma;
 	//Slice2d3dAny sl;
     };
+
+    class ContractionTest : public ::testing::Test
+    {
+    protected:
+	typedef MAT::SingleRange<char,MAT::RangeType::ANY> Range1dAny;
+	typedef MAT::MultiRange<Range1dAny,Range1dAny> Range2dAny;
+	typedef MAT::MultiRange<Range1dAny,Range1dAny,Range1dAny> Range3dAny;
+	typedef MAT::MultiArray<int,Range3dAny> MultiArray3dAny;
+	typedef MAT::Slice<int,Range2dAny,Range3dAny> Slice2d3dAny;
+	typedef MAT::MultiArray<int,Range2dAny> MultiArray2dAny;
+	
+	ContractionTest() : r1({'a','b','c'}), r2({'a','b','c','d'}), r3({'a','b'}),
+			    ra(r1,r3),
+			    rb(r1,r2),
+			    r3d(r1,r2,r3),
+			    ma(r3d, {-5,6,2,1,9,54,27,-7,-13,32,90,-67,
+					-10,16,-2,101,39,-64,81,-22,14,34,95,-62}) {}
+	
+	Range1dAny r1;
+	Range1dAny r2;
+	Range1dAny r3;
+	Range2dAny ra;
+	Range2dAny rb;
+	Range3dAny r3d;
+	MultiArray3dAny ma;
+    };
     
     TEST_F(OneDimTest, CorrectExtensions)
     {
@@ -454,6 +507,26 @@ namespace {
 	EXPECT_EQ(sl[j(j1 = 2, j2 = 0)], 14);
 	EXPECT_EQ(sl[j(j1 = 2, j2 = 1)], 34);
     }
+
+    TEST_F(ContractionTest, ContractionWorks)
+    {
+	MultiArray2dAny ma2(ra);
+
+	ma2("alpha","gamma") = ma("alpha","beta","gamma").contract<Range1dAny>(sum<int>(),"beta");
+
+	auto i = ma2.beginIndex();
+	auto i1 = i.template getIndex<0>();
+	auto i2 = i.template getIndex<1>();
+
+	EXPECT_EQ(ma2[i(i1 = 0, i2 = 0)], 33);
+	EXPECT_EQ(ma2[i(i1 = 0, i2 = 1)], 54);
+
+	EXPECT_EQ(ma2[i(i1 = 1, i2 = 0)], 65);
+	EXPECT_EQ(ma2[i(i1 = 1, i2 = 1)], 82);
+
+	EXPECT_EQ(ma2[i(i1 = 2, i2 = 0)], 229);
+	EXPECT_EQ(ma2[i(i1 = 2, i2 = 1)], -114);
+    }
     
 } // end namespace