commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r832577 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java site/xdoc/changes.xml test/java/org/apache/commons/math/linear/EigenDecompositionImplTest.java
Date Tue, 03 Nov 2009 22:04:13 GMT
Author: luc
Date: Tue Nov  3 22:04:08 2009
New Revision: 832577

URL: http://svn.apache.org/viewvc?rev=832577&view=rev
Log:
fixed an ArrayIndexOutOfBoundsException
Kudos to Dimitri who debugged this mess of fortran/java array indices translation
JIRA: MATH-308

Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java
    commons/proper/math/trunk/src/site/xdoc/changes.xml
    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=832577&r1=832576&r2=832577&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
Tue Nov  3 22:04:08 2009
@@ -595,8 +595,12 @@
         }
 
         final double dCurrent = main[m - 1];
-        work[lowerStart + m - 1] = dCurrent - eCurrent;
-        work[upperStart + m - 1] = dCurrent + eCurrent;
+        final double lower = dCurrent - eCurrent;
+        work[lowerStart + m - 1] = lower;
+        lowerSpectra = Math.min(lowerSpectra, lower);
+        final double upper = dCurrent + eCurrent;
+        work[upperStart + m - 1] = upper;
+        upperSpectra = Math.max(upperSpectra, upper);
         minPivot = MathUtils.SAFE_MIN * Math.max(1.0, eMax * eMax);
 
     }
@@ -647,7 +651,7 @@
 
                 tau = (range[1] - range[0]) * MathUtils.EPSILON * n + 2 * minPivot;
 
-                // decompose T&lambda;I as LDL<sup>T</sup>
+                // decompose T-&lambda;I as LDL<sup>T</sup>
                 ldlTDecomposition(lambda, begin, n);
 
                 // apply general dqd/dqds method
@@ -899,8 +903,8 @@
                     diagMax    = work[4 * i0];
                     offDiagMin = work[4 * i0 + 2];
                     double previousEMin = work[4 * i0 + 3];
-                    for (int i = 4 * i0; i < 4 * n0 - 11; i += 4) {
-                        if ((work[i + 3] <= TOLERANCE_2 * work[i]) &&
+                    for (int i = 4 * i0; i < 4 * n0 - 16; i += 4) {
+                        if ((work[i + 3] <= TOLERANCE_2 * work[i]) ||
                             (work[i + 2] <= TOLERANCE_2 * sigma)) {
                             // insert a split
                             work[i + 2]  = -sigma;
@@ -1537,7 +1541,7 @@
                 double a2 = (work[np - 8] / b2) * (1 + work[np - 4] / b1);
 
                 // approximate contribution to norm squared from i < nn-2.
-                if (end - start > 2) {
+                if (end - start > 3) {
                     b2 = work[nn - 13] / work[nn - 15];
                     a2 = a2 + b2;
                     for (int i4 = nn - 17; i4 >= 4 * start + 2 + pingPong; i4 -= 4) {

Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=832577&r1=832576&r2=832577&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Tue Nov  3 22:04:08 2009
@@ -39,6 +39,10 @@
   </properties>
   <body>
     <release version="2.1" date="TBD" description="TBD">
+      <action dev="luc" type="fix" issue="MATH-308" due-to="Dimitri Pourbaix">
+        Fixed an ArrayIndexOutOfBoundsException in eigen decomposition. Kudos to Dimitri
+        for debugging this, it was really difficult.
+      </action>
       <action dev="psteitz" type="fix" issue="MATH-309" due-to="Mikkel Meyer Andersen">
         Fixed parameter test in RandomDataImpl#nextExponential. The method now throws
         IllegalArgumentException for mean = 0.

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=832577&r1=832576&r2=832577&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
Tue Nov  3 22:04:08 2009
@@ -108,6 +108,40 @@
         assertEquals(0.1, ed.getRealEigenvalue(3), 1.0e-15);
     }
 
+    // the following test triggered an ArrayIndexOutOfBoundsException in commons-math 2.0
+    public void testMath308() {
+
+        double[] mainTridiagonal = {
+            22.330154644539597, 46.65485522478641, 17.393672330044705, 54.46687435351116,
80.17800767709437
+        };
+        double[] secondaryTridiagonal = {
+            13.04450406501361, -5.977590941539671, 2.9040909856707517, 7.1570352792841225
+        };
+
+        // the reference values have been computed using routine DSTEMR
+        // from the fortran library LAPACK version 3.2.1
+        double[] refEigenValues = {
+            82.044413207204002, 53.456697699894512, 52.536278520113882, 18.847969733754262,
14.138204224043099
+        };
+        RealVector[] refEigenVectors = {
+            new ArrayRealVector(new double[] { -0.000462690386766, -0.002118073109055,  0.011530080757413,
 0.252322434584915,  0.967572088232592 }),
+            new ArrayRealVector(new double[] {  0.314647769490148,  0.750806415553905, -0.167700312025760,
-0.537092972407375,  0.143854968127780 }),
+            new ArrayRealVector(new double[] {  0.222368839324646,  0.514921891363332, -0.021377019336614,
 0.801196801016305, -0.207446991247740 }),
+            new ArrayRealVector(new double[] {  0.713933751051495, -0.190582113553930,  0.671410443368332,
-0.056056055955050,  0.006541576993581 }),
+            new ArrayRealVector(new double[] {  0.584677060845929, -0.367177264979103, -0.721453187784497,
 0.052971054621812, -0.005740715188257 })
+        };
+
+        EigenDecomposition decomposition =
+            new EigenDecompositionImpl(mainTridiagonal, secondaryTridiagonal, MathUtils.SAFE_MIN);
+
+        double[] eigenValues = decomposition.getRealEigenvalues();
+        for (int i = 0; i < refEigenValues.length; ++i) {
+            assertEquals(refEigenValues[i], eigenValues[i], 1.0e-5);
+            assertEquals(0, refEigenVectors[i].subtract(decomposition.getEigenvector(i)).getNorm(),
2.0e-7);
+        }
+
+    }
+
     /** test a matrix already in tridiagonal form. */
     public void testTridiagonal() {
         Random r = new Random(4366663527842l);



Mime
View raw message