From 4477843bdd7555f634952357f6087021e2d8519b Mon Sep 17 00:00:00 2001 From: Adolfo Rodriguez Tsourouksdissian Date: Sun, 30 Oct 2011 23:55:20 -0400 Subject: [PATCH] bug #206 - part 4: Removes heap allocations from JacobiSVD and its preconditioners --- Eigen/src/QR/FullPivHouseholderQR.h | 2 +- Eigen/src/SVD/JacobiSVD.h | 269 ++++++++++++++++++++++------ 2 files changed, 212 insertions(+), 59 deletions(-) diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index 4992c2c3c..96cd4a1cf 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -562,7 +562,7 @@ public: // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] const Index rows = m_qr.rows(); const Index cols = m_qr.cols(); - const Index size = std::min(rows, cols); + const Index size = (std::min)(rows, cols); workspace.resize(rows); result.setIdentity(rows, rows); for (Index k = size-1; k >= 0; k--) diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index c2326d752..3daa29df5 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -61,9 +61,12 @@ template struct qr_preconditioner_impl {}; template -struct qr_preconditioner_impl +class qr_preconditioner_impl { - static bool run(JacobiSVD&, const MatrixType&) +public: + typedef typename MatrixType::Index Index; + void allocate(const JacobiSVD&) {} + bool run(JacobiSVD&, const MatrixType&) { return false; } @@ -72,134 +75,279 @@ struct qr_preconditioner_impl /*** preconditioner using FullPivHouseholderQR ***/ template -struct qr_preconditioner_impl +class qr_preconditioner_impl { - static bool run(JacobiSVD& svd, const MatrixType& matrix) +public: + typedef typename MatrixType::Index Index; + typedef typename MatrixType::Scalar Scalar; + enum + { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime + }; + typedef Matrix WorkspaceType; + + void allocate(const JacobiSVD& svd) + { + if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) + { + m_qr = FullPivHouseholderQR(svd.rows(), svd.cols()); + } + if (svd.m_computeFullU) m_workspace.resize(svd.rows()); + } + + bool run(JacobiSVD& svd, const MatrixType& matrix) { if(matrix.rows() > matrix.cols()) { - FullPivHouseholderQR qr(matrix); - svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView(); - if(svd.m_computeFullU) svd.m_matrixU = qr.matrixQ(); - if(svd.computeV()) svd.m_matrixV = qr.colsPermutation(); + m_qr.compute(matrix); + svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView(); + if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace); + if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation(); return true; } return false; } +private: + FullPivHouseholderQR m_qr; + WorkspaceType m_workspace; }; template -struct qr_preconditioner_impl +class qr_preconditioner_impl { - static bool run(JacobiSVD& svd, const MatrixType& matrix) +public: + typedef typename MatrixType::Index Index; + typedef typename MatrixType::Scalar Scalar; + enum + { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, + Options = MatrixType::Options + }; + typedef Matrix + TransposeTypeWithSameStorageOrder; + + void allocate(const JacobiSVD& svd) + { + if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) + { + m_qr = FullPivHouseholderQR(svd.cols(), svd.rows()); + } + m_adjoint.resize(svd.cols(), svd.rows()); + if (svd.m_computeFullV) m_workspace.resize(svd.cols()); + } + + bool run(JacobiSVD& svd, const MatrixType& matrix) { if(matrix.cols() > matrix.rows()) { - typedef Matrix - TransposeTypeWithSameStorageOrder; - FullPivHouseholderQR qr(matrix.adjoint()); - svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView().adjoint(); - if(svd.m_computeFullV) svd.m_matrixV = qr.matrixQ(); - if(svd.computeU()) svd.m_matrixU = qr.colsPermutation(); + m_adjoint = matrix.adjoint(); + m_qr.compute(m_adjoint); + svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView().adjoint(); + if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace); + if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation(); return true; } else return false; } +private: + FullPivHouseholderQR m_qr; + TransposeTypeWithSameStorageOrder m_adjoint; + typename internal::plain_row_type::type m_workspace; }; /*** preconditioner using ColPivHouseholderQR ***/ template -struct qr_preconditioner_impl +class qr_preconditioner_impl { - static bool run(JacobiSVD& svd, const MatrixType& matrix) +public: + typedef typename MatrixType::Index Index; + + void allocate(const JacobiSVD& svd) + { + if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) + { + m_qr = ColPivHouseholderQR(svd.rows(), svd.cols()); + } + if (svd.m_computeFullU) m_workspace.resize(svd.rows()); + else if (svd.m_computeThinU) m_workspace.resize(svd.cols()); + } + + bool run(JacobiSVD& svd, const MatrixType& matrix) { if(matrix.rows() > matrix.cols()) { - ColPivHouseholderQR qr(matrix); - svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView(); - if(svd.m_computeFullU) svd.m_matrixU = qr.householderQ(); - else if(svd.m_computeThinU) { + m_qr.compute(matrix); + svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView(); + if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace); + else if(svd.m_computeThinU) + { svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); - qr.householderQ().applyThisOnTheLeft(svd.m_matrixU); + m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace); } - if(svd.computeV()) svd.m_matrixV = qr.colsPermutation(); + if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation(); return true; } return false; } + +private: + ColPivHouseholderQR m_qr; + typename internal::plain_col_type::type m_workspace; }; template -struct qr_preconditioner_impl +class qr_preconditioner_impl { - static bool run(JacobiSVD& svd, const MatrixType& matrix) +public: + typedef typename MatrixType::Index Index; + typedef typename MatrixType::Scalar Scalar; + enum + { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, + Options = MatrixType::Options + }; + + typedef Matrix + TransposeTypeWithSameStorageOrder; + + void allocate(const JacobiSVD& svd) + { + if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) + { + m_qr = ColPivHouseholderQR(svd.cols(), svd.rows()); + } + if (svd.m_computeFullV) m_workspace.resize(svd.cols()); + else if (svd.m_computeThinV) m_workspace.resize(svd.rows()); + m_adjoint.resize(svd.cols(), svd.rows()); + } + + bool run(JacobiSVD& svd, const MatrixType& matrix) { if(matrix.cols() > matrix.rows()) { - typedef Matrix - TransposeTypeWithSameStorageOrder; - ColPivHouseholderQR qr(matrix.adjoint()); - svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView().adjoint(); - if(svd.m_computeFullV) svd.m_matrixV = qr.householderQ(); - else if(svd.m_computeThinV) { + m_adjoint = matrix.adjoint(); + m_qr.compute(m_adjoint); + + svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView().adjoint(); + if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace); + else if(svd.m_computeThinV) + { svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); - qr.householderQ().applyThisOnTheLeft(svd.m_matrixV); + m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace); } - if(svd.computeU()) svd.m_matrixU = qr.colsPermutation(); + if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation(); return true; } else return false; } + +private: + ColPivHouseholderQR m_qr; + TransposeTypeWithSameStorageOrder m_adjoint; + typename internal::plain_row_type::type m_workspace; }; /*** preconditioner using HouseholderQR ***/ template -struct qr_preconditioner_impl +class qr_preconditioner_impl { - static bool run(JacobiSVD& svd, const MatrixType& matrix) +public: + typedef typename MatrixType::Index Index; + + void allocate(const JacobiSVD& svd) + { + if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) + { + m_qr = HouseholderQR(svd.rows(), svd.cols()); + } + if (svd.m_computeFullU) m_workspace.resize(svd.rows()); + else if (svd.m_computeThinU) m_workspace.resize(svd.cols()); + } + + bool run(JacobiSVD& svd, const MatrixType& matrix) { if(matrix.rows() > matrix.cols()) { - HouseholderQR qr(matrix); - svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView(); - if(svd.m_computeFullU) svd.m_matrixU = qr.householderQ(); - else if(svd.m_computeThinU) { + m_qr.compute(matrix); + svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView(); + if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace); + else if(svd.m_computeThinU) + { svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); - qr.householderQ().applyThisOnTheLeft(svd.m_matrixU); + m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace); } if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols()); return true; } return false; } +private: + HouseholderQR m_qr; + typename internal::plain_col_type::type m_workspace; }; template -struct qr_preconditioner_impl +class qr_preconditioner_impl { - static bool run(JacobiSVD& svd, const MatrixType& matrix) +public: + typedef typename MatrixType::Index Index; + typedef typename MatrixType::Scalar Scalar; + enum + { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, + Options = MatrixType::Options + }; + + typedef Matrix + TransposeTypeWithSameStorageOrder; + + void allocate(const JacobiSVD& svd) + { + if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) + { + m_qr = HouseholderQR(svd.cols(), svd.rows()); + } + if (svd.m_computeFullV) m_workspace.resize(svd.cols()); + else if (svd.m_computeThinV) m_workspace.resize(svd.rows()); + m_adjoint.resize(svd.cols(), svd.rows()); + } + + bool run(JacobiSVD& svd, const MatrixType& matrix) { if(matrix.cols() > matrix.rows()) { - typedef Matrix - TransposeTypeWithSameStorageOrder; - HouseholderQR qr(matrix.adjoint()); - svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView().adjoint(); - if(svd.m_computeFullV) svd.m_matrixV = qr.householderQ(); - else if(svd.m_computeThinV) { + m_adjoint = matrix.adjoint(); + m_qr.compute(m_adjoint); + + svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView().adjoint(); + if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace); + else if(svd.m_computeThinV) + { svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); - qr.householderQ().applyThisOnTheLeft(svd.m_matrixV); + m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace); } if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows()); return true; } else return false; } + +private: + HouseholderQR m_qr; + TransposeTypeWithSameStorageOrder m_adjoint; + typename internal::plain_row_type::type m_workspace; }; /*** 2x2 SVD implementation @@ -316,7 +464,7 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, * Here's an example demonstrating basic usage: * \include JacobiSVD_basic.cpp * Output: \verbinclude JacobiSVD_basic.out - * + * * This JacobiSVD class is a two-sided Jacobi R-SVD decomposition, ensuring optimal reliability and accuracy. The downside is that it's slower than * bidiagonalizing SVD algorithms for large square matrices; however its complexity is still \f$ O(n^2p) \f$ where \a n is the smaller dimension and * \a p is the greater dimension, meaning that it is still of the same order of complexity as the faster bidiagonalizing R-SVD algorithms. @@ -324,7 +472,7 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, * * If the input matrix has inf or nan coefficients, the result of the computation is undefined, but the computation is guaranteed to * terminate in finite (and reasonable) time. - * + * * The possible values for QRPreconditioner are: * \li ColPivHouseholderQRPreconditioner is the default. In practice it's very safe. It uses column-pivoting QR. * \li FullPivHouseholderQRPreconditioner, is the safest and slowest. It uses full-pivoting QR. @@ -494,7 +642,7 @@ template class JacobiSVD * \param b the right-hand-side of the equation to solve. * * \note Solving requires both U and V to be computed. Thin U and V are enough, there is no need for full U or V. - * + * * \note SVD solving is implicitly least-squares. Thus, this method serves both purposes of exact solving and least-squares solving. * In other words, the returned solution is guaranteed to minimize the Euclidean norm \f$ \Vert A x - b \Vert \f$. */ @@ -535,6 +683,9 @@ template class JacobiSVD friend struct internal::svd_precondition_2x2_block_to_be_real; template friend struct internal::qr_preconditioner_impl; + + internal::qr_preconditioner_impl m_qr_precond_morecols; + internal::qr_preconditioner_impl m_qr_precond_morerows; }; template @@ -578,6 +729,9 @@ void JacobiSVD::allocate(Index rows, Index cols, u : m_computeThinV ? m_diagSize : 0); m_workMatrix.resize(m_diagSize, m_diagSize); + + m_qr_precond_morecols.allocate(*this); + m_qr_precond_morerows.allocate(*this); } template @@ -595,8 +749,7 @@ JacobiSVD::compute(const MatrixType& matrix, unsig /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ - if(!internal::qr_preconditioner_impl::run(*this, matrix) - && !internal::qr_preconditioner_impl::run(*this, matrix)) + if(!m_qr_precond_morecols.run(*this, matrix) && !m_qr_precond_morerows.run(*this, matrix)) { m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize); if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);