38 void cdiv(
double ar,
double ai,
double br,
double bi,
39 double &cr,
double &ci)
49 s = brs * brs + bis * bis;
51 cr = (ars * brs + ais * bis) / s;
52 ci = (ais * brs - ars * bis) / s;
66 lambda(M.lambda), vectors(M.vectors)
71 lambda(M.nrows()), vectors(M.nrows(), M.nrows(), 1.0)
82 throw SizeError(
"DoubleEigen::DoubleEigen()",
"Matrix is not square.");
95 elmhes(h, low, upp, index);
100 if(
hqr2(h, low, upp) != 0) {
102 "Unable to find all eigenvalues.");
107 if(
hqr(h, low, upp) != 0) {
109 "Unable to find all eigenvalues.");
132 for(
int i = 0; i <
n; i++) {
135 for(
int j = 0; j <
n; j++) {
136 R[j][i] = complex<double>(
vectors[j][i]);
140 for(
int j = 0; j <
n; j++) {
141 R[j][i] = complex<double>(
vectors[j][i], + vectors[j][i+1]);
142 R[j][i+1] = complex<double>(vectors[j][i], - vectors[j][i+1]);
203 for(
int j = upp; ; j--) {
204 for(
int i = 0; i <= upp; i++) {
205 if(i != j && copy[j][i] != 0.0)
goto next_column;
209 exchange(copy, j, upp, low, upp);
210 scale[upp] = double(j);
219 for(
int j = low; ; j++) {
220 for(
int i = low; i <= upp; i++) {
221 if(i != j && copy[i][j] != 0.0)
goto next_row;
225 exchange(copy, j, low, low, upp);
226 scale[low] = double(j);
235 for(
int i = low; i <= upp; i++) scale[i] = 1.0;
238 const double radix = 16.0;
239 const double b2 = radix * radix;
245 for(
int i = low; i <= upp; i++) {
249 for(
int j = low; j <= upp; j++) {
257 if(c != 0.0 && r != 0.0) {
258 double g = r / radix;
275 if((c + r) / f < s * .95) {
279 for(
int j = low; j <
n; j++) copy[i][j] *= g;
280 for(
int j = 0; j <= upp; j++) copy[j][i] *= f;
316 for(
int i = low; i <= upp; i++) {
318 for(
int j = 0; j <
n; j++)
vectors[i][j] *= s;
323 for(
int i = low; i-- > 0;) {
324 int k = int(scale[i]);
328 for(
int i = upp + 1; i <
n; i++) {
329 int k = int(scale[i]);
365 for(
int m = low + 1; m < upp; m++) {
369 for(
int j = m; j <= upp; j++) {
380 for(
int j = m - 1; j <
n; j++) swap(copy[i][j], copy[m][j]);
381 for(
int j = 0; j <= upp; j++) swap(copy[j][i], copy[j][m]);
385 for(i = m + 1; i <= upp; i++) {
386 double y = copy[i][m-1];
390 for(
int j = m; j <
n; j++) copy[i][j] -= y * copy[m][j];
391 for(
int j = 0; j <= upp; j++) copy[j][m] += y * copy[j][i];
430 for(
int mp = upp; --mp > low;) {
432 for(i = mp + 1; i <= upp; i++)
vectors[i][mp] = copy[i][mp-1];
436 for(
int j = mp; j <= upp; j++) {
448 int j,
int m,
int low,
int upp) {
452 for(
int i = 0; i <= upp; i++) swap(copy[i][j], copy[i][m]);
453 for(
int i = low; i <
n; i++) swap(copy[j][i], copy[m][i]);
498 for(
int i = 0; i < low; i++) lambda[i] = complex<double>(h[i][i]);
499 for(
int i = upp + 1; i < n; i++) lambda[i] = complex<double>(h[i][i]);
504 for(
int i = 0; i <
n; i++) {
505 for(
int j = k; j <
n; j++) norm +=
std::abs(h[i][j]);
511 double p = 0.0, q = 0.0, r = 0.0, s = 0.0;
512 double w = 0.0, x = 0.0, y = 0.0, z = 0.0;
514 for(
int en = upp + 1; en-- > low;) {
521 for(l = en; l > low; l--) {
523 if(s == 0.0) s = norm;
524 if((s +
std::abs(h[l][l-1])) == s)
break;
529 if(l == en)
goto one_root;
532 w = h[en][en-1] * h[en-1][en];
533 if(l == en - 1)
goto two_roots;
536 if(itn == 0)
return (en + 1);
539 if(its == 10 || its == 20) {
541 for(
int i = low; i <= en; i++) h[i][i] -= x;
552 for(m = en - 1; m-- > l;) {
556 p = (r * s - w) / h[m+1][m] + h[m][m+1];
557 q = h[m+1][m+1] - z - r - s;
566 if(tst2 == tst1)
break;
570 for(
int i = m + 3; i <= en; i++) h[i][i-2] = h[i][i-3] = 0.0;
573 for(k = m; k < en; k++) {
574 bool notLast = (k != en - 1);
579 r = notLast ? h[k+2][k-1] : 0.0;
581 if(x == 0.0)
continue;
587 s =
sqrt(p * p + q * q + r * r);
593 h[k][k-1] = - h[k][k-1];
603 for(
int j = k; j <= en; j++) {
604 p = h[k][j] + q * h[k+1][j];
614 int j = (en < k + 3) ? en : (k + 3);
615 for(
int i = l; i <= j; i++) {
616 p = x * h[i][k] + y * h[i][k+1];
628 lambda[en] = complex<double>(x + t);
639 z = (p > 0.0) ? (p + z) : (p - z);
640 lambda[en-1] = complex<double>(x + z);
641 lambda[en] = complex<double>((z != 0.0) ? (x - w / z) : x + z);
644 lambda[en-1] = complex<double>(x + p, + z);
645 lambda[en] = complex<double>(x + p, - z);
710 for(
int i = 0; i < low; i++) lambda[i] = complex<double>(h[i][i]);
711 for(
int i = upp + 1; i < n; i++) lambda[i] = complex<double>(h[i][i]);
716 for(
int i = 0; i <
n; i++) {
717 for(
int j = k; j <
n; j++) norm +=
std::abs(h[i][j]);
722 double p = 0.0, q = 0.0, r = 0.0, s = 0.0;
723 double w = 0.0, x = 0.0, y = 0.0, z = 0.0;
727 for(
int en = upp + 1; en-- > low;) {
733 for(l = en; l > low; l--) {
735 if(s == 0.0) s = norm;
736 if((s +
std::abs(h[l][l-1])) == s)
break;
741 if(l == en)
goto one_root;
744 w = h[en][en-1] * h[en-1][en];
745 if(l == en - 1)
goto two_roots;
748 if(itn == 0)
return (en + 1);
751 if(its == 10 || its == 20) {
753 for(
int i = low; i <= en; i++) h[i][i] -= x;
764 for(m = en - 1; m-- > l;) {
768 p = (r * s - w) / h[m+1][m] + h[m][m+1];
769 q = h[m+1][m+1] - z - r - s;
778 if(tst2 == tst1)
break;
782 for(
int i = m + 3; i <= en; i++) h[i][i-2] = h[i][i-3] = 0.0;
785 for(k = m; k < en; k++) {
786 bool notLast = (k != en - 1);
791 r = notLast ? h[k+2][k-1] : 0.0;
793 if(x == 0.0)
continue;
799 s =
sqrt(p * p + q * q + r * r);
805 h[k][k-1] = - h[k][k-1];
815 for(
int j = k; j <
n; j++) {
816 p = h[k][j] + q * h[k+1][j];
826 int j = (en < k + 3) ? en : (k + 3);
827 for(
int i = 0; i <= j; i++) {
828 p = x * h[i][k] + y * h[i][k+1];
838 for(
int i = low; i <= upp; i++) {
841 p += z * vectors[i][k+2];
842 vectors[i][k+2] -= p * r;
845 vectors[i][k+1] -= p * q;
851 lambda[en] = complex<double>(h[en][en] = x + t);
858 h[en][en] = x = x + t;
859 h[en-1][en-1] = y + t;
863 z = (p > 0.0) ? (p + z) : (p - z);
864 lambda[en-1] = complex<double>(x + z);
865 lambda[en] = complex<double>((z != 0.0) ? (x - w / z) : (x + z));
870 r =
sqrt(p * p + q * q);
875 for(
int j = en - 1; j <
n; j++) {
877 h[en-1][j] = q * z + p * h[en][j];
878 h[en][j] = q * h[en][j] - p * z;
882 for(
int i = 0; i <= en; i++) {
884 h[i][en-1] = q * z + p * h[i][en];
885 h[i][en] = q * h[i][en] - p * z;
889 for(
int i = low; i <= upp; i++) {
892 vectors[i][en] = q * vectors[i][en] - p * z;
896 lambda[en-1] = complex<double>(x + p, + z);
897 lambda[en] = complex<double>(x + p, - z);
905 if(norm == 0.0)
return 0;
907 for(
int en = n; en-- > 0;) {
918 h[en-1][en-1] = q / h[en][en-1];
919 h[en-1][en] = - (h[en][en] - p) / h[en][en-1];
921 cdiv(0.0, - h[en-1][en], h[en-1][en-1] - p, q,
922 h[en-1][en-1], h[en-1][en]);
929 for(
int i = en - 1; i-- > 0;) {
934 for(
int j = m; j <= en; j++) {
935 ra += h[i][j] * h[j][en-1];
936 sa += h[i][j] * h[j][en];
946 cdiv(- ra, -sa, w, q, h[i][en-1], h[i][en]);
955 if(vr == 0.0 && vi == 0.0) {
963 }
while(tst2 > tst1);
966 cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi,
967 h[i][en-1], h[i][en]);
970 h[i+1][en-1] = (- ra - w * h[i][en-1] + q * h[i][en]) / x;
971 h[i+1][en] = (- sa - w * h[i][en] - q * h[i][en-1]) / x;
973 cdiv(- r - y * h[i][en-1], - s - y * h[i][en], z, q,
974 h[i+1][en-1], h[i+1][en]);
981 if(t != 0.0 && (t + 1.0 / t) <= t) {
982 for(
int j = i; j <= en; j++) {
997 for(
int i = en; i-- > 0;) {
1000 for(
int j = m; j <= en; j++) r += h[i][j] * h[j][en];
1015 }
while(tst2 > tst1);
1025 t = (x * s - z * r) / q;
1028 (- r - w * t) / x : (- s - y * t) / z;
1033 if(t != 0.0 && (t + 1.0 / t) <= t) {
1034 for(
int j = i; j <= en; j++) h[j][en] /= t;
1044 for(
int i = 0; i < low; i++) {
1045 for(
int j = i; j <
n; j++)
vectors[i][j] = h[i][j];
1048 for(
int i = upp + 1; i <
n; i++) {
1049 for(
int j = i; j <
n; j++)
vectors[i][j] = h[i][j];
1054 for(
int j = n; j-- > low;) {
1055 int m = (j < upp) ? j : upp;
1057 for(
int i = low; i <= upp; i++) {
1059 for(k = low; k <= m; k++) z +=
vectors[i][k] * h[k][j];
static void balance(Matrix< double > &, int &low, int &high, Vector< double > &)
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
Double eigenvector routines.
void cdiv(double ar, double ai, double br, double bi, double &cr, double &ci)
static void exchange(Matrix< double > &, int j, int m, int low, int high)
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Matrix< double > packedEigenVectors() const
Get eigenvectors.
int ncols() const
Get number of columns.
void elmtran(Matrix< double > &, int low, int high, Array1D< int > &index)
Vector< complex< double > > lambda
FLieGenerator< T, N > real(const FLieGenerator< std::complex< T >, N > &)
Take real part of a complex generator.
Matrix< complex< double > > eigenVectors() const
Get eigenvectors.
constexpr double c
The velocity of light in m/s.
int hqr(Matrix< double > &, int low, int high)
int hqr2(Matrix< double > &, int low, int high)
Eigenvalue error exception.
Tps< T > sqrt(const Tps< T > &x)
Square root.
static void elmhes(Matrix< double > &, int low, int high, Array1D< int > &index)
FLieGenerator< T, N > imag(const FLieGenerator< std::complex< T >, N > &)
Take imaginary part of a complex generator.
int nrows() const
Get number of rows.
Vector< complex< double > > eigenValues() const
Get eigenvalues.
void balbak(int low, int high, Vector< double > &scale)
void swapRows(int row1, int row2)
Exchange rows.