From 1eae6d0fb9f433691d670cafa4f37ca43307a383 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 3 Feb 2011 11:25:34 +0100 Subject: an even more stable procedure --- bench/eig33.cpp | 37 ++++++++++++++++++++++++------------- 1 file changed, 24 insertions(+), 13 deletions(-) (limited to 'bench/eig33.cpp') 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::epsilon()) - evecs.col(1) = evecs.col(2).unitOrthogonal(); + // a more stable version: + if((evals(2)-evals(0))<=Eigen::NumTraits::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::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; -- cgit v1.2.3