[prev in list] [next in list] [prev in thread] [next in thread] 

List:       jakarta-commons-dev
Subject:    svn commit: r885268 -
From:       luc () apache ! org
Date:       2009-11-29 21:21:52
Message-ID: 20091129212152.891C223889CB () eris ! apache ! org
[Download RAW message or body]

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/c \
ommons/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;
     }


[prev in list] [next in list] [prev in thread] [next in thread]