diff --git a/src/include/functional_multi_array.h b/src/include/functional_multi_array.h
index c7149e1..4719b41 100644
--- a/src/include/functional_multi_array.h
+++ b/src/include/functional_multi_array.h
@@ -10,7 +10,7 @@ namespace MultiArrayTools
 
     template <class Index>
     auto indexToSlice(const std::shared_ptr<Index>& i)
-	-> Slice<typename Index::MetaType,typename Index::RangeType>;
+	-> ConstSlice<typename Index::MetaType,typename Index::RangeType>;
 
     
     template <typename T, class Function, class... SRanges>
@@ -19,15 +19,18 @@ namespace MultiArrayTools
     public:
 
 	typedef ContainerRange<T,SRanges...> CRange;
-	typedef MultiArrayBase<T,CRange> MAB;
+	typedef MultiArrayBase<T,SRanges...> MAB;
 	//typedef typename MultiArrayBase<T,CRange>::const_iterator const_iterator;
 	typedef ContainerIndex<T,typename SRanges::IndexType...> IndexType;
 	//typedef typename CRange::IndexType IndexType;
-
+	typedef MultiArray<T,SRanges...> MAType;
+	
     private:
 	mutable T mVal;
 	Function mFunc;
 
+	mutable std::shared_ptr<MAType> mMaPtr;
+
     public:
 	
 	DEFAULT_MEMBERS(FunctionalMultiArray);
@@ -45,6 +48,9 @@ namespace MultiArrayTools
 
 	auto exec(std::shared_ptr<typename SRanges::IndexType>&... inds) const
 	    -> decltype( mkOperation( mFunc, ConstOperationRoot<T,SRanges>( indexToSlice( inds ), inds) ... ) );
+
+	virtual ConstOperationRoot<T,SRanges...>
+	operator()(std::shared_ptr<typename SRanges::IndexType>&... inds) const override;
 	
     };
 
@@ -60,9 +66,9 @@ namespace MultiArrayTools
     
     template <class Index>
     auto indexToSlice(const std::shared_ptr<Index>& i)
-	-> Slice<typename Index::MetaType, typename Index::RangeType>
+	-> ConstSlice<typename Index::MetaType, typename Index::RangeType>
     {
-	return Slice<typename Index::MetaType, typename Index::RangeType>( i->range(), i->metaPtr() );
+	return ConstSlice<typename Index::MetaType, typename Index::RangeType>( i->range(), i->metaPtr() );
     }
 
 
@@ -136,6 +142,18 @@ namespace MultiArrayTools
 	return false;
     }
 
+    template <typename T, class Function, class... SRanges>
+    ConstOperationRoot<T,SRanges...> FunctionalMultiArray<T,Function,SRanges...>::
+    operator()(std::shared_ptr<typename SRanges::IndexType>&... inds) const
+    {
+	if(not mMaPtr){
+	    mMaPtr = std::make_shared<MAType>( MAB::mRange->space() );
+	    (*mMaPtr)(inds...) = exec(inds...);
+	}
+	return ConstOperationRoot<T,SRanges...>( *mMaPtr, inds... );
+    }
+
+    
     template <typename T, class Function, class... SRanges>
     auto FunctionalMultiArray<T,Function,SRanges...>::
     exec(std::shared_ptr<typename SRanges::IndexType>&... inds) const
diff --git a/src/include/multi_array.h b/src/include/multi_array.h
index 7fedb46..014976c 100644
--- a/src/include/multi_array.h
+++ b/src/include/multi_array.h
@@ -22,7 +22,8 @@ namespace MultiArrayTools
 	MultiArray(const std::shared_ptr<SRanges>&... ranges);
 	MultiArray(const std::shared_ptr<SRanges>&... ranges, const std::vector<T>& vec);
 	MultiArray(const std::shared_ptr<SRanges>&... ranges, std::vector<T>&& vec);
-
+	MultiArray(const typename CRange::SpaceType& space);
+	
 	// Only if ALL ranges have default extensions:
 	//MultiArray(const std::vector<T>& vec);
 	//MultiArray(std::vector<T>&& vec);
@@ -70,6 +71,14 @@ namespace MultiArrayTools
      *  MultiArray     *	     
      *******************/
 
+    template <typename T, class... SRanges>
+    MultiArray<T,SRanges...>::MultiArray(const typename CRange::SpaceType& space) :
+	MutableMultiArrayBase<T,SRanges...>(space),
+	mCont(MAB::mRange->size())
+    {
+	MAB::mInit = true;
+    }
+	
     template <typename T, class... SRanges>
     MultiArray<T,SRanges...>::MultiArray(const std::shared_ptr<SRanges>&... ranges) :
 	MutableMultiArrayBase<T,SRanges...>(ranges...),
diff --git a/src/include/ranges/single_range.h b/src/include/ranges/single_range.h
index 8ab53ce..7a4c0b9 100644
--- a/src/include/ranges/single_range.h
+++ b/src/include/ranges/single_range.h
@@ -48,7 +48,7 @@ namespace MultiArrayTools
 	int mm(std::intptr_t idxPtrNum);
 	
 	U meta();
-	U* metaPtr();
+	const U* metaPtr() const;
 	SingleIndex& at(const U& metaPos);
 	
 	size_t dim(); // = 1
@@ -213,7 +213,7 @@ namespace MultiArrayTools
     }
 
     template <typename U, SpaceType TYPE>
-    U* SingleIndex<U,TYPE>::metaPtr()
+    const U* SingleIndex<U,TYPE>::metaPtr() const
     {
 	return mMetaPtr;
     }
diff --git a/src/include/slice.h b/src/include/slice.h
index 43afdd0..2fc81ab 100644
--- a/src/include/slice.h
+++ b/src/include/slice.h
@@ -7,6 +7,41 @@
 
 namespace MultiArrayTools
 {
+    template <typename T, class... SRanges>
+    class ConstSlice : public MultiArrayBase<T,SRanges...>
+    {
+    public:
+
+	typedef ContainerRange<T,SRanges...> CRange;
+	typedef MultiArrayBase<T,SRanges...> MAB;
+	typedef ContainerIndex<T,typename SRanges::IndexType...> IType;
+	
+	DEFAULT_MEMBERS(ConstSlice);
+
+	ConstSlice(const std::shared_ptr<SRanges>&... ranges, const T* data = nullptr);
+	    
+	virtual const T& operator[](const IType& i) const override;
+	virtual const T& at(const typename IType::MetaType& meta) const override;
+	
+	virtual const T* data() const override;
+	
+	virtual bool isSlice() const override;
+
+	virtual auto begin() const -> IType override;
+	virtual auto end() const -> IType override;
+
+	auto define(const std::shared_ptr<typename SRanges::IndexType>&... inds)
+	    -> SliceDef<T,SRanges...>;
+	
+    private:
+	friend SliceDef<T,SRanges...>;
+
+	void format(const std::array<size_t,sizeof...(SRanges)+1>& blocks);
+	
+	const T* mData;
+    };
+
+    
     template <typename T, class... SRanges>
     class Slice : public MutableMultiArrayBase<T,SRanges...>
     {
@@ -53,10 +88,14 @@ namespace MultiArrayTools
 	
     private:
 	IType mIndex;
-	Slice<T,SRanges...>& mSl;
+	Slice<T,SRanges...>* mSlPtr = nullptr;
+	ConstSlice<T,SRanges...>* mCSlPtr = nullptr;
 	
 	SliceDef() = default;
     public:
+	SliceDef(ConstSlice<T,SRanges...>& csl,
+		 const std::shared_ptr<typename SRanges::IndexType>&... inds);
+
 	SliceDef(Slice<T,SRanges...>& sl,
 		 const std::shared_ptr<typename SRanges::IndexType>&... inds);
 
@@ -72,10 +111,79 @@ namespace MultiArrayTools
 
 namespace MultiArrayTools
 {
+
+    /*******************
+     *   ConstSlice    *
+     *******************/
+
+    template <typename T, class... SRanges>
+    void ConstSlice<T,SRanges...>::format(const std::array<size_t,sizeof...(SRanges)+1>& blocks)
+    {
+	MAB::mProtoI->format(blocks);
+    }
     
-    /*************
-     *  Slice    *
-     *************/
+    template <typename T, class... SRanges>
+    ConstSlice<T,SRanges...>::ConstSlice(const std::shared_ptr<SRanges>&... ranges, const T* data) :
+	MultiArrayBase<T,SRanges...>(ranges...),
+	mData(data) {}
+    
+    template <typename T, class... SRanges>
+    const T& ConstSlice<T,SRanges...>::operator[](const IType& i) const
+    {
+	//assert(i.sliceMode()); // -> compare objects !!!!!
+	assert(i.container() == reinterpret_cast<std::intptr_t>(this));
+	return mData[ i.pos() ];
+    }
+    
+    template <typename T, class... SRanges>
+    const T& ConstSlice<T,SRanges...>::at(const typename IType::MetaType& meta) const
+    {
+	//auto x = begin().at(meta);
+	//VCHECK(x.pos());
+	return mData[ begin().at(meta).pos() ];
+    }
+    
+    template <typename T, class... SRanges>
+    const T* ConstSlice<T,SRanges...>::data() const
+    {
+	return mData;
+    }
+
+    template <typename T, class... SRanges>
+    bool ConstSlice<T,SRanges...>::isSlice() const
+    {
+	return true;
+    }
+
+    template <typename T, class... SRanges>
+    auto ConstSlice<T,SRanges...>::begin() const -> ConstSlice<T,SRanges...>::IType 
+    {
+	IType i(*MAB::mProtoI);
+	i = 0;
+	//i = mStartPos;
+	return i.setData(data());
+    }
+
+    template <typename T, class... SRanges>
+    auto ConstSlice<T,SRanges...>::end() const -> ConstSlice<T,SRanges...>::IType 
+    {
+	IType i(*MAB::mProtoI);
+	i = i.max(); // CHECK !!!
+	//i = std::get<sizeof...(SRanges)>(mBlockSizes);
+	return i.setData(data());
+    }
+
+    template <typename T, class... SRanges>
+    auto ConstSlice<T,SRanges...>::define(const std::shared_ptr<typename SRanges::IndexType>&... inds)
+	-> SliceDef<T,SRanges...>
+    {
+	return SliceDef<T,SRanges...>(*this, inds...);
+    }
+
+    
+    /**************
+     *   Slice    *
+     **************/
 
     template <typename T, class... SRanges>
     void Slice<T,SRanges...>::format(const std::array<size_t,sizeof...(SRanges)+1>& blocks)
@@ -163,11 +271,20 @@ namespace MultiArrayTools
 	return SliceDef<T,SRanges...>(*this, inds...);
     }
     
+    template <typename T, class... SRanges>
+    SliceDef<T,SRanges...>::SliceDef(ConstSlice<T,SRanges...>& csl,
+				     const std::shared_ptr<typename SRanges::IndexType>&... inds) :
+	mIndex(csl.begin()),
+	mCSlPtr(&csl)
+    {
+	mIndex(inds...);
+    }
+
     template <typename T, class... SRanges>
     SliceDef<T,SRanges...>::SliceDef(Slice<T,SRanges...>& sl,
 				     const std::shared_ptr<typename SRanges::IndexType>&... inds) :
 	mIndex(sl.begin()),
-	mSl(sl)
+	mSlPtr(&sl)
     {
 	mIndex(inds...);
     }
@@ -179,10 +296,14 @@ namespace MultiArrayTools
 	std::array<size_t,sizeof...(SRanges)+1> blocks;
 	PackNum<sizeof...(SRanges)-1>::
 	    template mkSliceBlocks<T,OperationRoot<T,ORanges...>,SRanges...>(blocks, mIndex, op);
-	mSl.format(blocks);
-	//VCHECK(blocks[0]);
-	//VCHECK(blocks[1]);
-	mSl.mData = op.data();
+	if(mCSlPtr){
+	    mCSlPtr->format(blocks);
+	    mCSlPtr->mData = op.data();
+	}
+	else {
+	    mSlPtr->format(blocks);
+	    mSlPtr->mData = op.data();
+	}
 	return *this;
     }