mirror of
https://gitlab.com/libeigen/eigen.git
synced 2024-12-15 07:10:37 +08:00
fix two perf issues in product.
fix positive definite test in Cholesky. remove #include <cstring> in CoreDeclaration.
This commit is contained in:
parent
49ae3fca89
commit
f81dfcf00b
@ -30,7 +30,6 @@
|
|||||||
|
|
||||||
#include <cstdlib>
|
#include <cstdlib>
|
||||||
#include <cmath>
|
#include <cmath>
|
||||||
#include <cstring>
|
|
||||||
#include <complex>
|
#include <complex>
|
||||||
#include <cassert>
|
#include <cassert>
|
||||||
#include <functional>
|
#include <functional>
|
||||||
|
@ -96,7 +96,7 @@ void Cholesky<MatrixType>::compute(const MatrixType& a)
|
|||||||
|
|
||||||
RealScalar x;
|
RealScalar x;
|
||||||
x = ei_real(a.coeff(0,0));
|
x = ei_real(a.coeff(0,0));
|
||||||
m_isPositiveDefinite = x > precision<Scalar>() && ei_isMuchSmallerThan(ei_imag(m_matrix.coeff(0,0)), RealScalar(1));
|
m_isPositiveDefinite = x > precision<Scalar>() && ei_isMuchSmallerThan(ei_imag(a.coeff(0,0)), RealScalar(1));
|
||||||
m_matrix.coeffRef(0,0) = ei_sqrt(x);
|
m_matrix.coeffRef(0,0) = ei_sqrt(x);
|
||||||
m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0));
|
m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0));
|
||||||
for (int j = 1; j < size; ++j)
|
for (int j = 1; j < size; ++j)
|
||||||
@ -105,7 +105,7 @@ void Cholesky<MatrixType>::compute(const MatrixType& a)
|
|||||||
x = ei_real(tmp);
|
x = ei_real(tmp);
|
||||||
if (x < precision<Scalar>() || (!ei_isMuchSmallerThan(ei_imag(tmp), RealScalar(1))))
|
if (x < precision<Scalar>() || (!ei_isMuchSmallerThan(ei_imag(tmp), RealScalar(1))))
|
||||||
{
|
{
|
||||||
m_isPositiveDefinite = m_isPositiveDefinite;
|
m_isPositiveDefinite = false;
|
||||||
return;
|
return;
|
||||||
}
|
}
|
||||||
m_matrix.coeffRef(j,j) = x = ei_sqrt(x);
|
m_matrix.coeffRef(j,j) = x = ei_sqrt(x);
|
||||||
@ -117,6 +117,7 @@ void Cholesky<MatrixType>::compute(const MatrixType& a)
|
|||||||
m_matrix.col(j).end(endSize) =
|
m_matrix.col(j).end(endSize) =
|
||||||
(m_matrix.block(j+1, 0, endSize, j) * m_matrix.row(j).start(j).adjoint()).lazy();
|
(m_matrix.block(j+1, 0, endSize, j) * m_matrix.row(j).start(j).adjoint()).lazy();
|
||||||
|
|
||||||
|
// FIXME could use a.col instead of a.row
|
||||||
m_matrix.col(j).end(endSize) = (a.row(j).end(endSize).adjoint()
|
m_matrix.col(j).end(endSize) = (a.row(j).end(endSize).adjoint()
|
||||||
- m_matrix.col(j).end(endSize) ) / x;
|
- m_matrix.col(j).end(endSize) ) / x;
|
||||||
}
|
}
|
||||||
|
@ -402,8 +402,9 @@ EIGEN_DONT_INLINE static void ei_cache_friendly_product_colmajor_times_vector(
|
|||||||
{
|
{
|
||||||
ei_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0 || size<PacketSize);
|
ei_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0 || size<PacketSize);
|
||||||
|
|
||||||
for (; skipColumns<PacketSize && alignedStart != lhsAlignmentOffset + alignmentStep*skipColumns; ++skipColumns)
|
while (skipColumns<PacketSize &&
|
||||||
{}
|
alignedStart != ((lhsAlignmentOffset + alignmentStep*skipColumns)%PacketSize))
|
||||||
|
++skipColumns;
|
||||||
if (skipColumns==PacketSize)
|
if (skipColumns==PacketSize)
|
||||||
{
|
{
|
||||||
// nothing can be aligned, no need to skip any column
|
// nothing can be aligned, no need to skip any column
|
||||||
@ -568,7 +569,7 @@ EIGEN_DONT_INLINE static void ei_cache_friendly_product_rowmajor_times_vector(
|
|||||||
typedef typename ei_packet_traits<Scalar>::type Packet;
|
typedef typename ei_packet_traits<Scalar>::type Packet;
|
||||||
const int PacketSize = sizeof(Packet)/sizeof(Scalar);
|
const int PacketSize = sizeof(Packet)/sizeof(Scalar);
|
||||||
|
|
||||||
enum { AllAligned, EvenAligned, FirstAligned, NoneAligned };
|
enum { AllAligned=0, EvenAligned=1, FirstAligned=2, NoneAligned=3 };
|
||||||
const int rowsAtOnce = 4;
|
const int rowsAtOnce = 4;
|
||||||
const int peels = 2;
|
const int peels = 2;
|
||||||
const int PacketAlignedMask = PacketSize-1;
|
const int PacketAlignedMask = PacketSize-1;
|
||||||
@ -595,8 +596,9 @@ EIGEN_DONT_INLINE static void ei_cache_friendly_product_rowmajor_times_vector(
|
|||||||
{
|
{
|
||||||
ei_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0 || size<PacketSize);
|
ei_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0 || size<PacketSize);
|
||||||
|
|
||||||
for (; skipRows<PacketSize && alignedStart != lhsAlignmentOffset + alignmentStep*skipRows; ++skipRows)
|
while (skipRows<PacketSize &&
|
||||||
{}
|
alignedStart != ((lhsAlignmentOffset + alignmentStep*skipRows)%PacketSize))
|
||||||
|
++skipRows;
|
||||||
if (skipRows==PacketSize)
|
if (skipRows==PacketSize)
|
||||||
{
|
{
|
||||||
// nothing can be aligned, no need to skip any column
|
// nothing can be aligned, no need to skip any column
|
||||||
@ -611,7 +613,7 @@ EIGEN_DONT_INLINE static void ei_cache_friendly_product_rowmajor_times_vector(
|
|||||||
ei_internal_assert((alignmentPattern==NoneAligned) || PacketSize==1
|
ei_internal_assert((alignmentPattern==NoneAligned) || PacketSize==1
|
||||||
|| (size_t(lhs+alignedStart+lhsStride*skipRows)%sizeof(Packet))==0);
|
|| (size_t(lhs+alignedStart+lhsStride*skipRows)%sizeof(Packet))==0);
|
||||||
}
|
}
|
||||||
|
|
||||||
int rowBound = ((res.size()-skipRows)/rowsAtOnce)*rowsAtOnce + skipRows;
|
int rowBound = ((res.size()-skipRows)/rowsAtOnce)*rowsAtOnce + skipRows;
|
||||||
for (int i=skipRows; i<rowBound; i+=rowsAtOnce)
|
for (int i=skipRows; i<rowBound; i+=rowsAtOnce)
|
||||||
{
|
{
|
||||||
|
@ -637,7 +637,7 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo
|
|||||||
typedef typename ei_traits<ProductType>::_LhsNested Lhs;
|
typedef typename ei_traits<ProductType>::_LhsNested Lhs;
|
||||||
enum {
|
enum {
|
||||||
UseLhsDirectly = ((ei_packet_traits<Scalar>::size==1) || (Lhs::Flags&ActualPacketAccessBit))
|
UseLhsDirectly = ((ei_packet_traits<Scalar>::size==1) || (Lhs::Flags&ActualPacketAccessBit))
|
||||||
&& (!(Lhs::Flags & RowMajorBit)) };
|
&& (Lhs::Flags & RowMajorBit) };
|
||||||
|
|
||||||
template<typename DestDerived>
|
template<typename DestDerived>
|
||||||
inline static void run(DestDerived& res, const ProductType& product)
|
inline static void run(DestDerived& res, const ProductType& product)
|
||||||
|
@ -49,8 +49,9 @@ template<typename _MatrixType> class HessenbergDecomposition
|
|||||||
enum {
|
enum {
|
||||||
Size = MatrixType::RowsAtCompileTime,
|
Size = MatrixType::RowsAtCompileTime,
|
||||||
SizeMinusOne = MatrixType::RowsAtCompileTime==Dynamic
|
SizeMinusOne = MatrixType::RowsAtCompileTime==Dynamic
|
||||||
? Dynamic
|
? Dynamic
|
||||||
: MatrixType::RowsAtCompileTime-1};
|
: MatrixType::RowsAtCompileTime-1
|
||||||
|
};
|
||||||
|
|
||||||
typedef Matrix<Scalar, SizeMinusOne, 1> CoeffVectorType;
|
typedef Matrix<Scalar, SizeMinusOne, 1> CoeffVectorType;
|
||||||
typedef Matrix<RealScalar, Size, 1> DiagonalType;
|
typedef Matrix<RealScalar, Size, 1> DiagonalType;
|
||||||
@ -59,8 +60,7 @@ template<typename _MatrixType> class HessenbergDecomposition
|
|||||||
typedef typename NestByValue<DiagonalCoeffs<MatrixType> >::RealReturnType DiagonalReturnType;
|
typedef typename NestByValue<DiagonalCoeffs<MatrixType> >::RealReturnType DiagonalReturnType;
|
||||||
|
|
||||||
typedef typename NestByValue<DiagonalCoeffs<
|
typedef typename NestByValue<DiagonalCoeffs<
|
||||||
NestByValue<Block<
|
NestByValue<Block<MatrixType,SizeMinusOne,SizeMinusOne> > > >::RealReturnType SubDiagonalReturnType;
|
||||||
MatrixType,SizeMinusOne,SizeMinusOne> > > >::RealReturnType SubDiagonalReturnType;
|
|
||||||
|
|
||||||
/** This constructor initializes a HessenbergDecomposition object for
|
/** This constructor initializes a HessenbergDecomposition object for
|
||||||
* further use with HessenbergDecomposition::compute()
|
* further use with HessenbergDecomposition::compute()
|
||||||
@ -171,11 +171,11 @@ void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVector
|
|||||||
|
|
||||||
// first let's do A = H A
|
// first let's do A = H A
|
||||||
matA.corner(BottomRight,n-i-1,n-i-1) -= ((ei_conj(h) * matA.col(i).end(n-i-1)) *
|
matA.corner(BottomRight,n-i-1,n-i-1) -= ((ei_conj(h) * matA.col(i).end(n-i-1)) *
|
||||||
(matA.col(i).end(n-i-1).adjoint() * matA.corner(BottomRight,n-i-1,n-i-1)).lazy()).lazy();
|
(matA.col(i).end(n-i-1).adjoint() * matA.corner(BottomRight,n-i-1,n-i-1))).lazy();
|
||||||
|
|
||||||
// now let's do A = A H
|
// now let's do A = A H
|
||||||
matA.corner(BottomRight,n,n-i-1) -= ((matA.corner(BottomRight,n,n-i-1) * matA.col(i).end(n-i-1)).lazy() *
|
matA.corner(BottomRight,n,n-i-1) -= ((matA.corner(BottomRight,n,n-i-1) * matA.col(i).end(n-i-1))
|
||||||
(h * matA.col(i).end(n-i-1).adjoint())).lazy();
|
* (h * matA.col(i).end(n-i-1).adjoint())).lazy();
|
||||||
|
|
||||||
matA.col(i).coeffRef(i+1) = beta;
|
matA.col(i).coeffRef(i+1) = beta;
|
||||||
hCoeffs.coeffRef(i) = h;
|
hCoeffs.coeffRef(i) = h;
|
||||||
|
Loading…
Reference in New Issue
Block a user