From 64d10867dd251b4fca90a6eda0d526ef6c3b9d0e Mon Sep 17 00:00:00 2001
From: Christian Zimmermann <christian.zimmermann@ur.de>
Date: Mon, 18 Jan 2021 19:00:12 +0100
Subject: [PATCH] finish parts of auto-vectorization

---
 src/include/arith.h                    |   8 +-
 src/include/dynamic_operation.h        |   7 ++
 src/include/multi_array_operation.cc.h |  63 ++++++++-----
 src/include/multi_array_operation.h    | 117 ++++++++++++++++++++-----
 src/include/type_operations.h          |  31 ++++++-
 5 files changed, 180 insertions(+), 46 deletions(-)

diff --git a/src/include/arith.h b/src/include/arith.h
index db34d75..e9339d8 100644
--- a/src/include/arith.h
+++ b/src/include/arith.h
@@ -48,7 +48,8 @@ namespace MultiArrayTools
     {
 	static constexpr bool FISSTATIC = true;
 	typedef T value_type;
-	
+	typedef F function;
+
 	template <class... Ops>
 	static auto mk(const Ops&... ops)
 	    -> Operation<T,F,Ops...>
@@ -74,6 +75,11 @@ namespace MultiArrayTools
 	{
 	    return a;
 	}
+
+	static inline T selfApply(T& a1, const T& a2)
+	{
+	    return a1 = a2;
+	}
     };
 
     template <typename T, typename U>
diff --git a/src/include/dynamic_operation.h b/src/include/dynamic_operation.h
index 12c1936..3835019 100644
--- a/src/include/dynamic_operation.h
+++ b/src/include/dynamic_operation.h
@@ -17,6 +17,7 @@ namespace MultiArrayTools
 	
 	static constexpr size_t SIZE = 1;
 	static constexpr bool CONT = true;
+	static constexpr bool VABLE = false;
 	
 	DynamicOperationBase() = default;
 	DynamicOperationBase(const DynamicOperationBase& in) = default;
@@ -125,6 +126,7 @@ namespace MultiArrayTools
 	
 	static constexpr size_t SIZE = 1;
 	static constexpr bool CONT = true;
+	static constexpr bool VABLE = false;
 	
 	DynamicO() = default;
 	DynamicO(const DynamicO& in) : mOp(in.mOp ? in.mOp->deepCopy() : nullptr) {}
@@ -144,6 +146,11 @@ namespace MultiArrayTools
 	
         template <class X>
 	inline T get(const DExtTX<X>& pos) const { return mOp->get(pos.reduce()); }
+
+	template <typename V,class X>
+	inline auto vget(const DExtTX<X>& pos) const
+	{ return mOp->template vget<V>(pos.reduce()); }
+	
         template <class X>
 	inline DynamicO& set(const DExtTX<X>& pos) { mOp->set(pos.reduce()); return *this; }
 	inline DExtT rootSteps(std::intptr_t iPtrNum = 0) const { return mOp->rootSteps(iPtrNum); }
diff --git a/src/include/multi_array_operation.cc.h b/src/include/multi_array_operation.cc.h
index b9d0207..0f1fa41 100644
--- a/src/include/multi_array_operation.cc.h
+++ b/src/include/multi_array_operation.cc.h
@@ -555,55 +555,78 @@ namespace MultiArrayTools
     template <typename T, class... Ranges>
     template <class OpClass>
     auto OperationRoot<T,Ranges...>::assign(const OpClass& in) const
-        -> decltype(this->template asx<IAssign>(in))
+        -> decltype(this->template asx<IAssign<T>>(in))
     {
-        return this->template asx<IAssign>(in);
+        return this->template asx<IAssign<T>>(in);
     }
 
     template <typename T, class... Ranges>
     template <class OpClass>
     auto OperationRoot<T,Ranges...>::assignExpr(const OpClass& in) const
-        -> decltype(this->template asxExpr<IAssign>(in))
+        -> decltype(this->template asxExpr<IAssign<T>>(in))
     {
-        return this->template asxExpr<IAssign>(in);
+        return this->template asxExpr<IAssign<T>>(in);
     }
     
     template <typename T, class... Ranges>
     template <class OpClass, class Index>
     auto OperationRoot<T,Ranges...>::assign(const OpClass& in, const std::shared_ptr<Index>& i) const
-        -> decltype(this->template asx<IAssign>(in,i))
+        -> decltype(this->template asx<IAssign<T>>(in,i))
     {
-        return this->template asx<IAssign>(in,i);
+        return this->template asx<IAssign<T>>(in,i);
     }
 
     template <typename T, class... Ranges>
     template <class OpClass>
     auto OperationRoot<T,Ranges...>::plus(const OpClass& in) const
-        -> decltype(this->template asx<IPlus>(in))
+        -> decltype(this->template asx<IPlus<T>>(in))
     {
-        return this->template asx<IPlus>(in);
+        return this->template asx<IPlus<T>>(in);
     }
 
     template <typename T, class... Ranges>
     template <class OpClass, class Index>
     auto OperationRoot<T,Ranges...>::plus(const OpClass& in, const std::shared_ptr<Index>& i) const
-        -> decltype(this->template asx<IPlus>(in,i))
+        -> decltype(this->template asx<IPlus<T>>(in,i))
     {
-        return this->template asx<IPlus>(in,i);
+        return this->template asx<IPlus<T>>(in,i);
     }
-    
+
+    template <bool VABLE = false>
+    struct VExec
+    {
+	template <typename TarOp, class OpClass>
+	static inline void exec(TarOp& th, const OpClass& in)
+	{
+	    th.assign(in)();
+	}
+    };
+
+    template <>
+    struct VExec<true>
+    {
+	template <typename TarOp, class OpClass>
+	static inline void exec(TarOp& th, const OpClass& in)
+	{
+	    CHECK;
+	    typedef typename TarOp::value_type T;
+	    auto x = th.template asx<IVAssign<typename VType<T>::type,T>>(in);
+	    const size_t inum = x.vec(VType<T>::MULT);
+	    if(x.rootSteps(inum) == 1){
+		CHECK;
+		x();
+	    }
+	    else {
+		th.assign(in)();
+	    }
+	}
+    };
+
     template <typename T, class... Ranges>
     template <class OpClass>
     OperationRoot<T,Ranges...>& OperationRoot<T,Ranges...>::operator=(const OpClass& in)
     {
-        auto x = this->template asx<IVAssign<typename VType<T>::type>>(in);
-	const size_t inum = x.vec(VType<T>::MULT);
-	if(x.rootSteps(inum) == 1){
-	    x();
-	}
-	else {
-	    assign(in)();
-	}
+	VExec<OpClass::VABLE>::exec(*this,in);
         return *this;
     }
 
@@ -896,7 +919,7 @@ namespace MultiArrayTools
     {
 	typedef std::tuple<Ops...> OpTuple;
 	return PackNum<sizeof...(Ops)-1>::
-	    template mkVOpExpr<SIZE,V,ET,OpTuple,OpFunction>(mF, pos, mOps); // implement!!!
+	    template mkVOpExpr<SIZE,V,ET,OpTuple,VFunc<OpFunction>>(mkVFuncPtr(mF), pos, mOps); // implement!!!
     }
 
     template <typename T, class OpFunction, class... Ops>
diff --git a/src/include/multi_array_operation.h b/src/include/multi_array_operation.h
index 53e284b..bc35d0e 100644
--- a/src/include/multi_array_operation.h
+++ b/src/include/multi_array_operation.h
@@ -204,12 +204,74 @@ namespace MultiArrayTools
                            TARGET = 1
     };
 
+    template <class T>
+    struct VType
+    {
+	typedef T type;
+	static constexpr size_t MULT = sizeof(type)/sizeof(T);
+    };
+
+    template <>
+    struct VType<double>
+    {
+	typedef v256 type;
+	static constexpr size_t MULT = sizeof(type)/sizeof(double);
+    };
+
+    template <template <typename...> class F,typename... Ts>
+    inline auto mkVFuncPtr(const std::shared_ptr<F<Ts...>>& f)
+    {
+	return std::shared_ptr<F<typename VType<Ts>::type...>>();
+	// empty, implement corresponding constructors...!!!
+    }
+
+    template <template <typename...> class F,typename... Ts>
+    inline auto mkVFunc(const F<Ts...>& f)
+    {
+	return F<typename VType<Ts>::type...>();
+	// empty, implement corresponding constructors...!!!
+    }
+
+    template <class F>
+    using VFunc = decltype(mkVFunc(std::declval<F>()));
+
+    template <class F>
+    struct IAccess
+    {
+	template <typename T, typename Op, class ExtType>
+	static inline void f(T*& t, size_t pos, const Op& op, ExtType e)
+	{
+	    F::selfApply(t[pos],op.get(e));
+	}
+    };
+
+    template <typename V, class F>
+    struct IVAccess
+    {
+	template <typename T, typename Op, class ExtType>
+	static inline void f(T*& t, size_t pos, const Op& op, ExtType e)
+	{
+	    VFunc<F>::selfApply(reinterpret_cast<V*>(t)[pos],op.template vget<V>(e));
+	}
+    };
+
+    template <typename T>
+    using IAssign = IAccess<identity<T>>;
+
+    template <typename T>
+    using IPlus = IAccess<plus<T>>;
+
+    template <typename V, typename T>
+    using IVAssign = IVAccess<V,identity<T>>;
+
+    template <typename V, typename T>
+    using IVPlus = IVAccess<V,plus<T>>;
+    /*
     struct IAssign
     {
 	template <typename T, typename Op, class ExtType>
 	static inline void f(T*& t, size_t pos, const Op& op, ExtType e)
 	{
-	    //a = b;
 	    t[pos] = op.get(e);
 	}
     };
@@ -223,7 +285,7 @@ namespace MultiArrayTools
 	    reinterpret_cast<V*>(t)[pos] = op.template vget<V>(e);
 	}
     };
-
+    
     struct IPlus
     {
 	template <typename T, typename Op, class ExtType>
@@ -242,20 +304,8 @@ namespace MultiArrayTools
 	    reinterpret_cast<V*>(t)[pos] += op.template vget<V>(e);
 	}
     };
+    */
 
-    template <class T>
-    struct VType
-    {
-	typedef T type;
-	static constexpr size_t MULT = sizeof(type)/sizeof(T);
-    };
-
-    template <>
-    struct VType<double>
-    {
-	typedef v256 type;
-	static constexpr size_t MULT = sizeof(type)/sizeof(double);
-    };
     
     template <typename T, class IOp, class Target, class OpClass, OpIndexAff OIA=OpIndexAff::EXTERN>
     class AssignmentExpr : public ExpressionBase
@@ -295,10 +345,10 @@ namespace MultiArrayTools
     };
 
     template <typename T, class Target, class OpClass, OpIndexAff OIA=OpIndexAff::EXTERN>
-    using AssignmentExpr2 = AssignmentExpr<T,IAssign,Target,OpClass,OIA>;
+    using AssignmentExpr2 = AssignmentExpr<T,IAssign<T>,Target,OpClass,OIA>;
 
     template <typename T, class Target, class OpClass, OpIndexAff OIA=OpIndexAff::EXTERN>
-    using AddExpr = AssignmentExpr<T,IPlus,Target,OpClass,OIA>;
+    using AddExpr = AssignmentExpr<T,IPlus<T>,Target,OpClass,OIA>;
 
     template <typename T, class... Ops>
     class MOp
@@ -401,6 +451,7 @@ namespace MultiArrayTools
 
 	static constexpr size_t SIZE = 1;
         static constexpr bool CONT = true;
+        static constexpr bool VABLE = true;
 	
 	ConstOperationRoot(const MultiArrayBase<T,Ranges...>& ma,
 			   const std::shared_ptr<typename Ranges::IndexType>&... indices);
@@ -449,6 +500,7 @@ namespace MultiArrayTools
 
 	static constexpr size_t SIZE = Op::SIZE;
         static constexpr bool CONT = false;
+        static constexpr bool VABLE = false;
 	
 	StaticCast(const Op& op);
 
@@ -487,6 +539,7 @@ namespace MultiArrayTools
 	
 	static constexpr size_t SIZE = 1;
         static constexpr bool CONT = false;
+        static constexpr bool VABLE = false;
         
 	MetaOperationRoot(const std::shared_ptr<IndexType>& ind);
 
@@ -522,6 +575,7 @@ namespace MultiArrayTools
 
 	static constexpr size_t SIZE = 1;
         static constexpr bool CONT = true;
+        static constexpr bool VABLE = true;
 
     private:
 
@@ -554,23 +608,23 @@ namespace MultiArrayTools
             
         template <class OpClass>
         auto assign(const OpClass& in) const
-	    -> decltype(this->template asx<IAssign>(in));
+	    -> decltype(this->template asx<IAssign<T>>(in));
 
 	template <class OpClass>
         auto assignExpr(const OpClass& in) const
-	    -> decltype(this->template asxExpr<IAssign>(in));
+	    -> decltype(this->template asxExpr<IAssign<T>>(in));
 	    
         template <class OpClass, class Index>
         auto assign(const OpClass& in, const std::shared_ptr<Index>& i) const
-	    -> decltype(this->template asx<IAssign>(in,i));
+	    -> decltype(this->template asx<IAssign<T>>(in,i));
 	
 	template <class OpClass>
         auto plus(const OpClass& in) const
-	    -> decltype(this->template asx<IPlus>(in));
+	    -> decltype(this->template asx<IPlus<T>>(in));
 
         template <class OpClass, class Index>
         auto plus(const OpClass& in, const std::shared_ptr<Index>& i) const
-	    -> decltype(this->template asx<IPlus>(in,i));
+	    -> decltype(this->template asx<IPlus<T>>(in,i));
             
         template <class OpClass>
         OperationRoot& operator=(const OpClass& in);
@@ -616,6 +670,7 @@ namespace MultiArrayTools
 
 	static constexpr size_t SIZE = 1;
         static constexpr bool CONT = true;
+        static constexpr bool VABLE = true;
 
     private:
 
@@ -686,6 +741,7 @@ namespace MultiArrayTools
 
 	static constexpr size_t SIZE = 0;
         static constexpr bool CONT = true;
+        static constexpr bool VABLE = false;
 
 	OperationValue(const T& val);
 
@@ -706,8 +762,18 @@ namespace MultiArrayTools
     private:
 	T mVal;
     };
-    
 
+    template <class Op> 
+    inline constexpr bool isVAble()
+    {
+	return Op::VABLE;
+    }
+
+    template <class Op1, class Op2, class... Ops>
+    inline constexpr bool isVAble()
+    {
+	return Op1::VABLE and isVAble<Op2,Ops...>();
+    }
     
     template <typename T, class OpFunction, class... Ops>
     class Operation : public OperationTemplate<T,Operation<T,OpFunction,Ops...> >
@@ -721,6 +787,7 @@ namespace MultiArrayTools
 	static constexpr size_t SIZE = RootSum<Ops...>::SIZE;
 	static constexpr bool FISSTATIC = OpFunction::FISSTATIC;
         static constexpr bool CONT = false;
+        static constexpr bool VABLE = isVAble<Ops...>();
 
     private:
 	std::tuple<Ops...> mOps;
@@ -795,6 +862,7 @@ namespace MultiArrayTools
 
 	static constexpr size_t SIZE = Op::SIZE;
         static constexpr bool CONT = Op::CONT;
+        static constexpr bool VABLE = Op::VABLE;
         
     private:
 
@@ -839,7 +907,8 @@ namespace MultiArrayTools
 
 	static constexpr size_t SIZE = Op::SIZE;
         static constexpr bool CONT = false;
-        
+	static constexpr bool VABLE = false;
+
     private:
 
 	mutable Op mOp;
diff --git a/src/include/type_operations.h b/src/include/type_operations.h
index cb6903b..7b438f2 100644
--- a/src/include/type_operations.h
+++ b/src/include/type_operations.h
@@ -10,6 +10,7 @@
 #include "mbase_def.h"
 
 #include "pack_num.h"
+#include <cmath>
 
 namespace MultiArrayTools
 {
@@ -470,7 +471,35 @@ namespace MultiArrayTools
         xsdiv<4>( o._x, a._x );
 	return o;
     }
-   
+
+    inline double xpow(const double& b, const double& e)
+    {
+	return pow(b,e);
+    }
+
+    inline v256 pow(const v256& b, const v256& e)
+    {
+	v256 out;
+	for(int i = 0; i < 4; i++){
+	    out._x[i] = xpow(b._x[i],e._x[i]);
+	}
+	return out;
+    }
+
+    inline double xexp(const double& a)
+    {
+	return exp(a);
+    }
+
+    inline v256 exp(const v256& a)
+    {
+	v256 out;
+	for(int i = 0; i < 4; i++){
+	    out._x[i] = xexp(a._x[i]);
+	}
+	return out;
+    }
+
 } // namespace MultiArrayTools
 
 #endif