summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
authormarcandre <marcandre@b2dd03c8-39d4-4d8f-98ff-823fe69b080e>2011-07-01 06:13:35 +0000
committermarcandre <marcandre@b2dd03c8-39d4-4d8f-98ff-823fe69b080e>2011-07-01 06:13:35 +0000
commit004c34f9de58067e627eb714083ca99b17e3fcfd (patch)
treea9ae5d66790766d0966d591795d91f2ccb6b47cf /lib
parent376b8251271da10ff3571b9a15c25333bb4492bc (diff)
* lib/matrix: Add Eigenvalue Decomposition
git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@32351 b2dd03c8-39d4-4d8f-98ff-823fe69b080e
Diffstat (limited to 'lib')
-rw-r--r--lib/matrix.rb22
-rw-r--r--lib/matrix/eigenvalue_decomposition.rb886
2 files changed, 908 insertions, 0 deletions
diff --git a/lib/matrix.rb b/lib/matrix.rb
index 12f9992..8c6b3d6 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 0000000..55e41aa
--- /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