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

Fix possible underflow issues in SelfAdjointEigenSolver

Browse Source
This commit is contained in:
Gael Guennebaud 2012-07-10 09:51:26 +02:00
parent ff12a6cd43
commit a2c3003be2
1 changed files with 8 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -751,15 +751,15 @@ namespace internal {
template<int StorageOrder,typename RealScalar, typename Scalar, typename Index> template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n) static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
{ {
// NOTE this version avoids over & underflow, however since the matrix is prescaled, overflow cannot occur,
// and underflows should be meaningless anyway. So I don't any reason to enable this version, but I keep
// it here for reference:
// RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
// RealScalar e = subdiag[end-1];
// RealScalar mu = diag[end] - (e / (td + (td>0 ? 1 : -1))) * (e / hypot(td,e));
RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5); RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
RealScalar e2 = abs2(subdiag[end-1]); RealScalar e = subdiag[end-1];
RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2)); // 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 = 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] - (e / (td + (td>0 ? 1 : -1))) * (e / hypot(td,e));
RealScalar x = diag[start] - mu; RealScalar x = diag[start] - mu;
RealScalar z = subdiag[start]; RealScalar z = subdiag[start];
for (Index k = start; k < end; ++k) for (Index k = start; k < end; ++k)