corrections + partially pre-compile hl operations

.gitignore

@@ -4,3 +4,4 @@ build/
 install/
 */build/
 */install/
+src/lib/hl_ops/*.cc

src/include/high_level_operation.cc.h

@@ -176,8 +176,6 @@ namespace MultiArrayTools
         typename B::template RetT<Inds...> res;
         Create<N-1>::template cx<Inds...>::template ccx<ROP,OpF>::template cccx<N>
 	    (res,mIn,inds...);
-	//	assert(res.op.init());
-	//assert(res.outer.init());
         return res;
     }
 
@@ -262,26 +260,47 @@ namespace MultiArrayTools
         if(dim > 2){
 	    auto ci1 = di->getP(dim-2)->reduced();
 	    auto ci2 = di->getP(dim-1)->reduced();
-            //auto ci1 = std::dynamic_pointer_cast<IndexWrapper<CI>>(di->getP(dim-2));
-            //auto ci2 = std::dynamic_pointer_cast<IndexWrapper<CI>>(di->getP(dim-1));
             assert(ci1 != nullptr);
             assert(ci2 != nullptr);
             auto odi = mkSubSpaceX(di, dim-2);
             auto mi = mkMIndex(is..., odi);
-            //this->assign(in, mi, ci1->getIndex(), ci2->getIndex());
 	    this->assign(in, mi, ci1, ci2);
         }
         else {
             assert(dim == 2 or dim == 1);
-            //auto ci1 = std::dynamic_pointer_cast<IndexWrapper<CI>>(di->getP(dim-2));
 	    auto ci1 = di->getP(dim-1)->reduced();
             assert(ci1 != nullptr);
             auto odi = mkSubSpaceX(di, dim-1);
             auto mi = mkMIndex(is..., odi);
             this->assign(in, mi, ci1);
         }
-        //INDS<ROP,Indices...>::template CallHLOp<> call;
+        return *this;
-        //call.assign(*this, in, is..., di);
+    }
+
+    template <class ROP>
+    template <class... Indices>
+    HighLevelOpHolder<ROP>& HighLevelOpHolder<ROP>::xplus(const HighLevelOpHolder& in,
+							  const std::shared_ptr<DynamicIndex>& di,
+							  const std::shared_ptr<Indices>&... is)
+    {
+        const size_t dim = di->dim();
+        if(dim > 2){
+	    auto ci1 = di->getP(dim-2)->reduced();
+	    auto ci2 = di->getP(dim-1)->reduced();
+            assert(ci1 != nullptr);
+            assert(ci2 != nullptr);
+            auto odi = mkSubSpaceX(di, dim-2);
+            auto mi = mkMIndex(is..., odi);
+	    this->plus(in, mi, ci1, ci2);
+        }
+        else {
+            assert(dim == 2 or dim == 1);
+	    auto ci1 = di->getP(dim-1)->reduced();
+            assert(ci1 != nullptr);
+            auto odi = mkSubSpaceX(di, dim-1);
+            auto mi = mkMIndex(is..., odi);
+            this->plus(in, mi, ci1);
+        }
         return *this;
     }
 
@@ -305,12 +324,11 @@ namespace MultiArrayTools
                                                            const std::shared_ptr<MIndex>& mi,
                                                            const std::shared_ptr<Indices>&... inds)
     {
-        //VCHECK(printInd(inds...));
         auto xx = mkArrayPtr<double>(nullr());
         ROP& opr = *mOp->get();
 	if(in.root()){
 	    auto inx = in;
-	    opr.par() = *inx.get();
+	    opr.par().assign( *inx.get(), mkMIndex(mi,inds...) )();
 	    return *this;
 	}
         auto loop = mkPILoop
@@ -318,8 +336,6 @@ namespace MultiArrayTools
                 auto inx = in;
                 auto dop = inx.create(inds...);
 		DynamicO<size_t> gexp;
-		//VCHECK(dop.outer.init());
-		//VCHECK(dop.op.init());
 		if(dop.outer.init()){
 		    gexp = mkDynOp1<size_t>(mkMOp<size_t>(dop.outer,dop.op));
 		}
@@ -347,7 +363,7 @@ namespace MultiArrayTools
         ROP& opr = *mOp->get();
 	if(in.root()){
 	    auto inx = in;
-	    opr.par() += *inx.get();
+	    opr.par().plus( *inx.get(), mkMIndex(mi,inds...) )();
 	    return *this;
 	}
         auto loop = mkPILoop

src/include/high_level_operation.h

@@ -13,6 +13,11 @@ namespace MultiArrayTools
     typedef CR::IndexType CI;
     typedef std::shared_ptr<CI> CIP;
 
+    typedef OperationRoot<double,CR,DynamicRange> OpCD;
+    typedef OperationRoot<double,DynamicRange> OpD;
+    extern template class OperationRoot<double,CR,DynamicRange>;
+    extern template class OperationRoot<double,DynamicRange>;
+    
     template <typename T, class Op>
     DynamicO<T> mkDynOp1(const Op& op);
 
@@ -90,6 +95,11 @@ namespace MultiArrayTools
    
     };
 
+    extern template class HighLevelOpBase<OpCD>;
+    extern template class HighLevelOpBase<OpD>;
+    extern template class HighLevelOpRoot<OpCD>;
+    extern template class HighLevelOpRoot<OpD>;
+    
     template <class OpF, class... Ops>
     auto mkFOp(const Ops&... ops)
     {
@@ -125,6 +135,20 @@ namespace MultiArrayTools
 #undef reg_ind3
     };
 
+    extern template class HighLevelOp<OpCD,plusx<double,double>,2>;
+    extern template class HighLevelOp<OpCD,minusx<double,double>,2>;
+    extern template class HighLevelOp<OpCD,multipliesx<double,double>,2>;
+    extern template class HighLevelOp<OpCD,dividesx<double,double>,2>;
+    extern template class HighLevelOp<OpD,plusx<double,double>,2>;
+    extern template class HighLevelOp<OpD,minusx<double,double>,2>;
+    extern template class HighLevelOp<OpD,multipliesx<double,double>,2>;
+    extern template class HighLevelOp<OpD,dividesx<double,double>,2>;
+    
+#define regFunc1(fff) \
+    extern template class HighLevelOp<OpCD,x_##fff<double>,1>; \
+    extern template class HighLevelOp<OpD,x_##fff<double>,1>;
+#include "extensions/math.h"
+#undef regFunc1
     
     template <class ROP>
     class HighLevelOpHolder
@@ -161,6 +185,11 @@ namespace MultiArrayTools
                                    const std::shared_ptr<DynamicIndex>& di,
                                    const std::shared_ptr<Indices>&... is);
         
+        template <class... Indices>
+        HighLevelOpHolder& xplus(const HighLevelOpHolder& in,
+				 const std::shared_ptr<DynamicIndex>& di,
+				 const std::shared_ptr<Indices>&... is);
+
         template <class MIndex, class... Indices>
         HighLevelOpHolder& assign(const HighLevelOpHolder& in,
                                   const std::shared_ptr<MIndex>& mi,
@@ -172,6 +201,9 @@ namespace MultiArrayTools
                                 const std::shared_ptr<Indices>&... inds);
     };
 
+    extern template class HighLevelOpHolder<OpCD>;
+    extern template class HighLevelOpHolder<OpD>;
+    
     template <class F, class ROP, class... ROPs>
     HighLevelOpHolder<ROP> mkSFunc(const HighLevelOpHolder<ROP>& a, const HighLevelOpHolder<ROPs>&... as);
 
@@ -179,87 +211,19 @@ namespace MultiArrayTools
     template <class ROP>
     HighLevelOpHolder<ROP> mkHLO(const ROP& op);
 
-#define SP " "
 #define regFunc1(fff) template <class ROP> \
     HighLevelOpHolder<ROP> hl_##fff (const HighLevelOpHolder<ROP>& in);
 #include "extensions/math.h"
 #undef regFunc1
-#undef SP
-    /*
-    template <size_t N>
-    struct SetLInds
-    {
-	template <class ITuple>
-	static inline void mkLIT(const ITuple& itp, const std::shared_ptr<DynamicIndex>& di);
     
-	template <class Tar, class ITp, typename... Args>
+#define regFunc1(fff) template <class ROP>				\
-        struct xx
+    HighLevelOpHolder<ROP> hl_##fff (const HighLevelOpHolder<ROP>& in); \
-        {
+    extern template HighLevelOpHolder<OpCD> hl_##fff (const HighLevelOpHolder<OpCD>& in); \
-            template <class... Is>
+    extern template HighLevelOpHolder<OpD> hl_##fff (const HighLevelOpHolder<OpD>& in);
-            static inline void assign(Tar& tar, const Args&... args,
+#include "extensions/math.h"
-				      const ITp& itp, const std::shared_ptr<Is>&... is);
+#undef regFunc1
    
-            template <class... Is>
-            static inline void plus(Tar& tar, const Args&... args,
-				    const ITp& itp, const std::shared_ptr<Is>&... is);
-        };
-    };
 
-    template <>
-    struct SetLInds<0>
-    {
-	template <class ITuple>
-	static inline void mkLIT(const ITuple& itp, const std::shared_ptr<DynamicIndex>& di);
-
-        template <class Tar, class ITp, typename... Args>
-        struct xx
-        {
-            template <class... Is>
-            static inline void assign(Tar& tar, const Args&... args,
-				      const ITp& itp, const std::shared_ptr<Is>&... is);
-
-            template <class... Is>
-            static inline void plus(Tar& tar, const Args&... args,
-				    const ITp& itp, const std::shared_ptr<Is>&... is);
-        };
-    };
-
-    template <class ROP, class... Indices>
-    struct INDS
-    {
-        class CallHLOpBase
-        {
-        private:
-            size_t mDepth;
-        public:
-	    size_t depth() const;
-	    
-            void assign(const HighLevelOpHolder<ROP>& target, const HighLevelOpHolder<ROP>& source,
-                        const std::shared_ptr<Indices>&... is,
-                        const std::shared_ptr<DynamicIndex>& di) const = 0;
-
-            void plus(const HighLevelOpHolder<ROP>& target, const HighLevelOpHolder<ROP>& source,
-                      const std::shared_ptr<Indices>&... is,
-                      const std::shared_ptr<DynamicIndex>& di) const = 0;
-        };
-
-        template <class... LIndices>
-        class CallHLOp
-        {
-        private:
-            typedef std::tuple<std::shared_ptr<LIndices>...> ITuple;
-            static vector<std::shared_ptr<CallHLOpBase>> sNext;
-        public:
-            void assign(HighLevelOpHolder<ROP>& target, const HighLevelOpHolder<ROP>& source,
-                        const std::shared_ptr<Indices>&... is,
-                        const std::shared_ptr<DynamicIndex>& di) const override final;
-
-            void plus(HighLevelOpHolder<ROP>& target, const HighLevelOpHolder<ROP>& source,
-                      const std::shared_ptr<Indices>&... is,
-                      const std::shared_ptr<DynamicIndex>& di) const override final;
-        };
-    };
-    */
 }
 
 #endif

src/include/multi_array_operation.cc.h

@@ -711,6 +711,17 @@ namespace MultiArrayTools
                                       (mOrigDataPtr,*this,in)));
     }
 
+    template <typename T, class... Ranges>
+    template <class OpClass, class Index>
+    auto ParallelOperationRoot<T,Ranges...>::assign(const OpClass& in, const std::shared_ptr<Index>& i) const
+        -> decltype(i->pifor(1,in.loop(AssignmentExpr2<T,ParallelOperationRoot<T,Ranges...>,OpClass>
+				       (mOrigDataPtr,*this,in))))
+    {
+        static_assert( OpClass::SIZE == decltype(in.rootSteps())::SIZE, "Ext Size mismatch" );
+        return i->pifor(1,in.loop(AssignmentExpr2<T,ParallelOperationRoot<T,Ranges...>,OpClass>
+				  (mOrigDataPtr,*this,in)));
+    }
+
     template <typename T, class... Ranges>
     template <class OpClass>
     auto ParallelOperationRoot<T,Ranges...>::plus(const OpClass& in)
@@ -721,6 +732,17 @@ namespace MultiArrayTools
                                       (mOrigDataPtr,*this,in)));
     }
 
+    template <typename T, class... Ranges>
+    template <class OpClass, class Index>
+    auto ParallelOperationRoot<T,Ranges...>::plus(const OpClass& in, const std::shared_ptr<Index>& i) const
+        -> decltype(i->pifor(1,in.loop(AddExpr<T,ParallelOperationRoot<T,Ranges...>,OpClass>
+				       (mOrigDataPtr,*this,in))))
+    {
+        static_assert( OpClass::SIZE == decltype(in.rootSteps())::SIZE, "Ext Size mismatch" );
+        return i->pifor(1,in.loop(AddExpr<T,ParallelOperationRoot<T,Ranges...>,OpClass>
+				  (mOrigDataPtr,*this,in)));
+    }
+
     template <typename T, class... Ranges>
     template <class OpClass>
     ParallelOperationRoot<T,Ranges...>& ParallelOperationRoot<T,Ranges...>::operator=(const OpClass& in)

src/include/multi_array_operation.h

@@ -587,11 +587,21 @@ namespace MultiArrayTools
             -> decltype(mIndex.pifor(1,in.loop(AssignmentExpr2<T,ParallelOperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
                                                (mOrigDataPtr,*this,in))));
 
-        template <class OpClass>
+        template <class OpClass, class Index>
+        auto assign(const OpClass& in, const std::shared_ptr<Index>& i) const
+            -> decltype(i->pifor(1,in.loop(AssignmentExpr2<T,ParallelOperationRoot<T,Ranges...>,OpClass>
+                                          (mOrigDataPtr,*this,in))));
+
+	template <class OpClass>
         auto plus(const OpClass& in)
             -> decltype(mIndex.pifor(1,in.loop(AddExpr<T,ParallelOperationRoot<T,Ranges...>,OpClass,OpIndexAff::TARGET>
                                                (mOrigDataPtr,*this,in))));
 
+        template <class OpClass, class Index>
+        auto plus(const OpClass& in, const std::shared_ptr<Index>& i) const
+            -> decltype(i->pifor(1,in.loop(AddExpr<T,ParallelOperationRoot<T,Ranges...>,OpClass>
+					   (mOrigDataPtr,*this,in))));
+
         template <class OpClass>
         ParallelOperationRoot& operator=(const OpClass& in);
 

src/lib/CMakeLists.txt

@@ -1,4 +1,7 @@
 
+execute_process ( COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/mk_hl_op.sh
+  WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/ )
+
 set(libmultiarray_a_SOURCES
   ${CMAKE_SOURCE_DIR}/src/lib/ranges/range_base.cc
   ${CMAKE_SOURCE_DIR}/src/lib/ranges/anonymous_range.cc
@@ -16,6 +19,12 @@ foreach(ccfile ${cc_files})
     ${ccfile})
 endforeach(ccfile)
 
+file(GLOB cc_files "${CMAKE_SOURCE_DIR}/src/lib/hl_ops/*.cc")
+foreach(ccfile ${cc_files})
+  set(libmultiarray_a_SOURCES ${libmultiarray_a_SOURCES}
+    ${ccfile})
+endforeach(ccfile)
+
 add_library(multiarray_obj OBJECT
   ${libmultiarray_a_SOURCES}
   )

src/lib/high_level_operation.cc

@@ -21,4 +21,15 @@ namespace MultiArrayTools
 	return odi;
     }
 
+    template class OperationRoot<double,CR,DynamicRange>;
+    template class OperationRoot<double,DynamicRange>;
+
+    template class HighLevelOpHolder<OpCD>;
+    template class HighLevelOpHolder<OpD>;
+
+    template class HighLevelOpBase<OpCD>;
+    template class HighLevelOpBase<OpD>;
+    template class HighLevelOpRoot<OpCD>;
+    template class HighLevelOpRoot<OpD>;
+
 }

src/lib/mk_hl_op.sh

@@ -0,0 +1,34 @@
+#!/bin/bash
+
+for x in $(cat ../include/extensions/math.h) ; do
+    test "${x}" = "#ifdef" && continue
+    test "${x}" = "#endif" && continue
+    test "${x}" = "regFunc1" && continue
+    xx=${x#regFunc1\(}
+    fff=${xx%\)}
+    file=hl_ops/${fff}.cc
+    test -f ${file} && rm -f ${file}
+    echo "#include \"multi_array_header.h\"" >> ${file}
+    echo "#include \"high_level_operation.h\"" >> ${file}
+    echo "" >> ${file}
+    echo "namespace MultiArrayTools" >> ${file}
+    echo "{" >> ${file}
+    echo "    template class HighLevelOp<OpCD,x_${fff}<double>,1>;" >> ${file}
+    echo "    template class HighLevelOp<OpD,x_${fff}<double>,1>;" >> ${file}
+    echo "    template HighLevelOpHolder<OpCD> hl_${fff} (const HighLevelOpHolder<OpCD>& in);" >> ${file}
+    echo "    template HighLevelOpHolder<OpD> hl_${fff} (const HighLevelOpHolder<OpD>& in);" >> ${file}
+    echo "}" >> ${file}
+done
+
+for fff in plus minus multiplies divides ; do
+    file=hl_ops/${fff}.cc
+    test -f ${file} && rm -f ${fff}
+    echo "#include \"multi_array_header.h\"" >> ${file}
+    echo "#include \"high_level_operation.h\"" >> ${file}
+    echo "" >> ${file}
+    echo "namespace MultiArrayTools" >> ${file}
+    echo "{" >> ${file}
+    echo "    template class HighLevelOp<OpCD,${fff}x<double,double>,2>;" >> ${file}
+    echo "    template class HighLevelOp<OpD,${fff}x<double,double>,2>;" >> ${file}
+    echo "}" >> ${file}
+done