30 #ifndef CDPL_MATH_JACOBIDIAGONALIZATION_HPP
31 #define CDPL_MATH_JACOBIDIAGONALIZATION_HPP
66 template <
typename M1,
typename V,
typename M2>
67 bool jacobiDiagonalize(MatrixExpression<M1>& a, VectorExpression<V>& d, MatrixExpression<M2>& v, std::size_t max_iter = 50);
78 void jacobiRotate(
M& a,
typename M::SizeType i,
typename M::SizeType j,
79 typename M::SizeType k,
typename M::SizeType l,
80 typename M::ValueType tau,
typename M::ValueType
s)
82 typedef typename M::ValueType ValueType;
84 ValueType g = a(i, j);
85 ValueType h = a(k, l);
87 a(i, j) = g -
s * (h + g * tau);
88 a(k, l) = h +
s * (g - h * tau);
93 template <
typename M1,
typename V,
typename M2>
97 typename M2::ValueType>::Type ValueType;
100 typename M2::SizeType>::Type SizeType;
102 SizeType n = a().getSize1();
111 for (SizeType ip = 0; ip < n; ip++) {
112 b(ip) = d()(ip) = a()(ip, ip);
116 for (std::size_t i = 0; i < max_iter; i++) {
117 ValueType sm = ValueType();
119 for (SizeType ip = 0; ip < n - 1; ip++)
120 for (SizeType iq = ip + 1; iq < n; iq++)
123 if (sm == ValueType())
129 tresh = ValueType(0.2) * sm / (n * n);
133 for (SizeType ip = 0; ip < n - 1; ip++) {
134 for (SizeType iq = ip + 1; iq < n; iq++) {
139 a()(ip, iq) = ValueType();
142 ValueType h = d()(iq) - d()(ip);
146 t = (a()(ip, iq)) / h;
149 ValueType theta = 0.5 * h / a()(ip, iq);
152 if (theta < ValueType())
158 ValueType tau =
s / (1 + c);
164 a()(ip, iq) = ValueType();
166 for (SizeType j = 0; j < ip; j++)
167 jacobiRotate(a(), j, ip, j, iq, tau,
s);
169 for (SizeType j = ip + 1; j < iq; j++)
170 jacobiRotate(a(), ip, j, j, iq, tau,
s);
172 for (SizeType j = iq + 1; j < n; j++)
173 jacobiRotate(a(), ip, j, iq, j, tau,
s);
175 for (SizeType j = 0; j < n; j++)
176 jacobiRotate(v(), j, ip, j, iq, tau,
s);
181 for (SizeType ip = 0; ip < n; ip++) {
191 #endif // CDPL_MATH_JACOBIDIAGONALIZATION_HPP