GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-05-22 12:37:35 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix SparseQR for row-major inputs.

Browse Source
This commit is contained in:
Gael Guennebaud 2014-09-19 09:58:56 +02:00
parent 07c5500d70
commit 755e77266f
1 changed files with 27 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -284,9 +284,11 @@ template <typename MatrixType, typename OrderingType>
void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat) void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat)
{ {
eigen_assert(mat.isCompressed() && "SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR"); eigen_assert(mat.isCompressed() && "SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR");
// Copy to a column major matrix if the input is rowmajor
typename internal::conditional<MatrixType::IsRowMajor,QRMatrixType,const MatrixType&>::type matCpy(mat);
// Compute the column fill reducing ordering // Compute the column fill reducing ordering
OrderingType ord; OrderingType ord;
ord(mat, m_perm_c); ord(matCpy, m_perm_c);
Index n = mat.cols(); Index n = mat.cols();
Index m = mat.rows(); Index m = mat.rows();
Index diagSize = (std::min)(m,n); Index diagSize = (std::min)(m,n);
@ -299,7 +301,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(matCpy, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
m_isEtreeOk = true; m_isEtreeOk = true;
m_R.resize(m, n); m_R.resize(m, n);
@ -337,21 +339,35 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
m_R.setZero(); m_R.setZero();
m_Q.setZero(); m_Q.setZero();
m_pmat = mat;
if(!m_isEtreeOk) if(!m_isEtreeOk)
{ {
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(m_pmat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
m_isEtreeOk = true; m_isEtreeOk = true;
} }
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:
for (int i = 0; i < n; i++)
{ {
Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i; // If the input is row major, copy the original column indices,
m_pmat.outerIndexPtr()[p] = mat.outerIndexPtr()[i]; // otherwise directly use the input matrix
m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i]; //
IndexVector originalOuterIndicesCpy;
const Index *originalOuterIndices = mat.outerIndexPtr();
if(MatrixType::IsRowMajor)
{
originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1);
originalOuterIndices = originalOuterIndicesCpy.data();
}
for (int i = 0; i < n; i++)
{
Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
m_pmat.outerIndexPtr()[p] = originalOuterIndices[i];
m_pmat.innerNonZeroPtr()[p] = originalOuterIndices[i+1] - originalOuterIndices[i];
}
} }
/* Compute the default threshold as in MatLab, see: /* Compute the default threshold as in MatLab, see:
@ -386,7 +402,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
// all the nodes (with indexes lower than rank) reachable through the column elimination tree (etree) rooted at node k. // all the nodes (with indexes lower than rank) reachable through the column elimination tree (etree) rooted at node k.
// Note: if the diagonal entry does not exist, then its contribution must be explicitly added, // Note: if the diagonal entry does not exist, then its contribution must be explicitly added,
// thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found. // thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found.
for (typename MatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp) for (typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
{ {
Index curIdx = nonzeroCol; Index curIdx = nonzeroCol;
if(itp) curIdx = itp.row(); if(itp) curIdx = itp.row();
@ -544,7 +560,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
if(nonzeroCol<n) if(nonzeroCol<n)
{ {
// Permute the triangular factor to put the 'dead' columns to the end // Permute the triangular factor to put the 'dead' columns to the end
MatrixType tempR(m_R); QRMatrixType tempR(m_R);
m_R = tempR * m_pivotperm; m_R = tempR * m_pivotperm;
// Update the column permutation // Update the column permutation