commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r833740 - 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 Sat, 07 Nov 2009 19:57:03 GMT
Author: luc
Date: Sat Nov  7 19:57:02 2009
New Revision: 833740

URL: http://svn.apache.org/viewvc?rev=833740&view=rev
Log:
fixed yet another eigen decomposition bug identified, debugged and fixed by Dimitri!
Many thanks to him.

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=833740&r1=833739&r2=833740&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
Sat Nov  7 19:57:02 2009
@@ -837,7 +837,7 @@
         }
 
         // initial checks for splits (see Parlett & Marques section 3.3)
-        flipIfWarranted(n, 2);
+        flipEveryOtherIfWarranted(n);
 
         // two iterations with Li's test for initial splits
         initialSplits(n);
@@ -1051,7 +1051,7 @@
 
         // step 2: flip array if needed
         if ((dMin <= 0) || (deflatedEnd < end)) {
-            if (flipIfWarranted(deflatedEnd, 1)) {
+            if (flipAllIfWarranted(deflatedEnd)) {
                 dMin2 = Math.min(dMin2, work[l - 1]);
                 work[l - 1] =
                     Math.min(work[l - 1],
@@ -1123,27 +1123,59 @@
     }
 
     /**
-     * Flip qd array if warranted.
+     * Flip all elements of qd array if warranted.
      * @param n number of rows in the block
-     * @param step within the array (1 for flipping all elements, 2 for flipping
-     * only every other element)
      * @return true if qd array was flipped
      */
-    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);
-            for (int i = 0; i < j; i += 4) {
-                for (int k = 0; k < 4; k += step) {
-                    final double tmp = work[i + k];
-                    work[i + k] = work[j - k];
-                    work[j - k] = tmp;
-                }
-                j -= 4;
+    private boolean flipAllIfWarranted(final int n) {
+        if (1.5 * work[pingPong] >= work[4 * (n - 1) + pingPong]) {
+            return false;
+        }
+
+        int j = 4 * (n - 1);
+        for (int i = 0; i < j; i += 4) {
+            final double tmp1 = work[i];
+            work[i] = work[j];
+            work[j] = tmp1;
+            final double tmp2 = work[i+1];
+            work[i+1] = work[j+1];
+            work[j+1] = tmp2;
+            final double tmp3 = work[i+2];
+            work[i+2] = work[j-2];
+            work[j-2] = tmp3;
+            final double tmp4 = work[i+3];
+            work[i+3] = work[j-1];
+            work[j-1] = tmp4;
+            j -= 4;
+        }
+
+        return true;
+
+    }
+
+    /**
+     * Flip every other elements of qd array if warranted.
+     * @param n number of rows in the block
+     * @return true if qd array was flipped
+     */
+    private boolean flipEveryOtherIfWarranted(final int n) {
+        if (1.5 * work[pingPong] >= work[4 * (n - 1) + pingPong]) {
+            return false;
+        }
+
+        // flip array
+        int j = 4 * (n - 1);
+        for (int i = 0; i < j; i += 4) {
+            for (int k = 0; k < 4; k += 2) {
+                final double tmp = work[i + k];
+                work[i + k] = work[j - k];
+                work[j - k] = tmp;
             }
-            return true;
+            j -= 4;
         }
-        return false;
+
+        return true;
+
     }
 
     /**

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=833740&r1=833739&r2=833740&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
Sat Nov  7 19:57:02 2009
@@ -188,6 +188,50 @@
 
     }
 
+    public void testMathpbx03() {
+
+        double[] mainTridiagonal = {
+            1809.0978259647177,3395.4763425956166,1832.1894584712693,3804.364873592377,
+            806.0482458637571,2403.656427234185,28.48691431556015
+        };
+        double[] secondaryTridiagonal = {
+            -656.8932064545833,-469.30804108920734,-1021.7714889369421,
+            -1152.540497328983,-939.9765163817368,-12.885877015422391
+        };
+
+        // the reference values have been computed using routine DSTEMR
+        // from the fortran library LAPACK version 3.2.1
+        double[] refEigenValues = {
+            4603.121913685183245,3691.195818048970978,2743.442955402465032,1657.596442107321764,
+            1336.797819095331306,30.129865209677519,17.035352085224986
+        };
+
+        RealVector[] refEigenVectors = {
+            new ArrayRealVector(new double[] {-0.036249830202337,0.154184732411519,-0.346016328392363,0.867540105133093,-0.294483395433451,0.125854235969548,-0.000354507444044}),
+            new ArrayRealVector(new double[] {-0.318654191697157,0.912992309960507,-0.129270874079777,-0.184150038178035,0.096521712579439,-0.070468788536461,0.000247918177736}),
+            new ArrayRealVector(new double[] {-0.051394668681147,0.073102235876933,0.173502042943743,-0.188311980310942,-0.327158794289386,0.905206581432676,-0.004296342252659}),
+            new ArrayRealVector(new double[] {0.838150199198361,0.193305209055716,-0.457341242126146,-0.166933875895419,0.094512811358535,0.119062381338757,-0.000941755685226}),
+            new ArrayRealVector(new double[] {0.438071395458547,0.314969169786246,0.768480630802146,0.227919171600705,-0.193317045298647,-0.170305467485594,0.001677380536009}),
+            new ArrayRealVector(new double[] {-0.003726503878741,-0.010091946369146,-0.067152015137611,-0.113798146542187,-0.313123000097908,-0.118940107954918,0.932862311396062}),
+            new ArrayRealVector(new double[] {0.009373003194332,0.025570377559400,0.170955836081348,0.291954519805750,0.807824267665706,0.320108347088646,0.360202112392266}),
+        };
+
+        // 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-4);
+            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