Ft arm (cbdbb210) · Commits · Michael Tesch / cppduals

CMakeLists.txt

+6 −1

Original line number	Diff line number	Diff line
		@@ -3,7 +3,7 @@
		#
		cmake_minimum_required (VERSION 3.1)
		project (cppduals
		VERSION 0.4.3
		VERSION 0.5.0
		LANGUAGES C CXX
		)
		include (GNUInstallDirs)
		@@ -221,3 +221,8 @@ if (ETAGS)
		COMMAND ${ETAGS} --language=c++ --append `find ${PROJECT_BINARY_DIR}/thirdparty/eigenX/src/eigenX -type f`
		)
		endif (ETAGS)

		#
		# Packaging
		#
		include (CPack)

duals/arch/NEON/ComplexDual.h

0 → 100644

+228 −0

Original line number	Diff line number	Diff line
		// This file is part of Eigen, a lightweight C++ template library
		// for linear algebra.
		//
		// Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
		// Copyright (C) 2010 Konstantinos Margaritis <markos@freevec.org>
		//
		// This Source Code Form is subject to the terms of the Mozilla
		// Public License v. 2.0. If a copy of the MPL was not distributed
		// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

		#ifndef EIGEN_CDUAL_NEON_H
		#define EIGEN_CDUAL_NEON_H

		namespace Eigen {

		namespace internal {

		//---------- float ----------
		struct Packet1cdf
		{
		EIGEN_STRONG_INLINE Packet1cdf() {}
		EIGEN_STRONG_INLINE explicit Packet1cdf(const float32x4_t & a) : v(a) {}
		EIGEN_STRONG_INLINE explicit Packet1cdf(const Packet2df & a) : v(a) {}
		//__m128 v;
		Packet2df v;
		};

		// Use the packet_traits defined in AVX/ComplexDual.h instead if we're
		// going to leverage AVX instructions.
		#if !defined(CPPDUALS_DONT_VECTORIZE_CDUAL)
		#if !defined(EIGEN_VECTORIZE_AVX)
		template<> struct packet_traits<std::complex<duals::dual<float> > > : default_packet_traits
		{
		typedef Packet1cdf type;
		typedef Packet1cdf half;
		enum {
		Vectorizable = 1,
		AlignedOnScalar = 0,
		size = 1,
		HasHalfPacket = 0,

		HasAdd = 1,
		HasSub = 1,
		HasMul = 1,
		HasDiv = 1,
		HasNegate = 1,
		HasAbs = 0,
		HasAbs2 = 0,
		HasMin = 0,
		HasMax = 0,
		HasSetLinear = 0
		};
		};
		#endif
		#endif

		template<> struct unpacket_traits<Packet1cdf> {
		typedef std::complex<duals::dual<float> > type;
		enum {size=1, alignment=Aligned16, masked_load_available=false, masked_store_available=false, vectorizable=true};
		typedef Packet1cdf half;
		};

		template<> EIGEN_STRONG_INLINE Packet1cdf padd<Packet1cdf>(const Packet1cdf& a, const Packet1cdf& b)
		{ return Packet1cdf(padd<Packet4f>(a.v.v, b.v.v)); }
		template<> EIGEN_STRONG_INLINE Packet1cdf psub<Packet1cdf>(const Packet1cdf& a, const Packet1cdf& b)
		{ return Packet1cdf(psub<Packet4f>(a.v.v, b.v.v)); }
		template<> EIGEN_STRONG_INLINE Packet1cdf pnegate(const Packet1cdf& a) { return Packet1cdf(pnegate(a.v)); }
		template<> EIGEN_STRONG_INLINE Packet1cdf pconj(const Packet1cdf& a)
		{
		uint32x4_t mask = { 0x00000000, 0x00000000, 0x80000000, 0x80000000 };
		return Packet1cdf(vreinterpretq_f32_u32(veorq_u32(a.v.v, mask)));
		}

		template<> EIGEN_STRONG_INLINE Packet1cdf pmul<Packet1cdf>(const Packet1cdf& a, const Packet1cdf& b)
		{
		// can probably do this more efficiently...
		//= a0.b0 - a2.b2
		//= a0.b1 + a1.b0 - (a2.b3 + a3.b2)
		//= a0.b2 + a2.b0
		//= a0.b3 + a1.b2 + a2.b1 + a3.b0
		//
		uint32x4_t mask = {0x00000000, 0xffffffff, 0x00000000, 0xffffffff};
		uint32x4_t nega = {0x80000000, 0x80000000, 0x00000000, 0x00000000};
		return Packet1cdf(vaddq_f32(vaddq_f32(vmulq_f32(vdupq_lane_u32(vget_low_f32(a.v.v), 0),
		b.v.v),
		vandq_u32(mask,
		vmulq_f32(vdupq_lane_u64(vget_low_u64(a.v.v), 0),
		vcombine_u64(vdup_lane_u32(vget_low_u32(b.v.v), 0),
		vsli_n_u64(vdup_lane_u32(vget_low_u32(b.v.v), 0),
		vget_high_u64(b.v.v), 32))))),
		veorq_u32
		(nega,
		vaddq_f32(vmulq_f32(vdupq_lane_u32(vget_high_f32(a.v.v), 0),
		vcombine_u64(vget_high_u64(b.v.v),
		vget_low_u64(b.v.v))),
		vandq_u32(mask,
		vmulq_f32(vsriq_n_u64(vcombine_u64(vget_high_u64(a.v.v),
		vget_high_u64(a.v.v)),
		vshlq_n_u64(vcombine_u64(vget_low_u64(a.v.v),
		vget_low_u64(a.v.v)), 32), 32),
		vcombine_u64(vsli_n_u64(vget_low_u64(b.v.v),
		vget_high_u64(b.v.v), 32),
		vdup_lane_u32(vget_low_u32(b.v.v), 0))))))));
		}

		template<> EIGEN_STRONG_INLINE Packet1cdf pand <Packet1cdf>(const Packet1cdf& a, const Packet1cdf& b)
		{ return Packet1cdf(vandq_u32(a.v.v, b.v.v)); }
		template<> EIGEN_STRONG_INLINE Packet1cdf por <Packet1cdf>(const Packet1cdf& a, const Packet1cdf& b)
		{ return Packet1cdf(vorrq_u32(a.v.v, b.v.v)); }
		template<> EIGEN_STRONG_INLINE Packet1cdf pxor <Packet1cdf>(const Packet1cdf& a, const Packet1cdf& b)
		{ return Packet1cdf(veorq_u32(a.v.v, b.v.v)); }
		template<> EIGEN_STRONG_INLINE Packet1cdf pandnot<Packet1cdf>(const Packet1cdf& a, const Packet1cdf& b)
		{ return Packet1cdf(vbicq_s32(b.v.v, a.v.v)); }

		template<> EIGEN_STRONG_INLINE Packet1cdf pload <Packet1cdf>(const std::complex<duals::dual<float> >* from)
		{ EIGEN_DEBUG_ALIGNED_LOAD return Packet1cdf(ploadu<Packet4f>((const float*)from)); }
		template<> EIGEN_STRONG_INLINE Packet1cdf ploadu<Packet1cdf>(const std::complex<duals::dual<float> >* from)
		{ EIGEN_DEBUG_UNALIGNED_LOAD return Packet1cdf(ploadu<Packet4f>((const float*)from)); }
		template<> EIGEN_STRONG_INLINE Packet1cdf pset1<Packet1cdf>(const std::complex<duals::dual<float> >& from)
		{ /* here we really have to use unaligned loads :( */ return ploadu<Packet1cdf>(&from); }

		template<> EIGEN_STRONG_INLINE Packet1cdf ploaddup<Packet1cdf>(const std::complex<duals::dual<float> >* from)
		{ return pset1<Packet1cdf>(*from); }

		// FIXME force unaligned store, this is a temporary fix
		template<> EIGEN_STRONG_INLINE void
		pstore <std::complex<duals::dual<float> > >(std::complex<duals::dual<float> > * to, const Packet1cdf& from)
		{ EIGEN_DEBUG_ALIGNED_STORE pstoreu((float*)to, from.v.v); }
		template<> EIGEN_STRONG_INLINE void
		pstoreu<std::complex<duals::dual<float> > >(std::complex<duals::dual<float> > * to, const Packet1cdf& from)
		{ EIGEN_DEBUG_UNALIGNED_STORE pstoreu((float*)to, from.v.v); }

		template<> EIGEN_STRONG_INLINE void
		prefetch<std::complex<duals::dual<float> > >(const std::complex<duals::dual<float> > * addr)
		{ __builtin_prefetch(addr); }

		template<> EIGEN_STRONG_INLINE std::complex<duals::dual<float> > pfirst<Packet1cdf>(const Packet1cdf& a)
		{
		EIGEN_ALIGN16 float res[4];
		vst1q_f32(res, a.v.v);
		return std::complex<duals::dual<float> >(duals::dual<float>(res[0],res[1]),
		duals::dual<float>(res[2],res[3]));
		}

		template<> EIGEN_STRONG_INLINE Packet1cdf preverse(const Packet1cdf& a) { return a; }

		template<> EIGEN_STRONG_INLINE std::complex<duals::dual<float> > predux<Packet1cdf>(const Packet1cdf& a)
		{
		return pfirst(a);
		}

		template<> EIGEN_STRONG_INLINE Packet1cdf preduxp<Packet1cdf>(const Packet1cdf* vecs)
		{
		return vecs[0];
		}

		template<> EIGEN_STRONG_INLINE std::complex<duals::dual<float> > predux_mul<Packet1cdf>(const Packet1cdf& a)
		{
		return pfirst(a);
		}

		template<int Offset>
		struct palign_impl<Offset,Packet1cdf>
		{
		static EIGEN_STRONG_INLINE void run(Packet1cdf& /first/, const Packet1cdf& /second/)
		{
		// FIXME is it sure we never have to align a Packet1cdf?
		// Even though a std::complex<duals::dual<float> > has 16 bytes, it is not necessarily aligned on a 16 bytes boundary...
		}
		};

		template<> struct conj_helper<Packet1cdf, Packet1cdf, false,true>
		{
		EIGEN_STRONG_INLINE Packet1cdf pmadd(const Packet1cdf& x, const Packet1cdf& y, const Packet1cdf& c) const
		{ return padd(pmul(x,y),c); }

		EIGEN_STRONG_INLINE Packet1cdf pmul(const Packet1cdf& a, const Packet1cdf& b) const
		{
		return internal::pmul(a, pconj(b));
		}
		};

		template<> struct conj_helper<Packet1cdf, Packet1cdf, true,false>
		{
		EIGEN_STRONG_INLINE Packet1cdf pmadd(const Packet1cdf& x, const Packet1cdf& y, const Packet1cdf& c) const
		{ return padd(pmul(x,y),c); }

		EIGEN_STRONG_INLINE Packet1cdf pmul(const Packet1cdf& a, const Packet1cdf& b) const
		{
		return internal::pmul(pconj(a), b);
		}
		};

		template<> struct conj_helper<Packet1cdf, Packet1cdf, true,true>
		{
		EIGEN_STRONG_INLINE Packet1cdf pmadd(const Packet1cdf& x, const Packet1cdf& y, const Packet1cdf& c) const
		{ return padd(pmul(x,y),c); }

		EIGEN_STRONG_INLINE Packet1cdf pmul(const Packet1cdf& a, const Packet1cdf& b) const
		{
		return pconj(internal::pmul(a, b));
		}
		};

		//TODO
		EIGEN_MAKE_CONJ_HELPER_CPLX_REAL(Packet1cdf,Packet2df)

		template<> EIGEN_STRONG_INLINE Packet1cdf pdiv<Packet1cdf>(const Packet1cdf& a, const Packet1cdf& b)
		{
		Packet1cdf res = conj_helper<Packet1cdf,Packet1cdf,false,true>().pmul(a, b);
		Packet2df s = pmul<Packet2df>(b.v, b.v);
		Packet2df rev_s = preverse<Packet2df>(s);

		return Packet1cdf(pdiv(res.v, padd<Packet2df>(s,rev_s)));
		}

		EIGEN_STRONG_INLINE Packet1cdf pcplxflip/* <Packet1cdf> */(const Packet1cdf& x)
		{
		return Packet1cdf(vcombine_u64(vget_high_u64(x.v.v),
		vget_low_u64(x.v.v)));
		//vec4f_swizzle1(x.v.v, 2, 3, 0, 1));
		}

		} // end namespace internal

		} // end namespace Eigen

		#endif // EIGEN_CDUAL_NEON_H

duals/arch/NEON/Dual.h

0 → 100644

+477 −0

File added.

Preview size limit exceeded, changes collapsed.

duals/dual

+35 −12

Original line number	Diff line number	Diff line
		@@ -584,7 +584,7 @@ public:
		template<class X, CPPDUALS_ONLY_SAME_DEPTH_AS_T(T,X)>
		dual<T> & operator-=(const dual<X> & x) { _real -= x.rpart(); _dual -= x.dpart(); return *this; }

		/// Multiply this dual with a dual of same of lower depth.
		/// Multiply this dual with a dual of same depth.
		template<class X, CPPDUALS_ONLY_SAME_DEPTH_AS_T(T,X)>
		dual<T> & operator*=(const dual<X> & x) {
		_dual = _real * x.dpart() + _dual * x.rpart();
		@@ -835,13 +835,14 @@ template<class T> dual<T> log2(const dual<T> & x) {

		template<class T, class U, CPPDUALS_ENABLE_SAME_DEPTH_AND_COMMON_T(T,U)>
		common_t
		pow(const dual<T> & x, const dual<U> & y) {
		pow(const dual<T> & f, const dual<U> & g) {
		using std::pow;
		using std::log;
		T v = pow(x.rpart(), y.rpart());
		T v = pow(f.rpart(), g.rpart());
		return common_t(v,
		x.dpart() * y.rpart() * pow(x.rpart(), y.rpart() - T(1)) +
		y.dpart() * log(x.rpart()) * v);
		pow(f.rpart(), g.rpart() - T(1)) *
		(g.rpart() * f.dpart()
		+ f.rpart() * log(f.rpart()) * g.dpart()));
		}

		template<class T, class U, CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,U)>
		@@ -877,6 +878,7 @@ template<class T> dual<T> fabs(const dual<T> & x) {
		}

		#if 0
		// TODO?
		template<class T> dual<T> abs2(const dual<T> & x) {
		using std::abs;
		return dual<T>(x.rpart() * x.rpart(),
		@@ -984,13 +986,34 @@ template<class T> dual<T> atan(const dual<T> & x) {
		return dual<T>(v, x.dpart() / (1 + x.rpart() * x.rpart()));
		}

		template<class T> dual<T> atan2(const dual<T> & x, const dual<T> & y) {
		template<class T> dual<T> atan2(const dual<T> & y, const dual<T> & x) {
		using std::atan2;
		T v = atan2(x.rpart(), y.rpart());
		if (x.dpart() == T(0))
		return v;
		else
		return dual<T>(v, x.dpart() / (1 + x.rpart()*x.rpart()));
		T v = atan2(y.rpart(), x.rpart());
		return dual<T>(v,
		(x.rpart() * y.dpart() - y.rpart() * x.dpart()) /
		(y.rpart() * y.rpart() + x.rpart() * x.rpart()));
		}

		// TODO
		template<class T, class U, CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,U)>
		common_t
		atan2(const dual<T> & y, const U & x) {
		using std::atan2;
		T v = atan2(y.rpart(), x);
		return dual<T>(v,
		(x * y.dpart()) /
		(y.rpart() * y.rpart() + x * x));
		}

		// TODO
		template<class T, class U, CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,U)>
		common_t
		atan2(const U & y, const dual<T> & x) {
		using std::atan2;
		T v = atan2(y, x.rpart());
		return dual<T>(v,
		(-y * x.dpart()) /
		(y * y + x.rpart() * x.rpart()));
		}

		// TODO

duals/dual_eigen

+2 −1

Original line number	Diff line number	Diff line
		@@ -384,7 +384,8 @@ template<typename RealScalar,bool Conj> struct dconj_helper<RealScalar, duals::d
		#elif defined(EIGEN_VECTORIZE_ALTIVEC) \|\| defined(EIGEN_VECTORIZE_VSX)
		// #include "duals/arch/AltiVec/Dual.h" // TODO
		#elif defined EIGEN_VECTORIZE_NEON
		//#include "duals/arch/NEON/Dual.h" // TODO
		#include "duals/arch/NEON/Dual.h"
		#include "duals/arch/NEON/ComplexDual.h"
		#elif defined EIGEN_VECTORIZE_ZVECTOR
		// #include "duals/arch/ZVector/Dual.h" // TODO
		#endif