2
0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-02-11 18:00:51 +08:00

Fix SelfAdjoingEigenSolver ()

Adjust the relaxation step to use the condition
```
abs(subdiag[i]) <= epsilon * sqrt(abs(diag[i]) + abs(diag[i+1]))
```
for setting the subdiagonal entry to zero.

Also adjust Wilkinson shift for small `e = subdiag[end-1]` -
I couldn't find a reference for the original, and it was not
consistent with the Wilkinson definition.

Fixes .
This commit is contained in:
Antonio Sanchez 2021-03-25 11:08:19 +00:00
parent 3ddc0974ce
commit 90187a33e1

View File

@ -498,8 +498,6 @@ template<typename MatrixType, typename DiagType, typename SubDiagType>
EIGEN_DEVICE_FUNC
ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec)
{
EIGEN_USING_STD(abs);
ComputationInfo info;
typedef typename MatrixType::Scalar Scalar;
@ -510,15 +508,23 @@ ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag
typedef typename DiagType::RealScalar RealScalar;
const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon();
const RealScalar precision_inv = RealScalar(2)/NumTraits<RealScalar>::epsilon();
while (end>0)
{
for (Index i = start; i<end; ++i)
if (internal::isMuchSmallerThan(abs(subdiag[i]),(abs(diag[i])+abs(diag[i+1])),precision) || abs(subdiag[i]) <= considerAsZero)
subdiag[i] = 0;
for (Index i = start; i<end; ++i) {
if (numext::abs(subdiag[i]) < considerAsZero) {
subdiag[i] = RealScalar(0);
} else {
// abs(subdiag[i]) <= epsilon * sqrt(abs(diag[i]) + abs(diag[i+1]))
// Scaled to prevent underflows.
const RealScalar scaled_subdiag = precision_inv * subdiag[i];
if (scaled_subdiag * scaled_subdiag <= (numext::abs(diag[i])+numext::abs(diag[i+1]))) {
subdiag[i] = RealScalar(0);
}
}
}
// find the largest unreduced block
// find the largest unreduced block at the end of the matrix.
while (end>0 && subdiag[end-1]==RealScalar(0))
{
end--;
@ -821,32 +827,38 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
}
namespace internal {
// Francis implicit QR step.
template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
EIGEN_DEVICE_FUNC
static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
{
EIGEN_USING_STD(abs);
// Wilkinson Shift.
RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
RealScalar e = subdiag[end-1];
// Note that thanks to scaling, e^2 or td^2 cannot overflow, however they can still
// underflow thus leading to inf/NaN values when using the following commented code:
// RealScalar e2 = numext::abs2(subdiag[end-1]);
// RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
// RealScalar e2 = numext::abs2(subdiag[end-1]);
// RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
// This explain the following, somewhat more complicated, version:
RealScalar mu = diag[end];
if(td==RealScalar(0))
mu -= abs(e);
else
{
RealScalar e2 = numext::abs2(subdiag[end-1]);
RealScalar h = numext::hypot(td,e);
if(e2==RealScalar(0)) mu -= (e / (td + (td>RealScalar(0) ? RealScalar(1) : RealScalar(-1)))) * (e / h);
else mu -= e2 / (td + (td>RealScalar(0) ? h : -h));
if(td==RealScalar(0)) {
mu -= numext::abs(e);
} else if (e != RealScalar(0)) {
const RealScalar e2 = numext::abs2(e);
const RealScalar h = numext::hypot(td,e);
if(e2 == RealScalar(0)) {
mu -= e / ((td + (td>RealScalar(0) ? h : -h)) / e);
} else {
mu -= e2 / (td + (td>RealScalar(0) ? h : -h));
}
}
RealScalar x = diag[start] - mu;
RealScalar z = subdiag[start];
for (Index k = start; k < end; ++k)
// If z ever becomes zero, the Givens rotation will be the identity and
// z will stay zero for all future iterations.
for (Index k = start; k < end && z != RealScalar(0); ++k)
{
JacobiRotation<RealScalar> rot;
rot.makeGivens(x, z);
@ -859,12 +871,11 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta
diag[k+1] = rot.s() * sdk + rot.c() * dkp1;
subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
if (k > start)
subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
// "Chasing the bulge" to return to triangular form.
x = subdiag[k];
if (k < end - 1)
{
z = -rot.s() * subdiag[k+1];