* lib/matrix: Add Eigenvalue Decomposition

git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@32351 b2dd03c8-39d4-4d8f-98ff-823fe69b080e

2 files changed, 908 insertions, 0 deletions

diff --git a/lib/matrix.rb b/lib/matrix.rb
index 12f99925bc..8c6b3d6e87 100644
--- a/lib/matrix.rb
+++ b/lib/matrix.rb
@@ -95,6 +95,10 @@ end
 # * <tt> #transpose                     </tt>
 # * <tt> #t                             </tt>
 #
+# Matrix decompositions:
+# * <tt> #eigen                         </tt>
+# * <tt> #eigensystem                   </tt>
+#
 # Complex arithmetic:
 # * <tt> conj                           </tt>
 # * <tt> conjugate                      </tt>
@@ -117,6 +121,7 @@ end
 class Matrix
   include Enumerable
   include ExceptionForMatrix
+  autoload :EigenvalueDecomposition, "matrix/eigenvalue_decomposition"
 
   # instance creations
   private_class_method :new
@@ -1163,6 +1168,23 @@ class Matrix
   alias t transpose
 
   #--
+  # DECOMPOSITIONS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
+  #++
+
+  #
+  # Returns the Eigensystem of the matrix; see +EigenvalueDecomposition+.
+  #   m = Matrix[[1, 2], [3, 4]]
+  #   v, d, v_inv = m.eigensystem
+  #   d.diagonal? # => true
+  #   v.inv == v_inv # => true
+  #   (v * d * v_inv).round(5) == m # => true
+  #
+  def eigensystem
+    EigenvalueDecomposition.new(self)
+  end
+  alias eigen eigensystem
+
+  #--
   # COMPLEX ARITHMETIC -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
   #++
 
diff --git a/lib/matrix/eigenvalue_decomposition.rb b/lib/matrix/eigenvalue_decomposition.rb
new file mode 100644
index 0000000000..55e41aaf7c
--- /dev/null
+++ b/lib/matrix/eigenvalue_decomposition.rb
@@ -0,0 +1,886 @@
+class Matrix
+  # Adapted from JAMA: http://math.nist.gov/javanumerics/jama/
+
+  # Eigenvalues and eigenvectors of a real matrix.
+  #
+  # Computes the eigenvalues and eigenvectors of a matrix A.
+  #
+  # If A is diagonalizable, this provides matrices V and D
+  # such that A = V*D*V.inv, where D is the diagonal matrix with entries
+  # equal to the eigenvalues and V is formed by the eigenvectors.
+  #
+  # If A is symmetric, then V is orthogonal and thus A = V*D*V.t
+
+  class EigenvalueDecomposition
+
+    # Constructs the eigenvalue decomposition for a square matrix +A+
+    #
+    def initialize(a)
+      # @d, @e: Arrays for internal storage of eigenvalues.
+      # @v: Array for internal storage of eigenvectors.
+      # @h: Array for internal storage of nonsymmetric Hessenberg form.
+      raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix)
+      @size = a.row_size
+      @d = Array.new(@size, 0)
+      @e = Array.new(@size, 0)
+
+      if (@symmetric = a.symmetric?)
+        @v = a.to_a
+        tridiagonalize
+        diagonalize
+      else
+        @v = Array.new(@size) { Array.new(@size, 0) }
+        @h = a.to_a
+        @ort = Array.new(@size, 0)
+        reduce_to_hessenberg
+        hessenberg_to_real_schur
+      end
+    end
+
+    # Returns the eigenvector matrix +V+
+    #
+    def eigenvector_matrix
+      Matrix.send :new, build_eigenvectors.transpose
+    end
+    alias v eigenvector_matrix
+
+    # Returns the inverse of the eigenvector matrix +V+
+    #
+    def eigenvector_matrix_inv
+      r = Matrix.send :new, build_eigenvectors
+      r = r.transpose.inverse unless @symmetric
+      r
+    end
+    alias v_inv eigenvector_matrix_inv
+
+    # Returns the eigenvalues in an array
+    #
+    def eigenvalues
+      values = @d.dup
+      @e.each_with_index{|imag, i| values[i] = Complex(values[i], imag) unless imag == 0}
+      values
+    end
+
+    # Returns an array of the eigenvectors
+    #
+    def eigenvectors
+      build_eigenvectors.map{|ev| Vector.send :new, ev}
+    end
+
+    # Returns the block diagonal eigenvalue matrix +D+
+    #
+    def eigenvalue_matrix
+      Matrix.diagonal(*eigenvalues)
+    end
+    alias d eigenvalue_matrix
+
+    # Returns [eigenvector_matrix, eigenvalue_matrix, eigenvector_matrix_inv]
+    #
+    def to_ary
+      [v, d, v_inv]
+    end
+    alias_method :to_a, :to_ary
+
+  private
+    def build_eigenvectors
+      # JAMA stores complex eigenvectors in a strange way
+      # See http://cio.nist.gov/esd/emaildir/lists/jama/msg01021.html
+      @e.each_with_index.map do |imag, i|
+        if imag == 0
+          Array.new(@size){|j| @v[j][i]}
+        elsif imag > 0
+          Array.new(@size){|j| Complex(@v[j][i], @v[j][i+1])}
+        else
+          Array.new(@size){|j| Complex(@v[j][i], -@v[j][i-1])}
+        end
+      end
+    end
+    # Complex scalar division.
+
+    def cdiv(xr, xi, yr, yi)
+      if (yr.abs > yi.abs)
+        r = yi/yr
+        d = yr + r*yi
+        [(xr + r*xi)/d, (xi - r*xr)/d]
+      else
+        r = yr/yi
+        d = yi + r*yr
+        [(r*xr + xi)/d, (r*xi - xr)/d]
+      end
+    end
+
+
+    # Symmetric Householder reduction to tridiagonal form.
+
+    def tridiagonalize
+
+      #  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.
+
+        @size.times do |j|
+          @d[j] = @v[@size-1][j]
+        end
+
+        # Householder reduction to tridiagonal form.
+
+        (@size-1).downto(0+1) do |i|
+
+          # Scale to avoid under/overflow.
+
+          scale = 0.0
+          h = 0.0
+          i.times do |k|
+            scale = scale + @d[k].abs
+          end
+          if (scale == 0.0)
+            @e[i] = @d[i-1]
+            i.times do |j|
+              @d[j] = @v[i-1][j]
+              @v[i][j] = 0.0
+              @v[j][i] = 0.0
+            end
+          else
+
+            # Generate Householder vector.
+
+            i.times do |k|
+              @d[k] /= scale
+              h += @d[k] * @d[k]
+            end
+            f = @d[i-1]
+            g = Math.sqrt(h)
+            if (f > 0)
+              g = -g
+            end
+            @e[i] = scale * g
+            h -= f * g
+            @d[i-1] = f - g
+            i.times do |j|
+              @e[j] = 0.0
+            end
+
+            # Apply similarity transformation to remaining columns.
+
+            i.times do |j|
+              f = @d[j]
+              @v[j][i] = f
+              g = @e[j] + @v[j][j] * f
+              (j+1).upto(i-1) do |k|
+                g += @v[k][j] * @d[k]
+                @e[k] += @v[k][j] * f
+              end
+              @e[j] = g
+            end
+            f = 0.0
+            i.times do |j|
+              @e[j] /= h
+              f += @e[j] * @d[j]
+            end
+            hh = f / (h + h)
+            i.times do |j|
+              @e[j] -= hh * @d[j]
+            end
+            i.times do |j|
+              f = @d[j]
+              g = @e[j]
+              j.upto(i-1) do |k|
+                @v[k][j] -= (f * @e[k] + g * @d[k])
+              end
+              @d[j] = @v[i-1][j]
+              @v[i][j] = 0.0
+            end
+          end
+          @d[i] = h
+        end
+
+        # Accumulate transformations.
+
+        0.upto(@size-1-1) do |i|
+          @v[@size-1][i] = @v[i][i]
+          @v[i][i] = 1.0
+          h = @d[i+1]
+          if (h != 0.0)
+            0.upto(i) do |k|
+              @d[k] = @v[k][i+1] / h
+            end
+            0.upto(i) do |j|
+              g = 0.0
+              0.upto(i) do |k|
+                g += @v[k][i+1] * @v[k][j]
+              end
+              0.upto(i) do |k|
+                @v[k][j] -= g * @d[k]
+              end
+            end
+          end
+          0.upto(i) do |k|
+            @v[k][i+1] = 0.0
+          end
+        end
+        @size.times do |j|
+          @d[j] = @v[@size-1][j]
+          @v[@size-1][j] = 0.0
+        end
+        @v[@size-1][@size-1] = 1.0
+        @e[0] = 0.0
+      end
+
+
+    # Symmetric tridiagonal QL algorithm.
+
+    def diagonalize
+      #  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.
+
+      1.upto(@size-1) do |i|
+        @e[i-1] = @e[i]
+      end
+      @e[@size-1] = 0.0
+
+      f = 0.0
+      tst1 = 0.0
+      eps = Float::EPSILON
+      @size.times do |l|
+
+        # Find small subdiagonal element
+
+        tst1 = [tst1, @d[l].abs + @e[l].abs].max
+        m = l
+        while (m < @size) do
+          if (@e[m].abs <= eps*tst1)
+            break
+          end
+          m+=1
+        end
+
+        # If m == l, @d[l] is an eigenvalue,
+        # otherwise, iterate.
+
+        if (m > l)
+          iter = 0
+          begin
+            iter = iter + 1  # (Could check iteration count here.)
+
+            # Compute implicit shift
+
+            g = @d[l]
+            p = (@d[l+1] - g) / (2.0 * @e[l])
+            r = Math.hypot(p, 1.0)
+            if (p < 0)
+              r = -r
+            end
+            @d[l] = @e[l] / (p + r)
+            @d[l+1] = @e[l] * (p + r)
+            dl1 = @d[l+1]
+            h = g - @d[l]
+            (l+2).upto(@size-1) do |i|
+              @d[i] -= h
+            end
+            f += h
+
+            # Implicit QL transformation.
+
+            p = @d[m]
+            c = 1.0
+            c2 = c
+            c3 = c
+            el1 = @e[l+1]
+            s = 0.0
+            s2 = 0.0
+            (m-1).downto(l) do |i|
+              c3 = c2
+              c2 = c
+              s2 = s
+              g = c * @e[i]
+              h = c * p
+              r = Math.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.
+
+              @size.times do |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
+              end
+            end
+            p = -s * s2 * c3 * el1 * @e[l] / dl1
+            @e[l] = s * p
+            @d[l] = c * p
+
+            # Check for convergence.
+
+          end while (@e[l].abs > eps*tst1)
+        end
+        @d[l] = @d[l] + f
+        @e[l] = 0.0
+      end
+
+      # Sort eigenvalues and corresponding vectors.
+
+      0.upto(@size-2) do |i|
+        k = i
+        p = @d[i]
+        (i+1).upto(@size-1) do |j|
+          if (@d[j] < p)
+            k = j
+            p = @d[j]
+          end
+        end
+        if (k != i)
+          @d[k] = @d[i]
+          @d[i] = p
+          @size.times do |j|
+            p = @v[j][i]
+            @v[j][i] = @v[j][k]
+            @v[j][k] = p
+          end
+        end
+      end
+    end
+
+    # Nonsymmetric reduction to Hessenberg form.
+
+    def reduce_to_hessenberg
+      #  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.
+
+      low = 0
+      high = @size-1
+
+      (low+1).upto(high-1) do |m|
+
+        # Scale column.
+
+        scale = 0.0
+        m.upto(high) do |i|
+          scale = scale + @h[i][m-1].abs
+        end
+        if (scale != 0.0)
+
+          # Compute Householder transformation.
+
+          h = 0.0
+          high.downto(m) do |i|
+            @ort[i] = @h[i][m-1]/scale
+            h += @ort[i] * @ort[i]
+          end
+          g = Math.sqrt(h)
+          if (@ort[m] > 0)
+            g = -g
+          end
+          h -= @ort[m] * g
+          @ort[m] = @ort[m] - g
+
+          # Apply Householder similarity transformation
+          # @h = (I-u*u'/h)*@h*(I-u*u')/h)
+
+          m.upto(@size-1) do |j|
+            f = 0.0
+            high.downto(m) do |i|
+              f += @ort[i]*@h[i][j]
+            end
+            f = f/h
+            m.upto(high) do |i|
+              @h[i][j] -= f*@ort[i]
+            end
+          end
+
+          0.upto(high) do |i|
+            f = 0.0
+            high.downto(m) do |j|
+              f += @ort[j]*@h[i][j]
+            end
+            f = f/h
+            m.upto(high) do |j|
+              @h[i][j] -= f*@ort[j]
+            end
+          end
+          @ort[m] = scale*@ort[m]
+          @h[m][m-1] = scale*g
+        end
+      end
+
+      # Accumulate transformations (Algol's ortran).
+
+      @size.times do |i|
+        @size.times do |j|
+          @v[i][j] = (i == j ? 1.0 : 0.0)
+        end
+      end
+
+      (high-1).downto(low+1) do |m|
+        if (@h[m][m-1] != 0.0)
+          (m+1).upto(high) do |i|
+            @ort[i] = @h[i][m-1]
+          end
+          m.upto(high) do |j|
+            g = 0.0
+            m.upto(high) do |i|
+              g += @ort[i] * @v[i][j]
+            end
+            # Double division avoids possible underflow
+            g = (g / @ort[m]) / @h[m][m-1]
+            m.upto(high) do |i|
+              @v[i][j] += g * @ort[i]
+            end
+          end
+        end
+      end
+    end
+
+
+
+    # Nonsymmetric reduction from Hessenberg to real Schur form.
+
+    def hessenberg_to_real_schur
+
+      #  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
+
+      nn = @size
+      n = nn-1
+      low = 0
+      high = nn-1
+      eps = Float::EPSILON
+      exshift = 0.0
+      p=q=r=s=z=0
+
+      # Store roots isolated by balanc and compute matrix norm
+
+      norm = 0.0
+      nn.times do |i|
+        if (i < low || i > high)
+          @d[i] = @h[i][i]
+          @e[i] = 0.0
+        end
+        ([i-1, 0].max).upto(nn-1) do |j|
+          norm = norm + @h[i][j].abs
+        end
+      end
+
+      # Outer loop over eigenvalue index
+
+      iter = 0
+      while (n >= low) do
+
+        # Look for single small sub-diagonal element
+
+        l = n
+        while (l > low) do
+          s = @h[l-1][l-1].abs + @h[l][l].abs
+          if (s == 0.0)
+            s = norm
+          end
+          if (@h[l][l-1].abs < eps * s)
+            break
+          end
+          l-=1
+        end
+
+        # 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-=1
+          iter = 0
+
+        # Two roots found
+
+        elsif (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 = Math.sqrt(q.abs)
+          @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
+            end
+            @d[n-1] = x + z
+            @d[n] = @d[n-1]
+            if (z != 0.0)
+              @d[n] = x - w / z
+            end
+            @e[n-1] = 0.0
+            @e[n] = 0.0
+            x = @h[n][n-1]
+            s = x.abs + z.abs
+            p = x / s
+            q = z / s
+            r = Math.sqrt(p * p+q * q)
+            p /= r
+            q /= r
+
+            # Row modification
+
+            (n-1).upto(nn-1) do |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
+            end
+
+            # Column modification
+
+            0.upto(n) do |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
+            end
+
+            # Accumulate transformations
+
+            low.upto(high) do |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
+            end
+
+          # Complex pair
+
+          else
+            @d[n-1] = x + p
+            @d[n] = x + p
+            @e[n-1] = z
+            @e[n] = -z
+          end
+          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]
+          end
+
+          # Wilkinson's original ad hoc shift
+
+          if (iter == 10)
+            exshift += x
+            low.upto(n) do |i|
+              @h[i][i] -= x
+            end
+            s = @h[n][n-1].abs + @h[n-1][n-2].abs
+            x = y = 0.75 * s
+            w = -0.4375 * s * s
+          end
+
+          # MATLAB's new ad hoc shift
+
+          if (iter == 30)
+             s = (y - x) / 2.0
+             s *= s + w
+             if (s > 0)
+                s = Math.sqrt(s)
+                if (y < x)
+                  s = -s
+                end
+                s = x - w / ((y - x) / 2.0 + s)
+                low.upto(n) do |i|
+                  @h[i][i] -= s
+                end
+                exshift += s
+                x = y = w = 0.964
+             end
+          end
+
+          iter = iter + 1  # (Could check iteration count here.)
+
+          # Look for two consecutive small sub-diagonal elements
+
+          m = n-2
+          while (m >= l) do
+            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 = p.abs + q.abs + r.abs
+            p /= s
+            q /= s
+            r /= s
+            if (m == l)
+              break
+            end
+            if (@h[m][m-1].abs * (q.abs + r.abs) <
+              eps * (p.abs * (@h[m-1][m-1].abs + z.abs +
+              @h[m+1][m+1].abs)))
+                break
+            end
+            m-=1
+          end
+
+          (m+2).upto(n) do |i|
+            @h[i][i-2] = 0.0
+            if (i > m+2)
+              @h[i][i-3] = 0.0
+            end
+          end
+
+          # Double QR step involving rows l:n and columns m:n
+
+          m.upto(n-1) do |k|
+            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 = p.abs + q.abs + r.abs
+              if (x != 0.0)
+                p /= x
+                q /= x
+                r /= x
+              end
+            end
+            if (x == 0.0)
+              break
+            end
+            s = Math.sqrt(p * p + q * q + r * r)
+            if (p < 0)
+              s = -s
+            end
+            if (s != 0)
+              if (k != m)
+                @h[k][k-1] = -s * x
+              elsif (l != m)
+                @h[k][k-1] = -@h[k][k-1]
+              end
+              p += s
+              x = p / s
+              y = q / s
+              z = r / s
+              q /= p
+              r /= p
+
+              # Row modification
+
+              k.upto(nn-1) do |j|
+                p = @h[k][j] + q * @h[k+1][j]
+                if (notlast)
+                  p += r * @h[k+2][j]
+                  @h[k+2][j] = @h[k+2][j] - p * z
+                end
+                @h[k][j] = @h[k][j] - p * x
+                @h[k+1][j] = @h[k+1][j] - p * y
+              end
+
+              # Column modification
+
+              0.upto([n, k+3].min) do |i|
+                p = x * @h[i][k] + y * @h[i][k+1]
+                if (notlast)
+                  p += z * @h[i][k+2]
+                  @h[i][k+2] = @h[i][k+2] - p * r
+                end
+                @h[i][k] = @h[i][k] - p
+                @h[i][k+1] = @h[i][k+1] - p * q
+              end
+
+              # Accumulate transformations
+
+              low.upto(high) do |i|
+                p = x * @v[i][k] + y * @v[i][k+1]
+                if (notlast)
+                  p += z * @v[i][k+2]
+                  @v[i][k+2] = @v[i][k+2] - p * r
+                end
+                @v[i][k] = @v[i][k] - p
+                @v[i][k+1] = @v[i][k+1] - p * q
+              end
+            end  # (s != 0)
+          end  # k loop
+        end  # check convergence
+      end  # while (n >= low)
+
+      # Backsubstitute to find vectors of upper triangular form
+
+      if (norm == 0.0)
+        return
+      end
+
+      (nn-1).downto(0) do |n|
+        p = @d[n]
+        q = @e[n]
+
+        # Real vector
+
+        if (q == 0)
+          l = n
+          @h[n][n] = 1.0
+          (n-1).downto(0) do |i|
+            w = @h[i][i] - p
+            r = 0.0
+            l.upto(n) do |j|
+              r += @h[i][j] * @h[j][n]
+            end
+            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)
+                end
+
+              # 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 (x.abs > z.abs)
+                  @h[i+1][n] = (-r - w * t) / x
+                else
+                  @h[i+1][n] = (-s - y * t) / z
+                end
+              end
+
+              # Overflow control
+
+              t = @h[i][n].abs
+              if ((eps * t) * t > 1)
+                i.upto(n) do |j|
+                  @h[j][n] = @h[j][n] / t
+                end
+              end
+            end
+          end
+
+        # Complex vector
+
+        elsif (q < 0)
+          l = n-1
+
+          # Last vector component imaginary so matrix is triangular
+
+          if (@h[n][n-1].abs > @h[n-1][n].abs)
+            @h[n-1][n-1] = q / @h[n][n-1]
+            @h[n-1][n] = -(@h[n][n] - p) / @h[n][n-1]
+          else
+            cdivr, cdivi = 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
+          end
+          @h[n][n-1] = 0.0
+          @h[n][n] = 1.0
+          (n-2).downto(0) do |i|
+            ra = 0.0
+            sa = 0.0
+            l.upto(n) do |j|
+              ra = ra + @h[i][j] * @h[j][n-1]
+              sa = sa + @h[i][j] * @h[j][n]
+            end
+            w = @h[i][i] - p
+
+            if (@e[i] < 0.0)
+              z = w
+              r = ra
+              s = sa
+            else
+              l = i
+              if (@e[i] == 0)
+                cdivr, cdivi = 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 * (w.abs + q.abs +
+                  x.abs + y.abs + z.abs)
+                end
+                cdivr, cdivi = 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 (x.abs > (z.abs + q.abs))
+                  @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
+                  cdivr, cdivi = 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
+                end
+              end
+
+              # Overflow control
+
+              t = [@h[i][n-1].abs, @h[i][n].abs].max
+              if ((eps * t) * t > 1)
+                i.upto(n) do |j|
+                  @h[j][n-1] = @h[j][n-1] / t
+                  @h[j][n] = @h[j][n] / t
+                end
+              end
+            end
+          end
+        end
+      end
+
+      # Vectors of isolated roots
+
+      nn.times do |i|
+        if (i < low || i > high)
+          i.upto(nn-1) do |j|
+            @v[i][j] = @h[i][j]
+          end
+        end
+      end
+
+      # Back transformation to get eigenvectors of original matrix
+
+      (nn-1).downto(low) do |j|
+        low.upto(high) do |i|
+          z = 0.0
+          low.upto([j, high].min) do |k|
+            z += @v[i][k] * @h[k][j]
+          end
+          @v[i][j] = z
+        end
+      end
+    end
+
+  end
+end