42 static const double root3=1.73205080756887729352744634150587236694280525381038062805580;
63 template<
int Rows,
int Cols,
typename P,
typename B>
71 int N = evalues.size();
72 int lda = evalues.size();
78 syev_((
char*)
"V",(
char*)
"U",&N,&evectors[0][0],&lda,&evalues[0], &size,&lwork,&info);
83 syev_((
char*)
"V",(
char*)
"U",&N,&evectors[0][0],&lda,&evalues[0], &WORK[0],&lwork,&info);
86 std::cerr <<
"In SymEigen<"<<Size<<
">: " << info
87 <<
" off-diagonal elements of an intermediate tridiagonal form did not converge to zero." << std::endl
88 <<
"M = " << m << std::endl;
104 template<
typename P,
typename B>
106 double trace = m[0][0] + m[1][1];
107 double det = m[0][0]*m[1][1] - m[0][1]*m[1][0];
108 double disc = trace*trace - 4 * det;
111 double root_disc = sqrt(disc);
112 ev[0] = 0.5 * (trace - root_disc);
113 ev[1] = 0.5 * (trace + root_disc);
114 double a = m[0][0] - ev[0];
116 double magsq = a*a + b*b;
123 eig[0] *= 1.0/sqrt(magsq);
125 eig[1][0] = -eig[0][1];
126 eig[1][1] = eig[0][0];
141 template<
typename P,
typename B>
148 const double& a11 = m[0][0];
149 const double& a12 = m[0][1];
150 const double& a13 = m[0][2];
152 const double& a22 = m[1][1];
153 const double& a23 = m[1][2];
155 const double& a33 = m[2][2];
158 double a = -a11-a22-a33;
159 double b = a11*a22+a11*a33+a22*a33-a12*a12-a13*a13-a23*a23;
160 double c = a11*(a23*a23)+(a13*a13)*a22+(a12*a12)*a33-a12*a13*a23*2.0-a11*a22*a33;
163 double p = b - a*a/3;
164 double q = c + (2*a*a*a - 9*a*b)/27;
167 double beta_descriminant = q*q/4 + p*p*p/27;
170 double beta = sqrt(-beta_descriminant);
171 double r2 = alpha*alpha - beta_descriminant;
177 double cuberoot_r = pow(r2, 1./6);
178 double theta3 = atan2(beta, alpha)/3;
180 double A_plus_B = 2*cuberoot_r*cos(theta3);
181 double A_minus_B = -2*cuberoot_r*sin(theta3);
184 ev =
makeVector(A_plus_B, -A_plus_B/2 + A_minus_B * sqrt(3)/2, -A_plus_B/2 - A_minus_B * sqrt(3)/2) -
Ones * a/3;
194 eig[0][0]=a12 * a23 - a13 * (a22 - ev[0]);
195 eig[0][1]=a12 * a13 - a23 * (a11 - ev[0]);
196 eig[0][2]=(a11-ev[0])*(a22-ev[0]) - a12*a12;
198 eig[1][0]=a12 * a23 - a13 * (a22 - ev[1]);
199 eig[1][1]=a12 * a13 - a23 * (a11 - ev[1]);
200 eig[1][2]=(a11-ev[1])*(a22-ev[1]) - a12*a12;
202 eig[2][0]=a12 * a23 - a13 * (a22 - ev[2]);
203 eig[2][1]=a12 * a13 - a23 * (a11 - ev[2]);
204 eig[2][2]=(a11-ev[2])*(a22-ev[2]) - a12*a12;
259 template <
int Size=Dynamic,
typename Precision =
double>
272 template<
int R,
int C,
typename B>
278 template<
int R,
int C,
typename B>
289 template <
int S,
typename P,
typename B>
298 template <
int R,
int C,
typename P,
typename B>
322 if(fabs(
my_evalues[i]) * condition > max_diag) {
403 for (
int i = 0; i < my_evalues.size(); ++i) {
404 diag_isqrt[i] = std::sqrt(my_evalues[i]);
405 if(fabs(diag_isqrt[i]) * condition > max_diag) {
406 diag_isqrt[i] = 1/diag_isqrt[i];
Matrix< Size, Size, Precision > get_pinv(const double condition=Internal::symeigen_condition_no) const
SymEigen(const Matrix< R, C, Precision, B > &m)
void normalize(Vector< Size, Precision, Base > v)
Matrix< Size, Size, Precision > get_isqrtm(const double condition=Internal::symeigen_condition_no) const
void compute(const Matrix< R, C, Precision, B > &m)
Perform the eigen decomposition of a matrix.
Vector< Size, Precision > my_evalues
Matrix< Size, Size, Precision > get_sqrtm() const
Everything lives inside this namespace.
static void test(int s1, int s2)
static void compute(const Matrix< Rows, Cols, P, B > &m, Matrix< Size, Size, P > &evectors, Vector< Size, P > &evalues)
Precision get_determinant() const
Get the determinant of the matrix.
static const Operator< Internal::Scalars< Internal::One > > Ones
static void compute(const Matrix< 2, 2, P, B > &m, Matrix< 2, 2, P > &eig, Vector< 2, P > &ev)
Matrix< Size, Size, Precision > & get_evectors()
const Matrix< Size, Size, Precision > & get_evectors() const
Matrix< Size, C, Precision > backsub(const Matrix< R, C, P, B > &rhs) const
Precision trace(const Matrix< Rows, Cols, Precision, Base > &m)
void syev_(const char *JOBZ, const char *UPLO, int *N, double *A, int *lda, double *W, double *WORK, int *LWORK, int *INFO)
bool is_posdef() const
Is the matrix positive definite?
Matrix< Size, Size, Precision > my_evectors
bool is_negdef() const
Is the matrix negative definite?
Vector< 1 > makeVector(double x1)
const Vector< Size, Precision > & get_evalues() const
Vector< Size, Precision > backsub(const Vector< S, P, B > &rhs) const
static const double symeigen_condition_no
Default condition number for SymEigen::backsub, SymEigen::get_pinv and SymEigen::get_inv_diag.
static const double root3
Vector< Size, Precision > get_inv_diag(const double condition) const
Vector< Internal::Sizer< S1, S2 >::size, typename Internal::MultiplyType< P1, P2 >::type > diagmult(const Vector< S1, P1, B1 > &v1, const Vector< S2, P2, B2 > &v2)
static void compute(const Matrix< 3, 3, P, B > &m, Matrix< 3, 3, P > &eig, Vector< 3, P > &ev)
Vector< Size, Precision > & get_evalues()