diff options
author | Gael Guennebaud <g.gael@free.fr> | 2014-09-19 09:58:56 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2014-09-19 09:58:56 +0200 |
commit | 755e77266f7c1005688654f7e7b070894c43bd18 (patch) | |
tree | 33d5b045586cd4384bdc8ced00a58a3cd49d9747 /Eigen/src | |
parent | 07c5500d709eb7914bd46a1c4b3629bf42783f1d (diff) |
Fix SparseQR for row-major inputs.
Diffstat (limited to 'Eigen/src')
-rw-r--r-- | Eigen/src/SparseQR/SparseQR.h | 38 |
1 files changed, 27 insertions, 11 deletions
diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h index 1a6c84e00..6d85ea9be 100644 --- a/Eigen/src/SparseQR/SparseQR.h +++ b/Eigen/src/SparseQR/SparseQR.h @@ -284,9 +284,11 @@ template <typename MatrixType, typename OrderingType> 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"); + // 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 OrderingType ord; - ord(mat, m_perm_c); + ord(matCpy, m_perm_c); Index n = mat.cols(); Index m = mat.rows(); 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 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_R.resize(m, n); @@ -337,21 +339,35 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat) m_R.setZero(); m_Q.setZero(); + m_pmat = mat; if(!m_isEtreeOk) { 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_pmat = mat; + m_pmat.uncompress(); // To have the innerNonZeroPtr allocated + // 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; - m_pmat.outerIndexPtr()[p] = mat.outerIndexPtr()[i]; - m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i]; + // If the input is row major, copy the original column indices, + // otherwise directly use the input matrix + // + 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: @@ -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. // 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. - 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; if(itp) curIdx = itp.row(); @@ -544,7 +560,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat) if(nonzeroCol<n) { // 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; // Update the column permutation |