Tools/Math/LA_Eigenvalues.h

Go to the documentation of this file.

//------------------------------------------------------------------------------
/**
 * @file LA_Eigenvalues.h
 * Contains class LinAlg::Eigenvalues
 *
 * @author <a href="mailto:stefanuhrig@gmx.net">Stefan Uhrig</a>
 */
//------------------------------------------------------------------------------
#ifndef LINALG_EIGENVALUES_H_INCLUDED
#define LINALG_EIGENVALUES_H_INCLUDED
//------------------------------------------------------------------------------
#include "LA_Vector.h"
#include "LA_Matrix.h"
//------------------------------------------------------------------------------
namespace LA
{
//------------------------------------------------------------------------------
/** 
 * Computes eigenvalues and eigenvectors of a real (non-complex)
 * matrix. 
 *
 * <P>
 * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
 * diagonal and the eigenvector matrix V is orthogonal. That is,
 * the diagonal values of D are the eigenvalues, and
 * V*V' = I, where I is the identity matrix.  The columns of V 
 * represent the eigenvectors in the sense that A*V = V*D.
 *  
 * <P>
 * If A is not symmetric, then the eigenvalue matrix D is block diagonal
 * with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
 * a + i*b, in 2-by-2 blocks, [a, b; -b, a].  That is, if the complex
 * eigenvalues look like
 * <pre>
 *
 *        u + iv     .        .          .      .    .
 *          .      u - iv     .          .      .    .
 *          .        .      a + ib       .      .    .
 *          .        .        .        a - ib   .    .
 *          .        .        .          .      x    .
 *          .        .        .          .      .    y
 * </pre>
 * then D looks like
 * <pre>
 *
 *          u        v        .          .      .    .
 *         -v        u        .          .      .    . 
 *          .        .        a          b      .    .
 *          .        .       -b          a      .    .
 *          .        .        .          .      x    .
 *          .        .        .          .      .    y
 * </pre>
 * This keeps V a real matrix in both symmetric and non-symmetric
 * cases, and A*V = V*D.
 *   
 *   
 *   
 * <p>
 * The matrix V may be badly
 * conditioned, or even singular, so the validity of the equation
 * A = V*D*inverse(V) depends upon the condition number of V.
 *
 * <p>
 * (Adapted from JAMA, a Java Matrix Library, developed by jointly 
 *  by the Mathworks and NIST; see  http://math.nist.gov/javanumerics/jama).
 */
class Eigenvalues
{
public:
  /**
   * Check for symmetry, then construct the eigenvalue decomposition
   * @param A    Square real (non-complex) matrix
   */
  Eigenvalues(const Matrix& A);

private:
  Eigenvalues(const Eigenvalues& ev) {}
  Eigenvalues& operator=(const Eigenvalues& ev) { return *this; }

public:
  /**
   * Returns a reference to the matrix with the eigenvectors.
   * @return Reference to a matrix with the eigenvectors.
   */
  const Matrix& getV() const;
  
  /**
   * Returns a reference to the vector with the real eigenvalues.
   * @return Reference to the vector with the real eigenvalues.
   */
  const Vector& getRealEigenvalues() const;
  
  /**
   * Returns a reference to the vector with the real eigenvalues.
   * @return Reference to the vector with the real eigenvalues.
   */
  const Vector& getImagEigenvalues() const;

  /** 
   * Computes the block diagonal eigenvalue matrix.
   * If the original matrix A is not symmetric, then the eigenvalue 
   * matrix D is block diagonal with the real eigenvalues in 1-by-1 
   * blocks and any complex eigenvalues,
   * a + i*b, in 2-by-2 blocks, [a, b; -b, a].  That is, if the complex
   * eigenvalues look like
   * <pre>
   *
   *      u + iv     .        .          .      .    .
   *        .      u - iv     .          .      .    .
   *        .        .      a + ib       .      .    .
   *        .        .        .        a - ib   .    .
   *        .        .        .          .      x    .
   *        .        .        .          .      .    y
   * </pre>
   * then D looks like
   * <pre>
   *
   *        u        v        .          .      .    .
   *       -v        u        .          .      .    . 
   *        .        .        a          b      .    .
   *        .        .       -b          a      .    .
   *        .        .        .          .      x    .
   *        .        .        .          .      .    y
   * </pre>
   * This keeps V a real matrix in both symmetric and non-symmetric
   * cases, and A*V = V*D.
   *
   * @return upon return, the matrix is filled with the block diagonal 
   * eigenvalue matrix.
   */
   Matrix getD() const;

private:
  // Symmetric Householder reduction to tridiagonal form.
  void tred2();

  // Symmetric tridiagonal QL algorithm.
  void tql2();

  // Nonsymmetric reduction to Hessenberg form.
  void orthes();

  // Complex scalar division.
  void cdiv(double xr, double xi, double yr, double yi);

  // Nonsymmetric reduction from Hessenberg to real Schur form.
  void hqr2();

private:
  /** Row and column dimension (square matrix).  */
  int n;
  
  int issymmetric;
  
   /** Arrays for internal storage of eigenvalues. */
  Vector d;         /* real part */
  Vector e;         /* img part */

  /** Array for internal storage of eigenvectors. */
  Matrix V;

  /** Array for internal storage of nonsymmetric Hessenberg form.
  @serial internal storage of nonsymmetric Hessenberg form.
  */
  Matrix H;

  /** Working storage for nonsymmetric algorithm.
  @serial working storage for nonsymmetric algorithm.
  */
  Vector ort;

  double cdivr, cdivi;
};
//------------------------------------------------------------------------------
}
//------------------------------------------------------------------------------
#endif
//------------------------------------------------------------------------------

---