GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-30 18:25:11 +08:00
Code Issues Packages Projects Releases Wiki Activity

simplify a bit the 2x2 direct eigenvalue solver

Browse Source
This commit is contained in:
Gael Guennebaud 2011-07-22 11:21:43 +02:00
parent 47a2bca89f
commit e8313364c1
1 changed files with 26 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

70
Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
View File

@ -583,19 +583,19 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3
scaledMat = scaledMat.template selfadjointView<Lower>(); scaledMat = scaledMat.template selfadjointView<Lower>();
MatrixType tmp; MatrixType tmp;
tmp = scaledMat; tmp = scaledMat;
tmp.diagonal().array () -= eivals (2); tmp.diagonal().array () -= eivals(2);
eivecs.col (2) = tmp.row (0).cross (tmp.row (1)).normalized (); eivecs.col(2) = tmp.row(0).cross(tmp.row(1)).normalized();
tmp = scaledMat; tmp = scaledMat;
tmp.diagonal ().array () -= eivals (1); tmp.diagonal().array() -= eivals(1);
eivecs.col(1) = tmp.row (0).cross(tmp.row (1)); eivecs.col(1) = tmp.row(0).cross(tmp.row(1));
Scalar n1 = eivecs.col(1).norm(); Scalar n1 = eivecs.col(1).norm();
if(n1<=Eigen::NumTraits<Scalar>::epsilon()) if(n1<=Eigen::NumTraits<Scalar>::epsilon())
eivecs.col(1) = eivecs.col(2).unitOrthogonal(); eivecs.col(1) = eivecs.col(2).unitOrthogonal();
else else
eivecs.col(1) /= n1; eivecs.col(1) /= n1;
// make sure that evecs[1] is orthogonal to evecs[2] // make sure that eivecs[1] is orthogonal to eivecs[2]
eivecs.col(1) = eivecs.col(2).cross(eivecs.col(1).cross(eivecs.col(2))).normalized(); eivecs.col(1) = eivecs.col(2).cross(eivecs.col(1).cross(eivecs.col(2))).normalized();
eivecs.col(0) = eivecs.col(2).cross(eivecs.col(1)); eivecs.col(0) = eivecs.col(2).cross(eivecs.col(1));
} }
@ -617,6 +617,15 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2
typedef typename SolverType::RealVectorType VectorType; typedef typename SolverType::RealVectorType VectorType;
typedef typename SolverType::Scalar Scalar; typedef typename SolverType::Scalar Scalar;
inline static void computeRoots(const MatrixType& m, VectorType& roots)
{
using std::sqrt;
const Scalar t0 = Scalar(0.5) * sqrt( abs2(m(0,0)-m(1,1)) + Scalar(4)*m(1,0)*m(1,0));
const Scalar t1 = Scalar(0.5) * (m(0,0) + m(1,1));
roots(0) = t1 - t0;
roots(1) = t1 + t0;
}
inline static void run(SolverType& solver, const MatrixType& mat, int options) inline static void run(SolverType& solver, const MatrixType& mat, int options)
{ {
eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows()); eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows());
@ -634,56 +643,27 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2
MatrixType scaledMat = mat / scale; MatrixType scaledMat = mat / scale;
// Compute the eigenvalues // Compute the eigenvalues
const Scalar xx = scaledMat(0,0); computeRoots(scaledMat,eivals);
const Scalar yy = scaledMat(1,1);
const Scalar xy = scaledMat(1,0);
const Scalar root = std::sqrt((xx-yy)*(xx-yy)+Scalar(4.0)*xy*xy);
const Scalar eig_val1 = Scalar(0.5)*(xx+yy+root);
const Scalar eig_val2 = Scalar(0.5)*(xx+yy-root);
using std::sqrt;
using std::abs;
// compute the eigen vectors // compute the eigen vectors
if(computeEigenvectors) if(computeEigenvectors)
{ {
Scalar nx, ny; scaledMat.diagonal().array () -= eivals(1);
if (xx-yy > 0) Scalar a2 = abs2(scaledMat(0,0));
Scalar c2 = abs2(scaledMat(1,1));
Scalar b2 = abs2(scaledMat(1,0));
if(a2>c2)
{ {
Scalar len(1.0); eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0);
if (-yy+xx+root != Scalar(0.0)) eivecs.col(1) /= sqrt(a2+b2);
{
ny = -Scalar(2.0)*xy/(-yy+xx+root);
len = Scalar(1.0)/sqrt(1+ny*ny);
ny *= len;
}
else
ny = Scalar(0.0);
nx = -len;
} }
else else
{ {
Scalar len = Scalar(1.0); eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0);
if (-yy+xx-root != Scalar(0.0)) eivecs.col(1) /= sqrt(c2+b2);
{
nx = -Scalar(2.0)*xy/(-yy+xx-root);
len = Scalar(1.0)/sqrt(1+nx*nx);
nx *= len;
}
else
nx = Scalar(0.0);
ny = len;
} }
int id0(0), id1(1); eivecs.col(0) << eivecs.col(1).unitOrthogonal();
if(eig_val1>eig_val2)
std::swap(id0,id1);
eivecs.col(id0) << nx, ny;
eivecs.col(id1) << -ny, nx;
eivals(id0) = eig_val1;
eivals(id1) = eig_val2;
} }
// Rescale back to the original size. // Rescale back to the original size.