SPQR: fix fixed-ordering permutation and complete matrixQ()/matrixR() handling libeigen/eigen!2468 Closes #2376, #1377, #1121, #2049, and #2790
diff --git a/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h b/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h index 3132794..99598f5 100644 --- a/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h +++ b/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h
@@ -28,6 +28,11 @@ template <typename SPQRType> struct traits<SPQRMatrixQReturnType<SPQRType> > { typedef typename SPQRType::MatrixType ReturnType; + typedef typename ReturnType::Scalar Scalar; + typedef typename ReturnType::StorageIndex StorageIndex; + typedef typename ReturnType::StorageKind StorageKind; + static constexpr int RowsAtCompileTime = Dynamic; + static constexpr int ColsAtCompileTime = Dynamic; }; template <typename SPQRType> struct traits<SPQRMatrixQTransposeReturnType<SPQRType> > { @@ -120,7 +125,9 @@ cholmod_l_free_sparse(&m_cR, &m_cc); cholmod_l_free_dense(&m_HTau, &m_cc); std::free(m_E); + m_E = nullptr; std::free(m_HPinv); + m_HPinv = nullptr; } void compute(const MatrixType_& matrix) { @@ -150,6 +157,12 @@ m_isInitialized = false; return; } + if (!m_E && !initIdentityPermutation(m_cR->ncol)) { + SPQR_free(); + m_info = NumericalIssue; + m_isInitialized = false; + return; + } m_info = Success; m_isInitialized = true; m_isRUpToDate = false; @@ -193,7 +206,7 @@ /** \returns the sparse triangular factor R. It is a sparse matrix */ - const MatrixType matrixR() const { + const MatrixType& matrixR() const { eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); if (!m_isRUpToDate) { m_R = viewAsEigen<Scalar, StorageIndex>(*m_cR); @@ -255,6 +268,16 @@ mutable cholmod_common m_cc; // Workspace and parameters bool m_useDefaultThreshold; // Use default threshold Index m_rows; + + bool initIdentityPermutation(StorageIndex size) { + if (m_E || size == 0) return true; + // SuiteSparse can omit the permutation array when no column reordering is applied. + m_E = static_cast<StorageIndex*>(std::malloc(sizeof(StorageIndex) * size)); + if (!m_E) return false; + for (StorageIndex i = 0; i < size; ++i) m_E[i] = i; + return true; + } + template <typename, typename> friend struct SPQR_QProduct; }; @@ -267,33 +290,102 @@ SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr), m_other(other), m_transpose(transpose) {} - inline Index rows() const { return m_transpose ? m_spqr.rows() : m_spqr.cols(); } + inline Index rows() const { return m_spqr.rows(); } inline Index cols() const { return m_other.cols(); } // Assign to a vector template <typename ResType> void evalTo(ResType& res) const { + evalToImpl(res, m_other); + } + + private: + template <typename ResType, typename OtherDerived, + typename std::enable_if<(int(OtherDerived::Flags) & DirectAccessBit) == DirectAccessBit, int>::type = 0> + void evalToImpl(ResType& res, const MatrixBase<OtherDerived>& otherExpr) const { cholmod_dense y_cd; cholmod_dense* x_cd; int method = m_transpose ? SPQR_QTX : SPQR_QX; cholmod_common* cc = m_spqr.cholmodCommon(); - y_cd = viewAsCholmod(m_other.const_cast_derived()); + y_cd = viewAsCholmod(otherExpr.const_cast_derived()); x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc); res = Matrix<Scalar, ResType::RowsAtCompileTime, ResType::ColsAtCompileTime>::Map( reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol); cholmod_l_free_dense(&x_cd, cc); } + + template <typename ResType, typename OtherDerived, + typename std::enable_if<(int(OtherDerived::Flags) & DirectAccessBit) == 0, int>::type = 0> + void evalToImpl(ResType& res, const MatrixBase<OtherDerived>& otherExpr) const { + cholmod_dense y_cd; + cholmod_dense* x_cd; + int method = m_transpose ? SPQR_QTX : SPQR_QX; + cholmod_common* cc = m_spqr.cholmodCommon(); + typename OtherDerived::PlainObject other = otherExpr; + y_cd = viewAsCholmod(other); + x_cd = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cd, cc); + res = Matrix<Scalar, ResType::RowsAtCompileTime, ResType::ColsAtCompileTime>::Map( + reinterpret_cast<Scalar*>(x_cd->x), x_cd->nrow, x_cd->ncol); + cholmod_l_free_dense(&x_cd, cc); + } + + template <typename ResType, typename OtherDerived> + void evalToImpl(ResType& res, const SparseMatrixBase<OtherDerived>& otherExpr) const { + cholmod_sparse y_cs; + cholmod_sparse* x_cs; + int method = m_transpose ? SPQR_QTX : SPQR_QX; + cholmod_common* cc = m_spqr.cholmodCommon(); + typename OtherDerived::PlainObject other = otherExpr; + other.makeCompressed(); + y_cs = viewAsCholmod(other); + x_cs = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cs, cc); + res = viewAsEigen<Scalar, StorageIndex>(*x_cs); + cholmod_l_free_sparse(&x_cs, cc); + } + + template <typename ResType, typename OtherScalar, int OtherOptions, typename OtherStorageIndex, + typename std::enable_if<internal::is_same<OtherStorageIndex, StorageIndex>::value, int>::type = 0> + void evalToImpl(ResType& res, const SparseMatrix<OtherScalar, OtherOptions, OtherStorageIndex>& otherExpr) const { + cholmod_sparse y_cs; + cholmod_sparse* x_cs; + int method = m_transpose ? SPQR_QTX : SPQR_QX; + cholmod_common* cc = m_spqr.cholmodCommon(); + const SparseMatrix<OtherScalar, OtherOptions, OtherStorageIndex>* otherPtr = &otherExpr; + SparseMatrix<OtherScalar, OtherOptions, OtherStorageIndex> other; + + if (!otherExpr.isCompressed()) { + other = otherExpr; + other.makeCompressed(); + otherPtr = &other; + } + + y_cs = viewAsCholmod(*otherPtr); + x_cs = SuiteSparseQR_qmult<Scalar>(method, m_spqr.m_H, m_spqr.m_HTau, m_spqr.m_HPinv, &y_cs, cc); + res = viewAsEigen<Scalar, StorageIndex>(*x_cs); + cholmod_l_free_sparse(&x_cs, cc); + } + + public: const SPQRType& m_spqr; const Derived& m_other; bool m_transpose; }; template <typename SPQRType> -struct SPQRMatrixQReturnType { +struct SPQRMatrixQReturnType : public EigenBase<SPQRMatrixQReturnType<SPQRType> > { + typedef typename SPQRType::Scalar Scalar; + static constexpr int RowsAtCompileTime = Dynamic; + static constexpr int ColsAtCompileTime = Dynamic; SPQRMatrixQReturnType(const SPQRType& spqr) : m_spqr(spqr) {} template <typename Derived> SPQR_QProduct<SPQRType, Derived> operator*(const MatrixBase<Derived>& other) { return SPQR_QProduct<SPQRType, Derived>(m_spqr, other.derived(), false); } + template <typename Derived> + SPQR_QProduct<SPQRType, Derived> operator*(const SparseMatrixBase<Derived>& other) { + return SPQR_QProduct<SPQRType, Derived>(m_spqr, other.derived(), false); + } SPQRMatrixQTransposeReturnType<SPQRType> adjoint() const { return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr); } + inline Index rows() const { return m_spqr.rows(); } + inline Index cols() const { return m_spqr.rows(); } // To use for operations with the transpose of Q SPQRMatrixQTransposeReturnType<SPQRType> transpose() const { return SPQRMatrixQTransposeReturnType<SPQRType>(m_spqr); @@ -308,8 +400,47 @@ SPQR_QProduct<SPQRType, Derived> operator*(const MatrixBase<Derived>& other) { return SPQR_QProduct<SPQRType, Derived>(m_spqr, other.derived(), true); } + template <typename Derived> + SPQR_QProduct<SPQRType, Derived> operator*(const SparseMatrixBase<Derived>& other) { + return SPQR_QProduct<SPQRType, Derived>(m_spqr, other.derived(), true); + } const SPQRType& m_spqr; }; +namespace internal { + +template <typename SPQRType> +struct evaluator_traits<SPQRMatrixQReturnType<SPQRType> > { + typedef typename SPQRType::MatrixType MatrixType; + typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind; + typedef SparseShape Shape; +}; + +template <typename DstXprType, typename SPQRType> +struct Assignment<DstXprType, SPQRMatrixQReturnType<SPQRType>, + internal::assign_op<typename DstXprType::Scalar, typename DstXprType::Scalar>, Sparse2Sparse> { + typedef SPQRMatrixQReturnType<SPQRType> SrcXprType; + typedef typename DstXprType::Scalar Scalar; + + static void run(DstXprType& dst, const SrcXprType& src, const internal::assign_op<Scalar, Scalar>& /*func*/) { + typename DstXprType::PlainObject idMat(src.rows(), src.cols()); + idMat.setIdentity(); + dst = src.m_spqr.matrixQ() * idMat; + } +}; + +template <typename DstXprType, typename SPQRType> +struct Assignment<DstXprType, SPQRMatrixQReturnType<SPQRType>, + internal::assign_op<typename DstXprType::Scalar, typename DstXprType::Scalar>, Sparse2Dense> { + typedef SPQRMatrixQReturnType<SPQRType> SrcXprType; + typedef typename DstXprType::Scalar Scalar; + + static void run(DstXprType& dst, const SrcXprType& src, const internal::assign_op<Scalar, Scalar>& /*func*/) { + dst = src.m_spqr.matrixQ() * DstXprType::Identity(src.rows(), src.cols()); + } +}; + +} // namespace internal + } // End namespace Eigen #endif
diff --git a/test/spqr_support.cpp b/test/spqr_support.cpp index cd75825..8a2c047 100644 --- a/test/spqr_support.cpp +++ b/test/spqr_support.cpp
@@ -53,7 +53,104 @@ refX = dA.colPivHouseholderQr().solve(b); VERIFY(x.isApprox(refX, test_precision<Scalar>())); } + +void test_spqr_fixed_ordering_uses_identity_permutation() { + typedef SparseMatrix<double, ColMajor> MatrixType; + typedef Matrix<double, Dynamic, Dynamic> DenseMatrix; + typedef Matrix<double, Dynamic, 1> DenseVector; + + DenseMatrix dA(6, 4); + dA << 4.0, 1.0, 0.0, 0.0, // + 1.0, 0.0, 2.0, 0.0, // + -2.0, 3.0, 0.0, 6.0, // + 0.0, 5.0, -1.0, 0.0, // + 0.0, 0.0, 7.0, 2.0, // + 0.0, 0.0, 0.0, 3.0; + + MatrixType A = dA.sparseView(); + A.makeCompressed(); + + DenseVector b(6); + b << 1.0, -2.0, 0.5, 4.0, -1.0, 3.0; + + SPQR<MatrixType> solver; + solver.setSPQROrdering(SPQR_ORDERING_FIXED); + solver.setPivotThreshold(SPQR_NO_TOL); + solver.compute(A); + + VERIFY_IS_EQUAL(solver.info(), Success); + VERIFY_IS_EQUAL(solver.rank(), A.cols()); + + const auto permutation = solver.colsPermutation(); + VERIFY_IS_EQUAL(permutation.size(), A.cols()); + for (Index i = 0; i < permutation.size(); ++i) { + VERIFY_IS_EQUAL(permutation.indices()(i), i); + } + + const DenseVector refX = dA.colPivHouseholderQr().solve(b); + DenseVector x = solver.solve(b); + VERIFY_IS_EQUAL(solver.info(), Success); + VERIFY_IS_APPROX(x, refX); +} + +void test_spqr_matrix_q_times_identity_expression() { + typedef SparseMatrix<double, ColMajor> MatrixType; + typedef Matrix<double, Dynamic, Dynamic> DenseMatrix; + typedef SPQR<MatrixType> SolverType; + typedef typename SolverType::MatrixType SolverSparseMatrix; + typedef Matrix<double, Dynamic, 1> DenseVector; + + DenseMatrix dA(6, 4); + dA << 4.0, 1.0, 0.0, 0.0, // + 1.0, 0.0, 2.0, 0.0, // + -2.0, 3.0, 0.0, 6.0, // + 0.0, 5.0, -1.0, 0.0, // + 0.0, 0.0, 7.0, 2.0, // + 0.0, 0.0, 0.0, 3.0; + + MatrixType A = dA.sparseView(); + A.makeCompressed(); + + SolverType solver; + solver.compute(A); + + VERIFY_IS_EQUAL(solver.info(), Success); + + const DenseMatrix Q = solver.matrixQ() * DenseMatrix::Identity(A.rows(), A.rows()); + const DenseMatrix denseIdentity = DenseMatrix::Identity(A.rows(), A.rows()); + const auto qProduct = solver.matrixQ() * denseIdentity; + VERIFY_IS_EQUAL(qProduct.rows(), A.rows()); + VERIFY_IS_EQUAL(qProduct.cols(), A.rows()); + const DenseMatrix qFromDenseIdentity = qProduct; + DenseMatrix denseAssignedQ; + denseAssignedQ = solver.matrixQ(); + VERIFY_IS_EQUAL(Q.rows(), A.rows()); + VERIFY_IS_EQUAL(Q.cols(), A.rows()); + VERIFY_IS_APPROX(Q.transpose() * Q, DenseMatrix::Identity(A.rows(), A.rows())); + VERIFY_IS_APPROX(qFromDenseIdentity, Q); + VERIFY_IS_APPROX(denseAssignedQ, Q); + + const DenseMatrix R = DenseMatrix(solver.matrixR().template triangularView<Upper>()); + const auto sparseR = solver.matrixR().template triangularView<Upper>(); + SolverSparseMatrix sparseIdentity(A.rows(), A.rows()); + sparseIdentity.setIdentity(); + SolverSparseMatrix sparseAssignedQ(A.rows(), A.rows()); + sparseAssignedQ = solver.matrixQ(); + SolverSparseMatrix sparseProductQ(A.rows(), A.rows()); + sparseProductQ = solver.matrixQ() * sparseIdentity; + const DenseVector rhs = DenseVector::LinSpaced(A.cols(), 1.0, double(A.cols())); + const DenseVector x = sparseR.solve(rhs); + const DenseVector expected = R.template triangularView<Upper>().solve(rhs); + const DenseMatrix recoveredA = Q.leftCols(A.cols()) * R * solver.colsPermutation().transpose(); + VERIFY_IS_APPROX(x, expected); + VERIFY_IS_APPROX(DenseMatrix(sparseAssignedQ), denseAssignedQ); + VERIFY_IS_APPROX(DenseMatrix(sparseProductQ), denseAssignedQ); + VERIFY_IS_APPROX(recoveredA, dA); +} + EIGEN_DECLARE_TEST(spqr_support) { CALL_SUBTEST_1(test_spqr_scalar<double>()); CALL_SUBTEST_2(test_spqr_scalar<std::complex<double> >()); + CALL_SUBTEST_3(test_spqr_fixed_ordering_uses_identity_permutation()); + CALL_SUBTEST_3(test_spqr_matrix_q_times_identity_expression()); }