From 9a324dc5463ea6329e0b2e59bb46073117d68168 Mon Sep 17 00:00:00 2001
From: Luc Maisonobe <luc@apache.org>
Date: Sun, 29 Nov 2009 21:21:50 +0000
Subject: [PATCH] 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
---
 .../math/linear/EigenDecompositionImpl.java   | 42 ++++++++++++++++++-
 1 file changed, 41 insertions(+), 1 deletion(-)

diff --git a/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java b/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java
index 20a3114b5..db1c82716 100644
--- a/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java
+++ b/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java
@@ -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;
     }