All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
EigenvalueDecomposition< Real > Class Template Reference

#include <jama-eig.h>

Collaboration diagram for EigenvalueDecomposition< Real >:

Public Member Functions

 EigenvalueDecomposition (const MatrixBase< Real > &A)
 
 ~EigenvalueDecomposition ()
 
void GetV (MatrixBase< Real > *V_out)
 
void GetRealEigenvalues (VectorBase< Real > *r_out)
 
void GetImagEigenvalues (VectorBase< Real > *i_out)
 

Private Member Functions

Real & H (int r, int c)
 
Real & V (int r, int c)
 
void Hqr2 ()
 
void Tred2 ()
 
void Tql2 ()
 
void Orthes ()
 

Static Private Member Functions

static void cdiv (Real xr, Real xi, Real yr, Real yi, Real *cdivr, Real *cdivi)
 

Private Attributes

int n_
 
Real * d_
 
Real * e_
 
Real * V_
 
Real * H_
 
Real * ort_
 

Detailed Description

template<typename Real>
class kaldi::EigenvalueDecomposition< Real >

Definition at line 39 of file jama-eig.h.

Constructor & Destructor Documentation

EigenvalueDecomposition ( const MatrixBase< Real > &  A)

Definition at line 876 of file jama-eig.h.

References rnnlm::i, MatrixBase< Real >::IsSymmetric(), rnnlm::j, KALDI_ASSERT, MatrixBase< Real >::NumCols(), and MatrixBase< Real >::NumRows().

876  {
877  KALDI_ASSERT(A.NumCols() == A.NumRows() && A.NumCols() >= 1);
878  n_ = A.NumRows();
879  V_ = new Real[n_*n_];
880  d_ = new Real[n_];
881  e_ = new Real[n_];
882  H_ = NULL;
883  ort_ = NULL;
884  if (A.IsSymmetric(0.0)) {
885 
886  for (int i = 0; i < n_; i++)
887  for (int j = 0; j < n_; j++)
888  V(i, j) = A(i, j); // Note that V(i, j) is a member function; A(i, j) is an operator
889  // of the matrix A.
890  // Tridiagonalize.
891  Tred2();
892 
893  // Diagonalize.
894  Tql2();
895  } else {
896  H_ = new Real[n_*n_];
897  ort_ = new Real[n_];
898  for (int i = 0; i < n_; i++)
899  for (int j = 0; j < n_; j++)
900  H(i, j) = A(i, j); // as before: H is member function, A(i, j) is operator of matrix.
901 
902  // Reduce to Hessenberg form.
903  Orthes();
904 
905  // Reduce Hessenberg to real Schur form.
906  Hqr2();
907  }
908 }
Real & V(int r, int c)
Definition: jama-eig.h:70
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
Real & H(int r, int c)
Definition: jama-eig.h:69

Definition at line 911 of file jama-eig.h.

911  {
912  delete [] d_;
913  delete [] e_;
914  delete [] V_;
915  delete [] H_;
916  delete [] ort_;
917 }

Member Function Documentation

static void cdiv ( Real  xr,
Real  xi,
Real  yr,
Real  yi,
Real *  cdivr,
Real *  cdivi 
)
inlinestaticprivate

Definition at line 73 of file jama-eig.h.

References rnnlm::d.

73  {
74  Real r, d;
75  if (std::abs(yr) > std::abs(yi)) {
76  r = yi/yr;
77  d = yr + r*yi;
78  *cdivr = (xr + r*xi)/d;
79  *cdivi = (xi - r*xr)/d;
80  } else {
81  r = yr/yi;
82  d = yi + r*yr;
83  *cdivr = (r*xr + xi)/d;
84  *cdivi = (r*xi - xr)/d;
85  }
86  }
void GetImagEigenvalues ( VectorBase< Real > *  i_out)
inline

Definition at line 61 of file jama-eig.h.

References VectorBase< Real >::Dim(), EigenvalueDecomposition< Real >::e_, rnnlm::i, KALDI_ASSERT, and EigenvalueDecomposition< Real >::n_.

Referenced by MatrixBase< Real >::Eig().

61  {
62  // returns imaginary part of eigenvalues.
63  KALDI_ASSERT(i_out->Dim() == static_cast<MatrixIndexT>(n_));
64  for (int i = 0; i < n_; i++)
65  (*i_out)(i) = e_[i];
66  }
int32 MatrixIndexT
Definition: matrix-common.h:96
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
void GetRealEigenvalues ( VectorBase< Real > *  r_out)
inline

Definition at line 55 of file jama-eig.h.

References EigenvalueDecomposition< Real >::d_, VectorBase< Real >::Dim(), rnnlm::i, KALDI_ASSERT, and EigenvalueDecomposition< Real >::n_.

Referenced by MatrixBase< Real >::Eig().

55  {
56  // returns real part of eigenvalues.
57  KALDI_ASSERT(r_out->Dim() == static_cast<MatrixIndexT>(n_));
58  for (int i = 0; i < n_; i++)
59  (*r_out)(i) = d_[i];
60  }
int32 MatrixIndexT
Definition: matrix-common.h:96
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
void GetV ( MatrixBase< Real > *  V_out)
inline

Definition at line 47 of file jama-eig.h.

References rnnlm::i, rnnlm::j, KALDI_ASSERT, EigenvalueDecomposition< Real >::n_, MatrixBase< Real >::NumCols(), MatrixBase< Real >::NumRows(), and EigenvalueDecomposition< Real >::V().

Referenced by MatrixBase< Real >::Eig().

47  { // V is what we call P externally; it's the matrix of
48  // eigenvectors.
49  KALDI_ASSERT(V_out->NumRows() == static_cast<MatrixIndexT>(n_)
50  && V_out->NumCols() == static_cast<MatrixIndexT>(n_));
51  for (int i = 0; i < n_; i++)
52  for (int j = 0; j < n_; j++)
53  (*V_out)(i, j) = V(i, j); // V(i, j) is member function.
54  }
int32 MatrixIndexT
Definition: matrix-common.h:96
Real & V(int r, int c)
Definition: jama-eig.h:70
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
Real& H ( int  r,
int  c 
)
inlineprivate
void Hqr2 ( )
private

Definition at line 436 of file jama-eig.h.

References rnnlm::i, rnnlm::j, and rnnlm::n.

436  {
437  // This is derived from the Algol procedure hqr2,
438  // by Martin and Wilkinson, Handbook for Auto. Comp.,
439  // Vol.ii-Linear Algebra, and the corresponding
440  // Fortran subroutine in EISPACK.
441 
442  int nn = n_;
443  int n = nn-1;
444  int low = 0;
445  int high = nn-1;
446  Real eps = std::numeric_limits<Real>::epsilon();
447  Real exshift = 0.0;
448  Real p = 0, q = 0, r = 0, s = 0, z=0, t, w, x, y;
449 
450  // Store roots isolated by balanc and compute matrix norm
451 
452  Real norm = 0.0;
453  for (int i = 0; i < nn; i++) {
454  if (i < low || i > high) {
455  d_[i] = H(i, i);
456  e_[i] = 0.0;
457  }
458  for (int j = std::max(i-1, 0); j < nn; j++) {
459  norm = norm + std::abs(H(i, j));
460  }
461  }
462 
463  // Outer loop over eigenvalue index
464 
465  int iter = 0;
466  while (n >= low) {
467 
468  // Look for single small sub-diagonal element
469 
470  int l = n;
471  while (l > low) {
472  s = std::abs(H(l-1, l-1)) + std::abs(H(l, l));
473  if (s == 0.0) {
474  s = norm;
475  }
476  if (std::abs(H(l, l-1)) < eps * s) {
477  break;
478  }
479  l--;
480  }
481 
482  // Check for convergence
483  // One root found
484 
485  if (l == n) {
486  H(n, n) = H(n, n) + exshift;
487  d_[n] = H(n, n);
488  e_[n] = 0.0;
489  n--;
490  iter = 0;
491 
492  // Two roots found
493 
494  } else if (l == n-1) {
495  w = H(n, n-1) * H(n-1, n);
496  p = (H(n-1, n-1) - H(n, n)) / 2.0;
497  q = p * p + w;
498  z = std::sqrt(std::abs(q));
499  H(n, n) = H(n, n) + exshift;
500  H(n-1, n-1) = H(n-1, n-1) + exshift;
501  x = H(n, n);
502 
503  // Real pair
504 
505  if (q >= 0) {
506  if (p >= 0) {
507  z = p + z;
508  } else {
509  z = p - z;
510  }
511  d_[n-1] = x + z;
512  d_[n] = d_[n-1];
513  if (z != 0.0) {
514  d_[n] = x - w / z;
515  }
516  e_[n-1] = 0.0;
517  e_[n] = 0.0;
518  x = H(n, n-1);
519  s = std::abs(x) + std::abs(z);
520  p = x / s;
521  q = z / s;
522  r = std::sqrt(p * p+q * q);
523  p = p / r;
524  q = q / r;
525 
526  // Row modification
527 
528  for (int j = n-1; j < nn; j++) {
529  z = H(n-1, j);
530  H(n-1, j) = q * z + p * H(n, j);
531  H(n, j) = q * H(n, j) - p * z;
532  }
533 
534  // Column modification
535 
536  for (int i = 0; i <= n; i++) {
537  z = H(i, n-1);
538  H(i, n-1) = q * z + p * H(i, n);
539  H(i, n) = q * H(i, n) - p * z;
540  }
541 
542  // Accumulate transformations
543 
544  for (int i = low; i <= high; i++) {
545  z = V(i, n-1);
546  V(i, n-1) = q * z + p * V(i, n);
547  V(i, n) = q * V(i, n) - p * z;
548  }
549 
550  // Complex pair
551 
552  } else {
553  d_[n-1] = x + p;
554  d_[n] = x + p;
555  e_[n-1] = z;
556  e_[n] = -z;
557  }
558  n = n - 2;
559  iter = 0;
560 
561  // No convergence yet
562 
563  } else {
564 
565  // Form shift
566 
567  x = H(n, n);
568  y = 0.0;
569  w = 0.0;
570  if (l < n) {
571  y = H(n-1, n-1);
572  w = H(n, n-1) * H(n-1, n);
573  }
574 
575  // Wilkinson's original ad hoc shift
576 
577  if (iter == 10) {
578  exshift += x;
579  for (int i = low; i <= n; i++) {
580  H(i, i) -= x;
581  }
582  s = std::abs(H(n, n-1)) + std::abs(H(n-1, n-2));
583  x = y = 0.75 * s;
584  w = -0.4375 * s * s;
585  }
586 
587  // MATLAB's new ad hoc shift
588 
589  if (iter == 30) {
590  s = (y - x) / 2.0;
591  s = s * s + w;
592  if (s > 0) {
593  s = std::sqrt(s);
594  if (y < x) {
595  s = -s;
596  }
597  s = x - w / ((y - x) / 2.0 + s);
598  for (int i = low; i <= n; i++) {
599  H(i, i) -= s;
600  }
601  exshift += s;
602  x = y = w = 0.964;
603  }
604  }
605 
606  iter = iter + 1; // (Could check iteration count here.)
607 
608  // Look for two consecutive small sub-diagonal elements
609 
610  int m = n-2;
611  while (m >= l) {
612  z = H(m, m);
613  r = x - z;
614  s = y - z;
615  p = (r * s - w) / H(m+1, m) + H(m, m+1);
616  q = H(m+1, m+1) - z - r - s;
617  r = H(m+2, m+1);
618  s = std::abs(p) + std::abs(q) + std::abs(r);
619  p = p / s;
620  q = q / s;
621  r = r / s;
622  if (m == l) {
623  break;
624  }
625  if (std::abs(H(m, m-1)) * (std::abs(q) + std::abs(r)) <
626  eps * (std::abs(p) * (std::abs(H(m-1, m-1)) + std::abs(z) +
627  std::abs(H(m+1, m+1))))) {
628  break;
629  }
630  m--;
631  }
632 
633  for (int i = m+2; i <= n; i++) {
634  H(i, i-2) = 0.0;
635  if (i > m+2) {
636  H(i, i-3) = 0.0;
637  }
638  }
639 
640  // Double QR step involving rows l:n and columns m:n
641 
642  for (int k = m; k <= n-1; k++) {
643  bool notlast = (k != n-1);
644  if (k != m) {
645  p = H(k, k-1);
646  q = H(k+1, k-1);
647  r = (notlast ? H(k+2, k-1) : 0.0);
648  x = std::abs(p) + std::abs(q) + std::abs(r);
649  if (x != 0.0) {
650  p = p / x;
651  q = q / x;
652  r = r / x;
653  }
654  }
655  if (x == 0.0) {
656  break;
657  }
658  s = std::sqrt(p * p + q * q + r * r);
659  if (p < 0) {
660  s = -s;
661  }
662  if (s != 0) {
663  if (k != m) {
664  H(k, k-1) = -s * x;
665  } else if (l != m) {
666  H(k, k-1) = -H(k, k-1);
667  }
668  p = p + s;
669  x = p / s;
670  y = q / s;
671  z = r / s;
672  q = q / p;
673  r = r / p;
674 
675  // Row modification
676 
677  for (int j = k; j < nn; j++) {
678  p = H(k, j) + q * H(k+1, j);
679  if (notlast) {
680  p = p + r * H(k+2, j);
681  H(k+2, j) = H(k+2, j) - p * z;
682  }
683  H(k, j) = H(k, j) - p * x;
684  H(k+1, j) = H(k+1, j) - p * y;
685  }
686 
687  // Column modification
688 
689  for (int i = 0; i <= std::min(n, k+3); i++) {
690  p = x * H(i, k) + y * H(i, k+1);
691  if (notlast) {
692  p = p + z * H(i, k+2);
693  H(i, k+2) = H(i, k+2) - p * r;
694  }
695  H(i, k) = H(i, k) - p;
696  H(i, k+1) = H(i, k+1) - p * q;
697  }
698 
699  // Accumulate transformations
700 
701  for (int i = low; i <= high; i++) {
702  p = x * V(i, k) + y * V(i, k+1);
703  if (notlast) {
704  p = p + z * V(i, k+2);
705  V(i, k+2) = V(i, k+2) - p * r;
706  }
707  V(i, k) = V(i, k) - p;
708  V(i, k+1) = V(i, k+1) - p * q;
709  }
710  } // (s != 0)
711  } // k loop
712  } // check convergence
713  } // while (n >= low)
714 
715  // Backsubstitute to find vectors of upper triangular form
716 
717  if (norm == 0.0) {
718  return;
719  }
720 
721  for (n = nn-1; n >= 0; n--) {
722  p = d_[n];
723  q = e_[n];
724 
725  // Real vector
726 
727  if (q == 0) {
728  int l = n;
729  H(n, n) = 1.0;
730  for (int i = n-1; i >= 0; i--) {
731  w = H(i, i) - p;
732  r = 0.0;
733  for (int j = l; j <= n; j++) {
734  r = r + H(i, j) * H(j, n);
735  }
736  if (e_[i] < 0.0) {
737  z = w;
738  s = r;
739  } else {
740  l = i;
741  if (e_[i] == 0.0) {
742  if (w != 0.0) {
743  H(i, n) = -r / w;
744  } else {
745  H(i, n) = -r / (eps * norm);
746  }
747 
748  // Solve real equations
749 
750  } else {
751  x = H(i, i+1);
752  y = H(i+1, i);
753  q = (d_[i] - p) * (d_[i] - p) +e_[i] *e_[i];
754  t = (x * s - z * r) / q;
755  H(i, n) = t;
756  if (std::abs(x) > std::abs(z)) {
757  H(i+1, n) = (-r - w * t) / x;
758  } else {
759  H(i+1, n) = (-s - y * t) / z;
760  }
761  }
762 
763  // Overflow control
764 
765  t = std::abs(H(i, n));
766  if ((eps * t) * t > 1) {
767  for (int j = i; j <= n; j++) {
768  H(j, n) = H(j, n) / t;
769  }
770  }
771  }
772  }
773 
774  // Complex vector
775 
776  } else if (q < 0) {
777  int l = n-1;
778 
779  // Last vector component imaginary so matrix is triangular
780 
781  if (std::abs(H(n, n-1)) > std::abs(H(n-1, n))) {
782  H(n-1, n-1) = q / H(n, n-1);
783  H(n-1, n) = -(H(n, n) - p) / H(n, n-1);
784  } else {
785  Real cdivr, cdivi;
786  cdiv(0.0, -H(n-1, n), H(n-1, n-1)-p, q, &cdivr, &cdivi);
787  H(n-1, n-1) = cdivr;
788  H(n-1, n) = cdivi;
789  }
790  H(n, n-1) = 0.0;
791  H(n, n) = 1.0;
792  for (int i = n-2; i >= 0; i--) {
793  Real ra, sa, vr, vi;
794  ra = 0.0;
795  sa = 0.0;
796  for (int j = l; j <= n; j++) {
797  ra = ra + H(i, j) * H(j, n-1);
798  sa = sa + H(i, j) * H(j, n);
799  }
800  w = H(i, i) - p;
801 
802  if (e_[i] < 0.0) {
803  z = w;
804  r = ra;
805  s = sa;
806  } else {
807  l = i;
808  if (e_[i] == 0) {
809  Real cdivr, cdivi;
810  cdiv(-ra, -sa, w, q, &cdivr, &cdivi);
811  H(i, n-1) = cdivr;
812  H(i, n) = cdivi;
813  } else {
814  Real cdivr, cdivi;
815  // Solve complex equations
816 
817  x = H(i, i+1);
818  y = H(i+1, i);
819  vr = (d_[i] - p) * (d_[i] - p) +e_[i] *e_[i] - q * q;
820  vi = (d_[i] - p) * 2.0 * q;
821  if (vr == 0.0 && vi == 0.0) {
822  vr = eps * norm * (std::abs(w) + std::abs(q) +
823  std::abs(x) + std::abs(y) + std::abs(z));
824  }
825  cdiv(x*r-z*ra+q*sa, x*s-z*sa-q*ra, vr, vi, &cdivr, &cdivi);
826  H(i, n-1) = cdivr;
827  H(i, n) = cdivi;
828  if (std::abs(x) > (std::abs(z) + std::abs(q))) {
829  H(i+1, n-1) = (-ra - w * H(i, n-1) + q * H(i, n)) / x;
830  H(i+1, n) = (-sa - w * H(i, n) - q * H(i, n-1)) / x;
831  } else {
832  cdiv(-r-y*H(i, n-1), -s-y*H(i, n), z, q, &cdivr, &cdivi);
833  H(i+1, n-1) = cdivr;
834  H(i+1, n) = cdivi;
835  }
836  }
837 
838  // Overflow control
839 
840  t = std::max(std::abs(H(i, n-1)), std::abs(H(i, n)));
841  if ((eps * t) * t > 1) {
842  for (int j = i; j <= n; j++) {
843  H(j, n-1) = H(j, n-1) / t;
844  H(j, n) = H(j, n) / t;
845  }
846  }
847  }
848  }
849  }
850  }
851 
852  // Vectors of isolated roots
853 
854  for (int i = 0; i < nn; i++) {
855  if (i < low || i > high) {
856  for (int j = i; j < nn; j++) {
857  V(i, j) = H(i, j);
858  }
859  }
860  }
861 
862  // Back transformation to get eigenvectors of original matrix
863 
864  for (int j = nn-1; j >= low; j--) {
865  for (int i = low; i <= high; i++) {
866  z = 0.0;
867  for (int k = low; k <= std::min(j, high); k++) {
868  z = z + V(i, k) * H(k, j);
869  }
870  V(i, j) = z;
871  }
872  }
873 }
static void cdiv(Real xr, Real xi, Real yr, Real yi, Real *cdivr, Real *cdivi)
Definition: jama-eig.h:73
struct rnnlm::@11::@12 n
Real & V(int r, int c)
Definition: jama-eig.h:70
Real & H(int r, int c)
Definition: jama-eig.h:69
void Orthes ( )
private

Definition at line 345 of file jama-eig.h.

References rnnlm::i, and rnnlm::j.

345  {
346 
347  // This is derived from the Algol procedures orthes and ortran,
348  // by Martin and Wilkinson, Handbook for Auto. Comp.,
349  // Vol.ii-Linear Algebra, and the corresponding
350  // Fortran subroutines in EISPACK.
351 
352  int low = 0;
353  int high = n_-1;
354 
355  for (int m = low+1; m <= high-1; m++) {
356 
357  // Scale column.
358 
359  Real scale = 0.0;
360  for (int i = m; i <= high; i++) {
361  scale = scale + std::abs(H(i, m-1));
362  }
363  if (scale != 0.0) {
364 
365  // Compute Householder transformation.
366 
367  Real h = 0.0;
368  for (int i = high; i >= m; i--) {
369  ort_[i] = H(i, m-1)/scale;
370  h += ort_[i] * ort_[i];
371  }
372  Real g = std::sqrt(h);
373  if (ort_[m] > 0) {
374  g = -g;
375  }
376  h = h - ort_[m] * g;
377  ort_[m] = ort_[m] - g;
378 
379  // Apply Householder similarity transformation
380  // H = (I-u*u'/h)*H*(I-u*u')/h)
381 
382  for (int j = m; j < n_; j++) {
383  Real f = 0.0;
384  for (int i = high; i >= m; i--) {
385  f += ort_[i]*H(i, j);
386  }
387  f = f/h;
388  for (int i = m; i <= high; i++) {
389  H(i, j) -= f*ort_[i];
390  }
391  }
392 
393  for (int i = 0; i <= high; i++) {
394  Real f = 0.0;
395  for (int j = high; j >= m; j--) {
396  f += ort_[j]*H(i, j);
397  }
398  f = f/h;
399  for (int j = m; j <= high; j++) {
400  H(i, j) -= f*ort_[j];
401  }
402  }
403  ort_[m] = scale*ort_[m];
404  H(m, m-1) = scale*g;
405  }
406  }
407 
408  // Accumulate transformations (Algol's ortran).
409 
410  for (int i = 0; i < n_; i++) {
411  for (int j = 0; j < n_; j++) {
412  V(i, j) = (i == j ? 1.0 : 0.0);
413  }
414  }
415 
416  for (int m = high-1; m >= low+1; m--) {
417  if (H(m, m-1) != 0.0) {
418  for (int i = m+1; i <= high; i++) {
419  ort_[i] = H(i, m-1);
420  }
421  for (int j = m; j <= high; j++) {
422  Real g = 0.0;
423  for (int i = m; i <= high; i++) {
424  g += ort_[i] * V(i, j);
425  }
426  // Double division avoids possible underflow
427  g = (g / ort_[m]) / H(m, m-1);
428  for (int i = m; i <= high; i++) {
429  V(i, j) += g * ort_[i];
430  }
431  }
432  }
433  }
434 }
Real & V(int r, int c)
Definition: jama-eig.h:70
Real & H(int r, int c)
Definition: jama-eig.h:69
void Tql2 ( )
private

Definition at line 227 of file jama-eig.h.

References kaldi::Hypot(), rnnlm::i, and rnnlm::j.

227  {
228  // This is derived from the Algol procedures tql2, by
229  // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
230  // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
231  // Fortran subroutine in EISPACK.
232 
233  for (int i = 1; i < n_; i++) {
234  e_[i-1] = e_[i];
235  }
236  e_[n_-1] = 0.0;
237 
238  Real f = 0.0;
239  Real tst1 = 0.0;
240  Real eps = std::numeric_limits<Real>::epsilon();
241  for (int l = 0; l < n_; l++) {
242 
243  // Find small subdiagonal element
244 
245  tst1 = std::max(tst1, std::abs(d_[l]) + std::abs(e_[l]));
246  int m = l;
247  while (m < n_) {
248  if (std::abs(e_[m]) <= eps*tst1) {
249  break;
250  }
251  m++;
252  }
253 
254  // If m == l, d_[l] is an eigenvalue,
255  // otherwise, iterate.
256 
257  if (m > l) {
258  int iter = 0;
259  do {
260  iter = iter + 1; // (Could check iteration count here.)
261 
262  // Compute implicit shift
263 
264  Real g = d_[l];
265  Real p = (d_[l+1] - g) / (2.0 *e_[l]);
266  Real r = Hypot(p, static_cast<Real>(1.0)); // This is a Kaldi version of hypot that works with templates.
267  if (p < 0) {
268  r = -r;
269  }
270  d_[l] =e_[l] / (p + r);
271  d_[l+1] =e_[l] * (p + r);
272  Real dl1 = d_[l+1];
273  Real h = g - d_[l];
274  for (int i = l+2; i < n_; i++) {
275  d_[i] -= h;
276  }
277  f = f + h;
278 
279  // Implicit QL transformation.
280 
281  p = d_[m];
282  Real c = 1.0;
283  Real c2 = c;
284  Real c3 = c;
285  Real el1 =e_[l+1];
286  Real s = 0.0;
287  Real s2 = 0.0;
288  for (int i = m-1; i >= l; i--) {
289  c3 = c2;
290  c2 = c;
291  s2 = s;
292  g = c *e_[i];
293  h = c * p;
294  r = Hypot(p, e_[i]); // This is a Kaldi version of Hypot that works with templates.
295  e_[i+1] = s * r;
296  s =e_[i] / r;
297  c = p / r;
298  p = c * d_[i] - s * g;
299  d_[i+1] = h + s * (c * g + s * d_[i]);
300 
301  // Accumulate transformation.
302 
303  for (int k = 0; k < n_; k++) {
304  h = V(k, i+1);
305  V(k, i+1) = s * V(k, i) + c * h;
306  V(k, i) = c * V(k, i) - s * h;
307  }
308  }
309  p = -s * s2 * c3 * el1 *e_[l] / dl1;
310  e_[l] = s * p;
311  d_[l] = c * p;
312 
313  // Check for convergence.
314 
315  } while (std::abs(e_[l]) > eps*tst1);
316  }
317  d_[l] = d_[l] + f;
318  e_[l] = 0.0;
319  }
320 
321  // Sort eigenvalues and corresponding vectors.
322 
323  for (int i = 0; i < n_-1; i++) {
324  int k = i;
325  Real p = d_[i];
326  for (int j = i+1; j < n_; j++) {
327  if (d_[j] < p) {
328  k = j;
329  p = d_[j];
330  }
331  }
332  if (k != i) {
333  d_[k] = d_[i];
334  d_[i] = p;
335  for (int j = 0; j < n_; j++) {
336  p = V(j, i);
337  V(j, i) = V(j, k);
338  V(j, k) = p;
339  }
340  }
341  }
342 }
Real & V(int r, int c)
Definition: jama-eig.h:70
double Hypot(double x, double y)
Definition: kaldi-math.h:351
void Tred2 ( )
private

Definition at line 113 of file jama-eig.h.

References rnnlm::i, and rnnlm::j.

113  {
114  // This is derived from the Algol procedures tred2 by
115  // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
116  // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
117  // Fortran subroutine in EISPACK.
118 
119  for (int j = 0; j < n_; j++) {
120  d_[j] = V(n_-1, j);
121  }
122 
123  // Householder reduction to tridiagonal form.
124 
125  for (int i = n_-1; i > 0; i--) {
126 
127  // Scale to avoid under/overflow.
128 
129  Real scale = 0.0;
130  Real h = 0.0;
131  for (int k = 0; k < i; k++) {
132  scale = scale + std::abs(d_[k]);
133  }
134  if (scale == 0.0) {
135  e_[i] = d_[i-1];
136  for (int j = 0; j < i; j++) {
137  d_[j] = V(i-1, j);
138  V(i, j) = 0.0;
139  V(j, i) = 0.0;
140  }
141  } else {
142 
143  // Generate Householder vector.
144 
145  for (int k = 0; k < i; k++) {
146  d_[k] /= scale;
147  h += d_[k] * d_[k];
148  }
149  Real f = d_[i-1];
150  Real g = std::sqrt(h);
151  if (f > 0) {
152  g = -g;
153  }
154  e_[i] = scale * g;
155  h = h - f * g;
156  d_[i-1] = f - g;
157  for (int j = 0; j < i; j++) {
158  e_[j] = 0.0;
159  }
160 
161  // Apply similarity transformation to remaining columns.
162 
163  for (int j = 0; j < i; j++) {
164  f = d_[j];
165  V(j, i) = f;
166  g =e_[j] + V(j, j) * f;
167  for (int k = j+1; k <= i-1; k++) {
168  g += V(k, j) * d_[k];
169  e_[k] += V(k, j) * f;
170  }
171  e_[j] = g;
172  }
173  f = 0.0;
174  for (int j = 0; j < i; j++) {
175  e_[j] /= h;
176  f += e_[j] * d_[j];
177  }
178  Real hh = f / (h + h);
179  for (int j = 0; j < i; j++) {
180  e_[j] -= hh * d_[j];
181  }
182  for (int j = 0; j < i; j++) {
183  f = d_[j];
184  g = e_[j];
185  for (int k = j; k <= i-1; k++) {
186  V(k, j) -= (f * e_[k] + g * d_[k]);
187  }
188  d_[j] = V(i-1, j);
189  V(i, j) = 0.0;
190  }
191  }
192  d_[i] = h;
193  }
194 
195  // Accumulate transformations.
196 
197  for (int i = 0; i < n_-1; i++) {
198  V(n_-1, i) = V(i, i);
199  V(i, i) = 1.0;
200  Real h = d_[i+1];
201  if (h != 0.0) {
202  for (int k = 0; k <= i; k++) {
203  d_[k] = V(k, i+1) / h;
204  }
205  for (int j = 0; j <= i; j++) {
206  Real g = 0.0;
207  for (int k = 0; k <= i; k++) {
208  g += V(k, i+1) * V(k, j);
209  }
210  for (int k = 0; k <= i; k++) {
211  V(k, j) -= g * d_[k];
212  }
213  }
214  }
215  for (int k = 0; k <= i; k++) {
216  V(k, i+1) = 0.0;
217  }
218  }
219  for (int j = 0; j < n_; j++) {
220  d_[j] = V(n_-1, j);
221  V(n_-1, j) = 0.0;
222  }
223  V(n_-1, n_-1) = 1.0;
224  e_[0] = 0.0;
225 }
Real & V(int r, int c)
Definition: jama-eig.h:70
Real& V ( int  r,
int  c 
)
inlineprivate

Member Data Documentation

Real* d_
private

Definition at line 94 of file jama-eig.h.

Referenced by EigenvalueDecomposition< Real >::GetRealEigenvalues().

Real * e_
private

Definition at line 94 of file jama-eig.h.

Referenced by EigenvalueDecomposition< Real >::GetImagEigenvalues().

Real* H_
private

Definition at line 96 of file jama-eig.h.

Referenced by EigenvalueDecomposition< Real >::H().

Real* ort_
private

Definition at line 97 of file jama-eig.h.

Real* V_
private

Definition at line 95 of file jama-eig.h.

Referenced by EigenvalueDecomposition< Real >::V().


The documentation for this class was generated from the following file: