svn commit: r833740 - in /commons/proper/math/trunk/src
|
|
luc-17
|
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); |
||||||||||||||||
| Free Embeddable Forum Powered by Nabble | Help |
