41 #include <Eigen/Eigenvalues> 42 #include <Eigen/Geometry> 45 using namespace Eigen;
48 template<
typename Matrix,
typename Roots>
52 const Scalar s_inv3 = 1.0/3.0;
58 Scalar c0 =
m(0,0)*
m(1,1)*
m(2,2) +
Scalar(2)*
m(0,1)*
m(0,2)*
m(1,2) -
m(0,0)*
m(1,2)*
m(1,2) -
m(1,1)*
m(0,2)*
m(0,2) -
m(2,2)*
m(0,1)*
m(0,1);
59 Scalar c1 =
m(0,0)*
m(1,1) -
m(0,1)*
m(0,1) +
m(0,0)*
m(2,2) -
m(0,2)*
m(0,2) +
m(1,1)*
m(2,2) -
m(1,2)*
m(1,2);
60 Scalar c2 =
m(0,0) +
m(1,1) +
m(2,2);
64 Scalar c2_over_3 = c2*s_inv3;
65 Scalar a_over_3 = (c1 - c2*c2_over_3)*s_inv3;
69 Scalar half_b =
Scalar(0.5)*(c0 + c2_over_3*(
Scalar(2)*c2_over_3*c2_over_3 - c1));
71 Scalar
q = half_b*half_b + a_over_3*a_over_3*a_over_3;
80 roots(2) = c2_over_3 +
Scalar(2)*rho*cos_theta;
81 roots(0) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta);
82 roots(1) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta);
85 template<
typename Matrix,
typename Vector>
92 Scalar shift = mat.trace()/3;
94 scaledMat.diagonal().array() -= shift;
95 Scalar
scale = scaledMat.cwiseAbs().maxCoeff();
143 tmp.diagonal ().array () -= evals (2);
144 evecs.col (2) = tmp.row (0).cross (tmp.row (1)).normalized ();
147 tmp.diagonal ().array () -= evals (1);
148 evecs.col(1) = tmp.row (0).cross(tmp.row (1));
149 Scalar
n1 = evecs.col(1).norm();
151 evecs.col(1) = evecs.col(2).unitOrthogonal();
156 evecs.col(1) = evecs.col(2).cross(evecs.col(1).cross(evecs.col(2))).normalized();
157 evecs.col(0) = evecs.col(2).cross(evecs.col(1));
162 evals.array()+=shift;
170 typedef Matrix3d
Mat;
171 typedef Vector3d
Vec;
172 Mat
A = Mat::Random(3,3);
179 std::cout <<
"Eigen iterative: " << t.
best() <<
"s\n";
182 std::cout <<
"Eigen direct : " << t.
best() <<
"s\n";
187 std::cout <<
"Direct: " << t.
best() <<
"s\n\n";
EIGEN_DEVICE_FUNC SelfAdjointEigenSolver & computeDirect(const MatrixType &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix using a closed-form algorithm.
internal::traits< Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > >::Scalar Scalar
Jet< T, N > cos(const Jet< T, N > &f)
EIGEN_DEVICE_FUNC SelfAdjointEigenSolver & compute(const EigenBase< InputType > &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix.
Computes eigenvalues and eigenvectors of selfadjoint matrices.
Jet< T, N > sin(const Jet< T, N > &f)
void eigen33(const Matrix &mat, Matrix &evecs, Vector &evals)
Namespace containing all symbols from the Eigen library.
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
double best(int TIMER=CPU_TIMER) const
SelfAdjointEigenSolver< PlainMatrixType > eig(mat, computeVectors?ComputeEigenvectors:EigenvaluesOnly)
EIGEN_DEVICE_FUNC const Scalar & q
AnnoyingScalar atan2(const AnnoyingScalar &y, const AnnoyingScalar &x)
Matrix< Scalar, Dynamic, 1 > Vec
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy y set format x g set format y g set format x2 g set format y2 g set format z g set angles radians set nogrid set key title set key left top Right noreverse box linetype linewidth samplen spacing width set nolabel set noarrow set nologscale set logscale x set set pointsize set encoding default set nopolar set noparametric set set set set surface set nocontour set clabel set mapping cartesian set nohidden3d set cntrparam order set cntrparam linear set cntrparam levels auto set cntrparam points set size set set xzeroaxis lt lw set x2zeroaxis lt lw set yzeroaxis lt lw set y2zeroaxis lt lw set tics in set ticslevel set tics scale
Matrix< Scalar, Dynamic, Dynamic > Mat
void computeRoots(const Matrix &m, Roots &roots)
Jet< T, N > sqrt(const Jet< T, N > &f)
#define BENCH(TIMER, TRIES, REP, CODE)
The matrix class, also used for vectors and row-vectors.
Eigen::Matrix< double, Eigen::Dynamic, 1 > Vector