diff --git a/wpiutil/src/main/native/include/Eigen/Core b/wpiutil/src/main/native/include/Eigen/Core index b923b8c000..b1aa7887ec 100644 --- a/wpiutil/src/main/native/include/Eigen/Core +++ b/wpiutil/src/main/native/include/Eigen/Core @@ -97,7 +97,7 @@ // this include file manages BLAS and MKL related macros // and inclusion of their respective header files -#include "src/Core/util/MKL_support.h" +// #include "src/Core/util/MKL_support.h" // if alignment is disabled, then disable vectorization. Note: EIGEN_MAX_ALIGN_BYTES is the proper check, it takes into // account both the user's will (EIGEN_MAX_ALIGN_BYTES,EIGEN_DONT_ALIGN) and our own platform checks @@ -408,9 +408,9 @@ using std::ptrdiff_t; #endif // Half float support -#include "src/Core/arch/CUDA/Half.h" -#include "src/Core/arch/CUDA/PacketMathHalf.h" -#include "src/Core/arch/CUDA/TypeCasting.h" +// #include "src/Core/arch/CUDA/Half.h" +// #include "src/Core/arch/CUDA/PacketMathHalf.h" +// #include "src/Core/arch/CUDA/TypeCasting.h" #if defined EIGEN_VECTORIZE_CUDA #include "src/Core/arch/CUDA/PacketMath.h" @@ -428,7 +428,7 @@ using std::ptrdiff_t; // Specialized functors to enable the processing of complex numbers // on CUDA devices -#include "src/Core/arch/CUDA/Complex.h" +// #include "src/Core/arch/CUDA/Complex.h" #include "src/Core/IO.h" #include "src/Core/DenseCoeffsBase.h" diff --git a/wpiutil/src/main/native/include/Eigen/src/Core/util/BlasUtil.h b/wpiutil/src/main/native/include/Eigen/src/Core/util/BlasUtil.h new file mode 100644 index 0000000000..6e6ee119b6 --- /dev/null +++ b/wpiutil/src/main/native/include/Eigen/src/Core/util/BlasUtil.h @@ -0,0 +1,398 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009-2010 Gael Guennebaud +// +// 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_BLASUTIL_H +#define EIGEN_BLASUTIL_H + +// This file contains many lightweight helper classes used to +// implement and control fast level 2 and level 3 BLAS-like routines. + +namespace Eigen { + +namespace internal { + +// forward declarations +template +struct gebp_kernel; + +template +struct gemm_pack_rhs; + +template +struct gemm_pack_lhs; + +template< + typename Index, + typename LhsScalar, int LhsStorageOrder, bool ConjugateLhs, + typename RhsScalar, int RhsStorageOrder, bool ConjugateRhs, + int ResStorageOrder> +struct general_matrix_matrix_product; + +template +struct general_matrix_vector_product; + + +template struct conj_if; + +template<> struct conj_if { + template + inline T operator()(const T& x) const { return numext::conj(x); } + template + inline T pconj(const T& x) const { return internal::pconj(x); } +}; + +template<> struct conj_if { + template + inline const T& operator()(const T& x) const { return x; } + template + inline const T& pconj(const T& x) const { return x; } +}; + +// Generic implementation for custom complex types. +template +struct conj_helper +{ + typedef typename ScalarBinaryOpTraits::ReturnType Scalar; + + EIGEN_STRONG_INLINE Scalar pmadd(const LhsScalar& x, const RhsScalar& y, const Scalar& c) const + { return padd(c, pmul(x,y)); } + + EIGEN_STRONG_INLINE Scalar pmul(const LhsScalar& x, const RhsScalar& y) const + { return conj_if()(x) * conj_if()(y); } +}; + +template struct conj_helper +{ + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const { return internal::pmadd(x,y,c); } + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const { return internal::pmul(x,y); } +}; + +template struct conj_helper, std::complex, false,true> +{ + typedef std::complex Scalar; + EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const + { return c + pmul(x,y); } + + EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const + { return Scalar(numext::real(x)*numext::real(y) + numext::imag(x)*numext::imag(y), numext::imag(x)*numext::real(y) - numext::real(x)*numext::imag(y)); } +}; + +template struct conj_helper, std::complex, true,false> +{ + typedef std::complex Scalar; + EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const + { return c + pmul(x,y); } + + EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const + { return Scalar(numext::real(x)*numext::real(y) + numext::imag(x)*numext::imag(y), numext::real(x)*numext::imag(y) - numext::imag(x)*numext::real(y)); } +}; + +template struct conj_helper, std::complex, true,true> +{ + typedef std::complex Scalar; + EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const Scalar& y, const Scalar& c) const + { return c + pmul(x,y); } + + EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const + { return Scalar(numext::real(x)*numext::real(y) - numext::imag(x)*numext::imag(y), - numext::real(x)*numext::imag(y) - numext::imag(x)*numext::real(y)); } +}; + +template struct conj_helper, RealScalar, Conj,false> +{ + typedef std::complex Scalar; + EIGEN_STRONG_INLINE Scalar pmadd(const Scalar& x, const RealScalar& y, const Scalar& c) const + { return padd(c, pmul(x,y)); } + EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const RealScalar& y) const + { return conj_if()(x)*y; } +}; + +template struct conj_helper, false,Conj> +{ + typedef std::complex Scalar; + EIGEN_STRONG_INLINE Scalar pmadd(const RealScalar& x, const Scalar& y, const Scalar& c) const + { return padd(c, pmul(x,y)); } + EIGEN_STRONG_INLINE Scalar pmul(const RealScalar& x, const Scalar& y) const + { return x*conj_if()(y); } +}; + +template struct get_factor { + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE To run(const From& x) { return To(x); } +}; + +template struct get_factor::Real> { + EIGEN_DEVICE_FUNC + static EIGEN_STRONG_INLINE typename NumTraits::Real run(const Scalar& x) { return numext::real(x); } +}; + + +template +class BlasVectorMapper { + public: + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE BlasVectorMapper(Scalar *data) : m_data(data) {} + + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Scalar operator()(Index i) const { + return m_data[i]; + } + template + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Packet load(Index i) const { + return ploadt(m_data + i); + } + + template + EIGEN_DEVICE_FUNC bool aligned(Index i) const { + return (UIntPtr(m_data+i)%sizeof(Packet))==0; + } + + protected: + Scalar* m_data; +}; + +template +class BlasLinearMapper { + public: + typedef typename packet_traits::type Packet; + typedef typename packet_traits::half HalfPacket; + + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE BlasLinearMapper(Scalar *data) : m_data(data) {} + + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE void prefetch(int i) const { + internal::prefetch(&operator()(i)); + } + + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Scalar& operator()(Index i) const { + return m_data[i]; + } + + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Packet loadPacket(Index i) const { + return ploadt(m_data + i); + } + + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE HalfPacket loadHalfPacket(Index i) const { + return ploadt(m_data + i); + } + + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE void storePacket(Index i, const Packet &p) const { + pstoret(m_data + i, p); + } + + protected: + Scalar *m_data; +}; + +// Lightweight helper class to access matrix coefficients. +template +class blas_data_mapper { + public: + typedef typename packet_traits::type Packet; + typedef typename packet_traits::half HalfPacket; + + typedef BlasLinearMapper LinearMapper; + typedef BlasVectorMapper VectorMapper; + + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE blas_data_mapper(Scalar* data, Index stride) : m_data(data), m_stride(stride) {} + + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE blas_data_mapper + getSubMapper(Index i, Index j) const { + return blas_data_mapper(&operator()(i, j), m_stride); + } + + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE LinearMapper getLinearMapper(Index i, Index j) const { + return LinearMapper(&operator()(i, j)); + } + + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE VectorMapper getVectorMapper(Index i, Index j) const { + return VectorMapper(&operator()(i, j)); + } + + + EIGEN_DEVICE_FUNC + EIGEN_ALWAYS_INLINE Scalar& operator()(Index i, Index j) const { + return m_data[StorageOrder==RowMajor ? j + i*m_stride : i + j*m_stride]; + } + + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Packet loadPacket(Index i, Index j) const { + return ploadt(&operator()(i, j)); + } + + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE HalfPacket loadHalfPacket(Index i, Index j) const { + return ploadt(&operator()(i, j)); + } + + template + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE void scatterPacket(Index i, Index j, const SubPacket &p) const { + pscatter(&operator()(i, j), p, m_stride); + } + + template + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE SubPacket gatherPacket(Index i, Index j) const { + return pgather(&operator()(i, j), m_stride); + } + + EIGEN_DEVICE_FUNC const Index stride() const { return m_stride; } + EIGEN_DEVICE_FUNC const Scalar* data() const { return m_data; } + + EIGEN_DEVICE_FUNC Index firstAligned(Index size) const { + if (UIntPtr(m_data)%sizeof(Scalar)) { + return -1; + } + return internal::first_default_aligned(m_data, size); + } + + protected: + Scalar* EIGEN_RESTRICT m_data; + const Index m_stride; +}; + +// lightweight helper class to access matrix coefficients (const version) +template +class const_blas_data_mapper : public blas_data_mapper { + public: + EIGEN_ALWAYS_INLINE const_blas_data_mapper(const Scalar *data, Index stride) : blas_data_mapper(data, stride) {} + + EIGEN_ALWAYS_INLINE const_blas_data_mapper getSubMapper(Index i, Index j) const { + return const_blas_data_mapper(&(this->operator()(i, j)), this->m_stride); + } +}; + + +/* Helper class to analyze the factors of a Product expression. + * In particular it allows to pop out operator-, scalar multiples, + * and conjugate */ +template struct blas_traits +{ + typedef typename traits::Scalar Scalar; + typedef const XprType& ExtractType; + typedef XprType _ExtractType; + enum { + IsComplex = NumTraits::IsComplex, + IsTransposed = false, + NeedToConjugate = false, + HasUsableDirectAccess = ( (int(XprType::Flags)&DirectAccessBit) + && ( bool(XprType::IsVectorAtCompileTime) + || int(inner_stride_at_compile_time::ret) == 1) + ) ? 1 : 0 + }; + typedef typename conditional::type DirectLinearAccessType; + static inline ExtractType extract(const XprType& x) { return x; } + static inline const Scalar extractScalarFactor(const XprType&) { return Scalar(1); } +}; + +// pop conjugate +template +struct blas_traits, NestedXpr> > + : blas_traits +{ + typedef blas_traits Base; + typedef CwiseUnaryOp, NestedXpr> XprType; + typedef typename Base::ExtractType ExtractType; + + enum { + IsComplex = NumTraits::IsComplex, + NeedToConjugate = Base::NeedToConjugate ? 0 : IsComplex + }; + static inline ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); } + static inline Scalar extractScalarFactor(const XprType& x) { return conj(Base::extractScalarFactor(x.nestedExpression())); } +}; + +// pop scalar multiple +template +struct blas_traits, const CwiseNullaryOp,Plain>, NestedXpr> > + : blas_traits +{ + typedef blas_traits Base; + typedef CwiseBinaryOp, const CwiseNullaryOp,Plain>, NestedXpr> XprType; + typedef typename Base::ExtractType ExtractType; + static inline ExtractType extract(const XprType& x) { return Base::extract(x.rhs()); } + static inline Scalar extractScalarFactor(const XprType& x) + { return x.lhs().functor().m_other * Base::extractScalarFactor(x.rhs()); } +}; +template +struct blas_traits, NestedXpr, const CwiseNullaryOp,Plain> > > + : blas_traits +{ + typedef blas_traits Base; + typedef CwiseBinaryOp, NestedXpr, const CwiseNullaryOp,Plain> > XprType; + typedef typename Base::ExtractType ExtractType; + static inline ExtractType extract(const XprType& x) { return Base::extract(x.lhs()); } + static inline Scalar extractScalarFactor(const XprType& x) + { return Base::extractScalarFactor(x.lhs()) * x.rhs().functor().m_other; } +}; +template +struct blas_traits, const CwiseNullaryOp,Plain1>, + const CwiseNullaryOp,Plain2> > > + : blas_traits,Plain1> > +{}; + +// pop opposite +template +struct blas_traits, NestedXpr> > + : blas_traits +{ + typedef blas_traits Base; + typedef CwiseUnaryOp, NestedXpr> XprType; + typedef typename Base::ExtractType ExtractType; + static inline ExtractType extract(const XprType& x) { return Base::extract(x.nestedExpression()); } + static inline Scalar extractScalarFactor(const XprType& x) + { return - Base::extractScalarFactor(x.nestedExpression()); } +}; + +// pop/push transpose +template +struct blas_traits > + : blas_traits +{ + typedef typename NestedXpr::Scalar Scalar; + typedef blas_traits Base; + typedef Transpose XprType; + typedef Transpose ExtractType; // const to get rid of a compile error; anyway blas traits are only used on the RHS + typedef Transpose _ExtractType; + typedef typename conditional::type DirectLinearAccessType; + enum { + IsTransposed = Base::IsTransposed ? 0 : 1 + }; + static inline ExtractType extract(const XprType& x) { return ExtractType(Base::extract(x.nestedExpression())); } + static inline Scalar extractScalarFactor(const XprType& x) { return Base::extractScalarFactor(x.nestedExpression()); } +}; + +template +struct blas_traits + : blas_traits +{}; + +template::HasUsableDirectAccess> +struct extract_data_selector { + static const typename T::Scalar* run(const T& m) + { + return blas_traits::extract(m).data(); + } +}; + +template +struct extract_data_selector { + static typename T::Scalar* run(const T&) { return 0; } +}; + +template const typename T::Scalar* extract_data(const T& m) +{ + return extract_data_selector::run(m); +} + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_BLASUTIL_H diff --git a/wpiutil/src/test/native/cpp/EigenTest.cpp b/wpiutil/src/test/native/cpp/EigenTest.cpp new file mode 100644 index 0000000000..04a0bc2543 --- /dev/null +++ b/wpiutil/src/test/native/cpp/EigenTest.cpp @@ -0,0 +1,61 @@ +/*----------------------------------------------------------------------------*/ +/* Copyright (c) 2019 FIRST. All Rights Reserved. */ +/* Open Source Software - may be modified and shared by FRC teams. The code */ +/* must be accompanied by the FIRST BSD license file in the root directory of */ +/* the project. */ +/*----------------------------------------------------------------------------*/ + +#include +#include + +#include "gtest/gtest.h" + +TEST(EigenTest, MultiplicationTest) { + Eigen::Matrix m1; + m1 << 2, 1, 0, 1; + + Eigen::Matrix m2; + m2 << 3, 0, 0, 2.5; + + const auto result = m1 * m2; + + Eigen::Matrix expectedResult; + expectedResult << 6.0, 2.5, 0.0, 2.5; + + EXPECT_TRUE(expectedResult.isApprox(result)); + + Eigen::Matrix m3; + m3 << 1.0, 3.0, 0.5, 2.0, 4.3, 1.2; + + Eigen::Matrix m4; + m4 << 3.0, 1.5, 2.0, 4.5, 2.3, 1.0, 1.6, 3.1, 5.2, 2.1, 2.0, 1.0; + + const auto result2 = m3 * m4; + + Eigen::Matrix expectedResult2; + expectedResult2 << 12.5, 5.55, 7.8, 14.3, 22.13, 9.82, 13.28, 23.53; + + EXPECT_TRUE(expectedResult2.isApprox(result2)); +} + +TEST(EigenTest, TransposeTest) { + Eigen::Matrix vec; + vec << 1, 2, 3; + + const auto transpose = vec.transpose(); + + Eigen::Matrix expectedTranspose; + expectedTranspose << 1, 2, 3; + + EXPECT_TRUE(expectedTranspose.isApprox(transpose)); +} + +TEST(EigenTest, InverseTest) { + Eigen::Matrix mat; + mat << 1.0, 3.0, 2.0, 5.0, 2.0, 1.5, 0.0, 1.3, 2.5; + + const auto inverse = mat.inverse(); + const auto identity = Eigen::MatrixXd::Identity(3, 3); + + EXPECT_TRUE(identity.isApprox(mat * inverse)); +}