diff --git a/src/base_def.h b/src/base_def.h
index 0e9723c..f181565 100644
--- a/src/base_def.h
+++ b/src/base_def.h
@@ -146,6 +146,10 @@ namespace MultiArrayTools
     template <typename T, class OpFunction, class... Ops>
     class Operation;
 
+    // multi_array_operation.h
+    template <typename T, class Op, class IndexType>
+    class Contraction;
+
     /*
     // multi_array_operation.h
     template <typename T, class InRange, class TotalInRange, class OutRange, class TotalRange>
diff --git a/src/block.cc b/src/block.cc
index 113c5ba..a351545 100644
--- a/src/block.cc
+++ b/src/block.cc
@@ -30,6 +30,7 @@ namespace MultiArrayHelper
     void BlockBinaryOpSelf<T,OpFunc,BlockClass>::operator()(const BlockClass& arg)
     {
 	static OpFunc f;
+	if(mRes.size() == 0) { mRes.assign(arg.size(), static_cast<T>(0)); }
 	assert(mRes.size() == arg.size());
 	for(size_t i = 0; i != arg.size(); ++i){
 	    mRes[i] = f(mRes[i], arg[i]);
@@ -49,46 +50,7 @@ namespace MultiArrayHelper
     {
 	return mSize;
     }
-    /*
-    template <typename T>
-    template <class OpFunction>
-    BlockResult<T> BlockBase<T>::operate(const BlockBase<T>& in)
-    {
-	assert(mSize == in.size());
-	OpFunction f;
-	BlockResult<T> res(mSize);
-	CHECK;
-	for(size_t i = 0; i != mSize; ++i){
-	    res[i] = f((*this)[i], in[i]);
-	}
-	return res;
-    }
-    
-    
-    template <typename T>
-    BlockResult<T> BlockBase<T>::operator+(const BlockBase<T>& in)
-    {
-	return operate<std::plus<T> >(in);
-    }
 
-    template <typename T>
-    BlockResult<T> BlockBase<T>::operator-(const BlockBase<T>& in)
-    {
-	return operate<std::minus<T> >(in);
-    }
-    
-    template <typename T>
-    BlockResult<T> BlockBase<T>::operator*(const BlockBase<T>& in)
-    {
-	return operate<std::multiplies<T> >(in);
-    }
-    
-    template <typename T>
-    BlockResult<T> BlockBase<T>::operator/(const BlockBase<T>& in)
-    {
-	return operate<std::divides<T> >(in);
-    }
-    */
     /************************
      *   MutableBlockBase   *
      ************************/
@@ -244,49 +206,11 @@ namespace MultiArrayHelper
     }
 
     template <typename T>
-    BlockResult<T>& BlockResult<T>::assing(const T& val)
+    BlockResult<T>& BlockResult<T>::assign(size_t size, const T& val)
     {
-	mRes.assing(BB::mSize, val);
+	BB::mSize = size;
+	mRes.assign(BB::mSize, val);
 	return *this;
     }
     
-    /*
-    template <typename T>
-    BlockResult<T>& BlockResult<T>::operator+=(const BlockBase<T>& in)
-    {
-	return operateSelf<std::plus<T> >(in);
-    }
-
-    template <typename T>
-    BlockResult<T>& BlockResult<T>::operator-=(const BlockBase<T>& in)
-    {
-	return operateSelf<std::minus<T> >(in);
-    }
-
-    template <typename T>
-    BlockResult<T>& BlockResult<T>::operator*=(const BlockBase<T>& in)
-    {
-	return operateSelf<std::multiplies<T> >(in);
-    }
-
-    template <typename T>
-    BlockResult<T>& BlockResult<T>::operator/=(const BlockBase<T>& in)
-    {
-	return operateSelf<std::divides<T> >(in);
-    }
-
-    template <typename T>
-    template <class OpFunction>
-    BlockResult<T>& BlockResult<T>::operateSelf(const BlockBase<T>& in)
-    {
-	assert(BB::mSize == in.size());
-	OpFunction f;
-	//BlockResult<T> res(mSize);
-	for(size_t i = 0; i != BB::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 468e6b6..ffb6610 100644
--- a/src/block.h
+++ b/src/block.h
@@ -50,17 +50,9 @@ namespace MultiArrayHelper
 	BlockBase(size_t size);
 
 	size_t size() const;
-	/*
-	template <class OpFunction>
-	BlockResult<T> operate(const BlockBase& in);
-	
-	BlockResult<T> operator+(const BlockBase& in);
-	BlockResult<T> operator-(const BlockBase& in);
-	BlockResult<T> operator*(const BlockBase& in);
-	BlockResult<T> operator/(const BlockBase& in);
-	*/
+
     protected:
-	size_t mSize;
+	size_t mSize = 0;
     };
 
     template <typename T>
@@ -136,7 +128,7 @@ namespace MultiArrayHelper
 	template <class BlockClass>
 	BlockResult& operator=(const BlockClass& in);
 
-	BlockResult& assing(const T& val);
+	BlockResult& assign(size_t size, const T& val);
 	
 	BlockType type() const;
 	const T& operator[](size_t pos) const;
@@ -144,14 +136,6 @@ namespace MultiArrayHelper
 	BlockResult& set(size_t npos);
 	size_t stepSize() const;
 	
-	//BlockResult<T>& operator+=(const BlockBase<T>& in);
-	//BlockResult<T>& operator-=(const BlockBase<T>& in);
-	//BlockResult<T>& operator*=(const BlockBase<T>& in);
-	//BlockResult<T>& operator/=(const BlockBase<T>& in);
-
-	//template <class OpFunction>
-	//BlockResult<T>& operateSelf(const BlockBase<T>& in);
-	
     protected:
 	std::vector<T> mRes;
     };
diff --git a/src/multi_array_operation.cc b/src/multi_array_operation.cc
index 6e12a72..647e882 100644
--- a/src/multi_array_operation.cc
+++ b/src/multi_array_operation.cc
@@ -160,7 +160,7 @@ namespace MultiArrayTools
     {
     	return Operation<T,std::divides<T>,OperationClass,Second>(*mOc, in);
     }
-    /*
+
     template <typename T, class OperationClass>
     template <class IndexType>
     auto OperationTemplate<T,OperationClass>::c(std::shared_ptr<IndexType>& ind) const
@@ -168,7 +168,7 @@ namespace MultiArrayTools
     {
 	return Contraction<T,OperationClass,IndexType>(*mOc, ind);
     }
-    */
+
     
     /*************************
      *   OperationMaster     *
@@ -324,7 +324,6 @@ namespace MultiArrayTools
     const BlockResult<T>& Operation<T,OpFunction,Ops...>::get() const
     {
 	mRes = std::move( PackNum<sizeof...(Ops)-1>::template unpackArgs<T,OpFunction>(mOps) );
-	//CHECK;
 	return mRes;
     }
 
@@ -350,15 +349,14 @@ namespace MultiArrayTools
     template <typename T, class Op, class IndexType>
     Contraction<T,Op,IndexType>::Contraction(const Op& op, std::shared_ptr<IndexType> ind) :
 	OperationTemplate<T,Contraction<T,Op,IndexType> >(this),
-	mOp(op) {}
+	mOp(op),
+	mInd(ind) {}
 
     template <typename T, class Op, class IndexType>
     const BlockResult<T>& Contraction<T,Op,IndexType>::get() const
     {
 	BlockBinaryOpSelf<T,std::plus<T>,BlockResult<T> > f(mRes);
-	mRes.assign( static_cast<T>(0) );
 	for(*mInd = 0; mInd->pos() != mInd->max(); ++(*mInd)){
-	    //mRes += mOp.get();
 	    f(mOp.get());
 	}
 	return mRes;
@@ -367,9 +365,7 @@ namespace MultiArrayTools
     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<0>::makeBlockTypeVec(btv, std::make_tuple( mOp ), blockIndex);
-	return btv;
+	return mOp.block(blockIndex);
     }
 
     template <typename T, class Op, class IndexType>
diff --git a/src/multi_array_operation.h b/src/multi_array_operation.h
index 1a9bc39..6e17fa0 100644
--- a/src/multi_array_operation.h
+++ b/src/multi_array_operation.h
@@ -108,11 +108,11 @@ namespace MultiArrayTools
 	template <class Second>
 	auto operator/(const Second& in) const
 	    -> Operation<T,std::divides<T>,OperationClass,Second>;
-	/*
+
 	template <class IndexType>
 	auto c(std::shared_ptr<IndexType>& ind) const
 	    -> Contraction<T,OperationClass,IndexType>;
-	*/
+
     private:
 	OperationClass* mOc;
     };
@@ -244,7 +244,7 @@ namespace MultiArrayTools
 	
     protected:
 
-	const Op& mOp;
+	Op mOp;
 	std::shared_ptr<IndexType> mInd;
 	mutable BlockResult<T> mRes;
     };
diff --git a/src/op_unit_test.cc b/src/op_unit_test.cc
index 1ca0767..5c38d68 100644
--- a/src/op_unit_test.cc
+++ b/src/op_unit_test.cc
@@ -239,6 +239,22 @@ namespace {
 
     }
 
+    TEST_F(OpTest_MDim, ExecContract)
+    {
+	MultiArray<double,SRange> res(sr2ptr);
+	const MultiArray<double,SRange> ma1(sr2ptr, v1);
+	const MultiArray<double,SRange> ma2(sr4ptr, v2);
+
+	auto i1 = std::dynamic_pointer_cast<SRange::IndexType>( sr2ptr->index() );
+	auto i2 = std::dynamic_pointer_cast<SRange::IndexType>( sr4ptr->index() );
+
+	res(i1) = (ma1(i1) * ma2(i2)).c(i2);
+
+	EXPECT_EQ( xround( res.at('1') ), xround(2.917 * 8.870 + 2.917 * 4.790) );
+	EXPECT_EQ( xround( res.at('2') ), xround(9.436 * 8.870 + 9.436 * 4.790) );
+	EXPECT_EQ( xround( res.at('3') ), xround(0.373 * 8.870 + 0.373 * 4.790) );
+    }
+    
     TEST_F(OpTest_MDim, ExecOp2)
     {
 	MultiArray<double,MRange,SRange> res(mr1ptr,sr4ptr);