GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-04 18:54:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

In SparseQR, calling factorize() without analyzePattern() was broken.

Browse Source
This commit is contained in:
Gael Guennebaud 2014-08-26 23:32:32 +02:00
parent be3477e206
commit 25a3e65a68
2 changed files with 16 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
Eigen/src/SparseQR/SparseQR.h
View File

@ -75,7 +75,7 @@ class SparseQR
typedef Matrix<Scalar, Dynamic, 1> ScalarVector; typedef Matrix<Scalar, Dynamic, 1> ScalarVector;
typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
public: public:
SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false) SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
{ } { }
/** Construct a QR factorization of the matrix \a mat. /** Construct a QR factorization of the matrix \a mat.
@ -84,7 +84,7 @@ class SparseQR
* *
* \sa compute() * \sa compute()
*/ */
SparseQR(const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false) SparseQR(const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false)
{ {
compute(mat); compute(mat);
} }
@ -262,6 +262,7 @@ class SparseQR
IndexVector m_etree; // Column elimination tree IndexVector m_etree; // Column elimination tree
IndexVector m_firstRowElt; // First element in each row IndexVector m_firstRowElt; // First element in each row
bool m_isQSorted; // whether Q is sorted or not bool m_isQSorted; // whether Q is sorted or not
bool m_isEtreeOk; // whether the elimination tree match the initial input matrix
template <typename, typename > friend struct SparseQR_QProduct; template <typename, typename > friend struct SparseQR_QProduct;
template <typename > friend struct SparseQRMatrixQReturnType; template <typename > friend struct SparseQRMatrixQReturnType;
@ -297,6 +298,7 @@ void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat)
// Compute the column elimination tree of the permuted matrix // Compute the column elimination tree of the permuted matrix
m_outputPerm_c = m_perm_c.inverse(); m_outputPerm_c = m_perm_c.inverse();
internal::coletree(mat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data()); internal::coletree(mat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
m_isEtreeOk = true;
m_R.resize(m, n); m_R.resize(m, n);
m_Q.resize(m, diagSize); m_Q.resize(m, diagSize);
@ -331,6 +333,15 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
ScalarVector tval(m); // The dense vector used to compute the current column ScalarVector tval(m); // The dense vector used to compute the current column
RealScalar pivotThreshold = m_threshold; RealScalar pivotThreshold = m_threshold;
m_R.setZero();
m_Q.setZero();
if(!m_isEtreeOk)
{
m_outputPerm_c = m_perm_c.inverse();
internal::coletree(mat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
m_isEtreeOk = true;
}
m_pmat = mat; m_pmat = mat;
m_pmat.uncompress(); // To have the innerNonZeroPtr allocated m_pmat.uncompress(); // To have the innerNonZeroPtr allocated
// Apply the fill-in reducing permutation lazily: // Apply the fill-in reducing permutation lazily:
@ -513,6 +524,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
// Recompute the column elimination tree // Recompute the column elimination tree
internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data()); internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data());
m_isEtreeOk = false;
} }
} }

2
test/sparseqr.cpp
View File

@ -54,6 +54,8 @@ template<typename Scalar> void test_sparseqr_scalar()
b = dA * DenseVector::Random(A.cols()); b = dA * DenseVector::Random(A.cols());
solver.compute(A); solver.compute(A);
if(internal::random<float>(0,1)>0.5)
solver.factorize(A); // this checks that calling analyzePattern is not needed if the pattern do not change.
if (solver.info() != Success) if (solver.info() != Success)
{ {
std::cerr << "sparse QR factorization failed\n"; std::cerr << "sparse QR factorization failed\n";