From 5beb2f4f0d7b245853eb96a64f2367d91db056b9 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 3 Feb 2011 10:31:45 +0100 Subject: slightly more stable eigen vector computation --- bench/eig33.cpp | 35 +++++++++++++++++++++++++---------- 1 file changed, 25 insertions(+), 10 deletions(-) (limited to 'bench/eig33.cpp') diff --git a/bench/eig33.cpp b/bench/eig33.cpp index da8518012..4a2320ee1 100644 --- a/bench/eig33.cpp +++ b/bench/eig33.cpp @@ -139,19 +139,34 @@ void eigen33(const Matrix& mat, Matrix& evecs, Vector& evals) // } // evecs.col(2) = evecs.col(0).cross(evecs.col(1)).normalized(); - // naive version +// // naive version +// Matrix tmp; +// tmp = scaledMat; +// tmp.diagonal().array() -= evals(0); +// evecs.col(0) = 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)).normalized(); +// +// tmp = scaledMat; +// 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(0); - evecs.col(0) = 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)).normalized(); - + tmp.diagonal ().array () -= evals (2); + evecs.col (2) = tmp.row (0).cross (tmp.row (1)).normalized (); tmp = scaledMat; - tmp.diagonal().array() -= evals(2); - evecs.col(2) = tmp.row(0).cross(tmp.row(1)).normalized(); + 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; + evecs.col(0) = evecs.col(2).cross(evecs.col(1)); // Rescale back to the original size. evals *= scale; -- cgit v1.2.3