18 #include "EigenvalueDecomposition.h"
22 EigenvalueDecomposition::EigenvalueDecomposition(){
23 warningLog.setProceedingText(
"[WARNING EigenvalueDecomposition]");
26 EigenvalueDecomposition::~EigenvalueDecomposition(){
30 bool EigenvalueDecomposition::decompose(
const MatrixDouble &a){
34 realEigenvalues.resize(n);
35 complexEigenvalues.resize(n);
38 for(
int j = 0; (j < n) & issymmetric; j++) {
39 for(
int i = 0; (i < n) & issymmetric; i++) {
40 issymmetric = (a[i][j] == a[j][i]);
45 for(
int i = 0; i < n; i++) {
46 for(
int j = 0; j < n; j++) {
47 eigenvectors[i][j] = a[i][j];
61 for(
int j = 0; j < n; j++) {
62 for(
int i = 0; i < n; i++) {
79 for(
int j = 0; j < n; j++) {
80 realEigenvalues[j] = eigenvectors[n-1][j];
84 for(
int i = n-1; i > 0; i--) {
89 for (
int k = 0; k < i; k++) {
90 scale = scale + fabs(realEigenvalues[k]);
93 complexEigenvalues[i] = realEigenvalues[i-1];
94 for (
int j = 0; j < i; j++) {
95 realEigenvalues[j] = eigenvectors[i-1][j];
96 eigenvectors[i][j] = 0.0;
97 eigenvectors[j][i] = 0.0;
102 for (
int k = 0; k < i; k++) {
103 realEigenvalues[k] /= scale;
104 h += realEigenvalues[k] * realEigenvalues[k];
106 double f = realEigenvalues[i-1];
111 complexEigenvalues[i] = scale * g;
113 realEigenvalues[i-1] = f - g;
114 for (
int j = 0; j < i; j++) {
115 complexEigenvalues[j] = 0.0;
119 for (
int j = 0; j < i; j++) {
120 f = realEigenvalues[j];
121 eigenvectors[j][i] = f;
122 g = complexEigenvalues[j] + eigenvectors[j][j] * f;
123 for (
int k = j+1; k <= i-1; k++) {
124 g += eigenvectors[k][j] * realEigenvalues[k];
125 complexEigenvalues[k] += eigenvectors[k][j] * f;
127 complexEigenvalues[j] = g;
130 for (
int j = 0; j < i; j++) {
131 complexEigenvalues[j] /= h;
132 f += complexEigenvalues[j] * realEigenvalues[j];
134 double hh = f / (h + h);
135 for (
int j = 0; j < i; j++) {
136 complexEigenvalues[j] -= hh * realEigenvalues[j];
138 for (
int j = 0; j < i; j++) {
139 f = realEigenvalues[j];
140 g = complexEigenvalues[j];
141 for (
int k = j; k <= i-1; k++) {
142 eigenvectors[k][j] -= (f * complexEigenvalues[k] + g * realEigenvalues[k]);
144 realEigenvalues[j] = eigenvectors[i-1][j];
145 eigenvectors[i][j] = 0.0;
148 realEigenvalues[i] = h;
152 for(
int i = 0; i < n-1; i++) {
153 eigenvectors[n-1][i] = eigenvectors[i][i];
154 eigenvectors[i][i] = 1.0;
155 double h = realEigenvalues[i+1];
157 for (
int k = 0; k <= i; k++) {
158 realEigenvalues[k] = eigenvectors[k][i+1] / h;
160 for (
int j = 0; j <= i; j++) {
162 for (
int k = 0; k <= i; k++) {
163 g += eigenvectors[k][i+1] * eigenvectors[k][j];
165 for (
int k = 0; k <= i; k++) {
166 eigenvectors[k][j] -= g * realEigenvalues[k];
170 for(
int k = 0; k <= i; k++) {
171 eigenvectors[k][i+1] = 0.0;
174 for(
int j = 0; j < n; j++) {
175 realEigenvalues[j] = eigenvectors[n-1][j];
176 eigenvectors[n-1][j] = 0.0;
178 eigenvectors[n-1][n-1] = 1.0;
179 complexEigenvalues[0] = 0.0;
186 for(
int i = 1; i < n; i++) {
187 complexEigenvalues[i-1] = complexEigenvalues[i];
189 complexEigenvalues[n-1] = 0.0;
193 double eps = pow(2.0,-52.0);
194 for (
int l = 0; l < n; l++) {
197 tst1 = findMax(tst1,fabs(realEigenvalues[l]) + fabs(complexEigenvalues[l]));
200 if(fabs(complexEigenvalues[m]) <= eps*tst1) {
213 double g = realEigenvalues[l];
214 double p = (realEigenvalues[l+1] - g) / (2.0 * complexEigenvalues[l]);
215 double r = hypot(p,1.0);
219 realEigenvalues[l] = complexEigenvalues[l] / (p + r);
220 realEigenvalues[l+1] = complexEigenvalues[l] * (p + r);
221 double dl1 = realEigenvalues[l+1];
222 double h = g - realEigenvalues[l];
223 for (
int i = l+2; i < n; i++) {
224 realEigenvalues[i] -= h;
229 p = realEigenvalues[m];
233 double el1 = complexEigenvalues[l+1];
236 for (
int i = m-1; i >= l; i--) {
240 g = c * complexEigenvalues[i];
242 r = hypot(p,complexEigenvalues[i]);
243 complexEigenvalues[i+1] = s * r;
244 s = complexEigenvalues[i] / r;
246 p = c * realEigenvalues[i] - s * g;
247 realEigenvalues[i+1] = h + s * (c * g + s * realEigenvalues[i]);
250 for(
int k = 0; k < n; k++) {
251 h = eigenvectors[k][i+1];
252 eigenvectors[k][i+1] = s * eigenvectors[k][i] + c * h;
253 eigenvectors[k][i] = c * eigenvectors[k][i] - s * h;
256 p = -s * s2 * c3 * el1 * complexEigenvalues[l] / dl1;
257 complexEigenvalues[l] = s * p;
258 realEigenvalues[l] = c * p;
261 }
while (fabs(complexEigenvalues[l]) > eps*tst1);
263 realEigenvalues[l] = realEigenvalues[l] + f;
264 complexEigenvalues[l] = 0.0;
268 for(
int i = 0; i < n-1; i++) {
270 double p = realEigenvalues[i];
271 for (
int j = i+1; j < n; j++) {
272 if(realEigenvalues[j] < p) {
274 p = realEigenvalues[j];
278 realEigenvalues[k] = realEigenvalues[i];
279 realEigenvalues[i] = p;
280 for (
int j = 0; j < n; j++) {
281 p = eigenvectors[j][i];
282 eigenvectors[j][i] = eigenvectors[j][k];
283 eigenvectors[j][k] = p;
295 for(
int m = low+1; m <= high-1; m++) {
299 for (
int i = m; i <= high; i++) {
300 scale = scale + fabs(h[i][m-1]);
306 for(
int i = high; i >= m; i--) {
307 ort[i] = h[i][m-1]/scale;
308 ht += ort[i] * ort[i];
310 double g = sqrt( ht );
314 ht = ht - ort[m] * g;
319 for (
int j = m; j < n; j++) {
321 for (
int i = high; i >= m; i--) {
325 for (
int i = m; i <= high; i++) {
330 for(
int i = 0; i <= high; i++) {
332 for(
int j = high; j >= m; j--) {
336 for (
int j = m; j <= high; j++) {
340 ort[m] = scale*ort[m];
346 for (
int i = 0; i < n; i++) {
347 for (
int j = 0; j < n; j++) {
348 eigenvectors[i][j] = (i == j ? 1.0 : 0.0);
352 for (
int m = high-1; m >= low+1; m--) {
353 if (h[m][m-1] != 0.0) {
354 for (
int i = m+1; i <= high; i++) {
357 for (
int j = m; j <= high; j++) {
359 for (
int i = m; i <= high; i++) {
360 g += ort[i] * eigenvectors[i][j];
363 g = (g / ort[m]) / h[m][m-1];
364 for (
int i = m; i <= high; i++) {
365 eigenvectors[i][j] += g * ort[i];
380 double eps = pow(2.0,-52.0);
381 double exshift = 0.0;
382 double p=0,q=0,r=0,s=0,z=0,t,w,x,y;
386 for(
int i = 0; i < nn; i++) {
387 if( (i < low) | (i > high) ){
388 realEigenvalues[i] = h[i][i];
389 complexEigenvalues[i] = 0.0;
391 for (
int j = findMax(i-1,0); j < nn; j++) {
392 norm = norm + fabs(h[i][j]);
403 s = fabs(h[l-1][l-1]) + fabs(h[l][l]);
407 if(fabs(h[l][l-1]) < eps * s) {
416 h[n][n] = h[n][n] + exshift;
417 realEigenvalues[n] = h[n][n];
418 complexEigenvalues[n] = 0.0;
423 }
else if (l == n-1) {
424 w = h[n][n-1] * h[n-1][n];
425 p = (h[n-1][n-1] - h[n][n]) / 2.0;
428 h[n][n] = h[n][n] + exshift;
429 h[n-1][n-1] = h[n-1][n-1] + exshift;
439 realEigenvalues[n-1] = x + z;
440 realEigenvalues[n] = realEigenvalues[n-1];
442 realEigenvalues[n] = x - w / z;
444 complexEigenvalues[n-1] = 0.0;
445 complexEigenvalues[n] = 0.0;
447 s = fabs(x) + fabs(z);
450 r = sqrt(p * p+q * q);
455 for(
int j = n-1; j < nn; j++) {
457 h[n-1][j] = q * z + p * h[n][j];
458 h[n][j] = q * h[n][j] - p * z;
462 for(
int i = 0; i <= n; i++) {
464 h[i][n-1] = q * z + p * h[i][n];
465 h[i][n] = q * h[i][n] - p * z;
469 for(
int i = low; i <= high; i++) {
470 z = eigenvectors[i][n-1];
471 eigenvectors[i][n-1] = q * z + p * eigenvectors[i][n];
472 eigenvectors[i][n] = q * eigenvectors[i][n] - p * z;
477 realEigenvalues[n-1] = x + p;
478 realEigenvalues[n] = x + p;
479 complexEigenvalues[n-1] = z;
480 complexEigenvalues[n] = -z;
494 w = h[n][n-1] * h[n-1][n];
500 for (
int i = low; i <= n; i++) {
503 s = fabs(h[n][n-1]) + fabs(h[n-1][n-2]);
517 s = x - w / ((y - x) / 2.0 + s);
518 for(
int i = low; i <= n; i++) {
534 p = (r * s - w) / h[m+1][m] + h[m][m+1];
535 q = h[m+1][m+1] - z - r - s;
537 s = fabs(p) + fabs(q) + fabs(r);
544 if(fabs(h[m][m-1]) * (fabs(q) +fabs(r)) <
545 eps * (fabs(p) * (fabs(h[m-1][m-1]) + fabs(z) +
546 fabs(h[m+1][m+1])))) {
552 for(
int i = m+2; i <= n; i++) {
560 for(
int k = m; k <= n-1; k++) {
561 bool notlast = (k != n-1);
565 r = (notlast ? h[k+2][k-1] : 0.0);
566 x = fabs(p) + fabs(q) + fabs(r);
576 s = sqrt(p * p + q * q + r * r);
584 h[k][k-1] = -h[k][k-1];
594 for(
int j = k; j < nn; j++) {
595 p = h[k][j] + q * h[k+1][j];
597 p = p + r * h[k+2][j];
598 h[k+2][j] = h[k+2][j] - p * z;
600 h[k][j] = h[k][j] - p * x;
601 h[k+1][j] = h[k+1][j] - p * y;
605 for (
int i = 0; i <= findMin(n,k+3); i++) {
606 p = x * h[i][k] + y * h[i][k+1];
608 p = p + z * h[i][k+2];
609 h[i][k+2] = h[i][k+2] - p * r;
611 h[i][k] = h[i][k] - p;
612 h[i][k+1] = h[i][k+1] - p * q;
616 for(
int i = low; i <= high; i++) {
617 p = x * eigenvectors[i][k] + y * eigenvectors[i][k+1];
619 p = p + z * eigenvectors[i][k+2];
620 eigenvectors[i][k+2] = eigenvectors[i][k+2] - p * r;
622 eigenvectors[i][k] = eigenvectors[i][k] - p;
623 eigenvectors[i][k+1] = eigenvectors[i][k+1] - p * q;
635 for(n = nn-1; n >= 0; n--) {
636 p = realEigenvalues[n];
637 q = complexEigenvalues[n];
643 for(
int i = n-1; i >= 0; i--) {
646 for(
int j = l; j <= n; j++) {
647 r = r + h[i][j] * h[j][n];
649 if(complexEigenvalues[i] < 0.0) {
654 if (complexEigenvalues[i] == 0.0) {
658 h[i][n] = -r / (eps * norm);
665 q = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) + complexEigenvalues[i] * complexEigenvalues[i];
666 t = (x * s - z * r) / q;
668 if(fabs(x) > fabs(z)) {
669 h[i+1][n] = (-r - w * t) / x;
671 h[i+1][n] = (-s - y * t) / z;
677 if ((eps * t) * t > 1) {
678 for(
int j = i; j <= n; j++) {
679 h[j][n] = h[j][n] / t;
690 if (fabs(h[n][n-1]) > fabs(h[n-1][n])) {
691 h[n-1][n-1] = q / h[n][n-1];
692 h[n-1][n] = -(h[n][n] - p) / h[n][n-1];
694 cdiv(0.0,-h[n-1][n],h[n-1][n-1]-p,q);
700 for(
int i = n-2; i >= 0; i--) {
704 for (
int j = l; j <= n; j++) {
705 ra = ra + h[i][j] * h[j][n-1];
706 sa = sa + h[i][j] * h[j][n];
710 if(complexEigenvalues[i] < 0.0) {
716 if(complexEigenvalues[i] == 0) {
725 vr = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) + complexEigenvalues[i] * complexEigenvalues[i] - q * q;
726 vi = (realEigenvalues[i] - p) * 2.0 * q;
727 if((vr == 0.0) & (vi == 0.0)){
728 vr = eps * norm * (fabs(w) + fabs(q) + fabs(x) + fabs(y) + fabs(z));
730 cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi);
733 if (fabs(x) > (fabs(z) + fabs(q))) {
734 h[i+1][n-1] = (-ra - w * h[i][n-1] + q * h[i][n]) / x;
735 h[i+1][n] = (-sa - w * h[i][n] - q * h[i][n-1]) / x;
737 cdiv(-r-y*h[i][n-1],-s-y*h[i][n],z,q);
744 t = findMax(fabs(h[i][n-1]),fabs(h[i][n]));
745 if ((eps * t) * t > 1) {
746 for(
int j = i; j <= n; j++) {
747 h[j][n-1] = h[j][n-1] / t;
748 h[j][n] = h[j][n] / t;
757 for (
int i = 0; i < nn; i++) {
758 if((i < low) | (i > high)){
759 for (
int j = i; j < nn; j++) {
760 eigenvectors[i][j] = h[i][j];
766 for (
int j = nn-1; j >= low; j--) {
767 for (
int i = low; i <= high; i++) {
769 for (
int k = low; k <= findMin(j,high); k++) {
770 z = z + eigenvectors[i][k] * h[k][j];
772 eigenvectors[i][j] = z;
781 if(fabs(yr) > fabs(yi)){
784 cdivr = (xr + r*xi)/d;
785 cdivi = (xi - r*xr)/d;
789 cdivr = (r*xr + xi)/d;
790 cdivi = (r*xi - xr)/d;
798 for (
int i = 0; i < n; i++) {
799 for (
int j = 0; j < n; j++) {
802 x[i][i] = realEigenvalues[i];
803 if(complexEigenvalues[i] > 0) {
804 x[i][i+1] = complexEigenvalues[i];
805 }
else if(complexEigenvalues[i] < 0) {
806 x[i][i-1] = complexEigenvalues[i];
813 return realEigenvalues;
817 return complexEigenvalues;
VectorDouble getComplexEigenvalues()
VectorDouble getRealEigenvalues()
void cdiv(double xr, double xi, double yr, double yi)
MatrixDouble getDiagonalEigenvalueMatrix()
virtual bool resize(const unsigned int r, const unsigned int c)