Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Namespace Members | Class Members | File Members | Related Pages

# Tools/Math/LA_Eigenvalues.h

Go to the documentation of this file.

    00001 //------------------------------------------------------------------------------
    00002 /**
    00003  * @file LA_Eigenvalues.h
    00004  * Contains class LinAlg::Eigenvalues
    00005  *
    00006  * @author <a href="mailto:stefanuhrig@gmx.net">Stefan Uhrig</a>
    00007  */
    00008 //------------------------------------------------------------------------------
    00009 #ifndef LINALG_EIGENVALUES_H_INCLUDED
    00010 #define LINALG_EIGENVALUES_H_INCLUDED
    00011 //------------------------------------------------------------------------------
    00012 #include "LA_Vector.h"
    00013 #include "LA_Matrix.h"
    00014 //------------------------------------------------------------------------------
    00015 namespace LA
    00016 {
    00017 //------------------------------------------------------------------------------
    00018 /**
    00019  * Computes eigenvalues and eigenvectors of a real (non-complex)
    00020  * matrix.
    00021  *
    00022  * <P>
    00023  * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
    00024  * diagonal and the eigenvector matrix V is orthogonal. That is,
    00025  * the diagonal values of D are the eigenvalues, and
    00026  * V*V' = I, where I is the identity matrix.  The columns of V
    00027  * represent the eigenvectors in the sense that A*V = V*D.
    00028  *
    00029  * <P>
    00030  * If A is not symmetric, then the eigenvalue matrix D is block diagonal
    00031  * with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
    00032  * a + i*b, in 2-by-2 blocks, [a, b; -b, a].  That is, if the complex
    00033  * eigenvalues look like
    00034  * <pre>
    00035  *
    00036  *        u + iv     .        .          .      .    .
    00037  *          .      u - iv     .          .      .    .
    00038  *          .        .      a + ib       .      .    .
    00039  *          .        .        .        a - ib   .    .
    00040  *          .        .        .          .      x    .
    00041  *          .        .        .          .      .    y
    00042  * </pre>
    00043  * then D looks like
    00044  * <pre>
    00045  *
    00046  *          u        v        .          .      .    .
    00047  *         -v        u        .          .      .    .
    00048  *          .        .        a          b      .    .
    00049  *          .        .       -b          a      .    .
    00050  *          .        .        .          .      x    .
    00051  *          .        .        .          .      .    y
    00052  * </pre>
    00053  * This keeps V a real matrix in both symmetric and non-symmetric
    00054  * cases, and A*V = V*D.
    00055  *
    00056  *
    00057  *
    00058  * <p>
    00059  * The matrix V may be badly
    00060  * conditioned, or even singular, so the validity of the equation
    00061  * A = V*D*inverse(V) depends upon the condition number of V.
    00062  *
    00063  * <p>
    00064  * (Adapted from JAMA, a Java Matrix Library, developed by jointly
    00065  *  by the Mathworks and NIST; see  http://math.nist.gov/javanumerics/jama).
    00066  */
    00067 class Eigenvalues
    00068 {
    00069 public:
    00070   /**
    00071    * Check for symmetry, then construct the eigenvalue decomposition
    00072    * @param A    Square real (non-complex) matrix
    00073    */
    00074   Eigenvalues(const Matrix& A);
    00075
    00076 private:
    00077   Eigenvalues(const Eigenvalues& ev) {}
    00078   Eigenvalues& operator=(const Eigenvalues& ev) { return *this; }
    00079
    00080 public:
    00081   /**
    00082    * Returns a reference to the matrix with the eigenvectors.
    00083    * @return Reference to a matrix with the eigenvectors.
    00084    */
    00085   const Matrix& getV() const;
    00086
    00087   /**
    00088    * Returns a reference to the vector with the real eigenvalues.
    00089    * @return Reference to the vector with the real eigenvalues.
    00090    */
    00091   const Vector& getRealEigenvalues() const;
    00092
    00093   /**
    00094    * Returns a reference to the vector with the real eigenvalues.
    00095    * @return Reference to the vector with the real eigenvalues.
    00096    */
    00097   const Vector& getImagEigenvalues() const;
    00098
    00099   /**
    00100    * Computes the block diagonal eigenvalue matrix.
    00101    * If the original matrix A is not symmetric, then the eigenvalue
    00102    * matrix D is block diagonal with the real eigenvalues in 1-by-1
    00103    * blocks and any complex eigenvalues,
    00104    * a + i*b, in 2-by-2 blocks, [a, b; -b, a].  That is, if the complex
    00105    * eigenvalues look like
    00106    * <pre>
    00107    *
    00108    *      u + iv     .        .          .      .    .
    00109    *        .      u - iv     .          .      .    .
    00110    *        .        .      a + ib       .      .    .
    00111    *        .        .        .        a - ib   .    .
    00112    *        .        .        .          .      x    .
    00113    *        .        .        .          .      .    y
    00114    * </pre>
    00115    * then D looks like
    00116    * <pre>
    00117    *
    00118    *        u        v        .          .      .    .
    00119    *       -v        u        .          .      .    .
    00120    *        .        .        a          b      .    .
    00121    *        .        .       -b          a      .    .
    00122    *        .        .        .          .      x    .
    00123    *        .        .        .          .      .    y
    00124    * </pre>
    00125    * This keeps V a real matrix in both symmetric and non-symmetric
    00126    * cases, and A*V = V*D.
    00127    *
    00128    * @return upon return, the matrix is filled with the block diagonal
    00129    * eigenvalue matrix.
    00130    */
    00131    Matrix getD() const;
    00132
    00133 private:
    00134   // Symmetric Householder reduction to tridiagonal form.
    00135   void tred2();
    00136
    00137   // Symmetric tridiagonal QL algorithm.
    00138   void tql2();
    00139
    00140   // Nonsymmetric reduction to Hessenberg form.
    00141   void orthes();
    00142
    00143   // Complex scalar division.
    00144   void cdiv(double xr, double xi, double yr, double yi);
    00145
    00146   // Nonsymmetric reduction from Hessenberg to real Schur form.
    00147   void hqr2();
    00148
    00149 private:
    00150   /** Row and column dimension (square matrix).  */
    00151   int n;
    00152
    00153   int issymmetric;
    00154
    00155    /** Arrays for internal storage of eigenvalues. */
    00156   Vector d;         /* real part */
    00157   Vector e;         /* img part */
    00158
    00159   /** Array for internal storage of eigenvectors. */
    00160   Matrix V;
    00161
    00162   /** Array for internal storage of nonsymmetric Hessenberg form.
    00163   @serial internal storage of nonsymmetric Hessenberg form.
    00164   */
    00165   Matrix H;
    00166
    00167   /** Working storage for nonsymmetric algorithm.
    00168   @serial working storage for nonsymmetric algorithm.
    00169   */
    00170   Vector ort;
    00171
    00172   double cdivr, cdivi;
    00173 };
    00174 //------------------------------------------------------------------------------
    00175 }
    00176 //------------------------------------------------------------------------------
    00177 #endif
    00178 //------------------------------------------------------------------------------