| Re: [eigen] Allocation policy of eigen decompositions |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
- To: eigen@xxxxxxxxxxxxxxxxxxx
- Subject: Re: [eigen] Allocation policy of eigen decompositions
- From: Adolfo Rodríguez Tsouroukdissian <dofo79@xxxxxxxxx>
- Date: Wed, 21 Apr 2010 17:41:40 +0200
- Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:in-reply-to:references :date:received:message-id:subject:from:to:content-type; bh=Rh+z28ct9Rm69/rTyNgPvVWQ4tpJnmGoZXth4tg7+pA=; b=GuK8I6rxigtPCOkUeat3EzQIYp3M8V5Q/IPme6XIdKzTHd59q9NvIRsGmysJMwTQYJ TUg5AXdztLUlzGE7C6ypkD2ndJ1YCSxfUHbWzdvLKNdoLnPTwT5FNK/mWFRn9eExY2Cz QmDOCbGNkhc4NvuDoFDqQ8auiE5cvwQgHUGUk=
- Domainkey-signature: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type; b=xIC01xApJBdRXhwOAJgILNZBp5nkpayFNuCJAKZyt+vXTlUdav7RSp/9R2XRdVPDhM h9xQ4RJe7qhxAXa7wr1U0NjFfuKqulbuHDaiZ9A4LeyecgwubmIV4IYBhOJYhyT1pxUo 1K8DTPle6IKRQXaU+qms3uXJfhmc1m95NGDIU=
Please find attached a patch that:
- Adds problem size constructor to decompositions that did not have one. It preallocates member data structures.
- Updates 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 single 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.
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 (see previous thread posts).
Please let me know if something is not working as advertised ;)
Regards,
Adolfo
--
Adolfo Rodríguez Tsouroukdissian, Ph. D.
Robotics engineer
PAL ROBOTICS S.L
http://www.pal-robotics.com
Tel. +34.93.414.53.47
Fax.+34.93.209.11.09
AVISO DE CONFIDENCIALIDAD: Este mensaje y sus documentos adjuntos, pueden contener información privilegiada y/o confidencial que está dirigida exclusivamente a su destinatario. Si usted recibe este mensaje y no es el destinatario indicado, o el empleado encargado de su entrega a dicha persona, por favor, notifíquelo inmediatamente y remita el mensaje original a la dirección de correo electrónico indicada. Cualquier copia, uso o distribución no autorizados de esta comunicación queda estrictamente prohibida.
CONFIDENTIALITY NOTICE: This e-mail and the accompanying document(s) may contain confidential information which is privileged and intended only for the individual or entity to whom they are addressed. If you are not the intended recipient, you are hereby notified that any disclosure, copying, distribution or use of this e-mail and/or accompanying document(s) is strictly prohibited. If you have received this e-mail in error, please immediately notify the sender at the above e-mail address.
# HG changeset patch
# User Adolfo Rodriguez Tsouroukdissian <adolfo.rodriguez@xxxxxxxxxxxxxxxx>
# Date 1271862957 -7200
# Node ID 8372b74161e9b26c033ead7c7fe088c62b58f8de
# Parent 3dad067a0bc4f773dd36fc2e885dd94a0ff2d64f
- Added problem size constructor to decompositions that did not have one. It preallocates member data structures.
- 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.
diff -r 3dad067a0bc4 -r 8372b74161e9 Eigen/src/Cholesky/LDLT.h
--- a/Eigen/src/Cholesky/LDLT.h Wed Apr 21 00:10:03 2010 +0300
+++ b/Eigen/src/Cholesky/LDLT.h Wed Apr 21 17:15:57 2010 +0200
@@ -66,6 +66,7 @@
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef typename ei_plain_col_type<MatrixType, int>::type IntColVectorType;
+ typedef Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> TmpMatrixType;
/** \brief Default Constructor.
*
@@ -80,12 +81,17 @@
* according to the specified problem \a size.
* \sa LDLT()
*/
- LDLT(int size) : m_matrix(size,size), m_p(size), m_transpositions(size), m_isInitialized(false) {}
+ LDLT(int size) : m_matrix(size, size),
+ m_p(size),
+ m_transpositions(size),
+ m_temporary(size),
+ m_isInitialized(false) {}
LDLT(const MatrixType& matrix)
: m_matrix(matrix.rows(), matrix.cols()),
m_p(matrix.rows()),
m_transpositions(matrix.rows()),
+ m_temporary(matrix.rows()),
m_isInitialized(false)
{
compute(matrix);
@@ -175,6 +181,7 @@
MatrixType m_matrix;
IntColVectorType m_p;
IntColVectorType m_transpositions; // FIXME do we really need to store permanently the transpositions?
+ TmpMatrixType m_temporary;
int m_sign;
bool m_isInitialized;
};
@@ -206,7 +213,7 @@
// By using a temorary, packet-aligned products are guarenteed. In the LLT
// case this is unnecessary because the diagonal is included and will always
// have optimal alignment.
- Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> _temporary(size);
+ m_temporary.resize(size);
for (int j = 0; j < size; ++j)
{
@@ -251,11 +258,11 @@
int endSize = size - j - 1;
if (endSize > 0) {
- _temporary.tail(endSize).noalias() = m_matrix.block(j+1,0, endSize, j)
+ m_temporary.tail(endSize).noalias() = m_matrix.block(j+1,0, endSize, j)
* m_matrix.col(j).head(j).conjugate();
m_matrix.row(j).tail(endSize) = m_matrix.row(j).tail(endSize).conjugate()
- - _temporary.tail(endSize).transpose();
+ - m_temporary.tail(endSize).transpose();
if(ei_abs(Djj) > cutoff)
{
diff -r 3dad067a0bc4 -r 8372b74161e9 Eigen/src/Cholesky/LLT.h
--- a/Eigen/src/Cholesky/LLT.h Wed Apr 21 00:10:03 2010 +0300
+++ b/Eigen/src/Cholesky/LLT.h Wed Apr 21 17:15:57 2010 +0200
@@ -81,6 +81,15 @@
* perform decompositions via LLT::compute(const MatrixType&).
*/
LLT() : m_matrix(), 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 LLT()
+ */
+ LLT(int size) : m_matrix(size, size),
+ m_isInitialized(false) {}
LLT(const MatrixType& matrix)
: m_matrix(matrix.rows(), matrix.cols()),
diff -r 3dad067a0bc4 -r 8372b74161e9 Eigen/src/Eigenvalues/ComplexEigenSolver.h
--- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h Wed Apr 21 00:10:03 2010 +0300
+++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h Wed Apr 21 17:15:57 2010 +0200
@@ -95,7 +95,24 @@
* 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 @@
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 @@
protected:
EigenvectorType m_eivec;
EigenvalueType m_eivalues;
+ ComplexSchur<MatrixType> m_schur;
bool m_isInitialized;
};
@@ -193,8 +212,8 @@
// 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 @@
// 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 @@
}
// 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++)
{
diff -r 3dad067a0bc4 -r 8372b74161e9 Eigen/src/Eigenvalues/ComplexSchur.h
--- a/Eigen/src/Eigenvalues/ComplexSchur.h Wed Apr 21 00:10:03 2010 +0300
+++ b/Eigen/src/Eigenvalues/ComplexSchur.h Wed Apr 21 17:15:57 2010 +0200
@@ -96,7 +96,11 @@
* \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 @@
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 @@
protected:
ComplexMatrixType m_matT, m_matU;
+ HessenbergDecomposition<MatrixType> m_hess;
bool m_isInitialized;
bool m_matUisUptodate;
@@ -300,10 +306,10 @@
// 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.
diff -r 3dad067a0bc4 -r 8372b74161e9 Eigen/src/Eigenvalues/EigenSolver.h
--- a/Eigen/src/Eigenvalues/EigenSolver.h Wed Apr 21 00:10:03 2010 +0300
+++ b/Eigen/src/Eigenvalues/EigenSolver.h Wed Apr 21 17:15:57 2010 +0200
@@ -121,6 +121,17 @@
* \sa compute() for an example.
*/
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.
*
diff -r 3dad067a0bc4 -r 8372b74161e9 Eigen/src/Eigenvalues/HessenbergDecomposition.h
--- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h Wed Apr 21 00:10:03 2010 +0300
+++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h Wed Apr 21 17:15:57 2010 +0200
@@ -91,7 +91,8 @@
* \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 @@
* \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 @@
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 @@
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 @@
* \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
diff -r 3dad067a0bc4 -r 8372b74161e9 Eigen/src/Eigenvalues/RealSchur.h
--- a/Eigen/src/Eigenvalues/RealSchur.h Wed Apr 21 00:10:03 2010 +0300
+++ b/Eigen/src/Eigenvalues/RealSchur.h Wed Apr 21 17:15:57 2010 +0200
@@ -78,6 +78,7 @@
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 @@
*/
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 @@
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 @@
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 @@
// 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.
diff -r 3dad067a0bc4 -r 8372b74161e9 Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
--- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h Wed Apr 21 00:10:03 2010 +0300
+++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h Wed Apr 21 17:15:57 2010 +0200
@@ -58,14 +58,16 @@
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 @@
*/
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 @@
*/
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 @@
protected:
MatrixType m_eivec;
RealVectorType m_eivalues;
+ typename TridiagonalizationType::SubDiagonalType m_subdiag;
#ifndef NDEBUG
bool m_eigenvectorsOk;
#endif
@@ -187,27 +194,27 @@
// 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.
diff -r 3dad067a0bc4 -r 8372b74161e9 Eigen/src/Eigenvalues/Tridiagonalization.h
--- a/Eigen/src/Eigenvalues/Tridiagonalization.h Wed Apr 21 00:10:03 2010 +0300
+++ b/Eigen/src/Eigenvalues/Tridiagonalization.h Wed Apr 21 17:15:57 2010 +0200
@@ -210,8 +210,8 @@
// 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);
diff -r 3dad067a0bc4 -r 8372b74161e9 Eigen/src/LU/FullPivLU.h
--- a/Eigen/src/LU/FullPivLU.h Wed Apr 21 00:10:03 2010 +0300
+++ b/Eigen/src/LU/FullPivLU.h Wed Apr 21 17:15:57 2010 +0200
@@ -80,6 +80,14 @@
* perform decompositions via LU::compute(const MatrixType&).
*/
FullPivLU();
+
+ /** \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 FullPivLU()
+ */
+ FullPivLU(int rows, int cols);
/** Constructor.
*
@@ -377,6 +385,8 @@
MatrixType m_lu;
PermutationPType m_p;
PermutationQType m_q;
+ IntColVectorType m_rowsTranspositions;
+ IntRowVectorType m_colsTranspositions;
int m_det_pq, m_nonzero_pivots;
RealScalar m_maxpivot, m_prescribedThreshold;
bool m_isInitialized, m_usePrescribedThreshold;
@@ -389,8 +399,26 @@
}
template<typename MatrixType>
+FullPivLU<MatrixType>::FullPivLU(int rows, int cols)
+ : m_lu(rows, cols),
+ m_p(rows),
+ m_q(cols),
+ m_rowsTranspositions(rows),
+ m_colsTranspositions(cols),
+ m_isInitialized(false),
+ m_usePrescribedThreshold(false)
+{
+}
+
+template<typename MatrixType>
FullPivLU<MatrixType>::FullPivLU(const MatrixType& matrix)
- : m_isInitialized(false), m_usePrescribedThreshold(false)
+ : m_lu(matrix.rows(), matrix.cols()),
+ m_p(matrix.rows()),
+ m_q(matrix.cols()),
+ m_rowsTranspositions(matrix.rows()),
+ m_colsTranspositions(matrix.cols()),
+ m_isInitialized(false),
+ m_usePrescribedThreshold(false)
{
compute(matrix);
}
@@ -407,9 +435,9 @@
// will store the transpositions, before we accumulate them at the end.
// can't accumulate on-the-fly because that will be done in reverse order for the rows.
- IntColVectorType rows_transpositions(matrix.rows());
- IntRowVectorType cols_transpositions(matrix.cols());
- int number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. rows_transpositions[i]!=i
+ m_rowsTranspositions.resize(matrix.rows());
+ m_colsTranspositions.resize(matrix.cols());
+ int number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
m_maxpivot = RealScalar(0);
@@ -442,8 +470,8 @@
m_nonzero_pivots = k;
for(int i = k; i < size; ++i)
{
- rows_transpositions.coeffRef(i) = i;
- cols_transpositions.coeffRef(i) = i;
+ m_rowsTranspositions.coeffRef(i) = i;
+ m_colsTranspositions.coeffRef(i) = i;
}
break;
}
@@ -453,8 +481,8 @@
// Now that we've found the pivot, we need to apply the row/col swaps to
// bring it to the location (k,k).
- rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
- cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
+ m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
+ m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
if(k != row_of_biggest_in_corner) {
m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
++number_of_transpositions;
@@ -478,11 +506,11 @@
m_p.setIdentity(rows);
for(int k = size-1; k >= 0; --k)
- m_p.applyTranspositionOnTheRight(k, rows_transpositions.coeff(k));
+ m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
m_q.setIdentity(cols);
for(int k = 0; k < size; ++k)
- m_q.applyTranspositionOnTheRight(k, cols_transpositions.coeff(k));
+ m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
m_det_pq = (number_of_transpositions%2) ? -1 : 1;
return *this;
diff -r 3dad067a0bc4 -r 8372b74161e9 Eigen/src/LU/PartialPivLU.h
--- a/Eigen/src/LU/PartialPivLU.h Wed Apr 21 00:10:03 2010 +0300
+++ b/Eigen/src/LU/PartialPivLU.h Wed Apr 21 17:15:57 2010 +0200
@@ -82,6 +82,14 @@
* perform decompositions via PartialPivLU::compute(const MatrixType&).
*/
PartialPivLU();
+
+ /** \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 PartialPivLU()
+ */
+ PartialPivLU(int size);
/** Constructor.
*
@@ -176,6 +184,7 @@
protected:
MatrixType m_lu;
PermutationType m_p;
+ PermutationVectorType m_rowsTranspositions;
int m_det_p;
bool m_isInitialized;
};
@@ -184,6 +193,17 @@
PartialPivLU<MatrixType>::PartialPivLU()
: m_lu(),
m_p(),
+ m_rowsTranspositions(),
+ m_det_p(0),
+ m_isInitialized(false)
+{
+}
+
+template<typename MatrixType>
+PartialPivLU<MatrixType>::PartialPivLU(int size)
+ : m_lu(size, size),
+ m_p(size),
+ m_rowsTranspositions(size),
m_det_p(0),
m_isInitialized(false)
{
@@ -191,8 +211,9 @@
template<typename MatrixType>
PartialPivLU<MatrixType>::PartialPivLU(const MatrixType& matrix)
- : m_lu(),
- m_p(),
+ : m_lu(matrix.rows(), matrix.rows()),
+ m_p(matrix.rows()),
+ m_rowsTranspositions(matrix.rows()),
m_det_p(0),
m_isInitialized(false)
{
@@ -384,15 +405,15 @@
ei_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
const int size = matrix.rows();
- PermutationVectorType rows_transpositions(size);
+ m_rowsTranspositions.resize(size);
int nb_transpositions;
- ei_partial_lu_inplace(m_lu, rows_transpositions, nb_transpositions);
+ ei_partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
m_det_p = (nb_transpositions%2) ? -1 : 1;
m_p.setIdentity(size);
for(int k = size-1; k >= 0; --k)
- m_p.applyTranspositionOnTheRight(k, rows_transpositions.coeff(k));
+ m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
m_isInitialized = true;
return *this;
diff -r 3dad067a0bc4 -r 8372b74161e9 Eigen/src/QR/ColPivHouseholderQR.h
--- a/Eigen/src/QR/ColPivHouseholderQR.h Wed Apr 21 00:10:03 2010 +0300
+++ b/Eigen/src/QR/ColPivHouseholderQR.h Wed Apr 21 17:15:57 2010 +0200
@@ -70,11 +70,38 @@
* The default constructor is useful in cases in which the user intends to
* perform decompositions via ColPivHouseholderQR::compute(const MatrixType&).
*/
- ColPivHouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
+ ColPivHouseholderQR()
+ : m_qr(),
+ m_hCoeffs(),
+ m_colsPermutation(),
+ m_colsTranspositions(),
+ m_temp(),
+ m_colSqNorms(),
+ 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 ColPivHouseholderQR()
+ */
+ ColPivHouseholderQR(int rows, int cols)
+ : m_qr(rows, cols),
+ m_hCoeffs(std::min(rows,cols)),
+ m_colsPermutation(cols),
+ m_colsTranspositions(cols),
+ m_temp(cols),
+ m_colSqNorms(cols),
+ m_isInitialized(false),
+ m_usePrescribedThreshold(false) {}
ColPivHouseholderQR(const MatrixType& matrix)
: m_qr(matrix.rows(), matrix.cols()),
m_hCoeffs(std::min(matrix.rows(),matrix.cols())),
+ m_colsPermutation(matrix.cols()),
+ m_colsTranspositions(matrix.cols()),
+ m_temp(matrix.cols()),
+ m_colSqNorms(matrix.cols()),
m_isInitialized(false),
m_usePrescribedThreshold(false)
{
@@ -121,7 +148,7 @@
const PermutationType& colsPermutation() const
{
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
- return m_cols_permutation;
+ return m_colsPermutation;
}
/** \returns the absolute value of the determinant of the matrix of which
@@ -307,7 +334,10 @@
protected:
MatrixType m_qr;
HCoeffsType m_hCoeffs;
- PermutationType m_cols_permutation;
+ PermutationType m_colsPermutation;
+ IntRowVectorType m_colsTranspositions;
+ RowVectorType m_temp;
+ RealRowVectorType m_colSqNorms;
bool m_isInitialized, m_usePrescribedThreshold;
RealScalar m_prescribedThreshold, m_maxpivot;
int m_nonzero_pivots;
@@ -342,35 +372,35 @@
m_qr = matrix;
m_hCoeffs.resize(size);
- RowVectorType temp(cols);
+ m_temp.resize(cols);
- IntRowVectorType cols_transpositions(matrix.cols());
+ m_colsTranspositions.resize(matrix.cols());
int number_of_transpositions = 0;
- RealRowVectorType colSqNorms(cols);
+ m_colSqNorms.resize(cols);
for(int k = 0; k < cols; ++k)
- colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
+ m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
- RealScalar threshold_helper = colSqNorms.maxCoeff() * ei_abs2(NumTraits<Scalar>::epsilon()) / rows;
+ RealScalar threshold_helper = m_colSqNorms.maxCoeff() * ei_abs2(NumTraits<Scalar>::epsilon()) / rows;
m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
m_maxpivot = RealScalar(0);
for(int k = 0; k < size; ++k)
{
- // first, we look up in our table colSqNorms which column has the biggest squared norm
+ // first, we look up in our table m_colSqNorms which column has the biggest squared norm
int biggest_col_index;
- RealScalar biggest_col_sq_norm = colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
+ RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
biggest_col_index += k;
- // since our table colSqNorms accumulates imprecision at every step, we must now recompute
+ // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute
// the actual squared norm of the selected column.
// Note that not doing so does result in solve() sometimes returning inf/nan values
// when running the unit test with 1000 repetitions.
biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
// we store that back into our table: it can't hurt to correct our table.
- colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
+ m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
// if the current biggest column is smaller than epsilon times the initial biggest column,
// terminate to avoid generating nan/inf values.
@@ -388,10 +418,10 @@
}
// apply the transposition to the columns
- cols_transpositions.coeffRef(k) = biggest_col_index;
+ m_colsTranspositions.coeffRef(k) = biggest_col_index;
if(k != biggest_col_index) {
m_qr.col(k).swap(m_qr.col(biggest_col_index));
- std::swap(colSqNorms.coeffRef(k), colSqNorms.coeffRef(biggest_col_index));
+ std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
++number_of_transpositions;
}
@@ -407,15 +437,15 @@
// apply the householder transformation
m_qr.corner(BottomRight, rows-k, cols-k-1)
- .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1));
+ .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
// update our table of squared norms of the columns
- colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
+ m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
}
- m_cols_permutation.setIdentity(cols);
+ m_colsPermutation.setIdentity(cols);
for(int k = 0; k < m_nonzero_pivots; ++k)
- m_cols_permutation.applyTranspositionOnTheRight(k, cols_transpositions.coeff(k));
+ m_colsPermutation.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
m_det_pq = (number_of_transpositions%2) ? -1 : 1;
m_isInitialized = true;
diff -r 3dad067a0bc4 -r 8372b74161e9 Eigen/src/QR/FullPivHouseholderQR.h
--- a/Eigen/src/QR/FullPivHouseholderQR.h Wed Apr 21 00:10:03 2010 +0300
+++ b/Eigen/src/QR/FullPivHouseholderQR.h Wed Apr 21 17:15:57 2010 +0200
@@ -69,10 +69,38 @@
* The default constructor is useful in cases in which the user intends to
* perform decompositions via FullPivHouseholderQR::compute(const MatrixType&).
*/
- FullPivHouseholderQR() : m_isInitialized(false) {}
+ FullPivHouseholderQR()
+ : m_qr(),
+ m_hCoeffs(),
+ m_rows_transpositions(),
+ m_cols_transpositions(),
+ m_cols_permutation(),
+ m_temp(),
+ 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 FullPivHouseholderQR()
+ */
+ FullPivHouseholderQR(int rows, int cols)
+ : m_qr(rows, cols),
+ m_hCoeffs(std::min(rows,cols)),
+ m_rows_transpositions(rows),
+ m_cols_transpositions(cols),
+ m_cols_permutation(cols),
+ m_temp(std::min(rows,cols)),
+ m_isInitialized(false) {}
FullPivHouseholderQR(const MatrixType& matrix)
- : m_isInitialized(false)
+ : m_qr(matrix.rows(), matrix.cols()),
+ m_hCoeffs(std::min(matrix.rows(), matrix.cols())),
+ m_rows_transpositions(matrix.rows()),
+ m_cols_transpositions(matrix.cols()),
+ m_cols_permutation(matrix.cols()),
+ m_temp(std::min(matrix.rows(), matrix.cols())),
+ m_isInitialized(false)
{
compute(matrix);
}
@@ -233,7 +261,9 @@
MatrixType m_qr;
HCoeffsType m_hCoeffs;
IntColVectorType m_rows_transpositions;
+ IntRowVectorType m_cols_transpositions;
PermutationType m_cols_permutation;
+ RowVectorType m_temp;
bool m_isInitialized;
RealScalar m_precision;
int m_rank;
@@ -269,12 +299,12 @@
m_qr = matrix;
m_hCoeffs.resize(size);
- RowVectorType temp(cols);
+ m_temp.resize(cols);
m_precision = NumTraits<Scalar>::epsilon() * size;
m_rows_transpositions.resize(matrix.rows());
- IntRowVectorType cols_transpositions(matrix.cols());
+ m_cols_transpositions.resize(matrix.cols());
int number_of_transpositions = 0;
RealScalar biggest(0);
@@ -298,14 +328,14 @@
for(int i = k; i < size; i++)
{
m_rows_transpositions.coeffRef(i) = i;
- cols_transpositions.coeffRef(i) = i;
+ m_cols_transpositions.coeffRef(i) = i;
m_hCoeffs.coeffRef(i) = Scalar(0);
}
break;
}
m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
- cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
+ m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
if(k != row_of_biggest_in_corner) {
m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
++number_of_transpositions;
@@ -320,12 +350,12 @@
m_qr.coeffRef(k,k) = beta;
m_qr.corner(BottomRight, rows-k, cols-k-1)
- .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1));
+ .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
}
m_cols_permutation.setIdentity(cols);
for(int k = 0; k < size; ++k)
- m_cols_permutation.applyTranspositionOnTheRight(k, cols_transpositions.coeff(k));
+ m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
m_det_pq = (number_of_transpositions%2) ? -1 : 1;
m_isInitialized = true;
diff -r 3dad067a0bc4 -r 8372b74161e9 Eigen/src/QR/HouseholderQR.h
--- a/Eigen/src/QR/HouseholderQR.h Wed Apr 21 00:10:03 2010 +0300
+++ b/Eigen/src/QR/HouseholderQR.h Wed Apr 21 17:15:57 2010 +0200
@@ -71,11 +71,24 @@
* The default constructor is useful in cases in which the user intends to
* perform decompositions via HouseholderQR::compute(const MatrixType&).
*/
- HouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
+ HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), 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 HouseholderQR()
+ */
+ HouseholderQR(int rows, int cols)
+ : m_qr(rows, cols),
+ m_hCoeffs(std::min(rows,cols)),
+ m_temp(cols),
+ m_isInitialized(false) {}
HouseholderQR(const MatrixType& matrix)
: m_qr(matrix.rows(), matrix.cols()),
m_hCoeffs(std::min(matrix.rows(),matrix.cols())),
+ m_temp(matrix.cols()),
m_isInitialized(false)
{
compute(matrix);
@@ -159,6 +172,7 @@
protected:
MatrixType m_qr;
HCoeffsType m_hCoeffs;
+ RowVectorType m_temp;
bool m_isInitialized;
};
@@ -190,7 +204,7 @@
m_qr = matrix;
m_hCoeffs.resize(size);
- RowVectorType temp(cols);
+ m_temp.resize(cols);
for(int k = 0; k < size; ++k)
{
@@ -203,7 +217,7 @@
// apply H to remaining part of m_qr from the left
m_qr.corner(BottomRight, remainingRows, remainingCols)
- .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingRows-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1));
+ .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingRows-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
}
m_isInitialized = true;
return *this;
diff -r 3dad067a0bc4 -r 8372b74161e9 Eigen/src/SVD/JacobiSVD.h
--- a/Eigen/src/SVD/JacobiSVD.h Wed Apr 21 00:10:03 2010 +0300
+++ b/Eigen/src/SVD/JacobiSVD.h Wed Apr 21 17:15:57 2010 +0200
@@ -93,10 +93,35 @@
public:
+ /** \brief Default Constructor.
+ *
+ * The default constructor is useful in cases in which the user intends to
+ * perform decompositions via JacobiSVD::compute(const MatrixType&).
+ */
JacobiSVD() : m_isInitialized(false) {}
- JacobiSVD(const MatrixType& matrix) : 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 JacobiSVD()
+ */
+ JacobiSVD(int rows, int cols) : m_matrixU(rows, rows),
+ m_matrixV(cols, cols),
+ m_singularValues(std::min(rows, cols)),
+ m_workMatrix(rows, cols),
+ m_isInitialized(false) {}
+
+ JacobiSVD(const MatrixType& matrix) : m_matrixU(matrix.rows(), matrix.rows()),
+ m_matrixV(matrix.cols(), matrix.cols()),
+ m_singularValues(),
+ m_workMatrix(),
+ m_isInitialized(false)
{
+ const int minSize = std::min(matrix.rows(), matrix.cols());
+ m_singularValues.resize(minSize);
+ m_workMatrix.resize(minSize, minSize);
compute(matrix);
}
@@ -124,6 +149,7 @@
MatrixUType m_matrixU;
MatrixVType m_matrixV;
SingularValuesType m_singularValues;
+ WorkMatrixType m_workMatrix;
bool m_isInitialized;
template<typename _MatrixType, unsigned int _Options, bool _IsComplex>
@@ -289,17 +315,16 @@
template<typename MatrixType, unsigned int Options>
JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute(const MatrixType& matrix)
{
- WorkMatrixType work_matrix;
int rows = matrix.rows();
int cols = matrix.cols();
int diagSize = std::min(rows, cols);
m_singularValues.resize(diagSize);
const RealScalar precision = 2 * NumTraits<Scalar>::epsilon();
- if(!ei_svd_precondition_if_more_rows_than_cols<MatrixType, Options>::run(matrix, work_matrix, *this)
- && !ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options>::run(matrix, work_matrix, *this))
+ if(!ei_svd_precondition_if_more_rows_than_cols<MatrixType, Options>::run(matrix, m_workMatrix, *this)
+ && !ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options>::run(matrix, m_workMatrix, *this))
{
- work_matrix = matrix.block(0,0,diagSize,diagSize);
+ m_workMatrix = matrix.block(0,0,diagSize,diagSize);
if(ComputeU) m_matrixU.setIdentity(rows,rows);
if(ComputeV) m_matrixV.setIdentity(cols,cols);
}
@@ -312,19 +337,19 @@
{
for(int q = 0; q < p; ++q)
{
- if(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p)))
- > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q)))*precision)
+ if(std::max(ei_abs(m_workMatrix.coeff(p,q)),ei_abs(m_workMatrix.coeff(q,p)))
+ > std::max(ei_abs(m_workMatrix.coeff(p,p)),ei_abs(m_workMatrix.coeff(q,q)))*precision)
{
finished = false;
- ei_svd_precondition_2x2_block_to_be_real<MatrixType, Options>::run(work_matrix, *this, p, q);
+ ei_svd_precondition_2x2_block_to_be_real<MatrixType, Options>::run(m_workMatrix, *this, p, q);
PlanarRotation<RealScalar> j_left, j_right;
- ei_real_2x2_jacobi_svd(work_matrix, p, q, &j_left, &j_right);
+ ei_real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
- work_matrix.applyOnTheLeft(p,q,j_left);
+ m_workMatrix.applyOnTheLeft(p,q,j_left);
if(ComputeU) m_matrixU.applyOnTheRight(p,q,j_left.transpose());
- work_matrix.applyOnTheRight(p,q,j_right);
+ m_workMatrix.applyOnTheRight(p,q,j_right);
if(ComputeV) m_matrixV.applyOnTheRight(p,q,j_right);
}
}
@@ -333,9 +358,9 @@
for(int i = 0; i < diagSize; ++i)
{
- RealScalar a = ei_abs(work_matrix.coeff(i,i));
+ RealScalar a = ei_abs(m_workMatrix.coeff(i,i));
m_singularValues.coeffRef(i) = a;
- if(ComputeU && (a!=RealScalar(0))) m_matrixU.col(i) *= work_matrix.coeff(i,i)/a;
+ if(ComputeU && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
}
for(int i = 0; i < diagSize; i++)
diff -r 3dad067a0bc4 -r 8372b74161e9 Eigen/src/SVD/SVD.h
--- a/Eigen/src/SVD/SVD.h Wed Apr 21 00:10:03 2010 +0300
+++ b/Eigen/src/SVD/SVD.h Wed Apr 21 17:15:57 2010 +0200
@@ -73,11 +73,25 @@
*/
SVD() : m_matU(), m_matV(), m_sigma(), m_isInitialized(false) {}
- SVD(const MatrixType& matrix)
- : m_matU(matrix.rows(), matrix.rows()),
- m_matV(matrix.cols(),matrix.cols()),
- m_sigma(matrix.cols()),
- 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 JacobiSVD()
+ */
+ SVD(int rows, int cols) : m_matU(rows, rows),
+ m_matV(cols,cols),
+ m_sigma(std::min(rows, cols)),
+ m_workMatrix(rows, cols),
+ m_rv1(cols),
+ m_isInitialized(false) {}
+
+ SVD(const MatrixType& matrix) : m_matU(matrix.rows(), matrix.rows()),
+ m_matV(matrix.cols(),matrix.cols()),
+ m_sigma(std::min(matrix.rows(), matrix.cols())),
+ m_workMatrix(matrix.rows(), matrix.cols()),
+ m_rv1(matrix.cols()),
+ m_isInitialized(false)
{
compute(matrix);
}
@@ -165,6 +179,8 @@
MatrixVType m_matV;
/** \internal */
SingularValuesType m_sigma;
+ MatrixType m_workMatrix;
+ RowVector m_rv1;
bool m_isInitialized;
int m_rows, m_cols;
};
@@ -185,11 +201,12 @@
m_matU.setZero();
m_sigma.resize(n);
m_matV.resize(n,n);
+ m_workMatrix = matrix;
int max_iters = 30;
MatrixVType& V = m_matV;
- MatrixType A = matrix;
+ MatrixType& A = m_workMatrix;
SingularValuesType& W = m_sigma;
bool flag;
@@ -198,14 +215,14 @@
bool convergence = true;
Scalar eps = NumTraits<Scalar>::dummy_precision();
- RowVector rv1(n);
+ m_rv1.resize(n);
g = scale = anorm = 0;
// Householder reduction to bidiagonal form.
for (i=0; i<n; i++)
{
l = i+2;
- rv1[i] = scale*g;
+ m_rv1[i] = scale*g;
g = s = scale = 0.0;
if (i < m)
{
@@ -246,16 +263,16 @@
g = -sign(ei_sqrt(s),f);
h = f*g - s;
A(i,l-1) = f-g;
- rv1.tail(n-l+1) = A.row(i).tail(n-l+1)/h;
+ m_rv1.tail(n-l+1) = A.row(i).tail(n-l+1)/h;
for (j=l-1; j<m; j++)
{
s = A.row(i).tail(n-l+1).dot(A.row(j).tail(n-l+1));
- A.row(j).tail(n-l+1) += s*rv1.tail(n-l+1).transpose();
+ A.row(j).tail(n-l+1) += s*m_rv1.tail(n-l+1).transpose();
}
A.row(i).tail(n-l+1) *= scale;
}
}
- anorm = std::max( anorm, (ei_abs(W[i])+ei_abs(rv1[i])) );
+ anorm = std::max( anorm, (ei_abs(W[i])+ei_abs(m_rv1[i])) );
}
// Accumulation of right-hand transformations.
for (i=n-1; i>=0; i--)
@@ -277,7 +294,7 @@
V.col(i).tail(n-l).setZero();
}
V(i, i) = 1.0;
- g = rv1[i];
+ g = m_rv1[i];
l = i;
}
// Accumulation of left-hand transformations.
@@ -318,7 +335,7 @@
nm = l-1;
// Note that rv1[1] is always zero.
//if ((double)(ei_abs(rv1[l])+anorm) == anorm)
- if (l==0 || ei_abs(rv1[l]) <= eps*anorm)
+ if (l==0 || ei_abs(m_rv1[l]) <= eps*anorm)
{
flag = false;
break;
@@ -333,8 +350,8 @@
s = 1.0;
for (i=l ;i<k+1; i++)
{
- f = s*rv1[i];
- rv1[i] = c*rv1[i];
+ f = s*m_rv1[i];
+ m_rv1[i] = c*m_rv1[i];
//if ((double)(ei_abs(f)+anorm) == anorm)
if (ei_abs(f) <= eps*anorm)
break;
@@ -363,8 +380,8 @@
x = W[l]; // Shift from bottom 2-by-2 minor.
nm = k-1;
y = W[nm];
- g = rv1[nm];
- h = rv1[k];
+ g = m_rv1[nm];
+ h = m_rv1[k];
f = ((y-z)*(y+z) + (g-h)*(g+h))/(Scalar(2.0)*h*y);
g = pythag(f,1.0);
f = ((x-z)*(x+z) + h*((y/(f+sign(g,f)))-h))/x;
@@ -373,13 +390,13 @@
for (j=l; j<=nm; j++)
{
i = j+1;
- g = rv1[i];
+ g = m_rv1[i];
y = W[i];
h = s*g;
g = c*g;
z = pythag(f,h);
- rv1[j] = z;
+ m_rv1[j] = z;
c = f/z;
s = h/z;
f = x*c + g*s;
@@ -401,8 +418,8 @@
x = c*y - s*g;
A.applyOnTheRight(i,j,PlanarRotation<Scalar>(c,s));
}
- rv1[l] = 0.0;
- rv1[k] = f;
+ m_rv1[l] = 0.0;
+ m_rv1[k] = f;
W[k] = x;
}
}
diff -r 3dad067a0bc4 -r 8372b74161e9 test/cholesky.cpp
--- a/test/cholesky.cpp Wed Apr 21 00:10:03 2010 +0300
+++ b/test/cholesky.cpp Wed Apr 21 17:15:57 2010 +0200
@@ -163,4 +163,8 @@
CALL_SUBTEST_7( cholesky_verify_assert<Matrix3d>() );
CALL_SUBTEST_8( cholesky_verify_assert<MatrixXf>() );
CALL_SUBTEST_2( cholesky_verify_assert<MatrixXd>() );
+
+ // Test problem size constructors
+ CALL_SUBTEST_9( LLT<MatrixXf>(10) );
+ CALL_SUBTEST_9( LDLT<MatrixXf>(10) );
}
diff -r 3dad067a0bc4 -r 8372b74161e9 test/eigensolver_complex.cpp
--- a/test/eigensolver_complex.cpp Wed Apr 21 00:10:03 2010 +0300
+++ b/test/eigensolver_complex.cpp Wed Apr 21 17:15:57 2010 +0200
@@ -63,4 +63,7 @@
CALL_SUBTEST_3( eigensolver(Matrix<std::complex<float>, 1, 1>()) );
CALL_SUBTEST_4( eigensolver(Matrix3f()) );
}
+
+ // Test problem size constructors
+ CALL_SUBTEST_5(ComplexEigenSolver<MatrixXf>(10));
}
diff -r 3dad067a0bc4 -r 8372b74161e9 test/eigensolver_generic.cpp
--- a/test/eigensolver_generic.cpp Wed Apr 21 00:10:03 2010 +0300
+++ b/test/eigensolver_generic.cpp Wed Apr 21 17:15:57 2010 +0200
@@ -89,4 +89,7 @@
CALL_SUBTEST_2( eigensolver_verify_assert<MatrixXd>() );
CALL_SUBTEST_4( eigensolver_verify_assert<Matrix2d>() );
CALL_SUBTEST_5( eigensolver_verify_assert<MatrixXf>() );
+
+ // Test problem size constructors
+ CALL_SUBTEST_6(EigenSolver<MatrixXf>(10));
}
diff -r 3dad067a0bc4 -r 8372b74161e9 test/eigensolver_selfadjoint.cpp
--- a/test/eigensolver_selfadjoint.cpp Wed Apr 21 00:10:03 2010 +0300
+++ b/test/eigensolver_selfadjoint.cpp Wed Apr 21 17:15:57 2010 +0200
@@ -129,5 +129,9 @@
CALL_SUBTEST_6( selfadjointeigensolver(Matrix<double,1,1>()) );
CALL_SUBTEST_7( selfadjointeigensolver(Matrix<double,2,2>()) );
}
+
+ // Test problem size constructors
+ CALL_SUBTEST_8(SelfAdjointEigenSolver<MatrixXf>(10));
+ CALL_SUBTEST_8(Tridiagonalization<MatrixXf>(10));
}
diff -r 3dad067a0bc4 -r 8372b74161e9 test/hessenberg.cpp
--- a/test/hessenberg.cpp Wed Apr 21 00:10:03 2010 +0300
+++ b/test/hessenberg.cpp Wed Apr 21 17:15:57 2010 +0200
@@ -61,4 +61,7 @@
CALL_SUBTEST_3(( hessenberg<std::complex<float>,4>() ));
CALL_SUBTEST_4(( hessenberg<float,Dynamic>(ei_random<int>(1,320)) ));
CALL_SUBTEST_5(( hessenberg<std::complex<double>,Dynamic>(ei_random<int>(1,320)) ));
+
+ // Test problem size constructors
+ CALL_SUBTEST_6(HessenbergDecomposition<MatrixXf>(10));
}
diff -r 3dad067a0bc4 -r 8372b74161e9 test/jacobisvd.cpp
--- a/test/jacobisvd.cpp Wed Apr 21 00:10:03 2010 +0300
+++ b/test/jacobisvd.cpp Wed Apr 21 17:15:57 2010 +0200
@@ -106,4 +106,7 @@
CALL_SUBTEST_3(( svd_verify_assert<Matrix3d>() ));
CALL_SUBTEST_9(( svd_verify_assert<MatrixXf>() ));
CALL_SUBTEST_11(( svd_verify_assert<MatrixXd>() ));
+
+ // Test problem size constructors
+ CALL_SUBTEST_12( JacobiSVD<MatrixXf>(10, 20) );
}
diff -r 3dad067a0bc4 -r 8372b74161e9 test/lu.cpp
--- a/test/lu.cpp Wed Apr 21 00:10:03 2010 +0300
+++ b/test/lu.cpp Wed Apr 21 17:15:57 2010 +0200
@@ -212,5 +212,9 @@
CALL_SUBTEST_6( lu_verify_assert<MatrixXcd>() );
CALL_SUBTEST_7(( lu_non_invertible<Matrix<float,Dynamic,16> >() ));
+
+ // Test problem size constructors
+ CALL_SUBTEST_9( PartialPivLU<MatrixXf>(10) );
+ CALL_SUBTEST_9( FullPivLU<MatrixXf>(10, 20); );
}
}
diff -r 3dad067a0bc4 -r 8372b74161e9 test/qr.cpp
--- a/test/qr.cpp Wed Apr 21 00:10:03 2010 +0300
+++ b/test/qr.cpp Wed Apr 21 17:15:57 2010 +0200
@@ -133,4 +133,7 @@
CALL_SUBTEST_6(qr_verify_assert<MatrixXd>());
CALL_SUBTEST_7(qr_verify_assert<MatrixXcf>());
CALL_SUBTEST_8(qr_verify_assert<MatrixXcd>());
+
+ // Test problem size constructors
+ CALL_SUBTEST_12(HouseholderQR<MatrixXf>(10, 20));
}
diff -r 3dad067a0bc4 -r 8372b74161e9 test/qr_colpivoting.cpp
--- a/test/qr_colpivoting.cpp Wed Apr 21 00:10:03 2010 +0300
+++ b/test/qr_colpivoting.cpp Wed Apr 21 17:15:57 2010 +0200
@@ -157,4 +157,7 @@
CALL_SUBTEST_2(qr_verify_assert<MatrixXd>());
CALL_SUBTEST_6(qr_verify_assert<MatrixXcf>());
CALL_SUBTEST_3(qr_verify_assert<MatrixXcd>());
+
+ // Test problem size constructors
+ CALL_SUBTEST_9(ColPivHouseholderQR<MatrixXf>(10, 20));
}
diff -r 3dad067a0bc4 -r 8372b74161e9 test/qr_fullpivoting.cpp
--- a/test/qr_fullpivoting.cpp Wed Apr 21 00:10:03 2010 +0300
+++ b/test/qr_fullpivoting.cpp Wed Apr 21 17:15:57 2010 +0200
@@ -139,4 +139,7 @@
CALL_SUBTEST_2(qr_verify_assert<MatrixXd>());
CALL_SUBTEST_4(qr_verify_assert<MatrixXcf>());
CALL_SUBTEST_3(qr_verify_assert<MatrixXcd>());
+
+ // Test problem size constructors
+ CALL_SUBTEST_7(FullPivHouseholderQR<MatrixXf>(10, 20));
}
diff -r 3dad067a0bc4 -r 8372b74161e9 test/schur_complex.cpp
--- a/test/schur_complex.cpp Wed Apr 21 00:10:03 2010 +0300
+++ b/test/schur_complex.cpp Wed Apr 21 17:15:57 2010 +0200
@@ -64,4 +64,7 @@
CALL_SUBTEST_2(( schur<MatrixXcf>(ei_random<int>(1,50)) ));
CALL_SUBTEST_3(( schur<Matrix<std::complex<float>, 1, 1> >() ));
CALL_SUBTEST_4(( schur<Matrix<float, 3, 3, Eigen::RowMajor> >() ));
+
+ // Test problem size constructors
+ CALL_SUBTEST_5(ComplexSchur<MatrixXf>(10));
}
diff -r 3dad067a0bc4 -r 8372b74161e9 test/schur_real.cpp
--- a/test/schur_real.cpp Wed Apr 21 00:10:03 2010 +0300
+++ b/test/schur_real.cpp Wed Apr 21 17:15:57 2010 +0200
@@ -82,4 +82,7 @@
CALL_SUBTEST_2(( schur<MatrixXd>(ei_random<int>(1,50)) ));
CALL_SUBTEST_3(( schur<Matrix<float, 1, 1> >() ));
CALL_SUBTEST_4(( schur<Matrix<double, 3, 3, Eigen::RowMajor> >() ));
+
+ // Test problem size constructors
+ CALL_SUBTEST_5(RealSchur<MatrixXf>(10));
}
diff -r 3dad067a0bc4 -r 8372b74161e9 test/svd.cpp
--- a/test/svd.cpp Wed Apr 21 00:10:03 2010 +0300
+++ b/test/svd.cpp Wed Apr 21 17:15:57 2010 +0200
@@ -113,4 +113,7 @@
CALL_SUBTEST_2( svd_verify_assert<Matrix4d>() );
CALL_SUBTEST_3( svd_verify_assert<MatrixXf>() );
CALL_SUBTEST_4( svd_verify_assert<MatrixXd>() );
+
+ // Test problem size constructors
+ CALL_SUBTEST_9( SVD<MatrixXf>(10, 20) );
}