You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2009/11/03 23:04:13 UTC
svn commit: r832577 - 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
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);
Re: svn commit: r832577 - 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
Posted by Phil Steitz <ph...@gmail.com>.
Kudos indeed. Nice work!
Phil
luc@apache.org wrote:
> 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);
>
>
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org