Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

fixed some NaN appearing in eigenvectors when null pivots occurred in dstqds or dqds algorithms 

this is a partial fix for MATH-297 but not a complete one as for example computing the
eigendecomposition if identity leads to three times the same vector ...
JIRA: MATH-297

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@885268 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2009-11-29 21:21:50 +00:00
parent abb1b3fdd6
commit 9a324dc546
1 changed files with 41 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java
View File

@ -1832,14 +1832,35 @@ public class EigenDecompositionImpl implements EigenDecomposition {
for (int i = 0; i < nM1; ++i) {
final double di = d[i];
final double li = l[i];
final double ldi = li * di;
final double diP1 = di + si;
final double liP1 = li * di / diP1;
final double liP1 = ldi / diP1;
work[sixI] = si;
work[sixI + 1] = diP1;
work[sixI + 2] = liP1;
si = li * liP1 * si - lambda;
sixI += 6;
}
if (Double.isNaN(si)) {
// one of the pivot was null, use a slower but safer version of dstqds
si = -lambda;
sixI = 0;
for (int i = 0; i < nM1; ++i) {
final double di = d[i];
final double li = l[i];
final double ldi = li * di;
double diP1 = di + si;
if (Math.abs(diP1) < minPivot) {
diP1 = -minPivot;
}
final double liP1 = ldi / diP1;
work[sixI] = si;
work[sixI + 1] = diP1;
work[sixI + 2] = liP1;
si = li * ((liP1 == 0) ? li * di : liP1 * si) - lambda;
sixI += 6;
}
}
work[6 * nM1 + 1] = d[nM1] + si;
work[6 * nM1] = si;
}
@ -1868,6 +1889,25 @@ public class EigenDecompositionImpl implements EigenDecomposition {
pi = pi * t - lambda;
sixI -= 6;
}
if (Double.isNaN(pi)) {
// one of the pivot was null, use a slower but safer version of dqds
pi = d[nM1] - lambda;
sixI = 6 * (nM1 - 1);
for (int i = nM1 - 1; i >= 0; --i) {
final double di = d[i];
final double li = l[i];
double diP1 = di * li * li + pi;
if (Math.abs(diP1) < minPivot) {
diP1 = -minPivot;
}
final double t = di / diP1;
work[sixI + 9] = pi;
work[sixI + 10] = diP1;
work[sixI + 5] = li * t;
pi = ((t == 0) ? di : pi * t) - lambda;
sixI -= 6;
}
}
work[3] = pi;
work[4] = pi;
}