Author: luc
Date: Sun Nov 29 21:21:50 2009
New Revision: 885268
URL: http://svn.apache.org/viewvc?rev=885268&view=rev
Log:
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
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java?rev=885268&r1=885267&r2=885268&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java
Sun Nov 29 21:21:50 2009
@@ -1832,14 +1832,35 @@
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 @@
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;
}
