commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From billbar...@apache.org
Subject svn commit: r826550 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java test/java/org/apache/commons/math/linear/EigenDecompositionImplTest.java
Date Sun, 18 Oct 2009 23:29:37 GMT
Author: billbarker
Date: Sun Oct 18 23:29:37 2009
New Revision: 826550

URL: http://svn.apache.org/viewvc?rev=826550&view=rev
Log:
Fix for possible zero divide on an indefinite matrix

Fix for: MATH-297

Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/EigenDecompositionImplTest.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=826550&r1=826549&r2=826550&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 Oct 18 23:29:37 2009
@@ -1699,19 +1699,21 @@
         // perform an initial non-shifted LDLt decomposition
         final double[] d = new double[m];
         final double[] l = new double[m - 1];
-        double di = main[0];
+        // avoid zero divide on indefinite matrix
+        final double mu = realEigenvalues[m-1] <= 0 && realEigenvalues[0] >
0 ? 0.5-realEigenvalues[m-1] : 0;
+        double di = main[0]+mu;
         d[0] = di;
         for (int i = 1; i < m; ++i) {
             final double eiM1  = secondary[i - 1];
             final double ratio = eiM1 / di;
-            di       = main[i] - eiM1 * ratio;
+            di       = main[i] - eiM1 * ratio + mu;
             l[i - 1] = ratio;
             d[i]     = di;
         }
 
         // compute eigenvectors
         for (int i = 0; i < m; ++i) {
-            eigenvectors[i] = findEigenvector(realEigenvalues[i], d, l);
+            eigenvectors[i] = findEigenvector(realEigenvalues[i]+mu, d, l);
         }
 
     }

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/EigenDecompositionImplTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/EigenDecompositionImplTest.java?rev=826550&r1=826549&r2=826550&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/EigenDecompositionImplTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/EigenDecompositionImplTest.java
Sun Oct 18 23:29:37 2009
@@ -244,6 +244,24 @@
     }
 
     /**
+     * Verifies operation on indefinite matrix
+     */
+    public void testZeroDivide() {
+        RealMatrix indefinite = MatrixUtils.createRealMatrix(new double [][] {
+                { 0.0, 1.0, -1.0 }, 
+                { 1.0, 1.0, 0.0 }, 
+                { -1.0,0.0, 1.0 }        
+        });
+        EigenDecomposition ed = new EigenDecompositionImpl(indefinite, MathUtils.SAFE_MIN);
+        checkEigenValues((new double[] {2, 1, -1}), ed, 1E-12);
+        double isqrt3 = 1/Math.sqrt(3.0);
+        checkEigenVector((new double[] {isqrt3,isqrt3,-isqrt3}), ed, 1E-12);
+        double isqrt2 = 1/Math.sqrt(2.0);
+        checkEigenVector((new double[] {0.0,-isqrt2,-isqrt2}), ed, 1E-12);
+        double isqrt6 = 1/Math.sqrt(6.0);
+        checkEigenVector((new double[] {2*isqrt6,-isqrt6,isqrt6}), ed, 1E-12);
+    }
+    /**
      * Verifies that the given EigenDecomposition has eigenvalues equivalent to
      * the targetValues, ignoring the order of the values and allowing
      * values to differ by tolerance.
@@ -257,6 +275,7 @@
         }
     }
 
+    
     /**
      * Returns true iff there is an entry within tolerance of value in
      * searchArray.



Mime
View raw message