aboutsummaryrefslogtreecommitdiffhomepage
path: root/Eigen
diff options
context:
space:
mode:
Diffstat (limited to 'Eigen')
-rw-r--r--Eigen/src/Cholesky/LDLT.h15
-rw-r--r--Eigen/src/Cholesky/LLT.h9
-rw-r--r--Eigen/src/Eigenvalues/ComplexEigenSolver.h33
-rw-r--r--Eigen/src/Eigenvalues/ComplexSchur.h14
-rw-r--r--Eigen/src/Eigenvalues/EigenSolver.h11
-rw-r--r--Eigen/src/Eigenvalues/HessenbergDecomposition.h19
-rw-r--r--Eigen/src/Eigenvalues/RealSchur.h20
-rw-r--r--Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h29
-rw-r--r--Eigen/src/Eigenvalues/Tridiagonalization.h4
-rw-r--r--Eigen/src/LU/FullPivLU.h48
-rw-r--r--Eigen/src/LU/PartialPivLU.h31
-rw-r--r--Eigen/src/QR/ColPivHouseholderQR.h66
-rw-r--r--Eigen/src/QR/FullPivHouseholderQR.h46
-rw-r--r--Eigen/src/QR/HouseholderQR.h20
-rw-r--r--Eigen/src/SVD/JacobiSVD.h51
-rw-r--r--Eigen/src/SVD/SVD.h59
16 files changed, 354 insertions, 121 deletions
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index f92a72c7b..206ccef4d 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -66,6 +66,7 @@ template<typename _MatrixType> class LDLT
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 @@ template<typename _MatrixType> class LDLT
* 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 @@ template<typename _MatrixType> class LDLT
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 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
// 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 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
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 --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index 8b01e0ccb..22d0c91c8 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -82,6 +82,15 @@ template<typename _MatrixType, int _UpLo> class LLT
*/
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()),
m_isInitialized(false)
diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
index 0b2e0b0ec..8e5f1310a 100644
--- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h
+++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
@@ -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++)
{
diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h
index 7c9ab03c8..8874aa257 100644
--- a/Eigen/src/Eigenvalues/ComplexSchur.h
+++ b/Eigen/src/Eigenvalues/ComplexSchur.h
@@ -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.
diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h
index 32dee92e9..114bfab6f 100644
--- a/Eigen/src/Eigenvalues/EigenSolver.h
+++ b/Eigen/src/Eigenvalues/EigenSolver.h
@@ -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.
diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
index 3e8c4c64d..451604252 100644
--- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h
+++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
@@ -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
diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h
index a66ab8ace..dd4a1ab7e 100644
--- a/Eigen/src/Eigenvalues/RealSchur.h
+++ b/Eigen/src/Eigenvalues/RealSchur.h
@@ -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.
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
index b31f77e9b..41a600290 100644
--- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
+++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -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.
diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h
index a0b772888..c320242ef 100644
--- a/Eigen/src/Eigenvalues/Tridiagonalization.h
+++ b/Eigen/src/Eigenvalues/Tridiagonalization.h
@@ -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);
diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index 1fe15bfd2..241fb9d16 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -81,6 +81,14 @@ template<typename _MatrixType> class FullPivLU
*/
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.
*
* \param matrix the matrix of which to compute the LU decomposition.
@@ -377,6 +385,8 @@ template<typename _MatrixType> class FullPivLU
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 @@ FullPivLU<MatrixType>::FullPivLU()
}
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 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
// 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 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
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 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
// 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 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
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 --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h
index 5062fd0ce..42c74ee53 100644
--- a/Eigen/src/LU/PartialPivLU.h
+++ b/Eigen/src/LU/PartialPivLU.h
@@ -83,6 +83,14 @@ template<typename _MatrixType> class PartialPivLU
*/
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.
*
* \param matrix the matrix of which to compute the LU decomposition.
@@ -176,6 +184,7 @@ template<typename _MatrixType> class PartialPivLU
protected:
MatrixType m_lu;
PermutationType m_p;
+ PermutationVectorType m_rowsTranspositions;
int m_det_p;
bool m_isInitialized;
};
@@ -184,6 +193,17 @@ template<typename MatrixType>
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 @@ PartialPivLU<MatrixType>::PartialPivLU()
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 @@ PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& ma
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 --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h
index 27f80e046..5893125cd 100644
--- a/Eigen/src/QR/ColPivHouseholderQR.h
+++ b/Eigen/src/QR/ColPivHouseholderQR.h
@@ -70,11 +70,38 @@ template<typename _MatrixType> class ColPivHouseholderQR
* 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 @@ template<typename _MatrixType> class ColPivHouseholderQR
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 @@ template<typename _MatrixType> class ColPivHouseholderQR
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 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
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 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
}
// 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 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
// 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 --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h
index aeb551cc7..8df2ed1e0 100644
--- a/Eigen/src/QR/FullPivHouseholderQR.h
+++ b/Eigen/src/QR/FullPivHouseholderQR.h
@@ -69,10 +69,38 @@ template<typename _MatrixType> class FullPivHouseholderQR
* 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 @@ template<typename _MatrixType> class FullPivHouseholderQR
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 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
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 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
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 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
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 --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h
index c4abc16ff..10b9fb93f 100644
--- a/Eigen/src/QR/HouseholderQR.h
+++ b/Eigen/src/QR/HouseholderQR.h
@@ -71,11 +71,24 @@ template<typename _MatrixType> class HouseholderQR
* 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 @@ template<typename _MatrixType> class HouseholderQR
protected:
MatrixType m_qr;
HCoeffsType m_hCoeffs;
+ RowVectorType m_temp;
bool m_isInitialized;
};
@@ -190,7 +204,7 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
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 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
// 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 --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h
index b7cab023b..9323c0180 100644
--- a/Eigen/src/SVD/JacobiSVD.h
+++ b/Eigen/src/SVD/JacobiSVD.h
@@ -93,10 +93,35 @@ template<typename MatrixType, unsigned int Options> class JacobiSVD
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 @@ template<typename MatrixType, unsigned int Options> class JacobiSVD
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 @@ struct ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options, true>
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 @@ JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute(const Ma
{
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 @@ JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute(const Ma
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 --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h
index ed0e69f91..a9e22dfd4 100644
--- a/Eigen/src/SVD/SVD.h
+++ b/Eigen/src/SVD/SVD.h
@@ -73,11 +73,25 @@ template<typename _MatrixType> class SVD
*/
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 @@ template<typename _MatrixType> class SVD
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 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
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 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
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 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
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 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
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 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
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 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
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 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
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 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
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 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
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;
}
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-18 03:11:40 +0000