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);
