Add support for dynamic dispatch of MMA instructions for POWER 10
diff --git a/Eigen/src/Core/arch/AltiVec/MatrixProduct.h b/Eigen/src/Core/arch/AltiVec/MatrixProduct.h index 57227e2..b863675 100644 --- a/Eigen/src/Core/arch/AltiVec/MatrixProduct.h +++ b/Eigen/src/Core/arch/AltiVec/MatrixProduct.h
@@ -10,6 +10,19 @@ #ifndef EIGEN_MATRIX_PRODUCT_ALTIVEC_H #define EIGEN_MATRIX_PRODUCT_ALTIVEC_H +#include "MatrixProductCommon.h" + +#if __GNUC__ > 10 || \ + (__GNUC__ == 10 && (__GNUC_MINOR__ > 2 || \ + (__GNUC_MINOR__ == 2 && \ + __GNUC_PATCHLEVEL__ >= 1))) + #define ALTIVEC_MMA_SUPPORT +#endif + +#if defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA) + #include "MatrixProductMMA.h" +#endif + /************************************************************************************************** * TODO * * - Check StorageOrder on lhs_pack (the innermost second loop seems unvectorized when it could). * @@ -26,18 +39,6 @@ **************************/ const int QuadRegisterCount = 8; -#ifdef __MMA__ - -template<typename Packet> -union Packetx2u -{ - __vector_pair vectorpair; - PacketBlock<Packet, 2> pair; -}; - -#endif - - template<typename Scalar> struct quad_traits { @@ -82,17 +83,6 @@ 12, 13, 14, 15, 20, 21, 22, 23, 28, 29, 30, 31}; - -const static Packet16uc p16uc_SETCOMPLEX32_FIRST = { 0, 1, 2, 3, - 16, 17, 18, 19, - 4, 5, 6, 7, - 20, 21, 22, 23}; - -const static Packet16uc p16uc_SETCOMPLEX32_SECOND = { 8, 9, 10, 11, - 24, 25, 26, 27, - 12, 13, 14, 15, - 28, 29, 30, 31}; -//[a,ai],[b,bi] = [a,b] const static Packet16uc p16uc_GETREAL64 = { 0, 1, 2, 3, 4, 5, 6, 7, 16, 17, 18, 19, 20, 21, 22, 23}; @@ -100,14 +90,6 @@ const static Packet16uc p16uc_GETIMAG64 = { 8, 9, 10, 11, 12, 13, 14, 15, 24, 25, 26, 27, 28, 29, 30, 31}; -//[a,b],[ai,bi] = [a,ai] - This is equivalent to p16uc_GETREAL64 -const static Packet16uc p16uc_SETCOMPLEX64_FIRST = { 0, 1, 2, 3, 4, 5, 6, 7, - 16, 17, 18, 19, 20, 21, 22, 23}; - -//[a,b],[ai,bi] = [b,bi] - This is equivalent to p16uc_GETIMAG64 -const static Packet16uc p16uc_SETCOMPLEX64_SECOND = { 8, 9, 10, 11, 12, 13, 14, 15, - 24, 25, 26, 27, 28, 29, 30, 31}; - /********************************************* * Single precision real and complex packing * * *******************************************/ @@ -1316,154 +1298,6 @@ * GEMM utils * **************/ -// Grab two decouples real/imaginary PacketBlocks and return two coupled (real/imaginary pairs) PacketBlocks. -template<typename Packet, typename Packetc> -EIGEN_STRONG_INLINE void bcouple(PacketBlock<Packet,4>& taccReal, PacketBlock<Packet,4>& taccImag, PacketBlock<Packetc,8>& tRes, PacketBlock<Packetc, 4>& acc1, PacketBlock<Packetc, 4>& acc2) -{ - acc1.packet[0].v = vec_perm(taccReal.packet[0], taccImag.packet[0], p16uc_SETCOMPLEX32_FIRST); - acc1.packet[1].v = vec_perm(taccReal.packet[1], taccImag.packet[1], p16uc_SETCOMPLEX32_FIRST); - acc1.packet[2].v = vec_perm(taccReal.packet[2], taccImag.packet[2], p16uc_SETCOMPLEX32_FIRST); - acc1.packet[3].v = vec_perm(taccReal.packet[3], taccImag.packet[3], p16uc_SETCOMPLEX32_FIRST); - - acc2.packet[0].v = vec_perm(taccReal.packet[0], taccImag.packet[0], p16uc_SETCOMPLEX32_SECOND); - acc2.packet[1].v = vec_perm(taccReal.packet[1], taccImag.packet[1], p16uc_SETCOMPLEX32_SECOND); - acc2.packet[2].v = vec_perm(taccReal.packet[2], taccImag.packet[2], p16uc_SETCOMPLEX32_SECOND); - acc2.packet[3].v = vec_perm(taccReal.packet[3], taccImag.packet[3], p16uc_SETCOMPLEX32_SECOND); - - acc1.packet[0] = padd<Packetc>(tRes.packet[0], acc1.packet[0]); - acc1.packet[1] = padd<Packetc>(tRes.packet[1], acc1.packet[1]); - acc1.packet[2] = padd<Packetc>(tRes.packet[2], acc1.packet[2]); - acc1.packet[3] = padd<Packetc>(tRes.packet[3], acc1.packet[3]); - - acc2.packet[0] = padd<Packetc>(tRes.packet[4], acc2.packet[0]); - acc2.packet[1] = padd<Packetc>(tRes.packet[5], acc2.packet[1]); - acc2.packet[2] = padd<Packetc>(tRes.packet[6], acc2.packet[2]); - acc2.packet[3] = padd<Packetc>(tRes.packet[7], acc2.packet[3]); -} - -template<> -EIGEN_STRONG_INLINE void bcouple<Packet2d, Packet1cd>(PacketBlock<Packet2d,4>& taccReal, PacketBlock<Packet2d,4>& taccImag, PacketBlock<Packet1cd,8>& tRes, PacketBlock<Packet1cd, 4>& acc1, PacketBlock<Packet1cd, 4>& acc2) -{ - acc1.packet[0].v = vec_perm(taccReal.packet[0], taccImag.packet[0], p16uc_SETCOMPLEX64_FIRST); - acc1.packet[1].v = vec_perm(taccReal.packet[1], taccImag.packet[1], p16uc_SETCOMPLEX64_FIRST); - acc1.packet[2].v = vec_perm(taccReal.packet[2], taccImag.packet[2], p16uc_SETCOMPLEX64_FIRST); - acc1.packet[3].v = vec_perm(taccReal.packet[3], taccImag.packet[3], p16uc_SETCOMPLEX64_FIRST); - - acc2.packet[0].v = vec_perm(taccReal.packet[0], taccImag.packet[0], p16uc_SETCOMPLEX64_SECOND); - acc2.packet[1].v = vec_perm(taccReal.packet[1], taccImag.packet[1], p16uc_SETCOMPLEX64_SECOND); - acc2.packet[2].v = vec_perm(taccReal.packet[2], taccImag.packet[2], p16uc_SETCOMPLEX64_SECOND); - acc2.packet[3].v = vec_perm(taccReal.packet[3], taccImag.packet[3], p16uc_SETCOMPLEX64_SECOND); - - acc1.packet[0] = padd<Packet1cd>(tRes.packet[0], acc1.packet[0]); - acc1.packet[1] = padd<Packet1cd>(tRes.packet[1], acc1.packet[1]); - acc1.packet[2] = padd<Packet1cd>(tRes.packet[2], acc1.packet[2]); - acc1.packet[3] = padd<Packet1cd>(tRes.packet[3], acc1.packet[3]); - - acc2.packet[0] = padd<Packet1cd>(tRes.packet[4], acc2.packet[0]); - acc2.packet[1] = padd<Packet1cd>(tRes.packet[5], acc2.packet[1]); - acc2.packet[2] = padd<Packet1cd>(tRes.packet[6], acc2.packet[2]); - acc2.packet[3] = padd<Packet1cd>(tRes.packet[7], acc2.packet[3]); -} - -#ifdef __MMA__ -template<typename Packet> -EIGEN_STRONG_INLINE PacketBlock<Packet,2> pmul (const PacketBlock<Packet,2>& a, const Packet& b) -{ - PacketBlock<Packet,2> pb; - pb.packet[0] = a.packet[0]*b; - pb.packet[1] = a.packet[1]*b; - return pb; -} -template<typename DataMapper, typename Index, typename Packet> -EIGEN_STRONG_INLINE void storeAccumulator(Index i, Index j, const DataMapper& data, const Packet& alpha, __vector_quad *acc) -{ - PacketBlock<Packet, 4> result; - __builtin_mma_disassemble_acc(&result.packet, acc); - - PacketBlock<Packet, 4> block; - block.packet[0] = data.template loadPacket<Packet>(i, j + 0) + pmul<Packet>(alpha, result.packet[0]); - block.packet[1] = data.template loadPacket<Packet>(i, j + 1) + pmul<Packet>(alpha, result.packet[1]); - block.packet[2] = data.template loadPacket<Packet>(i, j + 2) + pmul<Packet>(alpha, result.packet[2]); - block.packet[3] = data.template loadPacket<Packet>(i, j + 3) + pmul<Packet>(alpha, result.packet[3]); - - data.template storePacketBlock<Packet, 4>(i, j, block); -} - -template<typename DataMapper, typename Index, typename Packet, typename Packetc, int N> -EIGEN_STRONG_INLINE void storeComplexAccumulator(Index i, Index j, const DataMapper& data, const Packet& alphaReal, const Packet& alphaImag, __vector_quad *accReal, __vector_quad *accImag, const int accColsC) -{ - PacketBlock<Packet, 4> resultReal, resultImag; - __builtin_mma_disassemble_acc(&resultReal.packet, accReal); - __builtin_mma_disassemble_acc(&resultImag.packet, accImag); - - PacketBlock<Packet,4> taccReal, taccImag; - taccReal.packet[0] = pmul<Packet>(resultReal.packet[0], alphaReal); - taccReal.packet[1] = pmul<Packet>(resultReal.packet[1], alphaReal); - taccReal.packet[2] = pmul<Packet>(resultReal.packet[2], alphaReal); - taccReal.packet[3] = pmul<Packet>(resultReal.packet[3], alphaReal); - - taccImag.packet[0] = pmul<Packet>(resultImag.packet[0], alphaReal); - taccImag.packet[1] = pmul<Packet>(resultImag.packet[1], alphaReal); - taccImag.packet[2] = pmul<Packet>(resultImag.packet[2], alphaReal); - taccImag.packet[3] = pmul<Packet>(resultImag.packet[3], alphaReal); - - taccReal.packet[0] = psub<Packet>(taccReal.packet[0], pmul<Packet>(resultImag.packet[0], alphaImag)); - taccReal.packet[1] = psub<Packet>(taccReal.packet[1], pmul<Packet>(resultImag.packet[1], alphaImag)); - taccReal.packet[2] = psub<Packet>(taccReal.packet[2], pmul<Packet>(resultImag.packet[2], alphaImag)); - taccReal.packet[3] = psub<Packet>(taccReal.packet[3], pmul<Packet>(resultImag.packet[3], alphaImag)); - - taccImag.packet[0] = pmadd<Packet>(resultReal.packet[0], alphaImag, taccImag.packet[0]); - taccImag.packet[1] = pmadd<Packet>(resultReal.packet[1], alphaImag, taccImag.packet[1]); - taccImag.packet[2] = pmadd<Packet>(resultReal.packet[2], alphaImag, taccImag.packet[2]); - taccImag.packet[3] = pmadd<Packet>(resultReal.packet[3], alphaImag, taccImag.packet[3]); - - PacketBlock<Packetc, 8> tRes; - tRes.packet[0] = data.template loadPacket<Packetc>(i + N*accColsC, j + 0); - tRes.packet[1] = data.template loadPacket<Packetc>(i + N*accColsC, j + 1); - tRes.packet[2] = data.template loadPacket<Packetc>(i + N*accColsC, j + 2); - tRes.packet[3] = data.template loadPacket<Packetc>(i + N*accColsC, j + 3); - - tRes.packet[4] = data.template loadPacket<Packetc>(i + (N+1)*accColsC, j + 0); - tRes.packet[5] = data.template loadPacket<Packetc>(i + (N+1)*accColsC, j + 1); - tRes.packet[6] = data.template loadPacket<Packetc>(i + (N+1)*accColsC, j + 2); - tRes.packet[7] = data.template loadPacket<Packetc>(i + (N+1)*accColsC, j + 3); - - PacketBlock<Packetc, 4> acc1, acc2; - bcouple<Packet, Packetc>(taccReal, taccImag, tRes, acc1, acc2); - - data.template storePacketBlock<Packetc, 4>(i + N*accColsC, j, acc1); - data.template storePacketBlock<Packetc, 4>(i + (N+1)*accColsC, j, acc2); -} - -// Defaults to float32, since Eigen still supports C++03 we can't use default template arguments -template<typename LhsPacket, typename RhsPacket, bool NegativeAccumulate> -EIGEN_STRONG_INLINE void pger(__vector_quad *acc, const RhsPacket& a, const LhsPacket& b) -{ - if(NegativeAccumulate) - { - __builtin_mma_xvf32gernp(acc, (__vector unsigned char)a, (__vector unsigned char)b); - } else { - __builtin_mma_xvf32gerpp(acc, (__vector unsigned char)a, (__vector unsigned char)b); - } -} - -template<> -EIGEN_STRONG_INLINE void pger<Packet2d, PacketBlock<Packet2d, 2>, false>(__vector_quad *acc, const PacketBlock<Packet2d,2>& a, const Packet2d& b) -{ - Packetx2u<Packet2d> p; - p.pair = a; - __builtin_mma_xvf64gerpp(acc, p.vectorpair, (__vector unsigned char)b); -} - -template<> -EIGEN_STRONG_INLINE void pger<Packet2d, PacketBlock<Packet2d, 2>, true>(__vector_quad *acc, const PacketBlock<Packet2d, 2>& a, const Packet2d& b) -{ - Packetx2u<Packet2d> p; - p.pair = a; - __builtin_mma_xvf64gernp(acc, p.vectorpair, (__vector unsigned char)b); -} -#else - // 512-bits rank1-update of acc. It can either positive or negative accumulate (useful for complex gemm). template<typename Scalar, typename Packet, bool NegativeAccumulate> EIGEN_STRONG_INLINE void pger(PacketBlock<Packet, 4> *acc, const Scalar* lhs, const Scalar* rhs) @@ -1561,25 +1395,6 @@ accImag.packet[3] = pmadd<Packet>(rhsV4, lhsVi, accImag.packet[3]); } } -#endif - -// This is necessary because ploadRhs for double returns a pair of vectors when MMA is enabled. -template<typename Scalar, typename Packet> -EIGEN_STRONG_INLINE Packet ploadRhs(const Scalar *rhs) -{ - return *((Packet *)rhs); -} - -#ifdef __MMA__ -template<> -EIGEN_STRONG_INLINE PacketBlock<Packet2d, 2> ploadRhs<double, PacketBlock<Packet2d, 2> >(const double *rhs) -{ - PacketBlock<Packet2d, 2> pair; - pair.packet[0] = *((Packet2d *)rhs ); - pair.packet[1] = *(((Packet2d *)rhs) + 1); - return pair; -} -#endif template<typename Scalar, typename Packet> EIGEN_STRONG_INLINE Packet ploadLhs(const Scalar *lhs) @@ -1587,7 +1402,6 @@ return *((Packet *)lhs); } -#ifndef __MMA__ // Zero the accumulator on PacketBlock. template<typename Scalar, typename Packet> EIGEN_STRONG_INLINE void bsetzero(PacketBlock<Packet,4>& acc) @@ -1656,7 +1470,7 @@ acc.packet[6] = res.template loadPacket<Packet>(row + (N+1)*accCols, col + 2); acc.packet[7] = res.template loadPacket<Packet>(row + (N+1)*accCols, col + 3); } -#endif + // PEEL loop factor. #define PEEL 10 @@ -1682,31 +1496,6 @@ const Scalar *lhs_base = blockA; Index row = 0; -#ifdef __MMA__ - for(; row + accCols <= rows; row += accCols) - { - const Scalar *rhs_ptr = rhs_base; - const Scalar *lhs_ptr1 = lhs_base + (row/accCols)*strideA*accCols; - - __vector_quad acc; - __builtin_mma_xxsetaccz(&acc); - - lhs_ptr1 += accCols*offsetA; - rhs_ptr += accRows*offsetB; - for(Index k = 0; k < depth; k++) - { - Packet lhsV = ploadLhs<Scalar, Packet>(lhs_ptr1); - RhsPacket rhsV = ploadRhs<Scalar, RhsPacket>(rhs_ptr); - - pger<Packet, RhsPacket, false>(&acc, rhsV, lhsV); - - lhs_ptr1 += accCols; - rhs_ptr += accRows; - } - - storeAccumulator<DataMapper, Index, Packet>(row, col, res, pAlpha, &acc); - } -#else for(; row + 6*accCols <= rows; row += 6*accCols) { #define MICRO() \ @@ -2135,7 +1924,6 @@ res.template storePacketBlock<Packet, 4>(row, col, acc1); #undef MICRO } -#endif if(remaining_rows > 0) { const Scalar *rhs_ptr = rhs_base; @@ -2239,60 +2027,6 @@ const Scalar *lhs_base = blockA; Index row = 0; -#ifdef __MMA__ - for(; row + accCols <= rows; row += accCols) - { - const Scalar *rhs_ptr = rhs_base; - const Scalar *rhs_ptr_imag = rhs_ptr + accRows*strideB; - const Scalar *lhs_ptr = lhs_base + ((advanceRows*row)/accCols)*strideA*accCols; - const Scalar *lhs_ptr_imag = lhs_ptr + accCols*strideA; - - __vector_quad accReal, accImag; - __builtin_mma_xxsetaccz(&accReal); - __builtin_mma_xxsetaccz(&accImag); - - lhs_ptr += accCols*offsetA; - if(!LhsIsReal) - lhs_ptr_imag += accCols*offsetA; - rhs_ptr += accRows*offsetB; - if(!RhsIsReal) - rhs_ptr_imag += accRows*offsetB; - for(Index k = 0; k < depth; k++) - { - Packet lhsV = ploadLhs<Scalar, Packet>(lhs_ptr); - RhsPacket rhsV = ploadRhs<Scalar, RhsPacket>(rhs_ptr); - - Packet lhsVi = ploadLhs<Scalar, Packet>(lhs_ptr_imag); - RhsPacket rhsVi = ploadRhs<Scalar, RhsPacket>(rhs_ptr_imag); - - if(ConjugateLhs && !LhsIsReal) lhsVi = pmul<Packet>(lhsVi, conj); - if(ConjugateRhs && !RhsIsReal) rhsVi = pmul<Packet>(rhsVi, conj); - - if(LhsIsReal) - { - pger<Packet, RhsPacket, false>(&accReal, rhsV, lhsV); - pger<Packet, RhsPacket, false>(&accImag, rhsVi, lhsV); - } else if(RhsIsReal) { - pger<Packet, RhsPacket, false>(&accReal, rhsV, lhsV); - pger<Packet, RhsPacket, false>(&accImag, rhsV, lhsVi); - } else { - pger<Packet, RhsPacket, false>(&accReal, rhsV, lhsV); - pger<Packet, RhsPacket, true>(&accReal, rhsVi, lhsVi); - pger<Packet, RhsPacket, false>(&accImag, rhsVi, lhsV); - pger<Packet, RhsPacket, false>(&accImag, rhsV, lhsVi); - } - - lhs_ptr += accCols; - rhs_ptr += accRows; - if(!LhsIsReal) - lhs_ptr_imag += accCols; - if(!RhsIsReal) - rhs_ptr_imag += accRows; - } - - storeComplexAccumulator<DataMapper, Index, Packet, Packetc, 0>(row, col, res, pAlphaReal, pAlphaImag, &accReal, &accImag, accColsC); - } -#else for(; row + accCols <= rows; row += accCols) { #define MICRO() \ @@ -2302,7 +2036,7 @@ if(!LhsIsReal) \ lhs_ptr_imag1 += accCols; \ if(!RhsIsReal) \ - rhs_ptr_imag += accRows; + rhs_ptr_imag += accRows; const Scalar *rhs_ptr = rhs_base; const Scalar *rhs_ptr_imag = rhs_ptr + accRows*strideB; @@ -2356,7 +2090,6 @@ res.template storePacketBlock<Packetc, 4>(row + accColsC, col, acc2); #undef MICRO } -#endif if(remaining_rows > 0) { const Scalar *rhs_ptr = rhs_base; @@ -2383,7 +2116,7 @@ lhsc.real(lhs_real); if(!LhsIsReal) { - if(ConjugateLhs) + if(ConjugateLhs) lhsc.imag(-lhs_imag); else lhsc.imag(lhs_imag); @@ -2457,7 +2190,7 @@ { Scalar lhs_real = lhs_ptr[arow]; Scalar lhs_imag; - if(!LhsIsReal) + if(!LhsIsReal) { lhs_imag = lhs_ptr_imag[arow]; @@ -2534,7 +2267,7 @@ lhsc.real(lhs_real); if(!LhsIsReal) { - if(ConjugateLhs) + if(ConjugateLhs) lhsc.imag(-lhs_imag); else lhsc.imag(lhs_imag); @@ -2819,8 +2552,22 @@ { const int accRows = quad_traits<float>::rows; const int accCols = quad_traits<float>::size; + void (*gemm_function)(const DataMapper&, const float*, const float*, Index, Index, Index, float, Index, Index, Index, Index, const int, const int); - gemm<float, Index, Packet, RhsPacket, DataMapper>(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols); + #ifdef EIGEN_ALTIVEC_MMA_ONLY + //generate with MMA only + gemm_function = &Eigen::internal::gemmMMA<float, Index, Packet, RhsPacket, DataMapper>; + #elif defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA) + if (__builtin_cpu_supports ("arch_3_1") && __builtin_cpu_supports ("mma")){ + gemm_function = &Eigen::internal::gemmMMA<float, Index, Packet, RhsPacket, DataMapper>; + } + else{ + gemm_function = &Eigen::internal::gemm<float, Index, Packet, RhsPacket, DataMapper>; + } + #else + gemm_function = &Eigen::internal::gemm<float, Index, Packet, RhsPacket, DataMapper>; + #endif + gemm_function(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols); } template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs> @@ -2843,8 +2590,23 @@ { const int accRows = quad_traits<float>::rows; const int accCols = quad_traits<float>::size; + void (*gemm_function)(const DataMapper&, const std::complex<float>*, const std::complex<float>*, + Index, Index, Index, std::complex<float>, Index, Index , Index, Index, const int, const int); - gemm_complex<std::complex<float>, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols); + #ifdef EIGEN_ALTIVEC_MMA_ONLY + //generate with MMA only + gemm_function = &Eigen::internal::gemm_complexMMA<std::complex<float>, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>; + #elif defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA) + if (__builtin_cpu_supports ("arch_3_1") && __builtin_cpu_supports ("mma")){ + gemm_function = &Eigen::internal::gemm_complexMMA<std::complex<float>, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>; + } + else{ + gemm_function = &Eigen::internal::gemm_complex<std::complex<float>, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>; + } + #else + gemm_function = &Eigen::internal::gemm_complex<std::complex<float>, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>; + #endif + gemm_function(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols); } template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs> @@ -2867,8 +2629,22 @@ { const int accRows = quad_traits<float>::rows; const int accCols = quad_traits<float>::size; - - gemm_complex<float, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols); + void (*gemm_function)(const DataMapper&, const float*, const std::complex<float>*, + Index, Index, Index, std::complex<float>, Index, Index , Index, Index, const int, const int); + #ifdef EIGEN_ALTIVEC_MMA_ONLY + //generate with MMA only + gemm_function = &Eigen::internal::gemm_complexMMA<float, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>; + #elif defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA) + if (__builtin_cpu_supports ("arch_3_1") && __builtin_cpu_supports ("mma")){ + gemm_function = &Eigen::internal::gemm_complexMMA<float, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>; + } + else{ + gemm_function = &Eigen::internal::gemm_complex<float, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>; + } + #else + gemm_function = &Eigen::internal::gemm_complex<float, std::complex<float>, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>; + #endif + gemm_function(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols); } template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs> @@ -2891,8 +2667,22 @@ { const int accRows = quad_traits<float>::rows; const int accCols = quad_traits<float>::size; - - gemm_complex<std::complex<float>, float, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols); + void (*gemm_function)(const DataMapper&, const std::complex<float>*, const float*, + Index, Index, Index, std::complex<float>, Index, Index , Index, Index, const int, const int); + #ifdef EIGEN_ALTIVEC_MMA_ONLY + //generate with MMA only + gemm_function = &Eigen::internal::gemm_complexMMA<std::complex<float>, float, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>; + #elif defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA) + if (__builtin_cpu_supports ("arch_3_1") && __builtin_cpu_supports ("mma")){ + gemm_function = &Eigen::internal::gemm_complexMMA<std::complex<float>, float, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>; + } + else{ + gemm_function = &Eigen::internal::gemm_complex<std::complex<float>, float, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>; + } + #else + gemm_function = &Eigen::internal::gemm_complex<std::complex<float>, float, std::complex<float>, float, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>; + #endif + gemm_function(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols); } template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs> @@ -2914,8 +2704,22 @@ { const int accRows = quad_traits<double>::rows; const int accCols = quad_traits<double>::size; + void (*gemm_function)(const DataMapper&, const double*, const double*, Index, Index, Index, double, Index, Index, Index, Index, const int, const int); - gemm<double, Index, Packet, RhsPacket, DataMapper>(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols); + #ifdef EIGEN_ALTIVEC_MMA_ONLY + //generate with MMA only + gemm_function = &Eigen::internal::gemmMMA<double, Index, Packet, RhsPacket, DataMapper>; + #elif defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA) + if (__builtin_cpu_supports ("arch_3_1") && __builtin_cpu_supports ("mma")){ + gemm_function = &Eigen::internal::gemmMMA<double, Index, Packet, RhsPacket, DataMapper>; + } + else{ + gemm_function = &Eigen::internal::gemm<double, Index, Packet, RhsPacket, DataMapper>; + } + #else + gemm_function = &Eigen::internal::gemm<double, Index, Packet, RhsPacket, DataMapper>; + #endif + gemm_function(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols); } template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs> @@ -2938,8 +2742,22 @@ { const int accRows = quad_traits<double>::rows; const int accCols = quad_traits<double>::size; - - gemm_complex<std::complex<double>, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols); + void (*gemm_function)(const DataMapper&, const std::complex<double>*, const std::complex<double>*, + Index, Index, Index, std::complex<double>, Index, Index , Index, Index, const int, const int); + #ifdef EIGEN_ALTIVEC_MMA_ONLY + //generate with MMA only + gemm_function = &Eigen::internal::gemm_complexMMA<std::complex<double>, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>; + #elif defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA) + if (__builtin_cpu_supports ("arch_3_1") && __builtin_cpu_supports ("mma")){ + gemm_function = &Eigen::internal::gemm_complexMMA<std::complex<double>, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>; + } + else{ + gemm_function = &Eigen::internal::gemm_complex<std::complex<double>, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>; + } + #else + gemm_function = &Eigen::internal::gemm_complex<std::complex<double>, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, false>; + #endif + gemm_function(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols); } template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs> @@ -2962,8 +2780,22 @@ { const int accRows = quad_traits<double>::rows; const int accCols = quad_traits<double>::size; - - gemm_complex<std::complex<double>, double, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols); + void (*gemm_function)(const DataMapper&, const std::complex<double>*, const double*, + Index, Index, Index, std::complex<double>, Index, Index , Index, Index, const int, const int); + #ifdef EIGEN_ALTIVEC_MMA_ONLY + //generate with MMA only + gemm_function = &Eigen::internal::gemm_complexMMA<std::complex<double>, double, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>; + #elif defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA) + if (__builtin_cpu_supports ("arch_3_1") && __builtin_cpu_supports ("mma")){ + gemm_function = &Eigen::internal::gemm_complexMMA<std::complex<double>, double, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>; + } + else{ + gemm_function = &Eigen::internal::gemm_complex<std::complex<double>, double, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>; + } + #else + gemm_function = &Eigen::internal::gemm_complex<std::complex<double>, double, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, false, true>; + #endif + gemm_function(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols); } template<typename Index, typename DataMapper, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs> @@ -2986,10 +2818,25 @@ { const int accRows = quad_traits<double>::rows; const int accCols = quad_traits<double>::size; - - gemm_complex<double, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols); + void (*gemm_function)(const DataMapper&, const double*, const std::complex<double>*, + Index, Index, Index, std::complex<double>, Index, Index , Index, Index, const int, const int); + #ifdef EIGEN_ALTIVEC_MMA_ONLY + //generate with MMA only + gemm_function = &Eigen::internal::gemm_complexMMA<double, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>; + #elif defined(ALTIVEC_MMA_SUPPORT) && !defined(EIGEN_ALTIVEC_DISABLE_MMA) + if (__builtin_cpu_supports ("arch_3_1") && __builtin_cpu_supports ("mma")){ + gemm_function = &Eigen::internal::gemm_complexMMA<double, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>; + } + else{ + gemm_function = &Eigen::internal::gemm_complex<double, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>; + } + #else + gemm_function = &Eigen::internal::gemm_complex<double, std::complex<double>, std::complex<double>, double, Index, Packet, Packetc, RhsPacket, DataMapper, ConjugateLhs, ConjugateRhs, true, false>; + #endif + gemm_function(res, blockA, blockB, rows, depth, cols, alpha, strideA, strideB, offsetA, offsetB, accRows, accCols); } } // end namespace internal } // end namespace Eigen -#endif // EIGEN_MATRIX_PRODUCT_ALTIVEC_H \ No newline at end of file + +#endif // EIGEN_MATRIX_PRODUCT_ALTIVEC_H
diff --git a/Eigen/src/Core/arch/AltiVec/MatrixProductCommon.h b/Eigen/src/Core/arch/AltiVec/MatrixProductCommon.h new file mode 100644 index 0000000..87b60c2 --- /dev/null +++ b/Eigen/src/Core/arch/AltiVec/MatrixProductCommon.h
@@ -0,0 +1,80 @@ +namespace Eigen { + +namespace internal { + +const static Packet16uc p16uc_SETCOMPLEX32_FIRST = { 0, 1, 2, 3, + 16, 17, 18, 19, + 4, 5, 6, 7, + 20, 21, 22, 23}; + +const static Packet16uc p16uc_SETCOMPLEX32_SECOND = { 8, 9, 10, 11, + 24, 25, 26, 27, + 12, 13, 14, 15, + 28, 29, 30, 31}; +//[a,b],[ai,bi] = [a,ai] - This is equivalent to p16uc_GETREAL64 +const static Packet16uc p16uc_SETCOMPLEX64_FIRST = { 0, 1, 2, 3, 4, 5, 6, 7, + 16, 17, 18, 19, 20, 21, 22, 23}; + +//[a,b],[ai,bi] = [b,bi] - This is equivalent to p16uc_GETIMAG64 +const static Packet16uc p16uc_SETCOMPLEX64_SECOND = { 8, 9, 10, 11, 12, 13, 14, 15, + 24, 25, 26, 27, 28, 29, 30, 31}; + + +// Grab two decouples real/imaginary PacketBlocks and return two coupled (real/imaginary pairs) PacketBlocks. +template<typename Packet, typename Packetc> +EIGEN_STRONG_INLINE void bcouple(PacketBlock<Packet,4>& taccReal, PacketBlock<Packet,4>& taccImag, PacketBlock<Packetc,8>& tRes, PacketBlock<Packetc, 4>& acc1, PacketBlock<Packetc, 4>& acc2) +{ + acc1.packet[0].v = vec_perm(taccReal.packet[0], taccImag.packet[0], p16uc_SETCOMPLEX32_FIRST); + acc1.packet[1].v = vec_perm(taccReal.packet[1], taccImag.packet[1], p16uc_SETCOMPLEX32_FIRST); + acc1.packet[2].v = vec_perm(taccReal.packet[2], taccImag.packet[2], p16uc_SETCOMPLEX32_FIRST); + acc1.packet[3].v = vec_perm(taccReal.packet[3], taccImag.packet[3], p16uc_SETCOMPLEX32_FIRST); + + acc2.packet[0].v = vec_perm(taccReal.packet[0], taccImag.packet[0], p16uc_SETCOMPLEX32_SECOND); + acc2.packet[1].v = vec_perm(taccReal.packet[1], taccImag.packet[1], p16uc_SETCOMPLEX32_SECOND); + acc2.packet[2].v = vec_perm(taccReal.packet[2], taccImag.packet[2], p16uc_SETCOMPLEX32_SECOND); + acc2.packet[3].v = vec_perm(taccReal.packet[3], taccImag.packet[3], p16uc_SETCOMPLEX32_SECOND); + + acc1.packet[0] = padd<Packetc>(tRes.packet[0], acc1.packet[0]); + acc1.packet[1] = padd<Packetc>(tRes.packet[1], acc1.packet[1]); + acc1.packet[2] = padd<Packetc>(tRes.packet[2], acc1.packet[2]); + acc1.packet[3] = padd<Packetc>(tRes.packet[3], acc1.packet[3]); + + acc2.packet[0] = padd<Packetc>(tRes.packet[4], acc2.packet[0]); + acc2.packet[1] = padd<Packetc>(tRes.packet[5], acc2.packet[1]); + acc2.packet[2] = padd<Packetc>(tRes.packet[6], acc2.packet[2]); + acc2.packet[3] = padd<Packetc>(tRes.packet[7], acc2.packet[3]); +} + +template<> +EIGEN_STRONG_INLINE void bcouple<Packet2d, Packet1cd>(PacketBlock<Packet2d,4>& taccReal, PacketBlock<Packet2d,4>& taccImag, PacketBlock<Packet1cd,8>& tRes, PacketBlock<Packet1cd, 4>& acc1, PacketBlock<Packet1cd, 4>& acc2) +{ + acc1.packet[0].v = vec_perm(taccReal.packet[0], taccImag.packet[0], p16uc_SETCOMPLEX64_FIRST); + acc1.packet[1].v = vec_perm(taccReal.packet[1], taccImag.packet[1], p16uc_SETCOMPLEX64_FIRST); + acc1.packet[2].v = vec_perm(taccReal.packet[2], taccImag.packet[2], p16uc_SETCOMPLEX64_FIRST); + acc1.packet[3].v = vec_perm(taccReal.packet[3], taccImag.packet[3], p16uc_SETCOMPLEX64_FIRST); + + acc2.packet[0].v = vec_perm(taccReal.packet[0], taccImag.packet[0], p16uc_SETCOMPLEX64_SECOND); + acc2.packet[1].v = vec_perm(taccReal.packet[1], taccImag.packet[1], p16uc_SETCOMPLEX64_SECOND); + acc2.packet[2].v = vec_perm(taccReal.packet[2], taccImag.packet[2], p16uc_SETCOMPLEX64_SECOND); + acc2.packet[3].v = vec_perm(taccReal.packet[3], taccImag.packet[3], p16uc_SETCOMPLEX64_SECOND); + + acc1.packet[0] = padd<Packet1cd>(tRes.packet[0], acc1.packet[0]); + acc1.packet[1] = padd<Packet1cd>(tRes.packet[1], acc1.packet[1]); + acc1.packet[2] = padd<Packet1cd>(tRes.packet[2], acc1.packet[2]); + acc1.packet[3] = padd<Packet1cd>(tRes.packet[3], acc1.packet[3]); + + acc2.packet[0] = padd<Packet1cd>(tRes.packet[4], acc2.packet[0]); + acc2.packet[1] = padd<Packet1cd>(tRes.packet[5], acc2.packet[1]); + acc2.packet[2] = padd<Packet1cd>(tRes.packet[6], acc2.packet[2]); + acc2.packet[3] = padd<Packet1cd>(tRes.packet[7], acc2.packet[3]); +} + +// This is necessary because ploadRhs for double returns a pair of vectors when MMA is enabled. +template<typename Scalar, typename Packet> +EIGEN_STRONG_INLINE Packet ploadRhs(const Scalar *rhs) +{ + return *((Packet *)rhs); +} + +} // end namespace internal +} // end namespace Eigen
diff --git a/Eigen/src/Core/arch/AltiVec/MatrixProductMMA.h b/Eigen/src/Core/arch/AltiVec/MatrixProductMMA.h new file mode 100644 index 0000000..1866a71 --- /dev/null +++ b/Eigen/src/Core/arch/AltiVec/MatrixProductMMA.h
@@ -0,0 +1,638 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2020 Everton Constantino (everton.constantino@ibm.com) +// +// 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_MATRIX_PRODUCT_MMA_ALTIVEC_H +#define EIGEN_MATRIX_PRODUCT_MMA_ALTIVEC_H + +#pragma GCC target("cpu=power10") + +namespace Eigen { + +namespace internal { + +template<typename Packet> +union Packetx2u +{ + __vector_pair vectorpair; + PacketBlock<Packet, 2> pair; +}; +const static Packet16uc MMA_p16uc_SETCOMPLEX32_FIRST = { 0, 1, 2, 3, + 16, 17, 18, 19, + 4, 5, 6, 7, + 20, 21, 22, 23}; + +const static Packet16uc MMA_p16uc_SETCOMPLEX32_SECOND = { 8, 9, 10, 11, + 24, 25, 26, 27, + 12, 13, 14, 15, + 28, 29, 30, 31}; +//[a,b],[ai,bi] = [a,ai] - This is equivalent to p16uc_GETREAL64 +const static Packet16uc MMA_p16uc_SETCOMPLEX64_FIRST = { 0, 1, 2, 3, 4, 5, 6, 7, + 16, 17, 18, 19, 20, 21, 22, 23}; + +//[a,b],[ai,bi] = [b,bi] - This is equivalent to p16uc_GETIMAG64 +const static Packet16uc MMA_p16uc_SETCOMPLEX64_SECOND = { 8, 9, 10, 11, 12, 13, 14, 15, + 24, 25, 26, 27, 28, 29, 30, 31}; + + + +// Grab two decouples real/imaginary PacketBlocks and return two coupled (real/imaginary pairs) PacketBlocks. +template<typename Packet, typename Packetc> +EIGEN_STRONG_INLINE void bcoupleMMA(PacketBlock<Packet,4>& taccReal, PacketBlock<Packet,4>& taccImag, PacketBlock<Packetc,8>& tRes, PacketBlock<Packetc, 4>& acc1, PacketBlock<Packetc, 4>& acc2) +{ + acc1.packet[0].v = vec_perm(taccReal.packet[0], taccImag.packet[0], MMA_p16uc_SETCOMPLEX32_FIRST); + acc1.packet[1].v = vec_perm(taccReal.packet[1], taccImag.packet[1], MMA_p16uc_SETCOMPLEX32_FIRST); + acc1.packet[2].v = vec_perm(taccReal.packet[2], taccImag.packet[2], MMA_p16uc_SETCOMPLEX32_FIRST); + acc1.packet[3].v = vec_perm(taccReal.packet[3], taccImag.packet[3], MMA_p16uc_SETCOMPLEX32_FIRST); + + acc2.packet[0].v = vec_perm(taccReal.packet[0], taccImag.packet[0], MMA_p16uc_SETCOMPLEX32_SECOND); + acc2.packet[1].v = vec_perm(taccReal.packet[1], taccImag.packet[1], MMA_p16uc_SETCOMPLEX32_SECOND); + acc2.packet[2].v = vec_perm(taccReal.packet[2], taccImag.packet[2], MMA_p16uc_SETCOMPLEX32_SECOND); + acc2.packet[3].v = vec_perm(taccReal.packet[3], taccImag.packet[3], MMA_p16uc_SETCOMPLEX32_SECOND); + + acc1.packet[0] = padd<Packetc>(tRes.packet[0], acc1.packet[0]); + acc1.packet[1] = padd<Packetc>(tRes.packet[1], acc1.packet[1]); + acc1.packet[2] = padd<Packetc>(tRes.packet[2], acc1.packet[2]); + acc1.packet[3] = padd<Packetc>(tRes.packet[3], acc1.packet[3]); + + acc2.packet[0] = padd<Packetc>(tRes.packet[4], acc2.packet[0]); + acc2.packet[1] = padd<Packetc>(tRes.packet[5], acc2.packet[1]); + acc2.packet[2] = padd<Packetc>(tRes.packet[6], acc2.packet[2]); + acc2.packet[3] = padd<Packetc>(tRes.packet[7], acc2.packet[3]); +} + +template<> +EIGEN_STRONG_INLINE void bcoupleMMA<Packet2d, Packet1cd>(PacketBlock<Packet2d,4>& taccReal, PacketBlock<Packet2d,4>& taccImag, PacketBlock<Packet1cd,8>& tRes, PacketBlock<Packet1cd, 4>& acc1, PacketBlock<Packet1cd, 4>& acc2) +{ + acc1.packet[0].v = vec_perm(taccReal.packet[0], taccImag.packet[0], MMA_p16uc_SETCOMPLEX64_FIRST); + acc1.packet[1].v = vec_perm(taccReal.packet[1], taccImag.packet[1], MMA_p16uc_SETCOMPLEX64_FIRST); + acc1.packet[2].v = vec_perm(taccReal.packet[2], taccImag.packet[2], MMA_p16uc_SETCOMPLEX64_FIRST); + acc1.packet[3].v = vec_perm(taccReal.packet[3], taccImag.packet[3], MMA_p16uc_SETCOMPLEX64_FIRST); + + acc2.packet[0].v = vec_perm(taccReal.packet[0], taccImag.packet[0], MMA_p16uc_SETCOMPLEX64_SECOND); + acc2.packet[1].v = vec_perm(taccReal.packet[1], taccImag.packet[1], MMA_p16uc_SETCOMPLEX64_SECOND); + acc2.packet[2].v = vec_perm(taccReal.packet[2], taccImag.packet[2], MMA_p16uc_SETCOMPLEX64_SECOND); + acc2.packet[3].v = vec_perm(taccReal.packet[3], taccImag.packet[3], MMA_p16uc_SETCOMPLEX64_SECOND); + + acc1.packet[0] = padd<Packet1cd>(tRes.packet[0], acc1.packet[0]); + acc1.packet[1] = padd<Packet1cd>(tRes.packet[1], acc1.packet[1]); + acc1.packet[2] = padd<Packet1cd>(tRes.packet[2], acc1.packet[2]); + acc1.packet[3] = padd<Packet1cd>(tRes.packet[3], acc1.packet[3]); + + acc2.packet[0] = padd<Packet1cd>(tRes.packet[4], acc2.packet[0]); + acc2.packet[1] = padd<Packet1cd>(tRes.packet[5], acc2.packet[1]); + acc2.packet[2] = padd<Packet1cd>(tRes.packet[6], acc2.packet[2]); + acc2.packet[3] = padd<Packet1cd>(tRes.packet[7], acc2.packet[3]); +} + +template<typename Scalar, typename Packet> +EIGEN_STRONG_INLINE Packet ploadLhsMMA(const Scalar *lhs) +{ + return *((Packet *)lhs); +} + +template<typename Packet> +EIGEN_STRONG_INLINE PacketBlock<Packet,2> pmul (const PacketBlock<Packet,2>& a, const Packet& b) +{ + PacketBlock<Packet,2> pb; + pb.packet[0] = a.packet[0]*b; + pb.packet[1] = a.packet[1]*b; + return pb; +} + +template<typename Scalar, typename Packet> +EIGEN_STRONG_INLINE void bsetzeroMMA(__vector_quad *acc) +{ + __builtin_mma_xxsetaccz(acc); +} + +template<typename DataMapper, typename Index, typename Packet> +EIGEN_STRONG_INLINE void storeAccumulator(Index i, Index j, const DataMapper& data, const Packet& alpha, __vector_quad *acc) +{ + PacketBlock<Packet, 4> result; + __builtin_mma_disassemble_acc(&result.packet, acc); + + PacketBlock<Packet, 4> block; + block.packet[0] = data.template loadPacket<Packet>(i, j + 0) + pmul<Packet>(alpha, result.packet[0]); + block.packet[1] = data.template loadPacket<Packet>(i, j + 1) + pmul<Packet>(alpha, result.packet[1]); + block.packet[2] = data.template loadPacket<Packet>(i, j + 2) + pmul<Packet>(alpha, result.packet[2]); + block.packet[3] = data.template loadPacket<Packet>(i, j + 3) + pmul<Packet>(alpha, result.packet[3]); + + data.template storePacketBlock<Packet, 4>(i, j, block); +} + +template<typename DataMapper, typename Index, typename Packet, typename Packetc, int N> +EIGEN_STRONG_INLINE void storeComplexAccumulator(Index i, Index j, const DataMapper& data, const Packet& alphaReal, const Packet& alphaImag, __vector_quad *accReal, __vector_quad *accImag, const int accColsC) +{ + PacketBlock<Packet, 4> resultReal, resultImag; + __builtin_mma_disassemble_acc(&resultReal.packet, accReal); + __builtin_mma_disassemble_acc(&resultImag.packet, accImag); + + PacketBlock<Packet,4> taccReal, taccImag; + taccReal.packet[0] = pmul<Packet>(resultReal.packet[0], alphaReal); + taccReal.packet[1] = pmul<Packet>(resultReal.packet[1], alphaReal); + taccReal.packet[2] = pmul<Packet>(resultReal.packet[2], alphaReal); + taccReal.packet[3] = pmul<Packet>(resultReal.packet[3], alphaReal); + + taccImag.packet[0] = pmul<Packet>(resultImag.packet[0], alphaReal); + taccImag.packet[1] = pmul<Packet>(resultImag.packet[1], alphaReal); + taccImag.packet[2] = pmul<Packet>(resultImag.packet[2], alphaReal); + taccImag.packet[3] = pmul<Packet>(resultImag.packet[3], alphaReal); + + taccReal.packet[0] = psub<Packet>(taccReal.packet[0], pmul<Packet>(resultImag.packet[0], alphaImag)); + taccReal.packet[1] = psub<Packet>(taccReal.packet[1], pmul<Packet>(resultImag.packet[1], alphaImag)); + taccReal.packet[2] = psub<Packet>(taccReal.packet[2], pmul<Packet>(resultImag.packet[2], alphaImag)); + taccReal.packet[3] = psub<Packet>(taccReal.packet[3], pmul<Packet>(resultImag.packet[3], alphaImag)); + + taccImag.packet[0] = pmadd<Packet>(resultReal.packet[0], alphaImag, taccImag.packet[0]); + taccImag.packet[1] = pmadd<Packet>(resultReal.packet[1], alphaImag, taccImag.packet[1]); + taccImag.packet[2] = pmadd<Packet>(resultReal.packet[2], alphaImag, taccImag.packet[2]); + taccImag.packet[3] = pmadd<Packet>(resultReal.packet[3], alphaImag, taccImag.packet[3]); + + PacketBlock<Packetc, 8> tRes; + tRes.packet[0] = data.template loadPacket<Packetc>(i + N*accColsC, j + 0); + tRes.packet[1] = data.template loadPacket<Packetc>(i + N*accColsC, j + 1); + tRes.packet[2] = data.template loadPacket<Packetc>(i + N*accColsC, j + 2); + tRes.packet[3] = data.template loadPacket<Packetc>(i + N*accColsC, j + 3); + + tRes.packet[4] = data.template loadPacket<Packetc>(i + (N+1)*accColsC, j + 0); + tRes.packet[5] = data.template loadPacket<Packetc>(i + (N+1)*accColsC, j + 1); + tRes.packet[6] = data.template loadPacket<Packetc>(i + (N+1)*accColsC, j + 2); + tRes.packet[7] = data.template loadPacket<Packetc>(i + (N+1)*accColsC, j + 3); + + PacketBlock<Packetc, 4> acc1, acc2; + bcoupleMMA<Packet, Packetc>(taccReal, taccImag, tRes, acc1, acc2); + + data.template storePacketBlock<Packetc, 4>(i + N*accColsC, j, acc1); + data.template storePacketBlock<Packetc, 4>(i + (N+1)*accColsC, j, acc2); +} + +// Defaults to float32, since Eigen still supports C++03 we can't use default template arguments +template<typename LhsPacket, typename RhsPacket, bool NegativeAccumulate> +EIGEN_STRONG_INLINE void pgerMMA(__vector_quad *acc, const RhsPacket& a, const LhsPacket& b) +{ + if(NegativeAccumulate) + { + __builtin_mma_xvf32gernp(acc, (__vector unsigned char)a, (__vector unsigned char)b); + } else { + __builtin_mma_xvf32gerpp(acc, (__vector unsigned char)a, (__vector unsigned char)b); + } +} + +template<> +EIGEN_STRONG_INLINE void pgerMMA<Packet2d, PacketBlock<Packet2d, 2>, false>(__vector_quad *acc, const PacketBlock<Packet2d,2>& a, const Packet2d& b) +{ + Packetx2u<Packet2d> p; + p.pair = a; + __builtin_mma_xvf64gerpp(acc, p.vectorpair, (__vector unsigned char)b); +} + +template<> +EIGEN_STRONG_INLINE void pgerMMA<Packet2d, PacketBlock<Packet2d, 2>, true>(__vector_quad *acc, const PacketBlock<Packet2d, 2>& a, const Packet2d& b) +{ + Packetx2u<Packet2d> p; + p.pair = a; + __builtin_mma_xvf64gernp(acc, p.vectorpair, (__vector unsigned char)b); +} + +// This is necessary because ploadRhs for double returns a pair of vectors when MMA is enabled. +template<typename Scalar, typename Packet> +EIGEN_STRONG_INLINE Packet ploadRhsMMA(const Scalar *rhs) +{ + return *((Packet *)rhs); +} + +template<> +EIGEN_STRONG_INLINE PacketBlock<Packet2d, 2> ploadRhsMMA<double, PacketBlock<Packet2d, 2> >(const double *rhs) +{ + PacketBlock<Packet2d, 2> pair; + pair.packet[0] = *((Packet2d *)rhs ); + pair.packet[1] = *(((Packet2d *)rhs) + 1); + return pair; +} + +template<typename Scalar, typename Index, typename Packet, typename RhsPacket, typename DataMapper> +void gemmMMA(const DataMapper& res, const Scalar* blockA, const Scalar* blockB, + Index rows, Index depth, Index cols, Scalar alpha, Index strideA, Index strideB, Index offsetA, Index offsetB, const int accRows, const int accCols) +{ + const Index remaining_rows = rows % accCols; + const Index remaining_cols = cols % accRows; + + if( strideA == -1 ) strideA = depth; + if( strideB == -1 ) strideB = depth; + + const Packet pAlpha = pset1<Packet>(alpha); + Index col = 0; + for(; col + accRows <= cols; col += accRows) + { + const Scalar *rhs_base = blockB + ( col/accRows )*strideB*accRows; + const Scalar *lhs_base = blockA; + + Index row = 0; + for(; row + accCols <= rows; row += accCols) + { + const Scalar *rhs_ptr = rhs_base; + const Scalar *lhs_ptr1 = lhs_base + (row/accCols)*strideA*accCols; + + __vector_quad acc; + bsetzeroMMA<Scalar, Packet>(&acc); + + lhs_ptr1 += accCols*offsetA; + rhs_ptr += accRows*offsetB; + for(Index k = 0; k < depth; k++) + { + Packet lhsV = ploadLhsMMA<Scalar, Packet>(lhs_ptr1); + RhsPacket rhsV = ploadRhsMMA<Scalar, RhsPacket>(rhs_ptr); + + pgerMMA<Packet, RhsPacket, false>(&acc, rhsV, lhsV); + + lhs_ptr1 += accCols; + rhs_ptr += accRows; + } + + storeAccumulator<DataMapper, Index, Packet>(row, col, res, pAlpha, &acc); + } + if(remaining_rows > 0) + { + const Scalar *rhs_ptr = rhs_base; + const Scalar *lhs_ptr = lhs_base + (row/accCols)*strideA*accCols; + + lhs_ptr += remaining_rows*offsetA; + rhs_ptr += accRows*offsetB; + for(Index k = 0; k < depth; k++) + { + for(Index arow = 0; arow < remaining_rows; arow++) + { + for(Index acol = 0; acol < accRows; acol++ ) + { + res(row + arow, col + acol) += alpha*lhs_ptr[arow]*rhs_ptr[acol]; + } + } + rhs_ptr += accRows; + lhs_ptr += remaining_rows; + } + } + } + + if(remaining_cols > 0) + { + const Scalar *rhs_base = blockB + (col/accRows)*strideB*accRows; + const Scalar *lhs_base = blockA; + + Index row = 0; + for(; row + accCols <= rows; row += accCols) + { + const Scalar *rhs_ptr = rhs_base; + const Scalar *lhs_ptr = lhs_base + (row/accCols)*strideA*accCols; + + lhs_ptr += accCols*offsetA; + rhs_ptr += remaining_cols*offsetB; + for(Index k = 0; k < depth; k++) + { + for(Index arow = 0; arow < accCols; arow++) + { + for(Index acol = 0; acol < remaining_cols; acol++ ) + { + res(row + arow, col + acol) += alpha*lhs_ptr[arow]*rhs_ptr[acol]; + } + } + rhs_ptr += remaining_cols; + lhs_ptr += accCols; + } + } + + if(remaining_rows > 0 ) + { + const Scalar *rhs_ptr = rhs_base; + const Scalar *lhs_ptr = lhs_base + (row/accCols)*strideA*accCols; + + lhs_ptr += remaining_rows*offsetA; + rhs_ptr += remaining_cols*offsetB; + for(Index k = 0; k < depth; k++) + { + for(Index arow = 0; arow < remaining_rows; arow++) + { + for(Index acol = 0; acol < remaining_cols; acol++ ) + { + res(row + arow, col + acol) += alpha*lhs_ptr[arow]*rhs_ptr[acol]; + } + } + rhs_ptr += remaining_cols; + lhs_ptr += remaining_rows; + } + } + } +} + +template<typename LhsScalar, typename RhsScalar, typename Scalarc, typename Scalar, typename Index, typename Packet, typename Packetc, typename RhsPacket, typename DataMapper, bool ConjugateLhs, bool ConjugateRhs, bool LhsIsReal, bool RhsIsReal> +void gemm_complexMMA(const DataMapper& res, const LhsScalar* blockAc, const RhsScalar* blockBc, + Index rows, Index depth, Index cols, Scalarc alpha, Index strideA, Index strideB, Index offsetA, Index offsetB, const int accRows, const int accCols) +{ + const int remaining_rows = rows % accCols; + const int remaining_cols = cols % accRows; + const int accColsC = accCols / 2; + int advanceCols = 2; + int advanceRows = 2; + + if(LhsIsReal) advanceRows = 1; + if(RhsIsReal) advanceCols = 1; + + if( strideA == -1 ) strideA = depth; + if( strideB == -1 ) strideB = depth; + + const Packet pAlphaReal = pset1<Packet>(alpha.real()); + const Packet pAlphaImag = pset1<Packet>(alpha.imag()); + + const Scalar *blockA = (Scalar *) blockAc; + const Scalar *blockB = (Scalar *) blockBc; + + Packet conj = pset1<Packet>((Scalar)-1.0f); + + Index col = 0; + for(; col + accRows <= cols; col += accRows) + { + const Scalar *rhs_base = blockB + ( (advanceCols*col)/accRows )*strideB*accRows; + const Scalar *lhs_base = blockA; + + Index row = 0; + + for(; row + accCols <= rows; row += accCols) + { + const Scalar *rhs_ptr = rhs_base; + const Scalar *rhs_ptr_imag = rhs_ptr + accRows*strideB; + const Scalar *lhs_ptr = lhs_base + ((advanceRows*row)/accCols)*strideA*accCols; + const Scalar *lhs_ptr_imag = lhs_ptr + accCols*strideA; + + __vector_quad accReal, accImag; + __builtin_mma_xxsetaccz(&accReal); + __builtin_mma_xxsetaccz(&accImag); + + lhs_ptr += accCols*offsetA; + if(!LhsIsReal) + lhs_ptr_imag += accCols*offsetA; + rhs_ptr += accRows*offsetB; + if(!RhsIsReal) + rhs_ptr_imag += accRows*offsetB; + for(Index k = 0; k < depth; k++) + { + Packet lhsV = ploadLhsMMA<Scalar, Packet>(lhs_ptr); + RhsPacket rhsV = ploadRhs<Scalar, RhsPacket>(rhs_ptr); + + Packet lhsVi = ploadLhsMMA<Scalar, Packet>(lhs_ptr_imag); + RhsPacket rhsVi = ploadRhs<Scalar, RhsPacket>(rhs_ptr_imag); + + if(ConjugateLhs && !LhsIsReal) lhsVi = pmul<Packet>(lhsVi, conj); + if(ConjugateRhs && !RhsIsReal) rhsVi = pmul<Packet>(rhsVi, conj); + + if(LhsIsReal) + { + pgerMMA<Packet, RhsPacket, false>(&accReal, rhsV, lhsV); + pgerMMA<Packet, RhsPacket, false>(&accImag, rhsVi, lhsV); + } else if(RhsIsReal) { + pgerMMA<Packet, RhsPacket, false>(&accReal, rhsV, lhsV); + pgerMMA<Packet, RhsPacket, false>(&accImag, rhsV, lhsVi); + } else { + pgerMMA<Packet, RhsPacket, false>(&accReal, rhsV, lhsV); + pgerMMA<Packet, RhsPacket, true>(&accReal, rhsVi, lhsVi); + pgerMMA<Packet, RhsPacket, false>(&accImag, rhsVi, lhsV); + pgerMMA<Packet, RhsPacket, false>(&accImag, rhsV, lhsVi); + } + + lhs_ptr += accCols; + rhs_ptr += accRows; + if(!LhsIsReal) + lhs_ptr_imag += accCols; + if(!RhsIsReal) + rhs_ptr_imag += accRows; + } + + storeComplexAccumulator<DataMapper, Index, Packet, Packetc, 0>(row, col, res, pAlphaReal, pAlphaImag, &accReal, &accImag, accColsC); + } + + if(remaining_rows > 0) + { + const Scalar *rhs_ptr = rhs_base; + const Scalar *rhs_ptr_imag = rhs_ptr + accRows*strideB; + const Scalar *lhs_ptr = lhs_base + ((advanceRows*row)/accCols)*strideA*accCols; + const Scalar *lhs_ptr_imag = lhs_ptr + remaining_rows*strideA; + + lhs_ptr += remaining_rows*offsetA; + if(!LhsIsReal) + lhs_ptr_imag += remaining_rows*offsetA; + rhs_ptr += accRows*offsetB; + if(!RhsIsReal) + rhs_ptr_imag += accRows*offsetB; + for(Index k = 0; k < depth; k++) + { + for(Index arow = 0; arow < remaining_rows; arow++) + { + Scalar lhs_real = lhs_ptr[arow]; + Scalar lhs_imag; + if(!LhsIsReal) lhs_imag = lhs_ptr_imag[arow]; + + Scalarc lhsc; + + lhsc.real(lhs_real); + if(!LhsIsReal) + { + if(ConjugateLhs) + lhsc.imag(-lhs_imag); + else + lhsc.imag(lhs_imag); + } else { + //Lazy approach for now + lhsc.imag((Scalar)0); + } + + for(int acol = 0; acol < accRows; acol++ ) + { + Scalar rhs_real = rhs_ptr[acol]; + Scalar rhs_imag; + if(!RhsIsReal) rhs_imag = rhs_ptr_imag[acol]; + Scalarc rhsc; + + rhsc.real(rhs_real); + if(!RhsIsReal) + { + if(ConjugateRhs) + rhsc.imag(-rhs_imag); + else + rhsc.imag(rhs_imag); + } else { + //Lazy approach for now + rhsc.imag((Scalar)0); + } + res(row + arow, col + acol) += alpha*lhsc*rhsc; + } + } + rhs_ptr += accRows; + lhs_ptr += remaining_rows; + if(!LhsIsReal) + lhs_ptr_imag += remaining_rows; + if(!RhsIsReal) + rhs_ptr_imag += accRows; + } + } + } + + if(remaining_cols > 0) + { + const Scalar *rhs_base = blockB + ( (advanceCols*col)/accRows )*strideB*accRows; + const Scalar *lhs_base = blockA; + Index row = 0; + + for(; row + accCols <= rows; row += accCols) + { + const Scalar *rhs_ptr = rhs_base; + const Scalar *rhs_ptr_imag = rhs_ptr + remaining_cols*strideB; + const Scalar *lhs_ptr = lhs_base + ((advanceRows*row)/accCols)*strideA*accCols; + const Scalar *lhs_ptr_imag = lhs_ptr + accCols*strideA; + + lhs_ptr += accCols*offsetA; + if(!LhsIsReal) + lhs_ptr_imag += accCols*offsetA; + rhs_ptr += remaining_cols*offsetB; + if(!RhsIsReal) + rhs_ptr_imag += remaining_cols*offsetB; + Scalarc scalarAcc[4][4]; + for(Index arow = 0; arow < 4; arow++ ) + { + for(Index acol = 0; acol < 4; acol++ ) + { + scalarAcc[arow][acol].real((Scalar)0.0f); + scalarAcc[arow][acol].imag((Scalar)0.0f); + } + } + for(Index k = 0; k < depth; k++) + { + for(Index arow = 0; arow < accCols; arow++) + { + Scalar lhs_real = lhs_ptr[arow]; + Scalar lhs_imag; + if(!LhsIsReal) + { + lhs_imag = lhs_ptr_imag[arow]; + + if(ConjugateLhs) + lhs_imag *= -1; + } else { + lhs_imag = (Scalar)0; + } + + for(int acol = 0; acol < remaining_cols; acol++ ) + { + Scalar rhs_real = rhs_ptr[acol]; + Scalar rhs_imag; + if(!RhsIsReal) + { + rhs_imag = rhs_ptr_imag[acol]; + + if(ConjugateRhs) + rhs_imag *= -1; + } else { + rhs_imag = (Scalar)0; + } + + scalarAcc[arow][acol].real(scalarAcc[arow][acol].real() + lhs_real*rhs_real - lhs_imag*rhs_imag); + scalarAcc[arow][acol].imag(scalarAcc[arow][acol].imag() + lhs_imag*rhs_real + lhs_real*rhs_imag); + } + } + rhs_ptr += remaining_cols; + lhs_ptr += accCols; + if(!RhsIsReal) + rhs_ptr_imag += remaining_cols; + if(!LhsIsReal) + lhs_ptr_imag += accCols; + } + for(int arow = 0; arow < accCols; arow++ ) + { + for(int acol = 0; acol < remaining_cols; acol++ ) + { + Scalar accR = scalarAcc[arow][acol].real(); + Scalar accI = scalarAcc[arow][acol].imag(); + Scalar aR = alpha.real(); + Scalar aI = alpha.imag(); + Scalar resR = res(row + arow, col + acol).real(); + Scalar resI = res(row + arow, col + acol).imag(); + + res(row + arow, col + acol).real(resR + accR*aR - accI*aI); + res(row + arow, col + acol).imag(resI + accR*aI + accI*aR); + } + } + } + + if(remaining_rows > 0) + { + const Scalar *rhs_ptr = rhs_base; + const Scalar *rhs_ptr_imag = rhs_ptr + remaining_cols*strideB; + const Scalar *lhs_ptr = lhs_base + ((advanceRows*row)/accCols)*strideA*accCols; + const Scalar *lhs_ptr_imag = lhs_ptr + remaining_rows*strideA; + + lhs_ptr += remaining_rows*offsetA; + if(!LhsIsReal) + lhs_ptr_imag += remaining_rows*offsetA; + rhs_ptr += remaining_cols*offsetB; + if(!RhsIsReal) + rhs_ptr_imag += remaining_cols*offsetB; + for(Index k = 0; k < depth; k++) + { + for(Index arow = 0; arow < remaining_rows; arow++) + { + Scalar lhs_real = lhs_ptr[arow]; + Scalar lhs_imag; + if(!LhsIsReal) lhs_imag = lhs_ptr_imag[arow]; + Scalarc lhsc; + + lhsc.real(lhs_real); + if(!LhsIsReal) + { + if(ConjugateLhs) + lhsc.imag(-lhs_imag); + else + lhsc.imag(lhs_imag); + } else { + lhsc.imag((Scalar)0); + } + + for(Index acol = 0; acol < remaining_cols; acol++ ) + { + Scalar rhs_real = rhs_ptr[acol]; + Scalar rhs_imag; + if(!RhsIsReal) rhs_imag = rhs_ptr_imag[acol]; + Scalarc rhsc; + + rhsc.real(rhs_real); + if(!RhsIsReal) + { + if(ConjugateRhs) + rhsc.imag(-rhs_imag); + else + rhsc.imag(rhs_imag); + } else { + rhsc.imag((Scalar)0); + } + res(row + arow, col + acol) += alpha*lhsc*rhsc; + } + } + rhs_ptr += remaining_cols; + lhs_ptr += remaining_rows; + if(!LhsIsReal) + lhs_ptr_imag += remaining_rows; + if(!RhsIsReal) + rhs_ptr_imag += remaining_cols; + } + } + } +} + +#pragma GCC reset_options +} // end namespace internal + +} // end namespace Eigen +#endif // EIGEN_MATRIX_PRODUCT_MMA_ALTIVEC_H +