commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r833433 - 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 Fri, 06 Nov 2009 15:11:57 GMT
Author: luc
Date: Fri Nov  6 15:11:57 2009
New Revision: 833433

URL: http://svn.apache.org/viewvc?rev=833433&view=rev
Log:
Fixed an index computation error in eigen decomposition.
Once again, kudos to Dimitri for debugging this.
JIRA: MATH-318

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=833433&r1=833432&r2=833433&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
Fri Nov  6 15:11:57 2009
@@ -1132,7 +1132,7 @@
     private boolean flipIfWarranted(final int n, final int step) {
         if (1.5 * work[pingPong] < work[4 * (n - 1) + pingPong]) {
             // flip array
-            int j = 4 * n - 1;
+            int j = 4 * (n - 1);
             for (int i = 0; i < j; i += 4) {
                 for (int k = 0; k < 4; k += step) {
                     final double tmp = work[i + k];

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=833433&r1=833432&r2=833433&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Fri Nov  6 15:11:57 2009
@@ -39,6 +39,10 @@
   </properties>
   <body>
     <release version="2.1" date="TBD" description="TBD">
+      <action dev="luc" type="fix" issue="MATH-318" due-to="Dimitri Pourbaix">
+        Fixed an index computation error in eigen decomposition. Once again, kudos to Dimitri
+        for debugging this.
+      </action>
       <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.

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=833433&r1=833432&r2=833433&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
Fri Nov  6 15:11:57 2009
@@ -142,6 +142,52 @@
 
     }
 
+    public void testMathpbx02() {
+
+        double[] mainTridiagonal = {
+        	  7484.860960227216, 18405.28129035345, 13855.225609560746,
+        	 10016.708722343366, 559.8117399576674, 6750.190788301587, 
+        	    71.21428769782159
+        };
+        double[] secondaryTridiagonal = {
+        	 -4175.088570476366,1975.7955858241994,5193.178422374075, 
+        	  1995.286659169179,75.34535882933804,-234.0808002076056
+        };
+
+        // the reference values have been computed using routine DSTEMR
+        // from the fortran library LAPACK version 3.2.1
+        double[] refEigenValues = {
+        		20654.744890306974412,16828.208208485466457,
+        		6893.155912634994820,6757.083016675340332,
+        		5887.799885688558788,64.309089923240379,
+        		57.992628792736340
+        };
+        RealVector[] refEigenVectors = {
+        		new ArrayRealVector(new double[] {-0.270356342026904, 0.852811091326997, 0.399639490702077,
0.198794657813990, 0.019739323307666, 0.000106983022327, -0.000001216636321}),
+        		new ArrayRealVector(new double[] {0.179995273578326,-0.402807848153042,0.701870993525734,0.555058211014888,0.068079148898236,0.000509139115227,-0.000007112235617}),
+        		new ArrayRealVector(new double[] {-0.399582721284727,-0.056629954519333,-0.514406488522827,0.711168164518580,0.225548081276367,0.125943999652923,-0.004321507456014}),
+        		new ArrayRealVector(new double[] {0.058515721572821,0.010200130057739,0.063516274916536,-0.090696087449378,-0.017148420432597,0.991318870265707,-0.034707338554096}),
+        		new ArrayRealVector(new double[] {0.855205995537564,0.327134656629775,-0.265382397060548,0.282690729026706,0.105736068025572,-0.009138126622039,0.000367751821196}),
+        		new ArrayRealVector(new double[] {-0.002913069901144,-0.005177515777101,0.041906334478672,-0.109315918416258,0.436192305456741,0.026307315639535,0.891797507436344}),
+        		new ArrayRealVector(new double[] {-0.005738311176435,-0.010207611670378,0.082662420517928,-0.215733886094368,0.861606487840411,-0.025478530652759,-0.451080697503958})
+        };
+
+        // the following line triggers the exception
+        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-3);
+            if (refEigenVectors[i].dotProduct(decomposition.getEigenvector(i)) < 0) {
+                assertEquals(0, refEigenVectors[i].add(decomposition.getEigenvector(i)).getNorm(),
1.0e-5);
+            } else {
+                assertEquals(0, refEigenVectors[i].subtract(decomposition.getEigenvector(i)).getNorm(),
1.0e-5);
+            }
+        }
+
+    }
+
     /** test a matrix already in tridiagonal form. */
     public void testTridiagonal() {
         Random r = new Random(4366663527842l);



Mime
View raw message