an even more stable procedure

1 files changed, 24 insertions, 13 deletions

diff --git a/bench/eig33.cpp b/bench/eig33.cpp
index 4a2320ee1..9f9e1f8a9 100644
--- a/bench/eig33.cpp
+++ b/bench/eig33.cpp
@@ -153,20 +153,31 @@ void eigen33(const Matrix& mat, Matrix& evecs, Vector& evals)
 //   tmp.diagonal().array() -= evals(2);
 //   evecs.col(2) = tmp.row(0).cross(tmp.row(1)).normalized();
   
-  // a slighlty more stable version:
-  Matrix tmp;
-  tmp = scaledMat;
-  tmp.diagonal ().array () -= evals (2);
-  evecs.col (2) = tmp.row (0).cross (tmp.row (1)).normalized ();
-  tmp = scaledMat;
-  tmp.diagonal ().array () -= evals (1);
-  evecs.col(1) = tmp.row (0).cross(tmp.row (1));
-  Scalar n1 = evecs.col(1).norm();
-  if(n1<=Eigen::NumTraits<Scalar>::epsilon())
-    evecs.col(1) = evecs.col(2).unitOrthogonal();
+  // a more stable version:
+  if((evals(2)-evals(0))<=Eigen::NumTraits<Scalar>::epsilon())
+  {
+    evecs.setIdentity();
+  }
   else
-    evecs.col(1) /= n1;
-  evecs.col(0) = evecs.col(2).cross(evecs.col(1));
+  {
+    Matrix tmp;
+    tmp = scaledMat;
+    tmp.diagonal ().array () -= evals (2);
+    evecs.col (2) = tmp.row (0).cross (tmp.row (1)).normalized ();
+    
+    tmp = scaledMat;
+    tmp.diagonal ().array () -= evals (1);
+    evecs.col(1) = tmp.row (0).cross(tmp.row (1));
+    Scalar n1 = evecs.col(1).norm();
+    if(n1<=Eigen::NumTraits<Scalar>::epsilon())
+      evecs.col(1) = evecs.col(2).unitOrthogonal();
+    else
+      evecs.col(1) /= n1;
+    
+    // make sure that evecs[1] is orthogonal to evecs[2]
+    evecs.col(1) = evecs.col(2).cross(evecs.col(1).cross(evecs.col(2))).normalized();
+    evecs.col(0) = evecs.col(2).cross(evecs.col(1));
+  }
   
   // Rescale back to the original size.
   evals *= scale;