diff --git a/src/include/helper_tools.h b/src/include/helper_tools.h
index 1932888..7a40a42 100644
--- a/src/include/helper_tools.h
+++ b/src/include/helper_tools.h
@@ -6,6 +6,7 @@
 #include "slice.h"
 #include <ostream>
 #include "pack_num.h"
+#include <functional>
 
 namespace MultiArrayTools
 {
@@ -48,6 +49,9 @@ namespace MultiArrayTools
     template <size_t N, class MArray>
     auto prtr(const MArray& ma)
 	-> decltype(ma.template getRangePtr<N>());
+
+    template <class IndexType>
+    inline void For(const std::shared_ptr<IndexType>& ind, const std::function<void(void)>& ll);
 }
 
 /* ========================= *
@@ -139,6 +143,14 @@ namespace MultiArrayTools
 	return ma.template getRangePtr<N>();
     }
 
+    template <class IndexType>
+    inline void For(const std::shared_ptr<IndexType>& ind, const std::function<void(void)>& ll)
+    {
+	for((*ind) = 0; ind->pos() != ind->max(); ++(*ind)){
+	    ll();
+	}
+    }
+
 }
 
 #endif
diff --git a/src/include/multi_array_operation.h b/src/include/multi_array_operation.h
index 337554e..d6f2de3 100644
--- a/src/include/multi_array_operation.h
+++ b/src/include/multi_array_operation.h
@@ -50,9 +50,13 @@ namespace MultiArrayTools
 	    -> Operation<T,divides<T>,OperationClass,Second>;
 
 	template <class IndexType>
-	auto c(std::shared_ptr<IndexType>& ind) const
+	auto c(const std::shared_ptr<IndexType>& ind) const
 	    -> Contraction<T,OperationClass,IndexType>;
 
+	template <class... Indices>
+	auto sl(const std::shared_ptr<Indices>&... inds) const
+	    -> ConstSlice<T,typename Indices::RangeType...>;
+
     private:		
 	friend OperationClass;
 	friend OperationTemplate<T,OperationClass>;
@@ -153,12 +157,14 @@ namespace MultiArrayTools
 
 	template <class Expr>
 	Expr loop(Expr exp) const;
-		
+
+	const T* data() const;
+	
     private:
 
 	//MultiArrayBase<T,Ranges...> const& mArrayRef;
 	const T* mDataPtr;
-	IndexType mIndex;
+	mutable IndexType mIndex;
 	//std::shared_ptr<MultiArrayBase<T,Ranges...> > mMaPtr;
     };
 
@@ -257,11 +263,15 @@ namespace MultiArrayTools
 
 	T* data() const;
 	
+	template <class... Indices>
+	auto sl(const std::shared_ptr<Indices>&... inds)
+	    -> Slice<T,typename Indices::RangeType...>;
+
     private:
 
 	//MutableMultiArrayBase<T,Ranges...>& mArrayRef;
 	T* mDataPtr;
-	IndexType mIndex;
+	mutable IndexType mIndex;
     };
 
     template <typename T>
@@ -454,12 +464,22 @@ namespace MultiArrayTools
 
     template <typename T, class OperationClass>
     template <class IndexType>
-    auto OperationBase<T,OperationClass>::c(std::shared_ptr<IndexType>& ind) const
+    auto OperationBase<T,OperationClass>::c(const std::shared_ptr<IndexType>& ind) const
 	-> Contraction<T,OperationClass,IndexType>
     {
 	return Contraction<T,OperationClass,IndexType>(THIS(), ind);
     }
 
+    template <typename T, class OperationClass>
+    template <class... Indices>
+    auto OperationBase<T,OperationClass>::sl(const std::shared_ptr<Indices>&... inds) const
+	-> ConstSlice<T,typename Indices::RangeType...>
+    {
+	ConstSlice<T,typename Indices::RangeType...> out(inds->range()...);
+	out.define(inds...) = *this;
+	return out;
+    }
+
     /*****************************************
      *   OperationMaster::AssignmentExpr     *
      *****************************************/
@@ -567,6 +587,12 @@ namespace MultiArrayTools
 	return mDataPtr[pos.val()];
     }
 
+    template <typename T, class... Ranges>
+    const T* ConstOperationRoot<T,Ranges...>::data() const
+    {
+	return mDataPtr + mIndex().pos();
+    }
+    
     template <typename T, class... Ranges>
     MExt<void> ConstOperationRoot<T,Ranges...>::rootSteps(std::intptr_t iPtrNum) const
     {
@@ -697,9 +723,23 @@ namespace MultiArrayTools
     template <typename T, class... Ranges>
     T* OperationRoot<T,Ranges...>::data() const
     {
-	return mDataPtr + mIndex.pos();
+	return mDataPtr + mIndex().pos();
     }
 
+    template <typename T, class... Ranges>
+    template <class... Indices>
+    auto OperationRoot<T,Ranges...>::sl(const std::shared_ptr<Indices>&... inds)
+	-> Slice<T,typename Indices::RangeType...>
+    {
+	Slice<T,typename Indices::RangeType...> out(inds->range()...);
+	out.define(inds...) = *this;
+	return out;
+    }
+
+
+    /************************
+     *   OperationValue     *
+     ************************/
 
     template <typename T>
     OperationValue<T>::OperationValue(const T& val) : mVal(val) {}