devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

- Added problem size constructor to decompositions that did not have one. It preallocates member data structures.

Browse Source 
- Updated unit tests to check above constructor.
- In the compute() method of decompositions: Made temporary matrices/vectors class members to avoid heap allocations during compute() (when dynamic matrices are used, of course).

These  changes can speed up decomposition computation time when a solver instance is used to solve multiple same-sized problems. An added benefit is that the compute() method can now be invoked in contexts were heap allocations are forbidden, such as in real-time control loops.

CAVEAT: Not all of the decompositions in the Eigenvalues module have a heap-allocation-free compute() method. A future patch may address this issue, but some required API changes need to be incorporated first.
This commit is contained in:
Adolfo Rodriguez Tsouroukdissian
2010-04-21 17:15:57 +02:00
parent faf8f7732d
commit 28dde19e40
29 changed files with 396 additions and 121 deletions
Download Patch File Download Diff File Expand all files Collapse all files

33
Eigen/src/Eigenvalues/ComplexEigenSolver.h
View File

@@ -95,7 +95,24 @@ template<typename _MatrixType> class ComplexEigenSolver
* The default constructor is useful in cases in which the user intends to
* perform decompositions via compute().
*/
ComplexEigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false)
ComplexEigenSolver()
: m_eivec(),
m_eivalues(),
m_schur(),
m_isInitialized(false)
{}
/** \brief Default Constructor with memory preallocation
*
* Like the default constructor but with preallocation of the internal data
* according to the specified problem \a size.
* \sa ComplexEigenSolver()
*/
ComplexEigenSolver(int size)
: m_eivec(size, size),
m_eivalues(size),
m_schur(size),
m_isInitialized(false)
{}
/** \brief Constructor; computes eigendecomposition of given matrix.
@@ -107,6 +124,7 @@ template<typename _MatrixType> class ComplexEigenSolver
ComplexEigenSolver(const MatrixType& matrix)
: m_eivec(matrix.rows(),matrix.cols()),
m_eivalues(matrix.cols()),
m_schur(matrix.rows()),
m_isInitialized(false)
{
compute(matrix);
@@ -179,6 +197,7 @@ template<typename _MatrixType> class ComplexEigenSolver
protected:
EigenvectorType m_eivec;
EigenvalueType m_eivalues;
ComplexSchur<MatrixType> m_schur;
bool m_isInitialized;
};
@@ -193,8 +212,8 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
// Step 1: Do a complex Schur decomposition, A = U T U^*
// The eigenvalues are on the diagonal of T.
ComplexSchur<MatrixType> schur(matrix);
m_eivalues = schur.matrixT().diagonal();
m_schur.compute(matrix);
m_eivalues = m_schur.matrixT().diagonal();
// Step 2: Compute X such that T = X D X^(-1), where D is the diagonal of T.
// The matrix X is unit triangular.
@@ -205,10 +224,10 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
// Compute X(i,k) using the (i,k) entry of the equation X T = D X
for(int i=k-1 ; i>=0 ; i--)
{
X.coeffRef(i,k) = -schur.matrixT().coeff(i,k);
X.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
if(k-i-1>0)
X.coeffRef(i,k) -= (schur.matrixT().row(i).segment(i+1,k-i-1) * X.col(k).segment(i+1,k-i-1)).value();
ComplexScalar z = schur.matrixT().coeff(i,i) - schur.matrixT().coeff(k,k);
X.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * X.col(k).segment(i+1,k-i-1)).value();
ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
if(z==ComplexScalar(0))
{
// If the i-th and k-th eigenvalue are equal, then z equals 0.
@@ -220,7 +239,7 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
}
// Step 3: Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
m_eivec = schur.matrixU() * X;
m_eivec = m_schur.matrixU() * X;
// .. and normalize the eigenvectors
for(int k=0 ; k<n ; k++)
{

14
Eigen/src/Eigenvalues/ComplexSchur.h
View File

@@ -96,7 +96,11 @@ template<typename _MatrixType> class ComplexSchur
* \sa compute() for an example.
*/
ComplexSchur(int size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
: m_matT(size,size), m_matU(size,size), m_isInitialized(false), m_matUisUptodate(false)
: m_matT(size,size),
m_matU(size,size),
m_hess(size),
m_isInitialized(false),
m_matUisUptodate(false)
{}
/** \brief Constructor; computes Schur decomposition of given matrix.
@@ -111,6 +115,7 @@ template<typename _MatrixType> class ComplexSchur
ComplexSchur(const MatrixType& matrix, bool skipU = false)
: m_matT(matrix.rows(),matrix.cols()),
m_matU(matrix.rows(),matrix.cols()),
m_hess(matrix.rows()),
m_isInitialized(false),
m_matUisUptodate(false)
{
@@ -182,6 +187,7 @@ template<typename _MatrixType> class ComplexSchur
protected:
ComplexMatrixType m_matT, m_matU;
HessenbergDecomposition<MatrixType> m_hess;
bool m_isInitialized;
bool m_matUisUptodate;
@@ -300,10 +306,10 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
// Reduce to Hessenberg form
// TODO skip Q if skipU = true
HessenbergDecomposition<MatrixType> hess(matrix);
m_hess.compute(matrix);
m_matT = hess.matrixH().template cast<ComplexScalar>();
if(!skipU) m_matU = hess.matrixQ().template cast<ComplexScalar>();
m_matT = m_hess.matrixH().template cast<ComplexScalar>();
if(!skipU) m_matU = m_hess.matrixQ().template cast<ComplexScalar>();
// Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.

11
Eigen/src/Eigenvalues/EigenSolver.h
View File

@@ -122,6 +122,17 @@ template<typename _MatrixType> class EigenSolver
*/
EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) {}
/** \brief Default Constructor with memory preallocation
*
* Like the default constructor but with preallocation of the internal data
* according to the specified problem \a size.
* \sa EigenSolver()
*/
EigenSolver(int size)
: m_eivec(size, size),
m_eivalues(size),
m_isInitialized(false) {}
/** \brief Constructor; computes eigendecomposition of given matrix.
*
* \param[in] matrix Square matrix whose eigendecomposition is to be computed.

19
Eigen/src/Eigenvalues/HessenbergDecomposition.h
View File

@@ -91,7 +91,8 @@ template<typename _MatrixType> class HessenbergDecomposition
* \sa compute() for an example.
*/
HessenbergDecomposition(int size = Size==Dynamic ? 2 : Size)
: m_matrix(size,size)
: m_matrix(size,size),
m_temp(size)
{
if(size>1)
m_hCoeffs.resize(size-1);
@@ -107,12 +108,13 @@ template<typename _MatrixType> class HessenbergDecomposition
* \sa matrixH() for an example.
*/
HessenbergDecomposition(const MatrixType& matrix)
: m_matrix(matrix)
: m_matrix(matrix),
m_temp(matrix.rows())
{
if(matrix.rows()<2)
return;
m_hCoeffs.resize(matrix.rows()-1,1);
_compute(m_matrix, m_hCoeffs);
_compute(m_matrix, m_hCoeffs, m_temp);
}
/** \brief Computes Hessenberg decomposition of given matrix.
@@ -137,7 +139,7 @@ template<typename _MatrixType> class HessenbergDecomposition
if(matrix.rows()<2)
return;
m_hCoeffs.resize(matrix.rows()-1,1);
_compute(m_matrix, m_hCoeffs);
_compute(m_matrix, m_hCoeffs, m_temp);
}
/** \brief Returns the Householder coefficients.
@@ -226,13 +228,14 @@ template<typename _MatrixType> class HessenbergDecomposition
private:
static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs);
typedef Matrix<Scalar, 1, Size, Options | RowMajor, 1, MaxSize> VectorType;
typedef typename NumTraits<Scalar>::Real RealScalar;
static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp);
protected:
MatrixType m_matrix;
CoeffVectorType m_hCoeffs;
VectorType m_temp;
};
#ifndef EIGEN_HIDE_HEAVY_CODE
@@ -250,11 +253,11 @@ template<typename _MatrixType> class HessenbergDecomposition
* \sa packedMatrix()
*/
template<typename MatrixType>
void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs)
void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp)
{
assert(matA.rows()==matA.cols());
int n = matA.rows();
VectorType temp(n);
temp.resize(n);
for (int i = 0; i<n-1; ++i)
{
// let's consider the vector v = i-th column starting at position i+1

20
Eigen/src/Eigenvalues/RealSchur.h
View File

@@ -78,6 +78,7 @@ template<typename _MatrixType> class RealSchur
typedef typename MatrixType::Scalar Scalar;
typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
/** \brief Default constructor.
*
@@ -92,7 +93,9 @@ template<typename _MatrixType> class RealSchur
*/
RealSchur(int size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
: m_matT(size, size),
m_matU(size, size),
m_matU(size, size),
m_workspaceVector(size),
m_hess(size),
m_isInitialized(false)
{ }
@@ -108,6 +111,8 @@ template<typename _MatrixType> class RealSchur
RealSchur(const MatrixType& matrix)
: m_matT(matrix.rows(),matrix.cols()),
m_matU(matrix.rows(),matrix.cols()),
m_workspaceVector(matrix.rows()),
m_hess(matrix.rows()),
m_isInitialized(false)
{
compute(matrix);
@@ -165,6 +170,8 @@ template<typename _MatrixType> class RealSchur
MatrixType m_matT;
MatrixType m_matU;
ColumnVectorType m_workspaceVector;
HessenbergDecomposition<MatrixType> m_hess;
bool m_isInitialized;
typedef Matrix<Scalar,3,1> Vector3s;
@@ -185,14 +192,13 @@ void RealSchur<MatrixType>::compute(const MatrixType& matrix)
// Step 1. Reduce to Hessenberg form
// TODO skip Q if skipU = true
HessenbergDecomposition<MatrixType> hess(matrix);
m_matT = hess.matrixH();
m_matU = hess.matrixQ();
m_hess.compute(matrix);
m_matT = m_hess.matrixH();
m_matU = m_hess.matrixQ();
// Step 2. Reduce to real Schur form
typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
ColumnVectorType workspaceVector(m_matU.cols());
Scalar* workspace = &workspaceVector.coeffRef(0);
m_workspaceVector.resize(m_matU.cols());
Scalar* workspace = &m_workspaceVector.coeffRef(0);
// The matrix m_matT is divided in three parts.
// Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.

29
Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
View File

@@ -58,14 +58,16 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
SelfAdjointEigenSolver()
: m_eivec(int(Size), int(Size)),
m_eivalues(int(Size))
m_eivalues(int(Size)),
m_subdiag(int(TridiagonalizationType::SizeMinusOne))
{
ei_assert(Size!=Dynamic);
}
SelfAdjointEigenSolver(int size)
: m_eivec(size, size),
m_eivalues(size)
m_eivalues(size),
m_subdiag(TridiagonalizationType::SizeMinusOne)
{}
/** Constructors computing the eigenvalues of the selfadjoint matrix \a matrix,
@@ -75,8 +77,10 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
*/
SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
: m_eivec(matrix.rows(), matrix.cols()),
m_eivalues(matrix.cols())
m_eivalues(matrix.cols()),
m_subdiag()
{
if (matrix.rows() > 1) m_subdiag.resize(matrix.rows() - 1);
compute(matrix, computeEigenvectors);
}
@@ -89,8 +93,10 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
*/
SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
: m_eivec(matA.rows(), matA.cols()),
m_eivalues(matA.cols())
m_eivalues(matA.cols()),
m_subdiag()
{
if (matA.rows() > 1) m_subdiag.resize(matA.rows() - 1);
compute(matA, matB, computeEigenvectors);
}
@@ -132,6 +138,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
protected:
MatrixType m_eivec;
RealVectorType m_eivalues;
typename TridiagonalizationType::SubDiagonalType m_subdiag;
#ifndef NDEBUG
bool m_eigenvectorsOk;
#endif
@@ -187,27 +194,27 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute(
// the latter avoids multiple memory allocation when the same SelfAdjointEigenSolver is used multiple times...
// (same for diag and subdiag)
RealVectorType& diag = m_eivalues;
typename TridiagonalizationType::SubDiagonalType subdiag(n-1);
TridiagonalizationType::decomposeInPlace(m_eivec, diag, subdiag, computeEigenvectors);
m_subdiag.resize(n-1);
TridiagonalizationType::decomposeInPlace(m_eivec, diag, m_subdiag, computeEigenvectors);
int end = n-1;
int start = 0;
while (end>0)
{
for (int i = start; i<end; ++i)
if (ei_isMuchSmallerThan(ei_abs(subdiag[i]),(ei_abs(diag[i])+ei_abs(diag[i+1]))))
subdiag[i] = 0;
if (ei_isMuchSmallerThan(ei_abs(m_subdiag[i]),(ei_abs(diag[i])+ei_abs(diag[i+1]))))
m_subdiag[i] = 0;
// find the largest unreduced block
while (end>0 && subdiag[end-1]==0)
while (end>0 && m_subdiag[end-1]==0)
end--;
if (end<=0)
break;
start = end - 1;
while (start>0 && subdiag[start-1]!=0)
while (start>0 && m_subdiag[start-1]!=0)
start--;
ei_tridiagonal_qr_step(diag.data(), subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
ei_tridiagonal_qr_step(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
}
// Sort eigenvalues and corresponding vectors.

4
Eigen/src/Eigenvalues/Tridiagonalization.h
View File

@@ -210,8 +210,8 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
// i.e., A = H A H' where H = I - h v v' and v = matA.col(i).tail(n-i-1)
matA.col(i).coeffRef(i+1) = 1;
hCoeffs.tail(n-i-1) = (matA.corner(BottomRight,remainingSize,remainingSize).template selfadjointView<Lower>()
* (ei_conj(h) * matA.col(i).tail(remainingSize)));
hCoeffs.tail(n-i-1).noalias() = (matA.corner(BottomRight,remainingSize,remainingSize).template selfadjointView<Lower>()
* (ei_conj(h) * matA.col(i).tail(remainingSize)));
hCoeffs.tail(n-i-1) += (ei_conj(h)*Scalar(-0.5)*(hCoeffs.tail(remainingSize).dot(matA.col(i).tail(remainingSize)))) * matA.col(i).tail(n-i-1);