Avoid memory allocation in tridiagonalization_inplace_selector::run.

(cherry picked from commit a5a7faeb455efd7f6edb1138eda2e37546039b7d)
This commit is contained in:
Rasmus Munk Larsen 2021-08-06 20:48:10 +00:00
parent 1e9f623f3e
commit 4e0357c6dd
2 changed files with 16 additions and 11 deletions

View File

@ -125,6 +125,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
: m_eivec(), : m_eivec(),
m_eivalues(), m_eivalues(),
m_subdiag(), m_subdiag(),
m_hcoeffs(),
m_info(InvalidInput), m_info(InvalidInput),
m_isInitialized(false), m_isInitialized(false),
m_eigenvectorsOk(false) m_eigenvectorsOk(false)
@ -147,6 +148,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
: m_eivec(size, size), : m_eivec(size, size),
m_eivalues(size), m_eivalues(size),
m_subdiag(size > 1 ? size - 1 : 1), m_subdiag(size > 1 ? size - 1 : 1),
m_hcoeffs(size > 1 ? size - 1 : 1),
m_isInitialized(false), m_isInitialized(false),
m_eigenvectorsOk(false) m_eigenvectorsOk(false)
{} {}
@ -172,6 +174,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
: m_eivec(matrix.rows(), matrix.cols()), : m_eivec(matrix.rows(), matrix.cols()),
m_eivalues(matrix.cols()), m_eivalues(matrix.cols()),
m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1), m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
m_hcoeffs(matrix.cols() > 1 ? matrix.cols() - 1 : 1),
m_isInitialized(false), m_isInitialized(false),
m_eigenvectorsOk(false) m_eigenvectorsOk(false)
{ {
@ -378,6 +381,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
EigenvectorsType m_eivec; EigenvectorsType m_eivec;
RealVectorType m_eivalues; RealVectorType m_eivalues;
typename TridiagonalizationType::SubDiagonalType m_subdiag; typename TridiagonalizationType::SubDiagonalType m_subdiag;
typename TridiagonalizationType::CoeffVectorType m_hcoeffs;
ComputationInfo m_info; ComputationInfo m_info;
bool m_isInitialized; bool m_isInitialized;
bool m_eigenvectorsOk; bool m_eigenvectorsOk;
@ -450,7 +454,8 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
if(scale==RealScalar(0)) scale = RealScalar(1); if(scale==RealScalar(0)) scale = RealScalar(1);
mat.template triangularView<Lower>() /= scale; mat.template triangularView<Lower>() /= scale;
m_subdiag.resize(n-1); m_subdiag.resize(n-1);
internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors); m_hcoeffs.resize(n-1);
internal::tridiagonalization_inplace(mat, diag, m_subdiag, m_hcoeffs, computeEigenvectors);
m_info = internal::computeFromTridiagonal_impl(diag, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec); m_info = internal::computeFromTridiagonal_impl(diag, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);

View File

@ -425,12 +425,13 @@ struct tridiagonalization_inplace_selector;
* *
* \sa class Tridiagonalization * \sa class Tridiagonalization
*/ */
template<typename MatrixType, typename DiagonalType, typename SubDiagonalType> template<typename MatrixType, typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType>
EIGEN_DEVICE_FUNC EIGEN_DEVICE_FUNC
void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag,
CoeffVectorType& hcoeffs, bool extractQ)
{ {
eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1); eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ); tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, hcoeffs, extractQ);
} }
/** \internal /** \internal
@ -443,10 +444,9 @@ struct tridiagonalization_inplace_selector
typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType; typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
template<typename DiagonalType, typename SubDiagonalType> template<typename DiagonalType, typename SubDiagonalType>
static EIGEN_DEVICE_FUNC static EIGEN_DEVICE_FUNC
void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, CoeffVectorType& hCoeffs, bool extractQ)
{ {
CoeffVectorType hCoeffs(mat.cols()-1); tridiagonalization_inplace(mat, hCoeffs);
tridiagonalization_inplace(mat,hCoeffs);
diag = mat.diagonal().real(); diag = mat.diagonal().real();
subdiag = mat.template diagonal<-1>().real(); subdiag = mat.template diagonal<-1>().real();
if(extractQ) if(extractQ)
@ -466,8 +466,8 @@ struct tridiagonalization_inplace_selector<MatrixType,3,false>
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::RealScalar RealScalar;
template<typename DiagonalType, typename SubDiagonalType> template<typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType>
static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, CoeffVectorType&, bool extractQ)
{ {
using std::sqrt; using std::sqrt;
const RealScalar tol = (std::numeric_limits<RealScalar>::min)(); const RealScalar tol = (std::numeric_limits<RealScalar>::min)();
@ -511,9 +511,9 @@ struct tridiagonalization_inplace_selector<MatrixType,1,IsComplex>
{ {
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
template<typename DiagonalType, typename SubDiagonalType> template<typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType>
static EIGEN_DEVICE_FUNC static EIGEN_DEVICE_FUNC
void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, bool extractQ) void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, CoeffVectorType&, bool extractQ)
{ {
diag(0,0) = numext::real(mat(0,0)); diag(0,0) = numext::real(mat(0,0));
if(extractQ) if(extractQ)