36 template<
class cmptType>
39 for (label j = 0; j < n_; ++j)
41 EValsRe_[j] = EVecs_[n_ - 1][j];
45 for (label i = n_ - 1; i > 0; --i)
48 cmptType scale =
Zero;
51 for (label
k = 0;
k < i; ++
k)
53 scale = scale +
mag(EValsRe_[
k]);
58 EValsIm_[i] = EValsRe_[i - 1];
60 for (label j = 0; j < i; ++j)
62 EValsRe_[j] = EVecs_[i - 1][j];
70 for (label
k = 0;
k < i; ++
k)
73 h += EValsRe_[
k]*EValsRe_[
k];
76 cmptType
f = EValsRe_[i - 1];
84 EValsIm_[i] = scale*
g;
86 EValsRe_[i - 1] =
f -
g;
88 for (label j = 0; j < i; ++j)
94 for (label j = 0; j < i; ++j)
98 g = EValsIm_[j] + EVecs_[j][j]*
f;
100 for (label
k = j + 1;
k <= i - 1; ++
k)
102 g += EVecs_[
k][j]*EValsRe_[
k];
103 EValsIm_[
k] += EVecs_[
k][j]*
f;
111 for (label j = 0; j < i; ++j)
114 f += EValsIm_[j]*EValsRe_[j];
117 cmptType hh =
f/(2.0*
h);
119 for (label j = 0; j < i; ++j)
121 EValsIm_[j] -= hh*EValsRe_[j];
124 for (label j = 0; j < i; ++j)
129 for (label
k = j;
k <= i - 1; ++
k)
132 (
f*EValsIm_[
k] +
g*EValsRe_[
k]);
135 EValsRe_[j] = EVecs_[i - 1][j];
144 for (label i = 0; i < n_ - 1; ++i)
146 EVecs_[n_ - 1][i] = EVecs_[i][i];
147 EVecs_[i][i] = cmptType(1);
148 cmptType
h = EValsRe_[i + 1];
152 for (label
k = 0;
k <= i; ++
k)
154 EValsRe_[
k] = EVecs_[
k][i + 1]/
h;
157 for (label j = 0; j <= i; ++j)
161 for (label
k = 0;
k <= i; ++
k)
163 g += EVecs_[
k][i + 1]*EVecs_[
k][j];
166 for (label
k = 0;
k <= i; ++
k)
168 EVecs_[
k][j] -=
g*EValsRe_[
k];
173 for (label
k = 0;
k <= i; ++
k)
175 EVecs_[
k][i + 1] =
Zero;
179 for (label j = 0; j < n_; ++j)
181 EValsRe_[j] = EVecs_[n_ - 1][j];
182 EVecs_[n_ - 1][j] =
Zero;
185 EVecs_[n_ - 1][n_ - 1] = cmptType(1);
190 template<
class cmptType>
193 for (label i = 1; i < n_; ++i)
195 EValsIm_[i - 1] = EValsIm_[i];
198 EValsIm_[n_-1] =
Zero;
201 cmptType tst1 =
Zero;
202 cmptType eps =
pow(2.0, -52.0);
204 for (label l = 0; l < n_; l++)
207 tst1 =
max(tst1,
mag(EValsRe_[l]) +
mag(EValsIm_[l]));
213 if (
mag(EValsIm_[m]) <= eps*tst1)
231 cmptType
g = EValsRe_[l];
232 cmptType
p = (EValsRe_[l + 1] -
g)/(2.0*EValsIm_[l]);
240 EValsRe_[l] = EValsIm_[l]/(
p + r);
241 EValsRe_[l + 1] = EValsIm_[l]*(
p + r);
242 cmptType dl1 = EValsRe_[l + 1];
243 cmptType
h =
g - EValsRe_[l];
245 for (label i = l + 2; i < n_; ++i)
254 cmptType
c = cmptType(1);
257 cmptType el1 = EValsIm_[l + 1];
261 for (label i = m - 1; i >= l; --i)
269 EValsIm_[i + 1] =
s*r;
272 p =
c*EValsRe_[i] -
s*
g;
273 EValsRe_[i + 1] =
h +
s*(
c*
g +
s*EValsRe_[i]);
276 for (label
k = 0;
k < n_; ++
k)
278 h = EVecs_[
k][i + 1];
279 EVecs_[
k][i + 1] =
s*EVecs_[
k][i] +
c*
h;
280 EVecs_[
k][i] =
c*EVecs_[
k][i] -
s*
h;
284 p = -
s*s2*c3*el1*EValsIm_[l]/dl1;
288 while (
mag(EValsIm_[l]) > eps*tst1);
291 EValsRe_[l] = EValsRe_[l] +
f;
296 for (label i = 0; i < n_ - 1; ++i)
299 cmptType
p = EValsRe_[i];
301 for (label j = i + 1; j < n_; ++j)
312 EValsRe_[
k] = EValsRe_[i];
315 for (label j = 0; j < n_; ++j)
318 EVecs_[j][i] = EVecs_[j][
k];
326 template<
class cmptType>
329 DiagonalMatrix<cmptType> orth_(n_);
334 for (label m = low + 1; m <= high - 1; ++m)
337 cmptType scale =
Zero;
339 for (label i = m; i <= high; ++i)
341 scale = scale +
mag(H_[i][m - 1]);
349 for (label i = high; i >= m; --i)
351 orth_[i] = H_[i][m - 1]/scale;
352 h += orth_[i]*orth_[i];
367 for (label j = m; j < n_; ++j)
371 for (label i = high; i >= m; --i)
373 f += orth_[i]*H_[i][j];
378 for (label i = m; i <= high; ++i)
380 H_[i][j] -=
f*orth_[i];
384 for (label i = 0; i <= high; ++i)
388 for (label j = high; j >= m; --j)
390 f += orth_[j]*H_[i][j];
395 for (label j = m; j <= high; ++j)
397 H_[i][j] -=
f*orth_[j];
401 orth_[m] = scale*orth_[m];
402 H_[m][m-1] = scale*
g;
407 for (label i = 0; i < n_; ++i)
409 for (label j = 0; j < n_; ++j)
411 EVecs_[i][j] = (i == j ? 1.0 : 0.0);
415 for (label m = high - 1; m >= low + 1; --m)
417 if (H_[m][m - 1] != 0.0)
419 for (label i = m + 1; i <= high; ++i)
421 orth_[i] = H_[i][m - 1];
424 for (label j = m; j <= high; ++j)
428 for (label i = m; i <= high; ++i)
430 g += orth_[i]*EVecs_[i][j];
434 g = (
g/orth_[m])/H_[m][m - 1];
436 for (label i = m; i <= high; ++i)
438 EVecs_[i][j] +=
g*orth_[i];
446 template<
class cmptType>
454 cmptType eps =
pow(2.0, -52.0);
455 cmptType exshift =
Zero;
456 cmptType
p = 0, q = 0, r = 0,
s = 0, z = 0, t, w,
x,
y;
459 cmptType norm =
Zero;
461 for (label i = 0; i < nn; ++i)
463 if ((i < low) || (i > high))
465 EValsRe_[i] = H_[i][i];
469 for (label j =
max(i - 1, 0); j < nn; ++j)
471 norm +=
mag(H_[i][j]);
484 s =
mag(H_[l - 1][l - 1]) +
mag(H_[l][l]);
491 if (
mag(H_[l][l - 1]) < eps*
s)
502 H_[
n][
n] = H_[
n][
n] + exshift;
503 EValsRe_[
n] = H_[
n][
n];
510 w = H_[
n][
n - 1]*H_[
n - 1][
n];
511 p = (H_[
n - 1][
n - 1] - H_[
n][
n])/2.0;
515 H_[
n - 1][
n - 1] += exshift;
530 EValsRe_[
n - 1] =
x + z;
531 EValsRe_[
n] = EValsRe_[
n - 1];
535 EValsRe_[
n] =
x - w/z;
538 EValsIm_[
n - 1] =
Zero;
549 for (label j =
n - 1; j < nn; ++j)
552 H_[
n - 1][j] = q*z +
p*H_[
n][j];
553 H_[
n][j] = q*H_[
n][j] -
p*z;
557 for (label i = 0; i <=
n; ++i)
560 H_[i][
n - 1] = q*z +
p*H_[i][
n];
561 H_[i][
n] = q*H_[i][
n] -
p*z;
565 for (label i = low; i <= high; ++i)
567 z = EVecs_[i][
n - 1];
568 EVecs_[i][
n - 1] = q*z +
p*EVecs_[i][
n];
569 EVecs_[i][
n] = q*EVecs_[i][
n] -
p*z;
574 EValsRe_[
n - 1] =
x +
p;
592 y = H_[
n - 1][
n - 1];
593 w = H_[
n][
n - 1]*H_[
n - 1][
n];
600 for (label i = low; i <=
n; ++i)
626 s =
x - w/((
y -
x)/2.0 +
s);
628 for (label i = low; i <=
n; ++i)
650 p = (r*
s - w)/H_[m + 1][m] + H_[m][m + 1];
651 q = H_[m + 1][m + 1] - z - r -
s;
652 r = H_[m + 2][m + 1];
666 < eps*(
mag(
p)*(
mag(H_[m - 1][m - 1])
667 +
mag(z) +
mag(H_[m + 1][m + 1])))
676 for (label i = m + 2; i <=
n; ++i)
687 for (label
k = m;
k <=
n - 1; ++
k)
689 label notlast = (
k !=
n - 1);
694 q = H_[
k + 1][
k - 1];
695 r = (notlast ? H_[
k + 2][
k - 1] : 0.0);
726 H_[
k][
k - 1] = -H_[
k][
k - 1];
737 for (label j =
k; j < nn; ++j)
739 p = H_[
k][j] + q*H_[
k + 1][j];
752 for (label i = 0; i <=
min(
n,
k + 3); ++i)
754 p =
x*H_[i][
k] +
y*H_[i][
k + 1];
767 for (label i = low; i <= high; ++i)
769 p =
x*EVecs_[i][
k] +
y*EVecs_[i][
k + 1];
773 p += z*EVecs_[i][
k + 2];
774 EVecs_[i][
k + 2] -=
p*r;
778 EVecs_[i][
k + 1] -=
p*q;
791 for (
n = nn-1;
n >= 0; --
n)
800 H_[
n][
n] = cmptType(1);
802 for (label i =
n-1; i >= 0; --i)
807 for (label j = l; j <=
n; ++j)
809 r += H_[i][j]*H_[j][
n];
812 if (EValsIm_[i] < 0.0)
821 if (EValsIm_[i] == 0.0)
829 H_[i][
n] = -r/(eps*norm);
836 q = (EValsRe_[i] -
p)*(EValsRe_[i] -
p)
837 + EValsIm_[i]*EValsIm_[i];
844 H_[i + 1][
n] = (-r - w*t)/
x;
848 H_[i + 1][
n] = (-
s -
y*t)/z;
857 for (label j = i; j <=
n; ++j)
872 H_[
n - 1][
n - 1] = q/H_[
n][
n - 1];
873 H_[
n - 1][
n] = -(H_[
n][
n] -
p)/H_[
n][
n - 1];
880 H_[
n - 1][
n - 1] = cDiv.Re();
881 H_[
n - 1][
n] = cDiv.Im();
885 H_[
n][
n] = cmptType(1);
887 for (label i =
n - 2; i >= 0; --i)
889 cmptType ra, sa, vr, vi;
893 for (label j = l; j <=
n; ++j)
895 ra += H_[i][j]*H_[j][
n-1];
896 sa += H_[i][j]*H_[j][
n];
901 if (EValsIm_[i] < 0.0)
911 if (EValsIm_[i] == 0)
914 H_[i][
n - 1] = cDiv.Re();
915 H_[i][
n] = cDiv.Im();
922 vr = (EValsRe_[i] -
p)*(EValsRe_[i] -
p)
923 + EValsIm_[i]*EValsIm_[i] - q*q;
925 vi = 2.0*(EValsRe_[i] -
p)*q;
927 if ((vr == 0.0) && (vi == 0.0))
937 H_[i][
n - 1] = cDiv.Re();
938 H_[i][
n] = cDiv.Im();
942 H_[i + 1][
n - 1] = (-ra - w*H_[i][
n - 1]
945 H_[i + 1][
n] = (-sa - w*H_[i][
n]
953 H_[i + 1][
n - 1] = cDiv.Re();
954 H_[i + 1][
n] = cDiv.Im();
963 for (label j = i; j <=
n; ++j)
975 for (label i = 0; i < nn; ++i)
977 if (i < low || i > high)
979 for (label j = i; j < nn; ++j)
981 EVecs_[i][j] = H_[i][j];
987 for (label j = nn - 1; j >= low; --j)
989 for (label i = low; i <= high; ++i)
993 for (label
k = low;
k <=
min(j, high); ++
k)
995 z = z + EVecs_[i][
k]*H_[
k][j];
1006 template<
class cmptType>
1018 <<
"Input matrix has zero size."
1025 tridiagonaliseSymmMatrix();
1026 symmTridiagonalQL();
1037 template<
class cmptType>
1053 <<
"Input matrix has zero size."
1060 tridiagonaliseSymmMatrix();
1061 symmTridiagonalQL();
1074 template<
class cmptType>
1085 [&](
const cmptType&
x){
return complex(
x); }
1088 for (label i = 0; i < EValsIm_.size(); ++i)
1090 if (
mag(EValsIm_[i]) > VSMALL)
1092 for (label j = 0; j < EVecs.m(); ++j)
1094 EVecs(j, i) =
complex(EVecs_(j, i), EVecs_(j, i+1));
1095 EVecs(j, i+1) = EVecs(j, i).conjugate();