16 static const double EPS = 0.00001;
22 double a_norm, a_normEPS, thr, thr_nn;
25 int i, j, k, ij, ik,
l, m, lm, mq, lq, ll, mm, imv, im, iq, ilv, il, nn;
27 double a_ij, a_lm, a_ll, a_mm, a_im, a_il;
31 double sinx, sinx_2, cosx, cosx_2, sincos;
35 nn = (n * (n + 1)) / 2;
40 for (ij = 0; ij < nn; ij++) {
48 v =
new double[n * n];
51 for (i = 0; i < n; i++) {
52 for (j = 0; j < n; j++) {
68 for (i = 1; i <= n; i++) {
69 for (j = 1; j <= i; j++) {
72 a_norm += a_ij * a_ij;
79 a_normEPS = a_norm *
EPS;
83 while (thr > a_normEPS && nb_iter <
MAX_ITER) {
87 for (
l = 1;
l < n;
l++) {
88 for (m =
l + 1; m <= n; m++) {
97 if (a_lm_2 < thr_nn) {
112 x = -
atan((a_lm + a_lm) / delta) / 2.0;
117 sinx_2 = sinx * sinx;
118 cosx_2 = cosx * cosx;
119 sincos = sinx * cosx;
125 for (i = 1; i <= n; i++) {
126 if (!
ELEM(i,
l, m)) {
127 iq = (i * i - i) / 2;
145 a[il] = a_il * cosx - a_im * sinx;
146 a[im] = a_il * sinx + a_im * cosx;
155 v[ilv] = cosx * v_ilv - sinx * v_imv;
156 v[imv] = sinx * v_ilv + cosx * v_imv;
162 a[ll] = a_ll * cosx_2 + a_mm * sinx_2 -
x;
163 a[mm] = a_ll * sinx_2 + a_mm * cosx_2 +
x;
166 thr =
fabs(thr - a_lm_2);
177 for (i = 0; i < n; i++) {
178 k = i + (i * (i + 1)) / 2;
187 for (i = 0; i < n; i++) {
191 for (i = 0; i < (n - 1); i++) {
195 for (j = i + 1; j < n; j++) {
196 if (
x < eigen_val[j]) {
202 eigen_val[k] = eigen_val[i];
216 for (k = 0; k < n; k++) {
218 for (i = 0; i < n; i++) {
219 eigen_vec[ij++] =
v[ik++];
ATTR_WARN_UNUSED_RESULT const BMLoop * l
ATTR_WARN_UNUSED_RESULT const BMVert * v
ccl_device_inline float2 fabs(const float2 &a)
void semi_definite_symmetric_eigen(const double *mat, int n, double *eigen_vec, double *eigen_val)
INLINE Rall1d< T, V, S > cos(const Rall1d< T, V, S > &arg)
INLINE Rall1d< T, V, S > atan(const Rall1d< T, V, S > &x)
INLINE Rall1d< T, V, S > sin(const Rall1d< T, V, S > &arg)