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

JacobiSVD: get rid of m_scaledMatrix, m_adjoint, hopefully fix some compiler warnings

Browse Source
This commit is contained in:
Charles Schlosser 2024-02-17 03:41:55 +00:00
parent 18a161bf17
commit 960892ca13
2 changed files with 26 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

51
Eigen/src/SVD/JacobiSVD.h
View File

@ -52,7 +52,10 @@ template <typename MatrixType, int Options, int QRPreconditioner, int Case>
class qr_preconditioner_impl<MatrixType, Options, QRPreconditioner, Case, false> { class qr_preconditioner_impl<MatrixType, Options, QRPreconditioner, Case, false> {
public: public:
void allocate(const JacobiSVD<MatrixType, Options>&) {} void allocate(const JacobiSVD<MatrixType, Options>&) {}
bool run(JacobiSVD<MatrixType, Options>&, const MatrixType&) { return false; } template <typename Xpr>
bool run(JacobiSVD<MatrixType, Options>&, const Xpr&) {
return false;
}
}; };
/*** preconditioner using FullPivHouseholderQR ***/ /*** preconditioner using FullPivHouseholderQR ***/
@ -75,8 +78,8 @@ class qr_preconditioner_impl<MatrixType, Options, FullPivHouseholderQRPreconditi
} }
if (svd.m_computeFullU) m_workspace.resize(svd.rows()); if (svd.m_computeFullU) m_workspace.resize(svd.rows());
} }
template <typename Xpr>
bool run(SVDType& svd, const MatrixType& matrix) { bool run(SVDType& svd, const Xpr& matrix) {
if (matrix.rows() > matrix.cols()) { if (matrix.rows() > matrix.cols()) {
m_qr.compute(matrix); m_qr.compute(matrix);
svd.m_workMatrix = m_qr.matrixQR().block(0, 0, matrix.cols(), matrix.cols()).template triangularView<Upper>(); svd.m_workMatrix = m_qr.matrixQR().block(0, 0, matrix.cols(), matrix.cols()).template triangularView<Upper>();
@ -117,14 +120,12 @@ class qr_preconditioner_impl<MatrixType, Options, FullPivHouseholderQRPreconditi
internal::destroy_at(&m_qr); internal::destroy_at(&m_qr);
internal::construct_at(&m_qr, svd.cols(), svd.rows()); internal::construct_at(&m_qr, svd.cols(), svd.rows());
} }
m_adjoint.resize(svd.cols(), svd.rows());
if (svd.m_computeFullV) m_workspace.resize(svd.cols()); if (svd.m_computeFullV) m_workspace.resize(svd.cols());
} }
template <typename Xpr>
bool run(SVDType& svd, const MatrixType& matrix) { bool run(SVDType& svd, const Xpr& matrix) {
if (matrix.cols() > matrix.rows()) { if (matrix.cols() > matrix.rows()) {
m_adjoint = matrix.adjoint(); m_qr.compute(matrix.adjoint());
m_qr.compute(m_adjoint);
svd.m_workMatrix = svd.m_workMatrix =
m_qr.matrixQR().block(0, 0, matrix.rows(), matrix.rows()).template triangularView<Upper>().adjoint(); m_qr.matrixQR().block(0, 0, matrix.rows(), matrix.rows()).template triangularView<Upper>().adjoint();
if (svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace); if (svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace);
@ -137,7 +138,6 @@ class qr_preconditioner_impl<MatrixType, Options, FullPivHouseholderQRPreconditi
private: private:
typedef FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType; typedef FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
QRType m_qr; QRType m_qr;
TransposeTypeWithSameStorageOrder m_adjoint;
typename plain_row_type<MatrixType>::type m_workspace; typename plain_row_type<MatrixType>::type m_workspace;
}; };
@ -167,8 +167,8 @@ class qr_preconditioner_impl<MatrixType, Options, ColPivHouseholderQRPreconditio
else if (svd.m_computeThinU) else if (svd.m_computeThinU)
m_workspace.resize(svd.cols()); m_workspace.resize(svd.cols());
} }
template <typename Xpr>
bool run(SVDType& svd, const MatrixType& matrix) { bool run(SVDType& svd, const Xpr& matrix) {
if (matrix.rows() > matrix.cols()) { if (matrix.rows() > matrix.cols()) {
m_qr.compute(matrix); m_qr.compute(matrix);
svd.m_workMatrix = m_qr.matrixQR().block(0, 0, matrix.cols(), matrix.cols()).template triangularView<Upper>(); svd.m_workMatrix = m_qr.matrixQR().block(0, 0, matrix.cols(), matrix.cols()).template triangularView<Upper>();
@ -222,13 +222,11 @@ class qr_preconditioner_impl<MatrixType, Options, ColPivHouseholderQRPreconditio
m_workspace.resize(svd.cols()); m_workspace.resize(svd.cols());
else if (svd.m_computeThinV) else if (svd.m_computeThinV)
m_workspace.resize(svd.rows()); m_workspace.resize(svd.rows());
m_adjoint.resize(svd.cols(), svd.rows());
} }
template <typename Xpr>
bool run(SVDType& svd, const MatrixType& matrix) { bool run(SVDType& svd, const Xpr& matrix) {
if (matrix.cols() > matrix.rows()) { if (matrix.cols() > matrix.rows()) {
m_adjoint = matrix.adjoint(); m_qr.compute(matrix.adjoint());
m_qr.compute(m_adjoint);
svd.m_workMatrix = svd.m_workMatrix =
m_qr.matrixQR().block(0, 0, matrix.rows(), matrix.rows()).template triangularView<Upper>().adjoint(); m_qr.matrixQR().block(0, 0, matrix.rows(), matrix.rows()).template triangularView<Upper>().adjoint();
@ -247,7 +245,6 @@ class qr_preconditioner_impl<MatrixType, Options, ColPivHouseholderQRPreconditio
private: private:
typedef ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType; typedef ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
QRType m_qr; QRType m_qr;
TransposeTypeWithSameStorageOrder m_adjoint;
WorkspaceType m_workspace; WorkspaceType m_workspace;
}; };
@ -276,8 +273,8 @@ class qr_preconditioner_impl<MatrixType, Options, HouseholderQRPreconditioner, P
else if (svd.m_computeThinU) else if (svd.m_computeThinU)
m_workspace.resize(svd.cols()); m_workspace.resize(svd.cols());
} }
template <typename Xpr>
bool run(SVDType& svd, const MatrixType& matrix) { bool run(SVDType& svd, const Xpr& matrix) {
if (matrix.rows() > matrix.cols()) { if (matrix.rows() > matrix.cols()) {
m_qr.compute(matrix); m_qr.compute(matrix);
svd.m_workMatrix = m_qr.matrixQR().block(0, 0, matrix.cols(), matrix.cols()).template triangularView<Upper>(); svd.m_workMatrix = m_qr.matrixQR().block(0, 0, matrix.cols(), matrix.cols()).template triangularView<Upper>();
@ -330,13 +327,12 @@ class qr_preconditioner_impl<MatrixType, Options, HouseholderQRPreconditioner, P
m_workspace.resize(svd.cols()); m_workspace.resize(svd.cols());
else if (svd.m_computeThinV) else if (svd.m_computeThinV)
m_workspace.resize(svd.rows()); m_workspace.resize(svd.rows());
m_adjoint.resize(svd.cols(), svd.rows());
} }
bool run(SVDType& svd, const MatrixType& matrix) { template <typename Xpr>
bool run(SVDType& svd, const Xpr& matrix) {
if (matrix.cols() > matrix.rows()) { if (matrix.cols() > matrix.rows()) {
m_adjoint = matrix.adjoint(); m_qr.compute(matrix.adjoint());
m_qr.compute(m_adjoint);
svd.m_workMatrix = svd.m_workMatrix =
m_qr.matrixQR().block(0, 0, matrix.rows(), matrix.rows()).template triangularView<Upper>().adjoint(); m_qr.matrixQR().block(0, 0, matrix.rows(), matrix.rows()).template triangularView<Upper>().adjoint();
@ -355,7 +351,6 @@ class qr_preconditioner_impl<MatrixType, Options, HouseholderQRPreconditioner, P
private: private:
typedef HouseholderQR<TransposeTypeWithSameStorageOrder> QRType; typedef HouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
QRType m_qr; QRType m_qr;
TransposeTypeWithSameStorageOrder m_adjoint;
WorkspaceType m_workspace; WorkspaceType m_workspace;
}; };
@ -509,7 +504,6 @@ class JacobiSVD : public SVDBase<JacobiSVD<MatrixType_, Options_> > {
typedef MatrixType_ MatrixType; typedef MatrixType_ MatrixType;
typedef typename Base::Scalar Scalar; typedef typename Base::Scalar Scalar;
typedef typename Base::RealScalar RealScalar; typedef typename Base::RealScalar RealScalar;
typedef typename Base::Index Index;
enum : int { enum : int {
Options = Options_, Options = Options_,
QRPreconditioner = internal::get_qr_preconditioner(Options), QRPreconditioner = internal::get_qr_preconditioner(Options),
@ -654,7 +648,6 @@ class JacobiSVD : public SVDBase<JacobiSVD<MatrixType_, Options_> > {
internal::qr_preconditioner_impl<MatrixType, Options, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> internal::qr_preconditioner_impl<MatrixType, Options, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols>
m_qr_precond_morerows; m_qr_precond_morerows;
WorkMatrixType m_workMatrix; WorkMatrixType m_workMatrix;
MatrixType m_scaledMatrix;
}; };
template <typename MatrixType, int Options> template <typename MatrixType, int Options>
@ -669,7 +662,6 @@ void JacobiSVD<MatrixType, Options>::allocate(Index rows_, Index cols_, unsigned
m_workMatrix.resize(diagSize(), diagSize()); m_workMatrix.resize(diagSize(), diagSize());
if (cols() > rows()) m_qr_precond_morecols.allocate(*this); if (cols() > rows()) m_qr_precond_morecols.allocate(*this);
if (rows() > cols()) m_qr_precond_morerows.allocate(*this); if (rows() > cols()) m_qr_precond_morerows.allocate(*this);
if (rows() != cols()) m_scaledMatrix.resize(rows(), cols());
} }
template <typename MatrixType, int Options> template <typename MatrixType, int Options>
@ -699,9 +691,8 @@ JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute_impl(con
/*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
if (rows() != cols()) { if (rows() != cols()) {
m_scaledMatrix = matrix / scale; m_qr_precond_morecols.run(*this, matrix / scale);
m_qr_precond_morecols.run(*this, m_scaledMatrix); m_qr_precond_morerows.run(*this, matrix / scale);
m_qr_precond_morerows.run(*this, m_scaledMatrix);
} else { } else {
m_workMatrix = m_workMatrix =
matrix.template topLeftCorner<DiagSizeAtCompileTime, DiagSizeAtCompileTime>(diagSize(), diagSize()) / scale; matrix.template topLeftCorner<DiagSizeAtCompileTime, DiagSizeAtCompileTime>(diagSize(), diagSize()) / scale;

11
Eigen/src/SVD/SVDBase.h
View File

@ -125,7 +125,6 @@ class SVDBase : public SolverBase<SVDBase<Derived> > {
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef typename Eigen::internal::traits<SVDBase>::StorageIndex StorageIndex; typedef typename Eigen::internal::traits<SVDBase>::StorageIndex StorageIndex;
typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
static constexpr bool ShouldComputeFullU = internal::traits<Derived>::ShouldComputeFullU; static constexpr bool ShouldComputeFullU = internal::traits<Derived>::ShouldComputeFullU;
static constexpr bool ShouldComputeThinU = internal::traits<Derived>::ShouldComputeThinU; static constexpr bool ShouldComputeThinU = internal::traits<Derived>::ShouldComputeThinU;
@ -355,11 +354,11 @@ class SVDBase : public SolverBase<SVDBase<Derived> > {
m_isInitialized(false), m_isInitialized(false),
m_isAllocated(false), m_isAllocated(false),
m_usePrescribedThreshold(false), m_usePrescribedThreshold(false),
m_computeFullU(false), m_computeFullU(ShouldComputeFullU),
m_computeThinU(false), m_computeThinU(ShouldComputeThinU),
m_computeFullV(false), m_computeFullV(ShouldComputeFullV),
m_computeThinV(false), m_computeThinV(ShouldComputeThinV),
m_computationOptions(0), m_computationOptions(internal::traits<Derived>::Options),
m_nonzeroSingularValues(0), m_nonzeroSingularValues(0),
m_rows(RowsAtCompileTime), m_rows(RowsAtCompileTime),
m_cols(ColsAtCompileTime), m_cols(ColsAtCompileTime),