src/math/Eigenvalue.hpp

#ifndef JAMA_EIG_H
#define JAMA_EIG_H

#include "math/DynamicRectMatrix.hpp"

#include <algorithm>
// for min(), max() below
#include <cmath>
// for abs() below

using namespace OpenMD;
using namespace std;

namespace JAMA
{
  
  /** 
      
  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).
  **/
  
  template <class Real>
  class Eigenvalue
  {
    
    
    /** Row and column dimension (square matrix).  */
    int n;
    
    int issymmetric; /* boolean*/
    
    /** Arrays for internal storage of eigenvalues. */
    
    DynamicVector<Real> d;         /* real part */
    DynamicVector<Real> e;         /* img part */
    
    /** Array for internal storage of eigenvectors. */
    DynamicRectMatrix<Real> V;
    
    /** Array for internal storage of nonsymmetric Hessenberg form.
        @serial internal storage of nonsymmetric Hessenberg form.
    */
    DynamicRectMatrix<Real> H;
    
    
    /** Working storage for nonsymmetric algorithm.
        @serial working storage for nonsymmetric algorithm.
    */
    DynamicVector<Real> ort;
    
    
    // Symmetric Householder reduction to tridiagonal form.
    
    void 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 + 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 += d(k) * d(k);
          }
          Real f = d(i-1);
          Real g = sqrt(h);
          if (f > 0) {
            g = -g;
          }
          e(i) = scale * g;
          h = h - f * g;
          d(i-1) = f - g;
          for (int j = 0; j < i; j++) {
            e(j) = 0.0;
          }
          
          // Apply similarity transformation to remaining columns.
          
          for (int j = 0; j < i; j++) {
            f = d(j);
            V(j,i) = f;
            g = e(j) + V(j,j) * f;
            for (int k = j+1; k <= i-1; k++) {
              g += V(k,j) * d(k);
              e(k) += V(k,j) * f;
            }
            e(j) = g;
          }
          f = 0.0;
          for (int j = 0; j < i; j++) {
            e(j) /= h;
            f += e(j) * d(j);
          }
          Real hh = f / (h + h);
          for (int j = 0; j < i; j++) {
            e(j) -= hh * d(j);
          }
          for (int j = 0; j < i; j++) {
            f = d(j);
            g = e(j);
            for (int k = j; k <= i-1; k++) {
              V(k,j) -= (f * e(k) + g * d(k));
            }
            d(j) = V(i-1,j);
            V(i,j) = 0.0;
          }
        }
        d(i) = h;
      }
      
      // Accumulate transformations.
      
      for (int i = 0; i < n-1; i++) {
        V(n-1,i) = V(i,i);
        V(i,i) = 1.0;
        Real h = d(i+1);
        if (h != 0.0) {
          for (int k = 0; k <= i; k++) {
            d(k) = V(k,i+1) / h;
          }
          for (int j = 0; j <= i; j++) {
            Real g = 0.0;
            for (int k = 0; k <= i; k++) {
              g += V(k,i+1) * V(k,j);
            }
            for (int k = 0; k <= i; k++) {
              V(k,j) -= g * d(k);
            }
          }
        }
        for (int k = 0; k <= i; k++) {
          V(k,i+1) = 0.0;
        }
      }
      for (int j = 0; j < n; j++) {
        d(j) = V(n-1,j);
        V(n-1,j) = 0.0;
      }
      V(n-1,n-1) = 1.0;
      e(0) = 0.0;
    } 
    
    // Symmetric tridiagonal QL algorithm.
    
    void tql2 () {
      
      //  This is derived from the Algol procedures tql2, by
      //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
      //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
      //  Fortran subroutine in EISPACK.
      
      for (int i = 1; i < n; i++) {
        e(i-1) = e(i);
      }
      e(n-1) = 0.0;
      
      Real f = 0.0;
      Real tst1 = 0.0;
      Real eps = pow(2.0,-52.0);
      for (int l = 0; l < n; l++) {
        
        // Find small subdiagonal element
        
        tst1 = max(tst1,abs(d(l)) + abs(e(l)));
        int m = l;
        
        // Original while-loop from Java code
        while (m < n) {
          if (abs(e(m)) <= eps*tst1) {
            break;
          }
          m++;
        }
        
        
        // If m == l, d(l) is an eigenvalue,
        // otherwise, iterate.
        
        if (m > l) {
          int iter = 0;
          do {
            iter = iter + 1;  // (Could check iteration count here.)
            
            // Compute implicit shift
            
            Real g = d(l);
            Real p = (d(l+1) - g) / (2.0 * e(l));
            Real r = hypot(p,1.0);
            if (p < 0) {
              r = -r;
            }
            d(l) = e(l) / (p + r);
            d(l+1) = e(l) * (p + r);
            Real dl1 = d(l+1);
            Real h = g - d(l);
            for (int i = l+2; i < n; i++) {
              d(i) -= h;
            }
            f = f + h;
            
            // Implicit QL transformation.
            
            p = d(m);
            Real c = 1.0;
            Real c2 = c;
            Real c3 = c;
            Real el1 = e(l+1);
            Real s = 0.0;
            Real s2 = 0.0;
            for (int i = m-1; i >= l; i--) {
              c3 = c2;
              c2 = c;
              s2 = s;
              g = c * e(i);
              h = c * p;
              r = hypot(p,e(i));
              e(i+1) = s * r;
              s = e(i) / r;
              c = p / r;
              p = c * d(i) - s * g;
              d(i+1) = h + s * (c * g + s * d(i));
              
              // Accumulate transformation.
              
              for (int k = 0; k < n; k++) {
                h = V(k,i+1);
                V(k,i+1) = s * V(k,i) + c * h;
                V(k,i) = c * V(k,i) - s * h;
              }
            }
            p = -s * s2 * c3 * el1 * e(l) / dl1;
            e(l) = s * p;
            d(l) = c * p;
            
            // Check for convergence.
            
          } while (abs(e(l)) > eps*tst1);
        }
        d(l) = d(l) + f;
        e(l) = 0.0;
      }
      
      // Sort eigenvalues and corresponding vectors.
      
      for (int i = 0; i < n-1; i++) {
        int k = i;
        Real p = d(i);
        for (int j = i+1; j < n; j++) {
          if (d(j) < p) {
            k = j;
            p = d(j);
          }
        }
        if (k != i) {
          d(k) = d(i);
          d(i) = p;
          for (int j = 0; j < n; j++) {
            p = V(j,i);
            V(j,i) = V(j,k);
            V(j,k) = p;
          }
        }
      }
    }
    
    // Nonsymmetric reduction to Hessenberg form.
    
    void orthes () {
      
      //  This is derived from the Algol procedures orthes and ortran,
      //  by Martin and Wilkinson, Handbook for Auto. Comp.,
      //  Vol.ii-Linear Algebra, and the corresponding
      //  Fortran subroutines in EISPACK.
   
      int low = 0;
      int high = n-1;
      
      for (int m = low+1; m <= high-1; m++) {
        
        // Scale column.
        
        Real scale = 0.0;
        for (int i = m; i <= high; i++) {
          scale = scale + abs(H(i,m-1));
        }
        if (scale != 0.0) {
          
          // Compute Householder transformation.
          
          Real h = 0.0;
          for (int i = high; i >= m; i--) {
            ort(i) = H(i,m-1)/scale;
            h += ort(i) * ort(i);
          }
          Real g = sqrt(h);
          if (ort(m) > 0) {
            g = -g;
          }
          h = h - ort(m) * g;
          ort(m) = ort(m) - g;
          
          // Apply Householder similarity transformation
          // H = (I-u*u'/h)*H*(I-u*u')/h)
          
          for (int j = m; j < n; j++) {
            Real f = 0.0;
            for (int i = high; i >= m; i--) {
              f += ort(i)*H(i,j);
            }
            f = f/h;
            for (int i = m; i <= high; i++) {
              H(i,j) -= f*ort(i);
            }
          }
          
          for (int i = 0; i <= high; i++) {
            Real f = 0.0;
            for (int j = high; j >= m; j--) {
              f += ort(j)*H(i,j);
            }
            f = f/h;
            for (int j = m; j <= high; j++) {
              H(i,j) -= f*ort(j);
            }
          }
          ort(m) = scale*ort(m);
          H(m,m-1) = scale*g;
        }
      }
      
      // Accumulate transformations (Algol's ortran).
      
      for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
          V(i,j) = (i == j ? 1.0 : 0.0);
        }
      }
      
      for (int m = high-1; m >= low+1; m--) {
        if (H(m,m-1) != 0.0) {
          for (int i = m+1; i <= high; i++) {
            ort(i) = H(i,m-1);
          }
          for (int j = m; j <= high; j++) {
            Real g = 0.0;
            for (int i = m; i <= high; i++) {
              g += ort(i) * V(i,j);
            }
            // Double division avoids possible underflow
            g = (g / ort(m)) / H(m,m-1);
            for (int i = m; i <= high; i++) {
              V(i,j) += g * ort(i);
            }
          }
        }
      }
    }
    
    
    // Complex scalar division.
    
    Real cdivr, cdivi;
    void cdiv(Real xr, Real xi, Real yr, Real yi) {
      Real r,d;
      if (abs(yr) > 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 () {
      
      //  This is derived from the Algol procedure hqr2,
      //  by Martin and Wilkinson, Handbook for Auto. Comp.,
      //  Vol.ii-Linear Algebra, and the corresponding
      //  Fortran subroutine in EISPACK.
      
      // Initialize
      
      int nn = this->n;
      int n = nn-1;
      int low = 0;
      int high = nn-1;
      Real eps = pow(2.0,-52.0);
      Real exshift = 0.0;
      Real p=0,q=0,r=0,s=0,z=0,t,w,x,y;
      
      // Store roots isolated by balanc and compute matrix norm
      
      Real norm = 0.0;
      for (int i = 0; i < nn; i++) {
        if ((i < low) || (i > high)) {
          d(i) = H(i,i);
          e(i) = 0.0;
        }
        for (int j = max(i-1,0); j < nn; j++) {
          norm = norm + abs(H(i,j));
        }
      }
      
      // Outer loop over eigenvalue index
   
      int iter = 0;
      while (n >= low) {
        
        // Look for single small sub-diagonal element
        
        int l = n;
        while (l > low) {
          s = abs(H(l-1,l-1)) + abs(H(l,l));
          if (s == 0.0) {
            s = norm;
          }
          if (abs(H(l,l-1)) < eps * s) {
            break;
          }
          l--;
        }
        
        // Check for convergence
        // One root found
        
        if (l == n) {
          H(n,n) = H(n,n) + exshift;
          d(n) = H(n,n);
          e(n) = 0.0;
          n--;
          iter = 0;
          
          // Two roots found
          
        } else if (l == n-1) {
          w = H(n,n-1) * H(n-1,n);
          p = (H(n-1,n-1) - H(n,n)) / 2.0;
          q = p * p + w;
          z = sqrt(abs(q));
          H(n,n) = H(n,n) + exshift;
          H(n-1,n-1) = H(n-1,n-1) + exshift;
          x = H(n,n);
          
          // Real pair
          
          if (q >= 0) {
            if (p >= 0) {
              z = p + z;
            } else {
              z = p - z;
            }
            d(n-1) = x + z;
            d(n) = d(n-1);
            if (z != 0.0) {
              d(n) = x - w / z;
            }
            e(n-1) = 0.0;
            e(n) = 0.0;
            x = H(n,n-1);
            s = abs(x) + abs(z);
            p = x / s;
            q = z / s;
            r = sqrt(p * p+q * q);
            p = p / r;
            q = q / r;
            
            // Row modification
            
            for (int j = n-1; j < nn; j++) {
              z = H(n-1,j);
              H(n-1,j) = q * z + p * H(n,j);
              H(n,j) = q * H(n,j) - p * z;
            }
            
            // Column modification
            
            for (int i = 0; i <= n; i++) {
              z = H(i,n-1);
              H(i,n-1) = q * z + p * H(i,n);
              H(i,n) = q * H(i,n) - p * z;
            }
            
            // Accumulate transformations
            
            for (int i = low; i <= high; i++) {
              z = V(i,n-1);
              V(i,n-1) = q * z + p * V(i,n);
              V(i,n) = q * V(i,n) - p * z;
            }
            
            // Complex pair
            
          } else {
            d(n-1) = x + p;
            d(n) = x + p;
            e(n-1) = z;
            e(n) = -z;
          }
          n = n - 2;
          iter = 0;
          
          // No convergence yet
          
        } else {
          
          // Form shift
          
          x = H(n,n);
          y = 0.0;
          w = 0.0;
          if (l < n) {
            y = H(n-1,n-1);
            w = H(n,n-1) * H(n-1,n);
          }
          
          // Wilkinson's original ad hoc shift
          
          if (iter == 10) {
            exshift += x;
            for (int i = low; i <= n; i++) {
              H(i,i) -= x;
            }
            s = abs(H(n,n-1)) + abs(H(n-1,n-2));
            x = y = 0.75 * s;
            w = -0.4375 * s * s;
          }
          
          // MATLAB's new ad hoc shift
          
          if (iter == 30) {
            s = (y - x) / 2.0;
            s = s * s + w;
            if (s > 0) {
              s = sqrt(s);
              if (y < x) {
                s = -s;
              }
              s = x - w / ((y - x) / 2.0 + s);
              for (int i = low; i <= n; i++) {
                H(i,i) -= s;
              }
              exshift += s;
              x = y = w = 0.964;
            }
          }
          
          iter = iter + 1;   // (Could check iteration count here.)
          
          // Look for two consecutive small sub-diagonal elements
          
          int m = n-2;
          while (m >= l) {
            z = H(m,m);
            r = x - z;
            s = y - z;
            p = (r * s - w) / H(m+1,m) + H(m,m+1);
            q = H(m+1,m+1) - z - r - s;
            r = H(m+2,m+1);
            s = abs(p) + abs(q) + abs(r);
            p = p / s;
            q = q / s;
            r = r / s;
            if (m == l) {
              break;
            }
            if (abs(H(m,m-1)) * (abs(q) + abs(r)) <
                eps * (abs(p) * (abs(H(m-1,m-1)) + abs(z) +
                                 abs(H(m+1,m+1))))) {
              break;
            }
            m--;
          }
          
          for (int i = m+2; i <= n; i++) {
            H(i,i-2) = 0.0;
            if (i > m+2) {
              H(i,i-3) = 0.0;
            }
          }
          
          // Double QR step involving rows l:n and columns m:n
          
          for (int k = m; k <= n-1; k++) {
            int notlast = (k != n-1);
            if (k != m) {
              p = H(k,k-1);
              q = H(k+1,k-1);
              r = (notlast ? H(k+2,k-1) : 0.0);
              x = abs(p) + abs(q) + abs(r);
              if (x != 0.0) {
                p = p / x;
                q = q / x;
                r = r / x;
              }
            }
            if (x == 0.0) {
              break;
            }
            s = sqrt(p * p + q * q + r * r);
            if (p < 0) {
              s = -s;
            }
            if (s != 0) {
              if (k != m) {
                H(k,k-1) = -s * x;
              } else if (l != m) {
                H(k,k-1) = -H(k,k-1);
              }
              p = p + s;
              x = p / s;
              y = q / s;
              z = r / s;
              q = q / p;
              r = r / p;
              
              // Row modification
              
              for (int j = k; j < nn; j++) {
                p = H(k,j) + q * H(k+1,j);
                if (notlast) {
                  p = p + r * H(k+2,j);
                  H(k+2,j) = H(k+2,j) - p * z;
                }
                H(k,j) = H(k,j) - p * x;
                H(k+1,j) = H(k+1,j) - p * y;
              }
              
              // Column modification
              
              for (int i = 0; i <= min(n,k+3); i++) {
                p = x * H(i,k) + y * H(i,k+1);
                if (notlast) {
                  p = p + z * H(i,k+2);
                  H(i,k+2) = H(i,k+2) - p * r;
                }
                H(i,k) = H(i,k) - p;
                H(i,k+1) = H(i,k+1) - p * q;
              }
              
              // Accumulate transformations
              
              for (int i = low; i <= high; i++) {
                p = x * V(i,k) + y * V(i,k+1);
                if (notlast) {
                  p = p + z * V(i,k+2);
                  V(i,k+2) = V(i,k+2) - p * r;
                }
                V(i,k) = V(i,k) - p;
                V(i,k+1) = V(i,k+1) - p * q;
              }
            }  // (s != 0)
          }  // k loop
        }  // check convergence
      }  // while (n >= low)
      
      // Backsubstitute to find vectors of upper triangular form
      
      if (norm == 0.0) {
        return;
      }
      
      for (n = nn-1; n >= 0; n--) {
        p = d(n);
        q = e(n);
        
        // Real vector
        
        if (q == 0) {
          int l = n;
          H(n,n) = 1.0;
          for (int i = n-1; i >= 0; i--) {
            w = H(i,i) - p;
            r = 0.0;
            for (int j = l; j <= n; j++) {
              r = r + H(i,j) * H(j,n);
            }
            if (e(i) < 0.0) {
              z = w;
              s = r;
            } else {
              l = i;
              if (e(i) == 0.0) {
                if (w != 0.0) {
                  H(i,n) = -r / w;
                } else {
                  H(i,n) = -r / (eps * norm);
                }
                
                // Solve real equations
                
              } else {
                x = H(i,i+1);
                y = H(i+1,i);
                q = (d(i) - p) * (d(i) - p) + e(i) * e(i);
                t = (x * s - z * r) / q;
                H(i,n) = t;
                if (abs(x) > abs(z)) {
                  H(i+1,n) = (-r - w * t) / x;
                } else {
                  H(i+1,n) = (-s - y * t) / z;
                }
              }
              
              // Overflow control
              
              t = abs(H(i,n));
              if ((eps * t) * t > 1) {
                for (int j = i; j <= n; j++) {
                  H(j,n) = H(j,n) / t;
                }
              }
            }
          }
          
          // Complex vector
          
        } else if (q < 0) {
          int l = n-1;
          
          // Last vector component imaginary so matrix is triangular
          
          if (abs(H(n,n-1)) > abs(H(n-1,n))) {
            H(n-1,n-1) = q / H(n,n-1);
            H(n-1,n) = -(H(n,n) - p) / H(n,n-1);
          } else {
            cdiv(0.0,-H(n-1,n),H(n-1,n-1)-p,q);
            H(n-1,n-1) = cdivr;
            H(n-1,n) = cdivi;
          }
          H(n,n-1) = 0.0;
          H(n,n) = 1.0;
          for (int i = n-2; i >= 0; i--) {
            Real ra,sa,vr,vi;
            ra = 0.0;
            sa = 0.0;
            for (int j = l; j <= n; j++) {
              ra = ra + H(i,j) * H(j,n-1);
              sa = sa + H(i,j) * H(j,n);
            }
            w = H(i,i) - p;
            
            if (e(i) < 0.0) {
              z = w;
              r = ra;
              s = sa;
            } else {
              l = i;
              if (e(i) == 0) {
                cdiv(-ra,-sa,w,q);
                H(i,n-1) = cdivr;
                H(i,n) = cdivi;
              } else {
                
                // Solve complex equations
                
                x = H(i,i+1);
                y = H(i+1,i);
                vr = (d(i) - p) * (d(i) - p) + e(i) * e(i) - q * q;
                vi = (d(i) - p) * 2.0 * q;
                if ((vr == 0.0) && (vi == 0.0)) {
                  vr = eps * norm * (abs(w) + abs(q) +
                                     abs(x) + abs(y) + abs(z));
                }
                cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi);
                H(i,n-1) = cdivr;
                H(i,n) = cdivi;
                if (abs(x) > (abs(z) + abs(q))) {
                  H(i+1,n-1) = (-ra - w * H(i,n-1) + q * H(i,n)) / x;
                  H(i+1,n) = (-sa - w * H(i,n) - q * H(i,n-1)) / x;
                } else {
                  cdiv(-r-y*H(i,n-1),-s-y*H(i,n),z,q);
                  H(i+1,n-1) = cdivr;
                  H(i+1,n) = cdivi;
                }
              }
              
              // Overflow control
              
              t = max(abs(H(i,n-1)),abs(H(i,n)));
              if ((eps * t) * t > 1) {
                for (int j = i; j <= n; j++) {
                  H(j,n-1) = H(j,n-1) / t;
                  H(j,n) = H(j,n) / t;
                }
              }
            }
          }
        }
      }
      
      // Vectors of isolated roots
      
      for (int i = 0; i < nn; i++) {
        if (i < low || i > high) {
          for (int j = i; j < nn; j++) {
            V(i,j) = H(i,j);
          }
        }
      }
      
      // Back transformation to get eigenvectors of original matrix
      
      for (int j = nn-1; j >= low; j--) {
        for (int i = low; i <= high; i++) {
          z = 0.0;
          for (int k = low; k <= min(j,high); k++) {
            z = z + V(i,k) * H(k,j);
            }
          V(i,j) = z;
        }
      }
    }
    
  public:

    
    /** Check for symmetry, then construct the eigenvalue decomposition
        @param A    Square real (non-complex) matrix
    */    
    Eigenvalue(const DynamicRectMatrix<Real> &A) {
      n = A.getNCol();
      V = DynamicRectMatrix<Real>(n,n);
      d = DynamicVector<Real>(n);
      e = DynamicVector<Real>(n);
      
      issymmetric = 1;
      for (int j = 0; (j < n) && issymmetric; j++) {
        for (int i = 0; (i < n) && issymmetric; i++) {
          issymmetric = (A(i,j) == A(j,i));
        }
      }
      
      if (issymmetric) {
        for (int i = 0; i < n; i++) {
          for (int j = 0; j < n; j++) {
            V(i,j) = A(i,j);
          }
        }
        
        // Tridiagonalize.
        tred2();
        
        // Diagonalize.
        tql2();
        
      } else {
        H = DynamicRectMatrix<Real>(n,n);
        ort = DynamicVector<Real>(n);
        
        for (int j = 0; j < n; j++) {
          for (int i = 0; i < n; i++) {
            H(i,j) = A(i,j);
          }
        }
        
        // Reduce to Hessenberg form.
        orthes();
        
        // Reduce Hessenberg to real Schur form.
        hqr2();
      }
    }
    
    
    /** Return the eigenvector matrix
        @return     V
    */    
    void getV (DynamicRectMatrix<Real> &V_) {
      V_ = V;
      return;
    }
    
    /** Return the real parts of the eigenvalues
        @return     real(diag(D))
    */    
    void getRealEigenvalues (DynamicVector<Real> &d_) {
      d_ = d;
      return ;
    }
    
    /** Return the imaginary parts of the eigenvalues
        in parameter e_.
        
        @pararm e_: new matrix with imaginary parts of the eigenvalues.
    */
    void getImagEigenvalues (DynamicVector<Real> &e_) {
      e_ = e;
      return;
    }
    
   
    /** 
        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.

        @param D: upon return, the matrix is filled with the block diagonal 
        eigenvalue matrix.
        
*/
    void getD (DynamicRectMatrix<Real> &D) {
      D = DynamicRectMatrix<Real>(n,n);
      for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
          D(i,j) = 0.0;
        }
        D(i,i) = d(i);
        if (e(i) > 0) {
          D(i,i+1) = e(i);
        } else if (e(i) < 0) {
          D(i,i-1) = e(i);
        }
      }
    }
  };
  
} //namespace JAMA


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