commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject 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
Date Fri, 29 Jan 2010 17:08:28 GMT
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[][] {



Mime
View raw message