diff options
| -rw-r--r-- | Eigen/src/Core/Product.h | 29 | ||||
| -rw-r--r-- | Eigen/src/Core/SolveTriangular.h | 148 | ||||
| -rw-r--r-- | Eigen/src/Core/products/GeneralMatrixVector.h | 11 | ||||
| -rw-r--r-- | test/product_extra.cpp | 57 | ||||
| -rw-r--r-- | test/triangular.cpp | 51 |
5 files changed, 181 insertions, 115 deletions
diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index d63a7aa95..44fde3dcf 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -69,7 +69,7 @@ struct ProductReturnType<Lhs,Rhs,CacheFriendlyProduct> typedef typename ei_nested<Rhs,1, typename ei_plain_matrix_type_column_major<Rhs>::type >::type RhsNested; - + typedef Product<LhsNested, RhsNested, CacheFriendlyProduct> Type; }; @@ -268,7 +268,9 @@ template<typename LhsNested, typename RhsNested, int ProductMode> class Product */ EIGEN_STRONG_INLINE bool _useCacheFriendlyProduct() const { - return m_lhs.cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD + #define EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD 16 + // TODO do something more accurate here + return m_lhs.cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD && ( rows()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD || cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD); } @@ -624,7 +626,7 @@ struct ei_cache_friendly_product_selector<ProductType,LhsRows,ColMajor,HasDirect typedef typename LhsProductTraits::ActualXprType ActualLhsType; typedef typename RhsProductTraits::ActualXprType ActualRhsType; - + template<typename DestDerived> inline static void run(DestDerived& res, const ProductType& product, typename ProductType::Scalar alpha) { @@ -633,7 +635,7 @@ struct ei_cache_friendly_product_selector<ProductType,LhsRows,ColMajor,HasDirect Scalar actualAlpha = alpha * LhsProductTraits::extractScalarFactor(product.lhs()) * RhsProductTraits::extractScalarFactor(product.rhs()); - + enum { EvalToRes = (ei_packet_traits<Scalar>::size==1) ||((DestDerived::Flags&ActualPacketAccessBit) && (!(DestDerived::Flags & RowMajorBit))) }; @@ -645,7 +647,7 @@ struct ei_cache_friendly_product_selector<ProductType,LhsRows,ColMajor,HasDirect _res = ei_aligned_stack_new(Scalar,res.size()); Map<Matrix<Scalar,DestDerived::RowsAtCompileTime,1> >(_res, res.size()) = res; } - +// std::cerr << "colmajor * vector " << EvalToRes << "\n"; ei_cache_friendly_product_colmajor_times_vector <LhsProductTraits::NeedToConjugate,RhsProductTraits::NeedToConjugate>( res.size(), @@ -706,7 +708,7 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo _res = ei_aligned_stack_new(Scalar, res.size()); Map<Matrix<Scalar,DestDerived::SizeAtCompileTime,1> >(_res, res.size()) = res; } - + ei_cache_friendly_product_colmajor_times_vector <RhsProductTraits::NeedToConjugate,LhsProductTraits::NeedToConjugate>(res.size(), &actualRhs.const_cast_derived().coeffRef(0,0), actualRhs.stride(), @@ -725,7 +727,7 @@ template<typename ProductType, int LhsRows, int RhsOrder, int RhsAccess> struct ei_cache_friendly_product_selector<ProductType,LhsRows,RowMajor,HasDirectAccess,1,RhsOrder,RhsAccess> { typedef typename ProductType::Scalar Scalar; - + typedef ei_product_factor_traits<typename ei_traits<ProductType>::_LhsNested> LhsProductTraits; typedef ei_product_factor_traits<typename ei_traits<ProductType>::_RhsNested> RhsProductTraits; @@ -753,7 +755,7 @@ struct ei_cache_friendly_product_selector<ProductType,LhsRows,RowMajor,HasDirect _rhs = ei_aligned_stack_new(Scalar, actualRhs.size()); Map<Matrix<Scalar,ActualRhsType::SizeAtCompileTime,1> >(_rhs, actualRhs.size()) = actualRhs; } - + ei_cache_friendly_product_rowmajor_times_vector <LhsProductTraits::NeedToConjugate,RhsProductTraits::NeedToConjugate>( &actualLhs.const_cast_derived().coeffRef(0,0), actualLhs.stride(), @@ -774,7 +776,7 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo typedef typename LhsProductTraits::ActualXprType ActualLhsType; typedef typename RhsProductTraits::ActualXprType ActualRhsType; - + enum { UseLhsDirectly = ((ei_packet_traits<Scalar>::size==1) || (ActualLhsType::Flags&ActualPacketAccessBit)) && (ActualLhsType::Flags & RowMajorBit) }; @@ -796,7 +798,7 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo _lhs = ei_aligned_stack_new(Scalar, actualLhs.size()); Map<Matrix<Scalar,ActualLhsType::SizeAtCompileTime,1> >(_lhs, actualLhs.size()) = actualLhs; } - + ei_cache_friendly_product_rowmajor_times_vector <RhsProductTraits::NeedToConjugate, LhsProductTraits::NeedToConjugate>( &actualRhs.const_cast_derived().coeffRef(0,0), actualRhs.stride(), @@ -825,11 +827,12 @@ template<typename Derived> template<typename Lhs,typename Rhs> inline Derived& MatrixBase<Derived>::operator+=(const Flagged<Product<Lhs,Rhs,CacheFriendlyProduct>, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other) -{ +{//std::cerr << "operator+=\n"; if (other._expression()._useCacheFriendlyProduct()) ei_cache_friendly_product_selector<Product<Lhs,Rhs,CacheFriendlyProduct> >::run(const_cast_derived(), other._expression(), Scalar(1)); - else + else { //std::cerr << "no cf\n"; lazyAssign(derived() + other._expression()); + } return derived(); } @@ -893,7 +896,7 @@ inline void Product<Lhs,Rhs,ProductMode>::_cacheFriendlyEvalAndAdd(DestDerived& typedef typename LhsProductTraits::ActualXprType ActualLhsType; typedef typename RhsProductTraits::ActualXprType ActualRhsType; - + const ActualLhsType& actualLhs = LhsProductTraits::extract(m_lhs); const ActualRhsType& actualRhs = RhsProductTraits::extract(m_rhs); diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h index 8a6bf3af3..f3567c96c 100644 --- a/Eigen/src/Core/SolveTriangular.h +++ b/Eigen/src/Core/SolveTriangular.h @@ -25,7 +25,9 @@ #ifndef EIGEN_SOLVETRIANGULAR_H #define EIGEN_SOLVETRIANGULAR_H -template<typename Lhs, typename Rhs, int Mode, +template<typename Lhs, typename Rhs, + int Mode, // Upper/Lower | UnitDiag +// bool ConjugateLhs, bool ConjugateRhs, int UpLo = (Mode & LowerTriangularBit) ? LowerTriangular : (Mode & UpperTriangularBit) @@ -33,15 +35,57 @@ template<typename Lhs, typename Rhs, int Mode, : -1, int StorageOrder = int(Lhs::Flags) & RowMajorBit > -struct ei_solve_triangular_selector; +struct ei_triangular_solver_selector; // forward substitution, row-major -template<typename Lhs, typename Rhs, int Mode, int UpLo> -struct ei_solve_triangular_selector<Lhs,Rhs,Mode,UpLo,RowMajor> +template<typename Lhs, typename Rhs, int Mode, /*bool ConjugateLhs, bool ConjugateRhs,*/ int UpLo> +struct ei_triangular_solver_selector<Lhs,Rhs,Mode,/*ConjugateLhs,ConjugateRhs,*/UpLo,RowMajor> { typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) - { + {std::cerr << "here\n"; + #if NOTDEF + const bool IsLowerTriangular = (UpLo==LowerTriangular); + const int size = lhs.cols(); + for(int c=0 ; c<other.cols() ; ++c) + { + const int PanelWidth = 4; + for(int pi=IsLowerTriangular ? 0 : size; + IsLowerTriangular ? pi<size : pi>0; + IsLowerTriangular ? pi+=PanelWidth : pi-=PanelWidth) + { + int actualPanelWidth = std::min(IsLowerTriangular ? size - pi : pi, PanelWidth); + int startBlock = IsLowerTriangular ? pi : pi-actualPanelWidth; + int endBlock = IsLowerTriangular ? pi + actualPanelWidth : 0; + + if (pi > 0) + { + int r = IsLowerTriangular ? size - endBlock : startBlock; // remaining size + ei_cache_friendly_product_colmajor_times_vector<false,false>( + r, + &(lhs.const_cast_derived().coeffRef(endBlock,startBlock)), lhs.stride(), + other.col(c).segment(startBlock, actualPanelWidth), + &(other.coeffRef(endBlock, c)), + Scalar(-1)); + } + + for(int k=0; k<actualPanelWidth; ++k) + { + int i = IsLowerTriangular ? pi+k : pi-k-1; + if(!(Mode & UnitDiagBit)) + other.coeffRef(i,c) /= lhs.coeff(i,i); + + int r = actualPanelWidth - k - 1; // remaining size + if (r>0) + { + other.col(c).segment((IsLowerTriangular ? i+1 : i-r), r) -= + other.coeffRef(i,c) + * Block<Lhs,Dynamic,1>(lhs, (IsLowerTriangular ? i+1 : i-r), i, r, 1); + } + } + } + } + #else const bool IsLowerTriangular = (UpLo==LowerTriangular); const int size = lhs.cols(); /* We perform the inverse product per block of 4 rows such that we perfectly match @@ -115,6 +159,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,Mode,UpLo,RowMajor> } } } + #endif } }; @@ -124,7 +169,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,Mode,UpLo,RowMajor> // - inv(UpperTriangular, ColMajor) * Column vector // - inv(UpperTriangular,UnitDiag,ColMajor) * Column vector template<typename Lhs, typename Rhs, int Mode, int UpLo> -struct ei_solve_triangular_selector<Lhs,Rhs,Mode,UpLo,ColMajor> +struct ei_triangular_solver_selector<Lhs,Rhs,Mode,UpLo,ColMajor> { typedef typename Rhs::Scalar Scalar; typedef typename ei_packet_traits<Scalar>::type Packet; @@ -132,77 +177,44 @@ struct ei_solve_triangular_selector<Lhs,Rhs,Mode,UpLo,ColMajor> static void run(const Lhs& lhs, Rhs& other) { + static const int PanelWidth = 4; // TODO make this a user definable constant static const bool IsLowerTriangular = (UpLo==LowerTriangular); const int size = lhs.cols(); for(int c=0 ; c<other.cols() ; ++c) { - /* let's perform the inverse product per block of 4 columns such that we perfectly match - * our optimized matrix * vector product. blockyEnd represents the number of rows - * we can process using the block version. - */ - int blockyEnd = (std::max(size-5,0)/4)*4; - if (!IsLowerTriangular) - blockyEnd = size-1 - blockyEnd; - for(int i=IsLowerTriangular ? 0 : size-1; IsLowerTriangular ? i<blockyEnd : i>blockyEnd;) + for(int pi=IsLowerTriangular ? 0 : size; + IsLowerTriangular ? pi<size : pi>0; + IsLowerTriangular ? pi+=PanelWidth : pi-=PanelWidth) { - /* Let's process the 4x4 sub-matrix as usual. - * btmp stores the diagonal coefficients used to update the remaining part of the result. - */ - int startBlock = i; - int endBlock = startBlock + (IsLowerTriangular ? 4 : -4); - Matrix<Scalar,4,1> btmp; - for (;IsLowerTriangular ? i<endBlock : i>endBlock; - i += IsLowerTriangular ? 1 : -1) + int actualPanelWidth = std::min(IsLowerTriangular ? size - pi : pi, PanelWidth); + int startBlock = IsLowerTriangular ? pi : pi-actualPanelWidth; + int endBlock = IsLowerTriangular ? pi + actualPanelWidth : 0; + + for(int k=0; k<actualPanelWidth; ++k) { + int i = IsLowerTriangular ? pi+k : pi-k-1; if(!(Mode & UnitDiagBit)) other.coeffRef(i,c) /= lhs.coeff(i,i); - int remainingSize = IsLowerTriangular ? endBlock-i-1 : i-endBlock-1; - if (remainingSize>0) - other.col(c).segment((IsLowerTriangular ? i : endBlock) + 1, remainingSize) -= - other.coeffRef(i,c) - * Block<Lhs,Dynamic,1>(lhs, (IsLowerTriangular ? i : endBlock) + 1, i, remainingSize, 1); - btmp.coeffRef(IsLowerTriangular ? i-startBlock : remainingSize) = -other.coeffRef(i,c); - } - - /* Now we can efficiently update the remaining part of the result as a matrix * vector product. - * NOTE in order to reduce both compilation time and binary size, let's directly call - * the fast product implementation. It is equivalent to the following code: - * other.col(c).end(size-endBlock) += (lhs.block(endBlock, startBlock, size-endBlock, endBlock-startBlock) - * * other.col(c).block(startBlock,endBlock-startBlock)).lazy(); - */ - // FIXME this is cool but what about conjugate/adjoint expressions ? do we want to evaluate them ? - // this is a more general problem though. - ei_cache_friendly_product_colmajor_times_vector( - IsLowerTriangular ? size-endBlock : endBlock+1, - &(lhs.const_cast_derived().coeffRef(IsLowerTriangular ? endBlock : 0, IsLowerTriangular ? startBlock : endBlock+1)), - lhs.stride(), - btmp, &(other.coeffRef(IsLowerTriangular ? endBlock : 0, c)), - Scalar(1)); -// if (IsLowerTriangular) -// other.col(c).end(size-endBlock) += (lhs.block(endBlock, startBlock, size-endBlock, endBlock-startBlock) -// * other.col(c).block(startBlock,endBlock-startBlock)).lazy(); -// else -// other.col(c).end(size-endBlock) += (lhs.block(endBlock, startBlock, size-endBlock, endBlock-startBlock) -// * other.col(c).block(startBlock,endBlock-startBlock)).lazy(); - } - - /* Now we have to process the remaining part as usual */ - int i; - for(i=blockyEnd; IsLowerTriangular ? i<size-1 : i>0; i += (IsLowerTriangular ? 1 : -1) ) - { - if(!(Mode & UnitDiagBit)) - other.coeffRef(i,c) /= lhs.coeff(i,i); - /* NOTE we cannot use lhs.col(i).end(size-i-1) because Part::coeffRef gets called by .col() to - * get the address of the start of the row - */ - if(IsLowerTriangular) - other.col(c).end(size-i-1) -= other.coeffRef(i,c) * Block<Lhs,Dynamic,1>(lhs, i+1,i, size-i-1,1); - else - other.col(c).start(i) -= other.coeffRef(i,c) * Block<Lhs,Dynamic,1>(lhs, 0,i, i, 1); + int r = actualPanelWidth - k - 1; // remaining size + if (r>0) + { + other.col(c).segment((IsLowerTriangular ? i+1 : i-r), r) -= + other.coeffRef(i,c) + * Block<Lhs,Dynamic,1>(lhs, (IsLowerTriangular ? i+1 : i-r), i, r, 1); + } + } + int r = IsLowerTriangular ? size - endBlock : startBlock; // remaining size + if (r > 0) + { + ei_cache_friendly_product_colmajor_times_vector<false,false>( + r, + &(lhs.const_cast_derived().coeffRef(endBlock,startBlock)), lhs.stride(), + other.col(c).segment(startBlock, actualPanelWidth), + &(other.coeffRef(endBlock, c)), + Scalar(-1)); + } } - if(!(Mode & UnitDiagBit)) - other.coeffRef(i,c) /= lhs.coeff(i,i); } } }; @@ -232,7 +244,7 @@ void TriangularView<MatrixType,Mode>::solveInPlace(const MatrixBase<RhsDerived>& typename ei_plain_matrix_type_column_major<RhsDerived>::type, RhsDerived&>::ret RhsCopy; RhsCopy rhsCopy(rhs); - ei_solve_triangular_selector<MatrixType, typename ei_unref<RhsCopy>::type, Mode>::run(_expression(), rhsCopy); + ei_triangular_solver_selector<MatrixType, typename ei_unref<RhsCopy>::type, Mode>::run(_expression(), rhsCopy); if (copy) rhs = rhsCopy; diff --git a/Eigen/src/Core/products/GeneralMatrixVector.h b/Eigen/src/Core/products/GeneralMatrixVector.h index 851bf808f..61b2cc67c 100644 --- a/Eigen/src/Core/products/GeneralMatrixVector.h +++ b/Eigen/src/Core/products/GeneralMatrixVector.h @@ -57,6 +57,8 @@ void ei_cache_friendly_product_colmajor_times_vector( if(ConjugateRhs) alpha = ei_conj(alpha); +// std::cerr << "prod " << size << " " << rhs.size() << "\n"; + typedef typename ei_packet_traits<Scalar>::type Packet; const int PacketSize = sizeof(Packet)/sizeof(Scalar); @@ -101,7 +103,10 @@ void ei_cache_friendly_product_colmajor_times_vector( // note that the skiped columns are processed later. } - ei_internal_assert((alignmentPattern==NoneAligned) || (size_t(lhs+alignedStart+lhsStride*skipColumns)%sizeof(Packet))==0); + ei_internal_assert( (alignmentPattern==NoneAligned) + || (skipColumns + columnsAtOnce >= rhs.size()) + || PacketSize > size + || (size_t(lhs+alignedStart+lhsStride*skipColumns)%sizeof(Packet))==0); } int offset1 = (FirstAligned && alignmentStep==1?3:1); @@ -127,7 +132,6 @@ void ei_cache_friendly_product_colmajor_times_vector( res[j] = cj.pmadd(lhs1[j], ei_pfirst(ptmp1), res[j]); res[j] = cj.pmadd(lhs2[j], ei_pfirst(ptmp2), res[j]); res[j] = cj.pmadd(lhs3[j], ei_pfirst(ptmp3), res[j]); -// res[j] += ei_pfirst(ptmp0)*lhs0[j] + ei_pfirst(ptmp1)*lhs1[j] + ei_pfirst(ptmp2)*lhs2[j] + ei_pfirst(ptmp3)*lhs3[j]; } if (alignedSize>alignedStart) @@ -193,7 +197,6 @@ void ei_cache_friendly_product_colmajor_times_vector( res[j] = cj.pmadd(lhs1[j], ei_pfirst(ptmp1), res[j]); res[j] = cj.pmadd(lhs2[j], ei_pfirst(ptmp2), res[j]); res[j] = cj.pmadd(lhs3[j], ei_pfirst(ptmp3), res[j]); -// res[j] += ei_pfirst(ptmp0)*lhs0[j] + ei_pfirst(ptmp1)*lhs1[j] + ei_pfirst(ptmp2)*lhs2[j] + ei_pfirst(ptmp3)*lhs3[j]; } } @@ -212,7 +215,7 @@ void ei_cache_friendly_product_colmajor_times_vector( /* explicit vectorization */ // process first unaligned result's coeffs for (int j=0; j<alignedStart; ++j) - res[j] = cj.pmul(lhs0[j], ei_pfirst(ptmp0)); + res[j] += cj.pmul(lhs0[j], ei_pfirst(ptmp0)); // process aligned result's coeffs if ((size_t(lhs0+alignedStart)%sizeof(Packet))==0) diff --git a/test/product_extra.cpp b/test/product_extra.cpp index 4fa4c23f5..d73974886 100644 --- a/test/product_extra.cpp +++ b/test/product_extra.cpp @@ -46,9 +46,9 @@ template<typename MatrixType> void product_extra(const MatrixType& m) res = MatrixType::Random(rows, rows), square2 = MatrixType::Random(cols, cols), res2 = MatrixType::Random(cols, cols); - RowVectorType v1 = RowVectorType::Random(rows), - v2 = RowVectorType::Random(rows), - vzero = RowVectorType::Zero(rows); + RowVectorType v1 = RowVectorType::Random(rows), vrres(rows); +// v2 = RowVectorType::Random(rows), +// vzero = RowVectorType::Zero(rows); ColVectorType vc2 = ColVectorType::Random(cols), vcres(cols); OtherMajorMatrixType tm1 = m1; @@ -56,9 +56,14 @@ template<typename MatrixType> void product_extra(const MatrixType& m) s2 = ei_random<Scalar>(), s3 = ei_random<Scalar>(); + int c0 = ei_random<int>(0,cols/2-1), + c1 = ei_random<int>(cols/2,cols), + r0 = ei_random<int>(0,rows/2-1), + r1 = ei_random<int>(rows/2,rows); + // all the expressions in this test should be compiled as a single matrix product // TODO: add internal checks to verify that - +/* VERIFY_IS_APPROX(m1 * m2.adjoint(), m1 * m2.adjoint().eval()); VERIFY_IS_APPROX(m1.adjoint() * square.adjoint(), m1.adjoint().eval() * square.adjoint().eval()); VERIFY_IS_APPROX(m1.adjoint() * m2, m1.adjoint().eval() * m2); @@ -71,7 +76,7 @@ template<typename MatrixType> void product_extra(const MatrixType& m) // test all possible conjugate combinations for the four matrix-vector product cases: - + // std::cerr << "a\n"; VERIFY_IS_APPROX((-m1.conjugate() * s2) * (s1 * vc2), (-m1.conjugate()*s2).eval() * (s1 * vc2).eval()); @@ -106,15 +111,49 @@ template<typename MatrixType> void product_extra(const MatrixType& m) VERIFY_IS_APPROX((-m1.adjoint() * s2) * (s1 * v1.adjoint()), (-m1.adjoint()*s2).eval() * (s1 * v1.adjoint()).eval()); + */ + // test with sub matrices + m2 = m1; + m3 = m1; + +// std::cerr << (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).rows() << " " << (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).cols() << " == " << vrres.segment(r0,r1-r0).rows() << " " << vrres.segment(r0,r1-r0).cols() << "\n"; +// m2.col(c0).segment(0,r1-r0) += (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).lazy(); +// m3.col(c0).segment(0,r1-r0) += (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).eval(); + Matrix<Scalar,Dynamic,1> a = m2.col(c0), b = a; + a.segment(5,r1-r0) += (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).lazy(); + b.segment(5,r1-r0) += (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).eval(); + +// m2.col(c0).segment(0,r1-r0) += (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).lazy(); +// m3.col(c0).segment(0,r1-r0) += (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).eval(); +// if (!m2.isApprox(m3)) + std::cerr << (a-b).cwise().abs().maxCoeff() << "\n"; + VERIFY_IS_APPROX(a,b); +// VERIFY_IS_APPROX( vrres.segment(0,r1-r0).transpose().eval(), +// v1.segment(0,r1-r0).transpose() + m1.block(r0,c0, r1-r0, c1-c0).eval() * (vc2.segment(c0,c1-c0)).eval()); } void test_product_extra() { -// for(int i = 0; i < g_repeat; i++) { + for(int i = 0; i < g_repeat; i++) { + int rows = ei_random<int>(2,10); + int cols = ei_random<int>(2,10); + int c0 = ei_random<int>(0,cols/2-1), + c1 = ei_random<int>(cols/2,cols), + r0 = ei_random<int>(0,rows/2-1), + r1 = ei_random<int>(rows/2,rows); + + MatrixXf m1 = MatrixXf::Random(rows,cols), m2 = m1; + Matrix<float,Dynamic,1> a = m2.col(c0), b = a; + Matrix<float,Dynamic,1> vc2 = Matrix<float,Dynamic,1>::Random(cols); + if (1+r1-r0<rows) { + a.segment(1,r1-r0) += (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).lazy(); + b.segment(1,r1-r0) += (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).eval(); + VERIFY_IS_APPROX(a,b); + } // CALL_SUBTEST( product_extra(MatrixXf(ei_random<int>(1,320), ei_random<int>(1,320))) ); -// CALL_SUBTEST( product(MatrixXd(ei_random<int>(1,320), ei_random<int>(1,320))) ); +// CALL_SUBTEST( product_extra(MatrixXd(ei_random<int>(1,320), ei_random<int>(1,320))) ); // CALL_SUBTEST( product(MatrixXi(ei_random<int>(1,320), ei_random<int>(1,320))) ); - CALL_SUBTEST( product_extra(MatrixXcf(ei_random<int>(50,50), ei_random<int>(50,50))) ); +// CALL_SUBTEST( product_extra(MatrixXcf(ei_random<int>(50,50), ei_random<int>(50,50))) ); // CALL_SUBTEST( product(Matrix<float,Dynamic,Dynamic,RowMajor>(ei_random<int>(1,320), ei_random<int>(1,320))) ); -// } + } } diff --git a/test/triangular.cpp b/test/triangular.cpp index 1829a2578..3550d1a74 100644 --- a/test/triangular.cpp +++ b/test/triangular.cpp @@ -90,31 +90,35 @@ template<typename MatrixType> void triangular(const MatrixType& m) Transpose<MatrixType> trm4(m4); // test back and forward subsitution m3 = m1.template triangularView<Eigen::LowerTriangular>(); - VERIFY(m3.template triangularView<Eigen::LowerTriangular>().solve(m3).cwise().abs().isIdentity(test_precision<RealScalar>())); - VERIFY(m3.transpose().template triangularView<Eigen::UpperTriangular>() - .solve(m3.transpose()).cwise().abs().isIdentity(test_precision<RealScalar>())); +// VERIFY(m3.template triangularView<Eigen::LowerTriangular>().solve(m3).cwise().abs().isIdentity(test_precision<RealScalar>())); +// VERIFY(m3.transpose().template triangularView<Eigen::UpperTriangular>() +// .solve(m3.transpose()).cwise().abs().isIdentity(test_precision<RealScalar>())); // check M * inv(L) using in place API m4 = m3; m3.transpose().template triangularView<Eigen::UpperTriangular>().solveInPlace(trm4); - VERIFY(m4.cwise().abs().isIdentity(test_precision<RealScalar>())); +// VERIFY(m4.cwise().abs().isIdentity(test_precision<RealScalar>())); - m3 = m1.template triangularView<Eigen::UpperTriangular>(); - VERIFY(m3.template triangularView<Eigen::UpperTriangular>().solve(m3).cwise().abs().isIdentity(test_precision<RealScalar>())); - VERIFY(m3.transpose().template triangularView<Eigen::LowerTriangular>() - .solve(m3.transpose()).cwise().abs().isIdentity(test_precision<RealScalar>())); - // check M * inv(U) using in place API +// m3 = m1.template triangularView<Eigen::UpperTriangular>(); +// VERIFY(m3.template triangularView<Eigen::UpperTriangular>().solve(m3).cwise().abs().isIdentity(test_precision<RealScalar>())); +// VERIFY(m3.transpose().template triangularView<Eigen::LowerTriangular>() +// .solve(m3.transpose()).cwise().abs().isIdentity(test_precision<RealScalar>())); +// // check M * inv(U) using in place API m4 = m3; m3.transpose().template triangularView<Eigen::LowerTriangular>().solveInPlace(trm4); VERIFY(m4.cwise().abs().isIdentity(test_precision<RealScalar>())); - m3 = m1.template triangularView<Eigen::UpperTriangular>(); - VERIFY(m2.isApprox(m3 * (m3.template triangularView<Eigen::UpperTriangular>().solve(m2)), largerEps)); - m3 = m1.template triangularView<Eigen::LowerTriangular>(); - VERIFY(m2.isApprox(m3 * (m3.template triangularView<Eigen::LowerTriangular>().solve(m2)), largerEps)); +// m3 = m1.template triangularView<Eigen::UpperTriangular>(); +// VERIFY(m2.isApprox(m3 * (m3.template triangularView<Eigen::UpperTriangular>().solve(m2)), largerEps)); +// m3 = m1.template triangularView<Eigen::LowerTriangular>(); + +// std::cerr << (m2 - +// (m3 * (m3.template triangularView<Eigen::LowerTriangular>().solve(m2)))).cwise().abs() /*.maxCoeff()*/ << "\n\n"; + +// VERIFY(m2.isApprox(m3 * (m3.template triangularView<Eigen::LowerTriangular>().solve(m2)), largerEps)); // check solve with unit diagonal - m3 = m1.template triangularView<Eigen::UnitUpperTriangular>(); - VERIFY(m2.isApprox(m3 * (m1.template triangularView<Eigen::UnitUpperTriangular>().solve(m2)), largerEps)); +// m3 = m1.template triangularView<Eigen::UnitUpperTriangular>(); +// VERIFY(m2.isApprox(m3 * (m1.template triangularView<Eigen::UnitUpperTriangular>().solve(m2)), largerEps)); // VERIFY(( m1.template triangularView<Eigen::UpperTriangular>() // * m2.template triangularView<Eigen::UpperTriangular>()).isUpperTriangular()); @@ -132,12 +136,17 @@ template<typename MatrixType> void triangular(const MatrixType& m) void test_triangular() { for(int i = 0; i < g_repeat ; i++) { - CALL_SUBTEST( triangular(Matrix<float, 1, 1>()) ); - CALL_SUBTEST( triangular(Matrix<float, 2, 2>()) ); - CALL_SUBTEST( triangular(Matrix3d()) ); - CALL_SUBTEST( triangular(MatrixXcf(4, 4)) ); - CALL_SUBTEST( triangular(Matrix<std::complex<float>,8, 8>()) ); +// CALL_SUBTEST( triangular(Matrix<float, 1, 1>()) ); +// CALL_SUBTEST( triangular(Matrix<float, 2, 2>()) ); +// CALL_SUBTEST( triangular(Matrix3d()) ); +// CALL_SUBTEST( triangular(MatrixXcf(4, 4)) ); +// CALL_SUBTEST( triangular(Matrix<std::complex<float>,8, 8>()) ); +// CALL_SUBTEST( triangular(MatrixXd(1,1)) ); +// CALL_SUBTEST( triangular(MatrixXd(2,2)) ); +// CALL_SUBTEST( triangular(MatrixXd(3,3)) ); +// CALL_SUBTEST( triangular(MatrixXd(5,5)) ); +// CALL_SUBTEST( triangular(MatrixXd(8,8)) ); CALL_SUBTEST( triangular(MatrixXd(17,17)) ); - CALL_SUBTEST( triangular(Matrix<float,Dynamic,Dynamic,RowMajor>(5, 5)) ); +// CALL_SUBTEST( triangular(Matrix<float,Dynamic,Dynamic,RowMajor>(5, 5)) ); } } |