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