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

Fix 2240, 2620

Browse Source
This commit is contained in:
Charles Schlosser 2023-03-06 23:11:06 +00:00 committed by Rasmus Munk Larsen
parent d670039309
commit cb8e6d4975
2 changed files with 21 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
View File

@ -108,6 +108,7 @@ template<typename MatrixType_> class SelfAdjointEigenSolver
* This is a column vector with entries of type #RealScalar. * This is a column vector with entries of type #RealScalar.
* The length of the vector is the size of \p MatrixType_. * The length of the vector is the size of \p MatrixType_.
*/ */
typedef typename internal::plain_col_type<MatrixType, Scalar>::type VectorType;
typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType; typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;
typedef Tridiagonalization<MatrixType> TridiagonalizationType; typedef Tridiagonalization<MatrixType> TridiagonalizationType;
typedef typename TridiagonalizationType::SubDiagonalType SubDiagonalType; typedef typename TridiagonalizationType::SubDiagonalType SubDiagonalType;
@ -125,6 +126,7 @@ template<typename MatrixType_> class SelfAdjointEigenSolver
EIGEN_DEVICE_FUNC EIGEN_DEVICE_FUNC
SelfAdjointEigenSolver() SelfAdjointEigenSolver()
: m_eivec(), : m_eivec(),
m_workspace(),
m_eivalues(), m_eivalues(),
m_subdiag(), m_subdiag(),
m_hcoeffs(), m_hcoeffs(),
@ -148,6 +150,7 @@ template<typename MatrixType_> class SelfAdjointEigenSolver
EIGEN_DEVICE_FUNC EIGEN_DEVICE_FUNC
explicit SelfAdjointEigenSolver(Index size) explicit SelfAdjointEigenSolver(Index size)
: m_eivec(size, size), : m_eivec(size, size),
m_workspace(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_hcoeffs(size > 1 ? size - 1 : 1),
@ -174,6 +177,7 @@ template<typename MatrixType_> class SelfAdjointEigenSolver
EIGEN_DEVICE_FUNC EIGEN_DEVICE_FUNC
explicit SelfAdjointEigenSolver(const EigenBase<InputType>& matrix, int options = ComputeEigenvectors) explicit SelfAdjointEigenSolver(const EigenBase<InputType>& matrix, int options = ComputeEigenvectors)
: m_eivec(matrix.rows(), matrix.cols()), : m_eivec(matrix.rows(), matrix.cols()),
m_workspace(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_hcoeffs(matrix.cols() > 1 ? matrix.cols() - 1 : 1),
@ -377,6 +381,7 @@ template<typename MatrixType_> class SelfAdjointEigenSolver
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar) EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
EigenvectorsType m_eivec; EigenvectorsType m_eivec;
VectorType m_workspace;
RealVectorType m_eivalues; RealVectorType m_eivalues;
typename TridiagonalizationType::SubDiagonalType m_subdiag; typename TridiagonalizationType::SubDiagonalType m_subdiag;
typename TridiagonalizationType::CoeffVectorType m_hcoeffs; typename TridiagonalizationType::CoeffVectorType m_hcoeffs;
@ -451,7 +456,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
mat.template triangularView<Lower>() /= scale; mat.template triangularView<Lower>() /= scale;
m_subdiag.resize(n-1); m_subdiag.resize(n-1);
m_hcoeffs.resize(n-1); m_hcoeffs.resize(n-1);
internal::tridiagonalization_inplace(mat, diag, m_subdiag, m_hcoeffs, computeEigenvectors); internal::tridiagonalization_inplace(mat, diag, m_subdiag, m_hcoeffs, m_workspace, 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);

28
Eigen/src/Eigenvalues/Tridiagonalization.h
View File

@ -427,13 +427,13 @@ struct tridiagonalization_inplace_selector;
* *
* \sa class Tridiagonalization * \sa class Tridiagonalization
*/ */
template<typename MatrixType, typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType> template<typename MatrixType, typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType, typename WorkSpaceType>
EIGEN_DEVICE_FUNC EIGEN_DEVICE_FUNC
void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag,
CoeffVectorType& hcoeffs, bool extractQ) CoeffVectorType& hcoeffs, WorkSpaceType& workspace, 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, hcoeffs, extractQ); tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, hcoeffs, workspace, extractQ);
} }
/** \internal /** \internal
@ -443,17 +443,19 @@ template<typename MatrixType, int Size, bool IsComplex>
struct tridiagonalization_inplace_selector struct tridiagonalization_inplace_selector
{ {
typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType; typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
template<typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType> template<typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType, typename WorkSpaceType>
static EIGEN_DEVICE_FUNC static EIGEN_DEVICE_FUNC
void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, CoeffVectorType& hCoeffs, bool extractQ) void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, CoeffVectorType& hCoeffs, WorkSpaceType& workspace, bool extractQ)
{ {
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) {
mat = HouseholderSequenceType(mat, hCoeffs.conjugate()) HouseholderSequenceType hh(mat, hCoeffs.conjugate());
.setLength(mat.rows() - 1) hh.setLength(mat.rows() - 1);
.setShift(1); hh.setShift(1);
hh.evalTo(mat, workspace);
}
} }
}; };
@ -467,8 +469,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, typename CoeffVectorType> template<typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType, typename WorkSpaceType>
static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, CoeffVectorType&, bool extractQ) static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, CoeffVectorType&, WorkSpaceType&, bool extractQ)
{ {
using std::sqrt; using std::sqrt;
const RealScalar tol = (std::numeric_limits<RealScalar>::min)(); const RealScalar tol = (std::numeric_limits<RealScalar>::min)();
@ -512,9 +514,9 @@ struct tridiagonalization_inplace_selector<MatrixType,1,IsComplex>
{ {
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
template<typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType> template<typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType, typename WorkSpaceType>
static EIGEN_DEVICE_FUNC static EIGEN_DEVICE_FUNC
void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, CoeffVectorType&, bool extractQ) void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, CoeffVectorType&, WorkSpaceType&, bool extractQ)
{ {
diag(0,0) = numext::real(mat(0,0)); diag(0,0) = numext::real(mat(0,0));
if(extractQ) if(extractQ)