class Matrix::EigenvalueDecomposition

# File lib/matrix/eigenvalue_decomposition.rb, line 18
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_count
 @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

# File lib/matrix/eigenvalue_decomposition.rb, line 72
def eigenvalue_matrix
 Matrix.diagonal(*eigenvalues)
end

# File lib/matrix/eigenvalue_decomposition.rb, line 58
def eigenvalues
 values = @d.dup
 @e.each_with_index{|imag, i| values[i] = Complex(values[i], imag) unless imag == 0}
 values
end

# File lib/matrix/eigenvalue_decomposition.rb, line 42
def eigenvector_matrix
 Matrix.send :new, build_eigenvectors.transpose
end

# File lib/matrix/eigenvalue_decomposition.rb, line 49
def eigenvector_matrix_inv
 r = Matrix.send :new, build_eigenvectors
 r = r.transpose.inverse unless @symmetric
 r
end

# File lib/matrix/eigenvalue_decomposition.rb, line 66
def eigenvectors
 build_eigenvectors.map{|ev| Vector.send :new, ev}
end

# File lib/matrix/eigenvalue_decomposition.rb, line 79
def to_ary
 [v, d, v_inv]
end

# File lib/matrix/eigenvalue_decomposition.rb, line 85
def build_eigenvectors
 # JAMA stores complex eigenvectors in a strange way
 # See http://web.archive.org/web/20111016032731/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-1], -@v[j][i])}
 end
 end
end

# File lib/matrix/eigenvalue_decomposition.rb, line 100
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

# File lib/matrix/eigenvalue_decomposition.rb, line 233
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

# File lib/matrix/eigenvalue_decomposition.rb, line 446
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
 next if x == 0
 p /= x
 q /= x
 r /= x
 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

# File lib/matrix/eigenvalue_decomposition.rb, line 352
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

# File lib/matrix/eigenvalue_decomposition.rb, line 115
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