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 2010/01/29 18:08:28 UTC
svn commit: r904561 - 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: Fri Jan 29 17:08:28 2010
New Revision: 904561
URL: http://svn.apache.org/viewvc?rev=904561&view=rev
Log:
fixed a spurious exception in EigenDecompositionImpl when a 3x3 block had two identical eigenvalues
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=904561&r1=904560&r2=904561&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 Fri Jan 29 17:08:28 2010
@@ -771,23 +771,33 @@
// solve cubic equation
final double b2 = b * b;
+ final double beta = b / 3;
final double q = (3 * c - b2) / 9;
final double r = ((9 * c - 2 * b2) * b - 27 * d) / 54;
final double delta = q * q * q + r * r;
- if (delta >= 0) {
- // in fact, there are solutions to the equation, but in the context
- // of symmetric realEigenvalues problem, there should be three distinct
- // real roots, so we throw an error if this condition is not met
+ double z0;
+ double z1;
+ double z2;
+ if (delta > 0) {
+ // there are two complex solutions, we cannot handle this
throw new InvalidMatrixException("cannot solve degree {0} equation", 3);
+ } else if (delta < 0) {
+ // three different real roots
+ final double sqrtMq = Math.sqrt(-q);
+ final double theta = Math.acos(r / (-q * sqrtMq));
+ final double alpha = 2 * sqrtMq;
+ z0 = alpha * Math.cos(theta / 3) - beta;
+ z1 = alpha * Math.cos((theta + 2 * Math.PI) / 3) - beta;
+ z2 = alpha * Math.cos((theta + 4 * Math.PI) / 3) - beta;
+ } else {
+ // three real roots, two of which are equal
+ final double cbrtR = Math.cbrt(r);
+ z0 = 2 * cbrtR - beta;
+ z1 = -cbrtR - beta;
+ z2 = z1;
}
- final double sqrtMq = Math.sqrt(-q);
- final double theta = Math.acos(r / (-q * sqrtMq));
- final double alpha = 2 * sqrtMq;
- final double beta = b / 3;
-
- double z0 = alpha * Math.cos(theta / 3) - beta;
- double z1 = alpha * Math.cos((theta + 2 * Math.PI) / 3) - beta;
- double z2 = alpha * Math.cos((theta + 4 * Math.PI) / 3) - beta;
+
+ // sort the eigenvalues
if (z0 < z1) {
final double t = z0;
z0 = z1;
@@ -803,6 +813,7 @@
z0 = z1;
z1 = t;
}
+
realEigenvalues[index] = z0;
realEigenvalues[index + 1] = z1;
realEigenvalues[index + 2] = z2;
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=904561&r1=904560&r2=904561&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Fri Jan 29 17:08:28 2010
@@ -1,4 +1,4 @@
-<?xml version="1.0"?>
+<?xml version="1.0"?>
<!--
Licensed to the Apache Software Foundation (ASF) under one or more
contributor license agreements. See the NOTICE file distributed with
@@ -39,6 +39,10 @@
</properties>
<body>
<release version="2.1" date="TBD" description="TBD">
+ <action dev="luc" type="fix" >
+ Fixed a spurious exception in EigenDecompositionImpl when a 3x3 block
+ had two identical eigenvalues.
+ </action>
<action dev="luc" type="add" issue="MATH-334" >
Added min/max getters for real vectors (not yet in the RealVector interface for
compatibility purposes, but in the AbstractRealVector abstract class).
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=904561&r1=904560&r2=904561&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 Fri Jan 29 17:08:28 2010
@@ -64,6 +64,19 @@
assertEquals( 3125.0, ed.getRealEigenvalue(2), 3.0e-11);
}
+ public void testDimension3MultipleRoot() {
+ RealMatrix matrix =
+ MatrixUtils.createRealMatrix(new double[][] {
+ { 5, 10, 15 },
+ { 10, 20, 30 },
+ { 15, 30, 45 }
+ });
+ EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
+ assertEquals(70.0, ed.getRealEigenvalue(0), 3.0e-11);
+ assertEquals(0.0, ed.getRealEigenvalue(1), 3.0e-11);
+ assertEquals(0.0, ed.getRealEigenvalue(2), 3.0e-11);
+ }
+
public void testDimension4WithSplit() {
RealMatrix matrix =
MatrixUtils.createRealMatrix(new double[][] {