+ IndexInfo

src/multi_array_operation.h

@@ -93,6 +93,10 @@ namespace MultiArrayTools
 	OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
 			std::shared_ptr<typename CRange::IndexType>& index);
 
+	OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
+			std::shared_ptr<typename CRange::IndexType>& index,
+			std::shared_ptr<VIWB> blockIndex);
+
 	MBlock<T>& get();
 	const Block<T>& get() const;
 
@@ -101,7 +105,7 @@ namespace MultiArrayTools
 	
     protected:
 	
-	//void performAssignment(const OperationBase<T>& in);
+	void performAssignment(std::intptr_t blockIndexNum);
 	OpClass const& mSecond;
 	MutableMultiArrayBase<T,Ranges...>& mArrayRef;
 	std::shared_ptr<IndexType> mIndex;
@@ -157,6 +161,8 @@ namespace MultiArrayTools
 	const MBlock<T>& get() const;
 	MBlock<T>& get();
 	
+	OperationRoot& set(std::shared_ptr<VIWB> blockIndex);
+	
 	std::vector<BTSS> block(const std::shared_ptr<VIWB> blockIndex, bool init = false) const;
 	const OperationRoot& block() const;
 	
@@ -165,6 +171,7 @@ namespace MultiArrayTools
 	MutableMultiArrayBase<T,Ranges...>& mArrayRef;
 	std::shared_ptr<IndexType> mIndex;
 	mutable MBlock<T> mBlock;
+	std::shared_ptr<VIWB> mBlockIndex; // predefine to save time
     };
     
     template <typename T, class OpFunction, class... Ops>
@@ -325,6 +332,33 @@ namespace MultiArrayTools
 	
 	block(blockIndex, true);
 	second.block(blockIndex, true);
+
+	performAssignment(blockIndexNum);
+    }
+
+    template <typename T, class OpClass, class... Ranges>
+    OperationMaster<T,OpClass,Ranges...>::
+    OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
+		    std::shared_ptr<typename CRange::IndexType>& index,
+		    std::shared_ptr<VIWB> blockIndex) :
+	mSecond(second), mArrayRef(ma), mIndex()
+    {
+	MultiRangeFactory<Ranges...> mrf( index->range() );
+	std::shared_ptr<MultiRange<Ranges...> > mr =
+	    std::dynamic_pointer_cast<MultiRange<Ranges...> >( mrf.create() );
+	mIndex = std::make_shared<IndexType>( mr->begin() );
+	(*mIndex) = *index;
+
+	std::intptr_t blockIndexNum = blockIndex->getPtrNum();
+	second.block(blockIndex, true);
+
+	performAssignment(blockIndexNum);
+    }
+
+    
+    template <typename T, class OpClass, class... Ranges>    
+    void OperationMaster<T,OpClass,Ranges...>::performAssignment(std::intptr_t blockIndexNum)
+    {
 	//size_t cnt = 0;
 	//std::clock_t cs = clock();
 	for(*mIndex = 0; mIndex->pos() != mIndex->max(); mIndex->pp(blockIndexNum) ){
@@ -417,7 +451,8 @@ namespace MultiArrayTools
     OperationRoot(MutableMultiArrayBase<T,Ranges...>& ma,
 		  const std::shared_ptr<typename Ranges::IndexType>&... indices) :
 	OperationTemplate<T,OperationRoot<T,Ranges...> >(this),
-	mArrayRef(ma), mIndex( std::make_shared<IndexType>( mArrayRef.range() ) )
+	mArrayRef(ma), mIndex( std::make_shared<IndexType>( mArrayRef.range() ) ),
+	mBlockIndex(nullptr)
     {
 	(*mIndex)(indices...);
     }
@@ -426,8 +461,13 @@ namespace MultiArrayTools
     template <class OpClass>
     OperationMaster<T,OpClass,Ranges...> OperationRoot<T,Ranges...>::operator=(const OpClass& in)
     {
+	if(mBlockIndex){
+	    return OperationMaster<T,OpClass,Ranges...>(mArrayRef, in, mIndex, mBlockIndex);
+	}
+	else {
 	    return OperationMaster<T,OpClass,Ranges...>(mArrayRef, in, mIndex);
 	}
+    }
 
     template <typename T, class... Ranges>
     const MBlock<T>& OperationRoot<T,Ranges...>::get() const
@@ -443,6 +483,14 @@ namespace MultiArrayTools
 	return mBlock;
     }
 
+    template <typename T, class... Ranges>
+    OperationRoot<T,Ranges...>&
+    OperationRoot<T,Ranges...>::set(std::shared_ptr<VIWB> blockIndex)
+    {
+	mBlockIndex = blockIndex;
+	return *this;
+    }
+    
     template <typename T, class... Ranges>
     std::vector<BTSS> OperationRoot<T,Ranges...>::block(const std::shared_ptr<VIWB> blockIndex, bool init) const
     {

src/ranges/container_range.h

@@ -78,6 +78,8 @@ namespace MultiArrayTools
 	std::shared_ptr<VIWB> getVPtr(size_t n);
 	size_t getStepSize(size_t n);
 
+	std::vector<IndexInfo> infoVec() const;
+	
 	std::string id();
 	void print(size_t offset);
 
@@ -332,6 +334,14 @@ namespace MultiArrayTools
 	return mBlockSizes[n+1];
     }
 
+    template <class... Indices>
+    std::vector<IndexInfo> ContainerIndex<Indices...>::infoVec() const
+    {
+	std::vector<IndexInfo> out(sizeof...(Indices));
+	RPackNum<sizeof...(Indices)-1>::buildInfoVec(out, mIPack);
+	return std::move( out );
+    }
+    
     template <class... Indices>
     std::string ContainerIndex<Indices...>::id()
     {

src/ranges/index_base.h

@@ -11,6 +11,7 @@
 #include "range_base.h"
 #include "index_type.h"
 #include "vindex_wrapper.h"
+#include "index_info.h"
 
 namespace MultiArrayTools
 {
@@ -57,6 +58,8 @@ namespace MultiArrayTools
 	*/
 	std::shared_ptr<VIWB> getVPtr(size_t n) const { return THIS().getVPtr(n); }
 
+	std::vector<IndexInfo> infoVec() const { return THIS().infoVec(); }
+	
 	size_t getStepSize(size_t n) const { return THIS().getStepSize(n); }
 	
 	operator size_t() const;

src/ranges/index_info.cc

@@ -0,0 +1,33 @@
+
+#include "index_info.h"
+
+namespace MultiArrayTools
+{
+
+    const IndexInfo* IndexInfo::getPtr(size_t inum) const
+    {
+	return &mNext[inum];
+    }
+    
+    std::intptr_t IndexInfo::getPtrNum() const
+    {
+	return mPtrNum;
+    }
+    
+    size_t IndexInfo::dim() const
+    {
+	return mDim;
+    }
+    
+    size_t IndexInfo::getStepSize(size_t inum) const
+    {
+	return mNext[inum]->getStepSzize();
+    }
+    
+    size_t IndexInfo::getStepSize() const
+    {
+	return mStepSize;
+    }
+
+    
+} // end namespace MultiArrayTools

src/ranges/index_info.h

@@ -0,0 +1,44 @@
+// -*- C++ -*-
+
+#ifndef __index_info_h__
+#define __index_info_h__
+
+#include <cstdlib>
+#include <vector>
+
+namespace MultiArrayTools
+{
+    class IndexInfo;
+    
+    class IndexInfo
+    {
+    public:
+	IndexInfo() = delete;
+
+	template <class IndexClass>
+	IndexInfo(const IndexClass& ind, size_t stepSize = 1);
+	
+	const IndexInfo* getPtr(size_t inum) const;
+	std::intptr_t getPtrNum() const;
+	size_t dim() const;
+	size_t getStepSize(size_t inum) const;
+	size_t getStepSize() const;
+
+    private:
+	std::vector<IndexInfo> mNext;
+	std::intptr_t mPtrNum;
+	size_t mDim;
+	size_t mStepSize;
+    };
+
+    template <class IndexClass>
+    IndexInfo::IndexInfo(const IndexClass& ind, size_t stepSize) :
+	mNext(ind.infoVec()),
+	mPtrNum( reinterpret_cast<std::intptr_t>( &ind ) ),
+	mDim(ind.rangePtr()->dim()),
+	mStepSize(stepSize)
+    {}
+    
+} // end namespace MultiArrayTools
+
+#endif

src/ranges/multi_range.h

@@ -94,6 +94,8 @@ namespace MultiArrayTools
 	std::shared_ptr<VIWB> getVPtr(size_t n);
 	size_t getStepSize(size_t n);
 
+	std::vector<IndexInfo> infoVec() const;
+	
 	std::string id();
 	void print(size_t offset);
     };
@@ -376,6 +378,14 @@ namespace MultiArrayTools
 	return mBlockSizes[n+1];
     }
 
+    template <class... Indices>
+    std::vector<IndexInfo> MultiIndex<Indices...>::infoVec() const
+    {
+	std::vector<IndexInfo> out(sizeof...(Indices));
+	RPackNum<sizeof...(Indices)-1>::buildInfoVec(out, mIPack);
+	return std::move( out );
+    }
+    
     template <class... Indices>
     std::string MultiIndex<Indices...>::id()
     {

src/ranges/rpack_num.h

@@ -4,6 +4,7 @@
 
 #include <memory>
 #include "vindex_wrapper.h"
+#include "index_info.h"
 
 namespace MultiArrayHelper
 {
@@ -207,6 +208,15 @@ namespace MultiArrayHelper
 	    static_assert( Range::defaultable, "not defaultable" );
 	    RPackNum<N-1>::template checkDefaultable<Ranges...>();
 	}
+
+	template <class... Indices>
+	static void buildInfoVec(std::vector<IndexInfo>& out,
+				 std::tuple<std::shared_ptr<Indices>...>& ip)
+	{
+	    out.emplace_back(*std::get<sizeof...(Indices)-N-1>(ip));
+	    RPackNum<N-1>::buildInfoVec(out, ip);
+	}
+
     };
 
     
@@ -361,6 +371,13 @@ namespace MultiArrayHelper
 	    static_assert( Range::defaultable, "not defaultable" );
 	}
 
+	template <class... Indices>
+	static void buildInfoVec(std::vector<IndexInfo>& out,
+				 std::tuple<std::shared_ptr<Indices>...>& ip)
+	{
+	    out.emplace_back(*std::get<sizeof...(Indices)-1>(ip));
+	}
+
     };
 	
 

src/ranges/single_range.h

@@ -53,6 +53,8 @@ namespace MultiArrayTools
 	std::shared_ptr<VIWB> getVPtr(size_t n);
 	size_t getStepSize(size_t n);
 
+	std::vector<IndexInfo> infoVec() const;
+	
 	std::string id();
 	void print(size_t offset);
     };
@@ -206,6 +208,12 @@ namespace MultiArrayTools
 	return std::shared_ptr<VIWB>();
     }
 
+    template <typename U, SpaceType TYPE>
+    std::vector<IndexInfo> SingleIndex<U,TYPE>::infoVec() const
+    {
+	return std::move( std::vector<IndexInfo>() );
+    }
+    
     template <typename U, SpaceType TYPE>
     size_t SingleIndex<U,TYPE>::getStepSize(size_t n)
     {

src/tests/op_perf_test.cc

@@ -9,6 +9,8 @@
 #include <ctime>
 #include <cmath>
 
+#define ONLY_SPIN
+
 namespace MAT = MultiArrayTools;
 
 namespace {
@@ -62,7 +64,7 @@ namespace {
     {
 	return std::make_tuple(static_cast<size_t>( ts )...);
     }
-    
+#ifndef ONLY_SPIN    
     class OpTest_Performance
     {
     public:
@@ -119,7 +121,7 @@ namespace {
 	std::vector<double> cv1;
 	std::vector<double> cv2;
     };
-
+#endif
     class OpTest_Spin
     {
     public:
@@ -160,11 +162,12 @@ namespace {
 	auto gamma = MAT::getIndex<SR>();
 	auto delta = MAT::getIndex<SR>();
 	auto deltap = MAT::getIndex<SR>();
+	auto vdeltap = MAT::make_viwb(deltap);
 	
 	auto mix = MAT::mkMIndex( alpha, beta, gamma );
 
 	std::clock_t begin = std::clock();
-	res1(delta, deltap) = ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
+	res1(delta, deltap).set(vdeltap) = ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
 	std::clock_t end = std::clock();
 	std::cout << "MultiArray time: " << static_cast<double>( end - begin ) / CLOCKS_PER_SEC
 		  << std::endl;
@@ -212,7 +215,7 @@ namespace {
 		  << std::endl;
 	std::cout << "ratio: " << static_cast<double>( end - begin ) / static_cast<double>( end2 - begin2 ) << std::endl;
     }    
-    
+#ifndef ONLY_SPIN    
     void OpTest_Performance::PCheck()
     {
 	MultiArray<double,MRange> ma2(mrptr, cv2);
@@ -251,16 +254,18 @@ namespace {
 	//assert( xround( res.at(mkt(700,900)) ) == xround(res2[700*vs1 + 900]) );
 	
     }
-    
+#endif    
 } // anonymous namspace
 
 int main(int argc, char** argv)
 {
+#ifndef ONLY_SPIN
     OpTest_Performance pt;
-    OpTest_Spin st;
-
-    st.contract();
     pt.PCheck();
+#endif
+    OpTest_Spin st;
+    st.contract();
+
     
     return 0;
 }

src/tests/op_unit_test.cc

@@ -215,11 +215,12 @@ namespace {
 	auto gamma = MAT::getIndex<SR>();
 	auto delta = MAT::getIndex<SR>();
 	auto deltap = MAT::getIndex<SR>();
+	auto vdeltap = MAT::make_viwb(deltap);
 	
 	auto mix = MAT::mkMIndex( alpha, beta, gamma );
 
 	std::clock_t begin = std::clock();
-	res1(delta, deltap) = ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
+	res1(delta, deltap).set(vdeltap) = ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
 	std::clock_t end = std::clock();
 	std::cout << "MultiArray time: " << static_cast<double>( end - begin ) / CLOCKS_PER_SEC
 		  << std::endl;