diff options
Diffstat (limited to 'Eigen/src/SparseQR/SparseQR.h')
| -rw-r--r-- | Eigen/src/SparseQR/SparseQR.h | 51 |
1 files changed, 35 insertions, 16 deletions
diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h index 2d4498b03..d1fb96f5c 100644 --- a/Eigen/src/SparseQR/SparseQR.h +++ b/Eigen/src/SparseQR/SparseQR.h @@ -41,18 +41,19 @@ namespace internal { /** * \ingroup SparseQR_Module * \class SparseQR - * \brief Sparse left-looking rank-revealing QR factorization + * \brief Sparse left-looking QR factorization with numerical column pivoting * - * This class implements a left-looking rank-revealing QR decomposition - * of sparse matrices. When a column has a norm less than a given tolerance + * This class implements a left-looking QR decomposition of sparse matrices + * with numerical column pivoting. + * When a column has a norm less than a given tolerance * it is implicitly permuted to the end. The QR factorization thus obtained is * given by A*P = Q*R where R is upper triangular or trapezoidal. * * P is the column permutation which is the product of the fill-reducing and the - * rank-revealing permutations. Use colsPermutation() to get it. + * numerical permutations. Use colsPermutation() to get it. * * Q is the orthogonal matrix represented as products of Householder reflectors. - * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose. + * Use matrixQ() to get an expression and matrixQ().adjoint() to get the adjoint. * You can then apply it to a vector. * * R is the sparse triangular or trapezoidal matrix. The later occurs when A is rank-deficient. @@ -64,7 +65,19 @@ namespace internal { * * \implsparsesolverconcept * + * The numerical pivoting strategy and default threshold are the same as in SuiteSparse QR, and + * detailed in the following paper: + * <i> + * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing + * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011. + * </i> + * Even though it is qualified as "rank-revealing", this strategy might fail for some + * rank deficient problems. When this class is used to solve linear or least-square problems + * it is thus strongly recommended to check the accuracy of the computed solution. If it + * failed, it usually helps to increase the threshold with setPivotThreshold. + * * \warning The input sparse matrix A must be in compressed mode (see SparseMatrix::makeCompressed()). + * \warning For complex matrices matrixQ().transpose() will actually return the adjoint matrix. * */ template<typename _MatrixType, typename _OrderingType> @@ -196,9 +209,9 @@ class SparseQR : public SparseSolverBase<SparseQR<_MatrixType,_OrderingType> > Index rank = this->rank(); - // Compute Q^T * b; + // Compute Q^* * b; typename Dest::PlainObject y, b; - y = this->matrixQ().transpose() * B; + y = this->matrixQ().adjoint() * B; b = y; // Solve with the triangular matrix R @@ -330,7 +343,7 @@ void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat) m_R.resize(m, n); m_Q.resize(m, diagSize); - // Allocate space for nonzero elements : rough estimation + // Allocate space for nonzero elements: rough estimation m_R.reserve(2*mat.nonZeros()); //FIXME Get a more accurate estimation through symbolic factorization with the etree m_Q.reserve(2*mat.nonZeros()); m_hcoeffs.resize(diagSize); @@ -604,7 +617,7 @@ struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived // Get the references SparseQR_QProduct(const SparseQRType& qr, const Derived& other, bool transpose) : m_qr(qr),m_other(other),m_transpose(transpose) {} - inline Index rows() const { return m_transpose ? m_qr.rows() : m_qr.cols(); } + inline Index rows() const { return m_qr.matrixQ().rows(); } inline Index cols() const { return m_other.cols(); } // Assign to a vector @@ -632,16 +645,20 @@ struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived } else { - eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes"); + eigen_assert(m_qr.matrixQ().cols() == m_other.rows() && "Non conforming object sizes"); + + res.conservativeResize(rows(), cols()); + // Compute res = Q * other column by column for(Index j = 0; j < res.cols(); j++) { - for (Index k = diagSize-1; k >=0; k--) + Index start_k = internal::is_identity<Derived>::value ? numext::mini(j,diagSize-1) : diagSize-1; + for (Index k = start_k; k >=0; k--) { Scalar tau = Scalar(0); tau = m_qr.m_Q.col(k).dot(res.col(j)); if(tau==Scalar(0)) continue; - tau = tau * m_qr.m_hcoeffs(k); + tau = tau * numext::conj(m_qr.m_hcoeffs(k)); res.col(j) -= tau * m_qr.m_Q.col(k); } } @@ -650,7 +667,7 @@ struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived const SparseQRType& m_qr; const Derived& m_other; - bool m_transpose; + bool m_transpose; // TODO this actually means adjoint }; template<typename SparseQRType> @@ -668,13 +685,14 @@ struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<Sp { return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),false); } + // To use for operations with the adjoint of Q SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint() const { return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr); } inline Index rows() const { return m_qr.rows(); } - inline Index cols() const { return (std::min)(m_qr.rows(),m_qr.cols()); } - // To use for operations with the transpose of Q + inline Index cols() const { return m_qr.rows(); } + // To use for operations with the transpose of Q FIXME this is the same as adjoint at the moment SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const { return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr); @@ -682,6 +700,7 @@ struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<Sp const SparseQRType& m_qr; }; +// TODO this actually represents the adjoint of Q template<typename SparseQRType> struct SparseQRMatrixQTransposeReturnType { @@ -712,7 +731,7 @@ struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal: typedef typename DstXprType::StorageIndex StorageIndex; static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &/*func*/) { - typename DstXprType::PlainObject idMat(src.m_qr.rows(), src.m_qr.rows()); + typename DstXprType::PlainObject idMat(src.rows(), src.cols()); idMat.setIdentity(); // Sort the sparse householder reflectors if needed const_cast<SparseQRType *>(&src.m_qr)->_sort_matrix_Q(); |