contraction code (no compile test)

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

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;
     };

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;
+    }
+    
 }