diff options
Diffstat (limited to 'bench')
-rw-r--r-- | bench/eig33.cpp | 37 |
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; |