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λI as LDL<sup>T</sup>
+ // decompose T-λ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);
