diff options
Diffstat (limited to 'Eigen/src/Core/products/SelfadjointMatrixMatrix.h')
| -rw-r--r-- | Eigen/src/Core/products/SelfadjointMatrixMatrix.h | 104 |
1 files changed, 55 insertions, 49 deletions
diff --git a/Eigen/src/Core/products/SelfadjointMatrixMatrix.h b/Eigen/src/Core/products/SelfadjointMatrixMatrix.h index ede8b77bf..741c5f630 100644 --- a/Eigen/src/Core/products/SelfadjointMatrixMatrix.h +++ b/Eigen/src/Core/products/SelfadjointMatrixMatrix.h @@ -25,12 +25,14 @@ #ifndef EIGEN_SELFADJOINT_MATRIX_MATRIX_H #define EIGEN_SELFADJOINT_MATRIX_MATRIX_H +namespace internal { + // pack a selfadjoint block diagonal for use with the gebp_kernel template<typename Scalar, typename Index, int Pack1, int Pack2, int StorageOrder> -struct ei_symm_pack_lhs +struct symm_pack_lhs { template<int BlockRows> inline - void pack(Scalar* blockA, const ei_const_blas_data_mapper<Scalar,Index,StorageOrder>& lhs, Index cols, Index i, Index& count) + void pack(Scalar* blockA, const const_blas_data_mapper<Scalar,Index,StorageOrder>& lhs, Index cols, Index i, Index& count) { // normal copy for(Index k=0; k<i; k++) @@ -41,9 +43,9 @@ struct ei_symm_pack_lhs for(Index k=i; k<i+BlockRows; k++) { for(Index w=0; w<h; w++) - blockA[count++] = ei_conj(lhs(k, i+w)); // transposed + blockA[count++] = conj(lhs(k, i+w)); // transposed - blockA[count++] = ei_real(lhs(k,k)); // real (diagonal) + blockA[count++] = real(lhs(k,k)); // real (diagonal) for(Index w=h+1; w<BlockRows; w++) blockA[count++] = lhs(i+w, k); // normal @@ -52,11 +54,11 @@ struct ei_symm_pack_lhs // transposed copy for(Index k=i+BlockRows; k<cols; k++) for(Index w=0; w<BlockRows; w++) - blockA[count++] = ei_conj(lhs(k, i+w)); // transposed + blockA[count++] = conj(lhs(k, i+w)); // transposed } void operator()(Scalar* blockA, const Scalar* _lhs, Index lhsStride, Index cols, Index rows) { - ei_const_blas_data_mapper<Scalar,Index,StorageOrder> lhs(_lhs,lhsStride); + const_blas_data_mapper<Scalar,Index,StorageOrder> lhs(_lhs,lhsStride); Index count = 0; Index peeled_mc = (rows/Pack1)*Pack1; for(Index i=0; i<peeled_mc; i+=Pack1) @@ -76,23 +78,23 @@ struct ei_symm_pack_lhs for(Index k=0; k<i; k++) blockA[count++] = lhs(i, k); // normal - blockA[count++] = ei_real(lhs(i, i)); // real (diagonal) + blockA[count++] = real(lhs(i, i)); // real (diagonal) for(Index k=i+1; k<cols; k++) - blockA[count++] = ei_conj(lhs(k, i)); // transposed + blockA[count++] = conj(lhs(k, i)); // transposed } } }; template<typename Scalar, typename Index, int nr, int StorageOrder> -struct ei_symm_pack_rhs +struct symm_pack_rhs { - enum { PacketSize = ei_packet_traits<Scalar>::size }; + enum { PacketSize = packet_traits<Scalar>::size }; void operator()(Scalar* blockB, const Scalar* _rhs, Index rhsStride, Index rows, Index cols, Index k2) { Index end_k = k2 + rows; Index count = 0; - ei_const_blas_data_mapper<Scalar,Index,StorageOrder> rhs(_rhs,rhsStride); + const_blas_data_mapper<Scalar,Index,StorageOrder> rhs(_rhs,rhsStride); Index packet_cols = (cols/nr)*nr; // first part: normal case @@ -118,12 +120,12 @@ struct ei_symm_pack_rhs // transpose for(Index k=k2; k<j2; k++) { - blockB[count+0] = ei_conj(rhs(j2+0,k)); - blockB[count+1] = ei_conj(rhs(j2+1,k)); + blockB[count+0] = conj(rhs(j2+0,k)); + blockB[count+1] = conj(rhs(j2+1,k)); if (nr==4) { - blockB[count+2] = ei_conj(rhs(j2+2,k)); - blockB[count+3] = ei_conj(rhs(j2+3,k)); + blockB[count+2] = conj(rhs(j2+2,k)); + blockB[count+3] = conj(rhs(j2+3,k)); } count += nr; } @@ -135,11 +137,11 @@ struct ei_symm_pack_rhs for (Index w=0 ; w<h; ++w) blockB[count+w] = rhs(k,j2+w); - blockB[count+h] = ei_real(rhs(k,k)); + blockB[count+h] = real(rhs(k,k)); // transpose for (Index w=h+1 ; w<nr; ++w) - blockB[count+w] = ei_conj(rhs(j2+w,k)); + blockB[count+w] = conj(rhs(j2+w,k)); count += nr; ++h; } @@ -162,12 +164,12 @@ struct ei_symm_pack_rhs { for(Index k=k2; k<end_k; k++) { - blockB[count+0] = ei_conj(rhs(j2+0,k)); - blockB[count+1] = ei_conj(rhs(j2+1,k)); + blockB[count+0] = conj(rhs(j2+0,k)); + blockB[count+1] = conj(rhs(j2+1,k)); if (nr==4) { - blockB[count+2] = ei_conj(rhs(j2+2,k)); - blockB[count+3] = ei_conj(rhs(j2+3,k)); + blockB[count+2] = conj(rhs(j2+2,k)); + blockB[count+3] = conj(rhs(j2+3,k)); } count += nr; } @@ -180,13 +182,13 @@ struct ei_symm_pack_rhs Index half = std::min(end_k,j2); for(Index k=k2; k<half; k++) { - blockB[count] = ei_conj(rhs(j2,k)); + blockB[count] = conj(rhs(j2,k)); count += 1; } if(half==j2 && half<k2+rows) { - blockB[count] = ei_real(rhs(j2,j2)); + blockB[count] = real(rhs(j2,j2)); count += 1; } else @@ -209,12 +211,12 @@ template <typename Scalar, typename Index, int LhsStorageOrder, bool LhsSelfAdjoint, bool ConjugateLhs, int RhsStorageOrder, bool RhsSelfAdjoint, bool ConjugateRhs, int ResStorageOrder> -struct ei_product_selfadjoint_matrix; +struct product_selfadjoint_matrix; template <typename Scalar, typename Index, int LhsStorageOrder, bool LhsSelfAdjoint, bool ConjugateLhs, int RhsStorageOrder, bool RhsSelfAdjoint, bool ConjugateRhs> -struct ei_product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,LhsSelfAdjoint,ConjugateLhs, RhsStorageOrder,RhsSelfAdjoint,ConjugateRhs,RowMajor> +struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,LhsSelfAdjoint,ConjugateLhs, RhsStorageOrder,RhsSelfAdjoint,ConjugateRhs,RowMajor> { static EIGEN_STRONG_INLINE void run( @@ -224,7 +226,7 @@ struct ei_product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,LhsSelfAdjoint Scalar* res, Index resStride, Scalar alpha) { - ei_product_selfadjoint_matrix<Scalar, Index, + product_selfadjoint_matrix<Scalar, Index, EIGEN_LOGICAL_XOR(RhsSelfAdjoint,RhsStorageOrder==RowMajor) ? ColMajor : RowMajor, RhsSelfAdjoint, NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(RhsSelfAdjoint,ConjugateRhs), EIGEN_LOGICAL_XOR(LhsSelfAdjoint,LhsStorageOrder==RowMajor) ? ColMajor : RowMajor, @@ -237,7 +239,7 @@ struct ei_product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,LhsSelfAdjoint template <typename Scalar, typename Index, int LhsStorageOrder, bool ConjugateLhs, int RhsStorageOrder, bool ConjugateRhs> -struct ei_product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,true,ConjugateLhs, RhsStorageOrder,false,ConjugateRhs,ColMajor> +struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,true,ConjugateLhs, RhsStorageOrder,false,ConjugateRhs,ColMajor> { static EIGEN_DONT_INLINE void run( @@ -249,10 +251,10 @@ struct ei_product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,true,Conjugate { Index size = rows; - ei_const_blas_data_mapper<Scalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride); - ei_const_blas_data_mapper<Scalar, Index, RhsStorageOrder> rhs(_rhs,rhsStride); + const_blas_data_mapper<Scalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride); + const_blas_data_mapper<Scalar, Index, RhsStorageOrder> rhs(_rhs,rhsStride); - typedef ei_gebp_traits<Scalar,Scalar> Traits; + typedef gebp_traits<Scalar,Scalar> Traits; Index kc = size; // cache block size along the K direction Index mc = rows; // cache block size along the M direction @@ -267,10 +269,10 @@ struct ei_product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,true,Conjugate Scalar* allocatedBlockB = ei_aligned_stack_new(Scalar, sizeB); Scalar* blockB = allocatedBlockB + sizeW; - ei_gebp_kernel<Scalar, Scalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp_kernel; - ei_symm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs; - ei_gemm_pack_rhs<Scalar, Index, Traits::nr,RhsStorageOrder> pack_rhs; - ei_gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder==RowMajor?ColMajor:RowMajor, true> pack_lhs_transposed; + gebp_kernel<Scalar, Scalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp_kernel; + symm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs; + gemm_pack_rhs<Scalar, Index, Traits::nr,RhsStorageOrder> pack_rhs; + gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder==RowMajor?ColMajor:RowMajor, true> pack_lhs_transposed; for(Index k2=0; k2<size; k2+=kc) { @@ -305,7 +307,7 @@ struct ei_product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,true,Conjugate for(Index i2=k2+kc; i2<size; i2+=mc) { const Index actual_mc = std::min(i2+mc,size)-i2; - ei_gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder,false>() + gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder,false>() (blockA, &lhs(i2, k2), lhsStride, actual_kc, actual_mc); gebp_kernel(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, cols, alpha); @@ -321,7 +323,7 @@ struct ei_product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,true,Conjugate template <typename Scalar, typename Index, int LhsStorageOrder, bool ConjugateLhs, int RhsStorageOrder, bool ConjugateRhs> -struct ei_product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,false,ConjugateLhs, RhsStorageOrder,true,ConjugateRhs,ColMajor> +struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,false,ConjugateLhs, RhsStorageOrder,true,ConjugateRhs,ColMajor> { static EIGEN_DONT_INLINE void run( @@ -333,9 +335,9 @@ struct ei_product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,false,Conjugat { Index size = cols; - ei_const_blas_data_mapper<Scalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride); + const_blas_data_mapper<Scalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride); - typedef ei_gebp_traits<Scalar,Scalar> Traits; + typedef gebp_traits<Scalar,Scalar> Traits; Index kc = size; // cache block size along the K direction Index mc = rows; // cache block size along the M direction @@ -348,9 +350,9 @@ struct ei_product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,false,Conjugat Scalar* allocatedBlockB = ei_aligned_stack_new(Scalar, sizeB); Scalar* blockB = allocatedBlockB + sizeW; - ei_gebp_kernel<Scalar, Scalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp_kernel; - ei_gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs; - ei_symm_pack_rhs<Scalar, Index, Traits::nr,RhsStorageOrder> pack_rhs; + gebp_kernel<Scalar, Scalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp_kernel; + gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs; + symm_pack_rhs<Scalar, Index, Traits::nr,RhsStorageOrder> pack_rhs; for(Index k2=0; k2<size; k2+=kc) { @@ -373,14 +375,18 @@ struct ei_product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,false,Conjugat } }; +} // end namespace internal + /*************************************************************************** -* Wrapper to ei_product_selfadjoint_matrix +* Wrapper to product_selfadjoint_matrix ***************************************************************************/ +namespace internal { template<typename Lhs, int LhsMode, typename Rhs, int RhsMode> -struct ei_traits<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false> > - : ei_traits<ProductBase<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false>, Lhs, Rhs> > +struct traits<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false> > + : traits<ProductBase<SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false>, Lhs, Rhs> > {}; +} template<typename Lhs, int LhsMode, typename Rhs, int RhsMode> struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false> @@ -399,7 +405,7 @@ struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false> template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const { - ei_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols()); + eigen_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols()); const ActualLhsType lhs = LhsBlasTraits::extract(m_lhs); const ActualRhsType rhs = RhsBlasTraits::extract(m_rhs); @@ -407,14 +413,14 @@ struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,RhsMode,false> Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs) * RhsBlasTraits::extractScalarFactor(m_rhs); - ei_product_selfadjoint_matrix<Scalar, Index, + internal::product_selfadjoint_matrix<Scalar, Index, EIGEN_LOGICAL_XOR(LhsIsUpper, - ei_traits<Lhs>::Flags &RowMajorBit) ? RowMajor : ColMajor, LhsIsSelfAdjoint, + internal::traits<Lhs>::Flags &RowMajorBit) ? RowMajor : ColMajor, LhsIsSelfAdjoint, NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(LhsIsUpper,bool(LhsBlasTraits::NeedToConjugate)), EIGEN_LOGICAL_XOR(RhsIsUpper, - ei_traits<Rhs>::Flags &RowMajorBit) ? RowMajor : ColMajor, RhsIsSelfAdjoint, + internal::traits<Rhs>::Flags &RowMajorBit) ? RowMajor : ColMajor, RhsIsSelfAdjoint, NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(RhsIsUpper,bool(RhsBlasTraits::NeedToConjugate)), - ei_traits<Dest>::Flags&RowMajorBit ? RowMajor : ColMajor> + internal::traits<Dest>::Flags&RowMajorBit ? RowMajor : ColMajor> ::run( lhs.rows(), rhs.cols(), // sizes &lhs.coeff(0,0), lhs.outerStride(), // lhs info |