diff --git a/src/block.cc b/src/block.cc
index ea91bc2..ce5ea44 100644
--- a/src/block.cc
+++ b/src/block.cc
@@ -226,5 +226,43 @@ namespace MultiArrayHelper
     {
 	return 1;
     }
+
+    template <typename T>
+    BlockResult<T>& BlockResult<T>::operator+=(const BlockBase& in)
+    {
+	return operateSelf<std::plus<T> >(in);
+    }
+
+    template <typename T>
+    BlockResult<T>& BlockResult<T>::operator-=(const BlockBase& in)
+    {
+	return operateSelf<std::minus<T> >(in);
+    }
+
+    template <typename T>
+    BlockResult<T>& BlockResult<T>::operator*=(const BlockBase& in)
+    {
+	return operateSelf<std::multiplies<T> >(in);
+    }
+
+    template <typename T>
+    BlockResult<T>& BlockResult<T>::operator/=(const BlockBase& in)
+    {
+	return operateSelf<std::divides<T> >(in);
+    }
+
+    template <typename T>
+    template <class OpFunction>
+    BlockResult<T>& BlockResult<T>::operateSelf(const BlockBase<T>& in)
+    {
+	assert(mSize == in.size());
+	OpFunction f;
+	//BlockResult<T> res(mSize);
+	for(size_t i = 0; i != mSize; ++i){
+	    (*this)[i] = f((*this)[i], in[i]);
+	}
+	return *this;
+    }
+
     
 } // end namespace MultiArrayHelper
diff --git a/src/block.h b/src/block.h
index 8efdfe0..09e9768 100644
--- a/src/block.h
+++ b/src/block.h
@@ -129,6 +129,14 @@ namespace MultiArrayHelper
 	BlockResult& set(size_t npos);
 	size_t stepSize() const;
 	
+	BlockResult<T>& operator+=(const BlockBase& in);
+	BlockResult<T>& operator-=(const BlockBase& in);
+	BlockResult<T>& operator*=(const BlockBase& in);
+	BlockResult<T>& operator/=(const BlockBase& in);
+
+	template <class OpFunction>
+	BlockResult<T>& operateSelf(const BlockBase& in);
+	
     protected:
 	std::vector<T> mRes;
     };
diff --git a/src/multi_array_operation.cc b/src/multi_array_operation.cc
index 9d4d275..3aaefaf 100644
--- a/src/multi_array_operation.cc
+++ b/src/multi_array_operation.cc
@@ -342,4 +342,36 @@ namespace MultiArrayTools
 	return *this;
     }
 
+    /*********************
+     *   Contraction     *
+     *********************/
+
+    template <typename T, class Op, class IndexType>
+    Contraction(const Op& op, std::shared_ptr<IndexType> ind) :
+	OperationTemplate<T,Contraction<T,Op,IndexType> >(this),
+	mOp(op) {}
+
+    const BlockResult<T>& get() const
+    {
+	// set mRes = 0 !!!
+	for(*mIndex = 0; mIndex->pos() != mIndex->max(); ++(*mIndex)){
+	    mRes += mOp.get();
+	}
+	return mRes;
+    }
+    
+    template <typename T, class Op, class IndexType>
+    std::vector<BTSS> Contraction<T,Op,IndexType>::block(const std::shared_ptr<IndexBase> blockIndex) const
+    {
+	std::vector<BTSS> btv;
+	PackNum<sizeof...(Ops)-1>::makeBlockTypeVec(btv, mOp, blockIndex);
+	return btv;
+    }
+
+    template <typename T, class Op, class IndexType>
+    const Contraction& block() const
+    {
+	return *this;
+    }
+    
 }