OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
Eigenvalue.hpp
1#ifndef JAMA_EIG_H
2#define JAMA_EIG_H
3
4#include <algorithm>
5
7// for min(), max() below
8#include <cmath>
9// for abs() below
10
11using namespace OpenMD;
12using namespace std;
13
14namespace JAMA {
15
16 /**
17
18 Computes eigenvalues and eigenvectors of a real (non-complex)
19 matrix.
20 <P>
21 If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
22 diagonal and the eigenvector matrix V is orthogonal. That is,
23 the diagonal values of D are the eigenvalues, and
24 V*V' = I, where I is the identity matrix. The columns of V
25 represent the eigenvectors in the sense that A*V = V*D.
26
27 <P>
28 If A is not symmetric, then the eigenvalue matrix D is block diagonal
29 with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
30 a + i*b, in 2-by-2 blocks, [a, b; -b, a]. That is, if the complex
31 eigenvalues look like
32<pre>
33
34 u + iv . . . . .
35 . u - iv . . . .
36 . . a + ib . . .
37 . . . a - ib . .
38 . . . . x .
39 . . . . . y
40</pre>
41 then D looks like
42<pre>
43
44 u v . . . .
45 -v u . . . .
46 . . a b . .
47 . . -b a . .
48 . . . . x .
49 . . . . . y
50</pre>
51 This keeps V a real matrix in both symmetric and non-symmetric
52 cases, and A*V = V*D.
53
54 <p>
55 The matrix V may be badly
56 conditioned, or even singular, so the validity of the equation
57 A = V*D*inverse(V) depends upon the condition number of V.
58
59 <p>
60 (Adapted from JAMA, a Java Matrix Library, developed by jointly
61 by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama).
62 **/
63
64 template<class Real>
65 class Eigenvalue {
66 /** Row and column dimension (square matrix). */
67 int n;
68
69 int issymmetric; /* boolean*/
70
71 /** Arrays for internal storage of eigenvalues. */
72
73 DynamicVector<Real> d; /* real part */
74 DynamicVector<Real> e; /* img part */
75
76 /** Array for internal storage of eigenvectors. */
77 DynamicRectMatrix<Real> V;
78
79 /** Array for internal storage of nonsymmetric Hessenberg form.
80 @serial internal storage of nonsymmetric Hessenberg form.
81 */
82 DynamicRectMatrix<Real> H;
83
84 /** Working storage for nonsymmetric algorithm.
85 @serial working storage for nonsymmetric algorithm.
86 */
87 DynamicVector<Real> ort;
88
89 // Symmetric Householder reduction to tridiagonal form.
90
91 void tred2() {
92 // This is derived from the Algol procedures tred2 by
93 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
94 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
95 // Fortran subroutine in EISPACK.
96
97 for (int j = 0; j < n; j++) {
98 d(j) = V(n - 1, j);
99 }
100
101 // Householder reduction to tridiagonal form.
102
103 for (int i = n - 1; i > 0; i--) {
104 // Scale to avoid under/overflow.
105
106 Real scale = 0.0;
107 Real h = 0.0;
108 for (int k = 0; k < i; k++) {
109 scale = scale + abs(d(k));
110 }
111 if (scale == 0.0) {
112 e(i) = d(i - 1);
113 for (int j = 0; j < i; j++) {
114 d(j) = V(i - 1, j);
115 V(i, j) = 0.0;
116 V(j, i) = 0.0;
117 }
118 } else {
119 // Generate Householder vector.
120
121 for (int k = 0; k < i; k++) {
122 d(k) /= scale;
123 h += d(k) * d(k);
124 }
125 Real f = d(i - 1);
126 Real g = sqrt(h);
127 if (f > 0) { g = -g; }
128 e(i) = scale * g;
129 h = h - f * g;
130 d(i - 1) = f - g;
131 for (int j = 0; j < i; j++) {
132 e(j) = 0.0;
133 }
134
135 // Apply similarity transformation to remaining columns.
136
137 for (int j = 0; j < i; j++) {
138 f = d(j);
139 V(j, i) = f;
140 g = e(j) + V(j, j) * f;
141 for (int k = j + 1; k <= i - 1; k++) {
142 g += V(k, j) * d(k);
143 e(k) += V(k, j) * f;
144 }
145 e(j) = g;
146 }
147 f = 0.0;
148 for (int j = 0; j < i; j++) {
149 e(j) /= h;
150 f += e(j) * d(j);
151 }
152 Real hh = f / (h + h);
153 for (int j = 0; j < i; j++) {
154 e(j) -= hh * d(j);
155 }
156 for (int j = 0; j < i; j++) {
157 f = d(j);
158 g = e(j);
159 for (int k = j; k <= i - 1; k++) {
160 V(k, j) -= (f * e(k) + g * d(k));
161 }
162 d(j) = V(i - 1, j);
163 V(i, j) = 0.0;
164 }
165 }
166 d(i) = h;
167 }
168
169 // Accumulate transformations.
170
171 for (int i = 0; i < n - 1; i++) {
172 V(n - 1, i) = V(i, i);
173 V(i, i) = 1.0;
174 Real h = d(i + 1);
175 if (h != 0.0) {
176 for (int k = 0; k <= i; k++) {
177 d(k) = V(k, i + 1) / h;
178 }
179 for (int j = 0; j <= i; j++) {
180 Real g = 0.0;
181 for (int k = 0; k <= i; k++) {
182 g += V(k, i + 1) * V(k, j);
183 }
184 for (int k = 0; k <= i; k++) {
185 V(k, j) -= g * d(k);
186 }
187 }
188 }
189 for (int k = 0; k <= i; k++) {
190 V(k, i + 1) = 0.0;
191 }
192 }
193 for (int j = 0; j < n; j++) {
194 d(j) = V(n - 1, j);
195 V(n - 1, j) = 0.0;
196 }
197 V(n - 1, n - 1) = 1.0;
198 e(0) = 0.0;
199 }
200
201 // Symmetric tridiagonal QL algorithm.
202
203 void tql2() {
204 // This is derived from the Algol procedures tql2, by
205 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
206 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
207 // Fortran subroutine in EISPACK.
208
209 for (int i = 1; i < n; i++) {
210 e(i - 1) = e(i);
211 }
212 e(n - 1) = 0.0;
213
214 Real f = 0.0;
215 Real tst1 = 0.0;
216 Real eps = pow(2.0, -52.0);
217 for (int l = 0; l < n; l++) {
218 // Find small subdiagonal element
219
220 tst1 = max(tst1, abs(d(l)) + abs(e(l)));
221 int m = l;
222
223 // Original while-loop from Java code
224 while (m < n) {
225 if (abs(e(m)) <= eps * tst1) { break; }
226 m++;
227 }
228
229 // If m == l, d(l) is an eigenvalue,
230 // otherwise, iterate.
231
232 if (m > l) {
233 int iter = 0;
234 do {
235 iter = iter + 1; // (Could check iteration count here.)
236
237 // Compute implicit shift
238
239 Real g = d(l);
240 Real p = (d(l + 1) - g) / (2.0 * e(l));
241 Real r = hypot(p, 1.0);
242 if (p < 0) { r = -r; }
243 d(l) = e(l) / (p + r);
244 d(l + 1) = e(l) * (p + r);
245 Real dl1 = d(l + 1);
246 Real h = g - d(l);
247 for (int i = l + 2; i < n; i++) {
248 d(i) -= h;
249 }
250 f = f + h;
251
252 // Implicit QL transformation.
253
254 p = d(m);
255 Real c = 1.0;
256 Real c2 = c;
257 Real c3 = c;
258 Real el1 = e(l + 1);
259 Real s = 0.0;
260 Real s2 = 0.0;
261 for (int i = m - 1; i >= l; i--) {
262 c3 = c2;
263 c2 = c;
264 s2 = s;
265 g = c * e(i);
266 h = c * p;
267 r = hypot(p, e(i));
268 e(i + 1) = s * r;
269 s = e(i) / r;
270 c = p / r;
271 p = c * d(i) - s * g;
272 d(i + 1) = h + s * (c * g + s * d(i));
273
274 // Accumulate transformation.
275
276 for (int k = 0; k < n; k++) {
277 h = V(k, i + 1);
278 V(k, i + 1) = s * V(k, i) + c * h;
279 V(k, i) = c * V(k, i) - s * h;
280 }
281 }
282 p = -s * s2 * c3 * el1 * e(l) / dl1;
283 e(l) = s * p;
284 d(l) = c * p;
285
286 // Check for convergence.
287
288 } while (abs(e(l)) > eps * tst1);
289 }
290 d(l) = d(l) + f;
291 e(l) = 0.0;
292 }
293
294 // Sort eigenvalues and corresponding vectors.
295
296 for (int i = 0; i < n - 1; i++) {
297 int k = i;
298 Real p = d(i);
299 for (int j = i + 1; j < n; j++) {
300 if (d(j) < p) {
301 k = j;
302 p = d(j);
303 }
304 }
305 if (k != i) {
306 d(k) = d(i);
307 d(i) = p;
308 for (int j = 0; j < n; j++) {
309 p = V(j, i);
310 V(j, i) = V(j, k);
311 V(j, k) = p;
312 }
313 }
314 }
315 }
316
317 // Nonsymmetric reduction to Hessenberg form.
318
319 void orthes() {
320 // This is derived from the Algol procedures orthes and ortran,
321 // by Martin and Wilkinson, Handbook for Auto. Comp.,
322 // Vol.ii-Linear Algebra, and the corresponding
323 // Fortran subroutines in EISPACK.
324
325 int low = 0;
326 int high = n - 1;
327
328 for (int m = low + 1; m <= high - 1; m++) {
329 // Scale column.
330
331 Real scale = 0.0;
332 for (int i = m; i <= high; i++) {
333 scale = scale + abs(H(i, m - 1));
334 }
335 if (scale != 0.0) {
336 // Compute Householder transformation.
337
338 Real h = 0.0;
339 for (int i = high; i >= m; i--) {
340 ort(i) = H(i, m - 1) / scale;
341 h += ort(i) * ort(i);
342 }
343 Real g = sqrt(h);
344 if (ort(m) > 0) { g = -g; }
345 h = h - ort(m) * g;
346 ort(m) = ort(m) - g;
347
348 // Apply Householder similarity transformation
349 // H = (I-u*u'/h)*H*(I-u*u')/h)
350
351 for (int j = m; j < n; j++) {
352 Real f = 0.0;
353 for (int i = high; i >= m; i--) {
354 f += ort(i) * H(i, j);
355 }
356 f = f / h;
357 for (int i = m; i <= high; i++) {
358 H(i, j) -= f * ort(i);
359 }
360 }
361
362 for (int i = 0; i <= high; i++) {
363 Real f = 0.0;
364 for (int j = high; j >= m; j--) {
365 f += ort(j) * H(i, j);
366 }
367 f = f / h;
368 for (int j = m; j <= high; j++) {
369 H(i, j) -= f * ort(j);
370 }
371 }
372 ort(m) = scale * ort(m);
373 H(m, m - 1) = scale * g;
374 }
375 }
376
377 // Accumulate transformations (Algol's ortran).
378
379 for (int i = 0; i < n; i++) {
380 for (int j = 0; j < n; j++) {
381 V(i, j) = (i == j ? 1.0 : 0.0);
382 }
383 }
384
385 for (int m = high - 1; m >= low + 1; m--) {
386 if (H(m, m - 1) != 0.0) {
387 for (int i = m + 1; i <= high; i++) {
388 ort(i) = H(i, m - 1);
389 }
390 for (int j = m; j <= high; j++) {
391 Real g = 0.0;
392 for (int i = m; i <= high; i++) {
393 g += ort(i) * V(i, j);
394 }
395 // Double division avoids possible underflow
396 g = (g / ort(m)) / H(m, m - 1);
397 for (int i = m; i <= high; i++) {
398 V(i, j) += g * ort(i);
399 }
400 }
401 }
402 }
403 }
404
405 // Complex scalar division.
406
407 Real cdivr, cdivi;
408 void cdiv(Real xr, Real xi, Real yr, Real yi) {
409 Real r, d;
410 if (abs(yr) > abs(yi)) {
411 r = yi / yr;
412 d = yr + r * yi;
413 cdivr = (xr + r * xi) / d;
414 cdivi = (xi - r * xr) / d;
415 } else {
416 r = yr / yi;
417 d = yi + r * yr;
418 cdivr = (r * xr + xi) / d;
419 cdivi = (r * xi - xr) / d;
420 }
421 }
422
423 // Nonsymmetric reduction from Hessenberg to real Schur form.
424
425 void hqr2() {
426 // This is derived from the Algol procedure hqr2,
427 // by Martin and Wilkinson, Handbook for Auto. Comp.,
428 // Vol.ii-Linear Algebra, and the corresponding
429 // Fortran subroutine in EISPACK.
430
431 // Initialize
432
433 int nn = this->n;
434 int n = nn - 1;
435 int low = 0;
436 int high = nn - 1;
437 Real eps = pow(2.0, -52.0);
438 Real exshift = 0.0;
439 Real p = 0, q = 0, r = 0, s = 0, z = 0, t, w, x, y;
440
441 // Store roots isolated by balanc and compute matrix norm
442
443 Real norm = 0.0;
444 for (int i = 0; i < nn; i++) {
445 if ((i < low) || (i > high)) {
446 d(i) = H(i, i);
447 e(i) = 0.0;
448 }
449 for (int j = max(i - 1, 0); j < nn; j++) {
450 norm = norm + abs(H(i, j));
451 }
452 }
453
454 // Outer loop over eigenvalue index
455
456 int iter = 0;
457 while (n >= low) {
458 // Look for single small sub-diagonal element
459
460 int l = n;
461 while (l > low) {
462 s = abs(H(l - 1, l - 1)) + abs(H(l, l));
463 if (s == 0.0) { s = norm; }
464 if (abs(H(l, l - 1)) < eps * s) { break; }
465 l--;
466 }
467
468 // Check for convergence
469 // One root found
470
471 if (l == n) {
472 H(n, n) = H(n, n) + exshift;
473 d(n) = H(n, n);
474 e(n) = 0.0;
475 n--;
476 iter = 0;
477
478 // Two roots found
479
480 } else if (l == n - 1) {
481 w = H(n, n - 1) * H(n - 1, n);
482 p = (H(n - 1, n - 1) - H(n, n)) / 2.0;
483 q = p * p + w;
484 z = sqrt(abs(q));
485 H(n, n) = H(n, n) + exshift;
486 H(n - 1, n - 1) = H(n - 1, n - 1) + exshift;
487 x = H(n, n);
488
489 // Real pair
490
491 if (q >= 0) {
492 if (p >= 0) {
493 z = p + z;
494 } else {
495 z = p - z;
496 }
497 d(n - 1) = x + z;
498 d(n) = d(n - 1);
499 if (z != 0.0) { d(n) = x - w / z; }
500 e(n - 1) = 0.0;
501 e(n) = 0.0;
502 x = H(n, n - 1);
503 s = abs(x) + abs(z);
504 p = x / s;
505 q = z / s;
506 r = sqrt(p * p + q * q);
507 p = p / r;
508 q = q / r;
509
510 // Row modification
511
512 for (int j = n - 1; j < nn; j++) {
513 z = H(n - 1, j);
514 H(n - 1, j) = q * z + p * H(n, j);
515 H(n, j) = q * H(n, j) - p * z;
516 }
517
518 // Column modification
519
520 for (int i = 0; i <= n; i++) {
521 z = H(i, n - 1);
522 H(i, n - 1) = q * z + p * H(i, n);
523 H(i, n) = q * H(i, n) - p * z;
524 }
525
526 // Accumulate transformations
527
528 for (int i = low; i <= high; i++) {
529 z = V(i, n - 1);
530 V(i, n - 1) = q * z + p * V(i, n);
531 V(i, n) = q * V(i, n) - p * z;
532 }
533
534 // Complex pair
535
536 } else {
537 d(n - 1) = x + p;
538 d(n) = x + p;
539 e(n - 1) = z;
540 e(n) = -z;
541 }
542 n = n - 2;
543 iter = 0;
544
545 // No convergence yet
546
547 } else {
548 // Form shift
549
550 x = H(n, n);
551 y = 0.0;
552 w = 0.0;
553 if (l < n) {
554 y = H(n - 1, n - 1);
555 w = H(n, n - 1) * H(n - 1, n);
556 }
557
558 // Wilkinson's original ad hoc shift
559
560 if (iter == 10) {
561 exshift += x;
562 for (int i = low; i <= n; i++) {
563 H(i, i) -= x;
564 }
565 s = abs(H(n, n - 1)) + abs(H(n - 1, n - 2));
566 x = y = 0.75 * s;
567 w = -0.4375 * s * s;
568 }
569
570 // MATLAB's new ad hoc shift
571
572 if (iter == 30) {
573 s = (y - x) / 2.0;
574 s = s * s + w;
575 if (s > 0) {
576 s = sqrt(s);
577 if (y < x) { s = -s; }
578 s = x - w / ((y - x) / 2.0 + s);
579 for (int i = low; i <= n; i++) {
580 H(i, i) -= s;
581 }
582 exshift += s;
583 x = y = w = 0.964;
584 }
585 }
586
587 iter = iter + 1; // (Could check iteration count here.)
588
589 // Look for two consecutive small sub-diagonal elements
590
591 int m = n - 2;
592 while (m >= l) {
593 z = H(m, m);
594 r = x - z;
595 s = y - z;
596 p = (r * s - w) / H(m + 1, m) + H(m, m + 1);
597 q = H(m + 1, m + 1) - z - r - s;
598 r = H(m + 2, m + 1);
599 s = abs(p) + abs(q) + abs(r);
600 p = p / s;
601 q = q / s;
602 r = r / s;
603 if (m == l) { break; }
604 if (abs(H(m, m - 1)) * (abs(q) + abs(r)) <
605 eps * (abs(p) * (abs(H(m - 1, m - 1)) + abs(z) +
606 abs(H(m + 1, m + 1))))) {
607 break;
608 }
609 m--;
610 }
611
612 for (int i = m + 2; i <= n; i++) {
613 H(i, i - 2) = 0.0;
614 if (i > m + 2) { H(i, i - 3) = 0.0; }
615 }
616
617 // Double QR step involving rows l:n and columns m:n
618
619 for (int k = m; k <= n - 1; k++) {
620 int notlast = (k != n - 1);
621 if (k != m) {
622 p = H(k, k - 1);
623 q = H(k + 1, k - 1);
624 r = (notlast ? H(k + 2, k - 1) : 0.0);
625 x = abs(p) + abs(q) + abs(r);
626 if (x != 0.0) {
627 p = p / x;
628 q = q / x;
629 r = r / x;
630 }
631 }
632 if (x == 0.0) { break; }
633 s = sqrt(p * p + q * q + r * r);
634 if (p < 0) { s = -s; }
635 if (s != 0) {
636 if (k != m) {
637 H(k, k - 1) = -s * x;
638 } else if (l != m) {
639 H(k, k - 1) = -H(k, k - 1);
640 }
641 p = p + s;
642 x = p / s;
643 y = q / s;
644 z = r / s;
645 q = q / p;
646 r = r / p;
647
648 // Row modification
649
650 for (int j = k; j < nn; j++) {
651 p = H(k, j) + q * H(k + 1, j);
652 if (notlast) {
653 p = p + r * H(k + 2, j);
654 H(k + 2, j) = H(k + 2, j) - p * z;
655 }
656 H(k, j) = H(k, j) - p * x;
657 H(k + 1, j) = H(k + 1, j) - p * y;
658 }
659
660 // Column modification
661
662 for (int i = 0; i <= min(n, k + 3); i++) {
663 p = x * H(i, k) + y * H(i, k + 1);
664 if (notlast) {
665 p = p + z * H(i, k + 2);
666 H(i, k + 2) = H(i, k + 2) - p * r;
667 }
668 H(i, k) = H(i, k) - p;
669 H(i, k + 1) = H(i, k + 1) - p * q;
670 }
671
672 // Accumulate transformations
673
674 for (int i = low; i <= high; i++) {
675 p = x * V(i, k) + y * V(i, k + 1);
676 if (notlast) {
677 p = p + z * V(i, k + 2);
678 V(i, k + 2) = V(i, k + 2) - p * r;
679 }
680 V(i, k) = V(i, k) - p;
681 V(i, k + 1) = V(i, k + 1) - p * q;
682 }
683 } // (s != 0)
684 } // k loop
685 } // check convergence
686 } // while (n >= low)
687
688 // Backsubstitute to find vectors of upper triangular form
689
690 if (norm == 0.0) { return; }
691
692 for (n = nn - 1; n >= 0; n--) {
693 p = d(n);
694 q = e(n);
695
696 // Real vector
697
698 if (q == 0) {
699 int l = n;
700 H(n, n) = 1.0;
701 for (int i = n - 1; i >= 0; i--) {
702 w = H(i, i) - p;
703 r = 0.0;
704 for (int j = l; j <= n; j++) {
705 r = r + H(i, j) * H(j, n);
706 }
707 if (e(i) < 0.0) {
708 z = w;
709 s = r;
710 } else {
711 l = i;
712 if (e(i) == 0.0) {
713 if (w != 0.0) {
714 H(i, n) = -r / w;
715 } else {
716 H(i, n) = -r / (eps * norm);
717 }
718
719 // Solve real equations
720
721 } else {
722 x = H(i, i + 1);
723 y = H(i + 1, i);
724 q = (d(i) - p) * (d(i) - p) + e(i) * e(i);
725 t = (x * s - z * r) / q;
726 H(i, n) = t;
727 if (abs(x) > abs(z)) {
728 H(i + 1, n) = (-r - w * t) / x;
729 } else {
730 H(i + 1, n) = (-s - y * t) / z;
731 }
732 }
733
734 // Overflow control
735
736 t = abs(H(i, n));
737 if ((eps * t) * t > 1) {
738 for (int j = i; j <= n; j++) {
739 H(j, n) = H(j, n) / t;
740 }
741 }
742 }
743 }
744
745 // Complex vector
746
747 } else if (q < 0) {
748 int l = n - 1;
749
750 // Last vector component imaginary so matrix is triangular
751
752 if (abs(H(n, n - 1)) > abs(H(n - 1, n))) {
753 H(n - 1, n - 1) = q / H(n, n - 1);
754 H(n - 1, n) = -(H(n, n) - p) / H(n, n - 1);
755 } else {
756 cdiv(0.0, -H(n - 1, n), H(n - 1, n - 1) - p, q);
757 H(n - 1, n - 1) = cdivr;
758 H(n - 1, n) = cdivi;
759 }
760 H(n, n - 1) = 0.0;
761 H(n, n) = 1.0;
762 for (int i = n - 2; i >= 0; i--) {
763 Real ra, sa, vr, vi;
764 ra = 0.0;
765 sa = 0.0;
766 for (int j = l; j <= n; j++) {
767 ra = ra + H(i, j) * H(j, n - 1);
768 sa = sa + H(i, j) * H(j, n);
769 }
770 w = H(i, i) - p;
771
772 if (e(i) < 0.0) {
773 z = w;
774 r = ra;
775 s = sa;
776 } else {
777 l = i;
778 if (e(i) == 0) {
779 cdiv(-ra, -sa, w, q);
780 H(i, n - 1) = cdivr;
781 H(i, n) = cdivi;
782 } else {
783 // Solve complex equations
784
785 x = H(i, i + 1);
786 y = H(i + 1, i);
787 vr = (d(i) - p) * (d(i) - p) + e(i) * e(i) - q * q;
788 vi = (d(i) - p) * 2.0 * q;
789 if ((vr == 0.0) && (vi == 0.0)) {
790 vr =
791 eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(z));
792 }
793 cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi);
794 H(i, n - 1) = cdivr;
795 H(i, n) = cdivi;
796 if (abs(x) > (abs(z) + abs(q))) {
797 H(i + 1, n - 1) = (-ra - w * H(i, n - 1) + q * H(i, n)) / x;
798 H(i + 1, n) = (-sa - w * H(i, n) - q * H(i, n - 1)) / x;
799 } else {
800 cdiv(-r - y * H(i, n - 1), -s - y * H(i, n), z, q);
801 H(i + 1, n - 1) = cdivr;
802 H(i + 1, n) = cdivi;
803 }
804 }
805
806 // Overflow control
807
808 t = max(abs(H(i, n - 1)), abs(H(i, n)));
809 if ((eps * t) * t > 1) {
810 for (int j = i; j <= n; j++) {
811 H(j, n - 1) = H(j, n - 1) / t;
812 H(j, n) = H(j, n) / t;
813 }
814 }
815 }
816 }
817 }
818 }
819
820 // Vectors of isolated roots
821
822 for (int i = 0; i < nn; i++) {
823 if (i < low || i > high) {
824 for (int j = i; j < nn; j++) {
825 V(i, j) = H(i, j);
826 }
827 }
828 }
829
830 // Back transformation to get eigenvectors of original matrix
831
832 for (int j = nn - 1; j >= low; j--) {
833 for (int i = low; i <= high; i++) {
834 z = 0.0;
835 for (int k = low; k <= min(j, high); k++) {
836 z = z + V(i, k) * H(k, j);
837 }
838 V(i, j) = z;
839 }
840 }
841 }
842
843 public:
844 /** Check for symmetry, then construct the eigenvalue decomposition
845 @param A Square real (non-complex) matrix
846 */
847 Eigenvalue(const DynamicRectMatrix<Real>& A) {
848 n = A.getNCol();
849 V = DynamicRectMatrix<Real>(n, n);
850 d = DynamicVector<Real>(n);
851 e = DynamicVector<Real>(n);
852
853 issymmetric = 1;
854 for (int j = 0; (j < n) && issymmetric; j++) {
855 for (int i = 0; (i < n) && issymmetric; i++) {
856 issymmetric = (A(i, j) == A(j, i));
857 }
858 }
859
860 if (issymmetric) {
861 for (int i = 0; i < n; i++) {
862 for (int j = 0; j < n; j++) {
863 V(i, j) = A(i, j);
864 }
865 }
866
867 // Tridiagonalize.
868 tred2();
869
870 // Diagonalize.
871 tql2();
872
873 } else {
874 H = DynamicRectMatrix<Real>(n, n);
875 ort = DynamicVector<Real>(n);
876
877 for (int j = 0; j < n; j++) {
878 for (int i = 0; i < n; i++) {
879 H(i, j) = A(i, j);
880 }
881 }
882
883 // Reduce to Hessenberg form.
884 orthes();
885
886 // Reduce Hessenberg to real Schur form.
887 hqr2();
888 }
889 }
890
891 /** Return the eigenvector matrix
892 @return V
893 */
894 void getV(DynamicRectMatrix<Real>& V_) {
895 V_ = V;
896 return;
897 }
898
899 /** Return the real parts of the eigenvalues
900 @return real(diag(D))
901 */
902 void getRealEigenvalues(DynamicVector<Real>& d_) {
903 d_.reserve(d.size());
904 d_ = d;
905 return;
906 }
907
908 /** Return the imaginary parts of the eigenvalues
909 in parameter e_.
910
911 @param e_: new matrix with imaginary parts of the eigenvalues.
912 */
913 void getImagEigenvalues(DynamicVector<Real>& e_) {
914 e_ = e;
915 return;
916 }
917
918 /**
919 Computes the block diagonal eigenvalue matrix.
920 If the original matrix A is not symmetric, then the eigenvalue
921 matrix D is block diagonal with the real eigenvalues in 1-by-1
922 blocks and any complex eigenvalues,
923 a + i*b, in 2-by-2 blocks, (a, b; -b, a). That is, if the complex
924 eigenvalues look like
925<pre>
926
927 u + iv . . . . .
928 . u - iv . . . .
929 . . a + ib . . .
930 . . . a - ib . .
931 . . . . x .
932 . . . . . y
933</pre>
934 then D looks like
935<pre>
936
937 u v . . . .
938 -v u . . . .
939 . . a b . .
940 . . -b a . .
941 . . . . x .
942 . . . . . y
943</pre>
944 This keeps V a real matrix in both symmetric and non-symmetric
945 cases, and A*V = V*D.
946
947 @param D: upon return, the matrix is filled with the block diagonal
948 eigenvalue matrix.
949
950*/
951 void getD(DynamicRectMatrix<Real>& D) {
952 D = DynamicRectMatrix<Real>(n, n);
953 for (int i = 0; i < n; i++) {
954 for (int j = 0; j < n; j++) {
955 D(i, j) = 0.0;
956 }
957 D(i, i) = d(i);
958 if (e(i) > 0) {
959 D(i, i + 1) = e(i);
960 } else if (e(i) < 0) {
961 D(i, i - 1) = e(i);
962 }
963 }
964 }
965 };
966
967} // namespace JAMA
968
969#endif
970// JAMA_EIG_H
Computes eigenvalues and eigenvectors of a real (non-complex) matrix.
void getD(DynamicRectMatrix< Real > &D)
void getImagEigenvalues(DynamicVector< Real > &e_)
Return the imaginary parts of the eigenvalues in parameter e_.
void getRealEigenvalues(DynamicVector< Real > &d_)
Return the real parts of the eigenvalues.
Eigenvalue(const DynamicRectMatrix< Real > &A)
Check for symmetry, then construct the eigenvalue decomposition.
void getV(DynamicRectMatrix< Real > &V_)
Return the eigenvector matrix.
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.