// matrix/jama-eig.h
// Copyright 2009-2011 Microsoft Corporation
// See ../../COPYING for clarification regarding multiple authors
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
// http://www.apache.org/licenses/LICENSE-2.0
// THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
// KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
// WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
// MERCHANTABLITY OR NON-INFRINGEMENT.
// See the Apache 2 License for the specific language governing permissions and
// limitations under the License.
// This file consists of a port and modification of materials from
// JAMA: A Java Matrix Package
// under the following notice: This software is a cooperative product of
// The MathWorks and the National Institute of Standards and Technology (NIST)
// which has been released to the public. This notice and the original code are
// available at http://math.nist.gov/javanumerics/jama/domain.notice
#ifndef KALDI_MATRIX_JAMA_EIG_H_
#define KALDI_MATRIX_JAMA_EIG_H_ 1
#include "matrix/kaldi-matrix.h"
namespace kaldi {
// This class is not to be used externally. See the Eig function in the Matrix
// class in kaldi-matrix.h. This is the external interface.
template<typename Real> class EigenvalueDecomposition {
// This class is based on the EigenvalueDecomposition class from the JAMA
// library (version 1.0.2).
public:
EigenvalueDecomposition(const MatrixBase<Real> &A);
~EigenvalueDecomposition(); // free memory.
void GetV(MatrixBase<Real> *V_out) { // V is what we call P externally; it's the matrix of
// eigenvectors.
KALDI_ASSERT(V_out->NumRows() == static_cast<MatrixIndexT>(n_)
&& V_out->NumCols() == static_cast<MatrixIndexT>(n_));
for (int i = 0; i < n_; i++)
for (int j = 0; j < n_; j++)
(*V_out)(i, j) = V(i, j); // V(i, j) is member function.
}
void GetRealEigenvalues(VectorBase<Real> *r_out) {
// returns real part of eigenvalues.
KALDI_ASSERT(r_out->Dim() == static_cast<MatrixIndexT>(n_));
for (int i = 0; i < n_; i++)
(*r_out)(i) = d_[i];
}
void GetImagEigenvalues(VectorBase<Real> *i_out) {
// returns imaginary part of eigenvalues.
KALDI_ASSERT(i_out->Dim() == static_cast<MatrixIndexT>(n_));
for (int i = 0; i < n_; i++)
(*i_out)(i) = e_[i];
}
private:
inline Real &H(int r, int c) { return H_[r*n_ + c]; }
inline Real &V(int r, int c) { return V_[r*n_ + c]; }
// complex division
inline static void cdiv(Real xr, Real xi, Real yr, Real yi, Real *cdivr, Real *cdivi) {
Real r, d;
if (std::abs(yr) > std::abs(yi)) {
r = yi/yr;
d = yr + r*yi;
*cdivr = (xr + r*xi)/d;
*cdivi = (xi - r*xr)/d;
} else {
r = yr/yi;
d = yi + r*yr;
*cdivr = (r*xr + xi)/d;
*cdivi = (r*xi - xr)/d;
}
}
// Nonsymmetric reduction from Hessenberg to real Schur form.
void Hqr2 ();
int n_; // matrix dimension.
Real *d_, *e_; // real and imaginary parts of eigenvalues.
Real *V_; // the eigenvectors (P in our external notation)
Real *H_; // the nonsymmetric Hessenberg form.
Real *ort_; // working storage for nonsymmetric algorithm.
// Symmetric Householder reduction to tridiagonal form.
void Tred2 ();
// Symmetric tridiagonal QL algorithm.
void Tql2 ();
// Nonsymmetric reduction to Hessenberg form.
void Orthes ();
};
template class EigenvalueDecomposition<float>; // force instantiation.
template class EigenvalueDecomposition<double>; // force instantiation.
template<typename Real> void EigenvalueDecomposition<Real>::Tred2() {
// This is derived from the Algol procedures tred2 by
// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
// Fortran subroutine in EISPACK.
for (int j = 0; j < n_; j++) {
d_[j] = V(n_-1, j);
}
// Householder reduction to tridiagonal form.
for (int i = n_-1; i > 0; i--) {
// Scale to avoid under/overflow.
Real scale = 0.0;
Real h = 0.0;
for (int k = 0; k < i; k++) {
scale = scale + std::abs(d_[k]);
}
if (scale == 0.0) {
e_[i] = d_[i-1];
for (int j = 0; j < i; j++) {
d_[j] = V(i-1, j);
V(i, j) = 0.0;
V(j, i) = 0.0;
}
} else {
// Generate Householder vector.
for (int k = 0; k < i; k++) {
d_[k] /= scale;
h +=