diff options
| -rw-r--r-- | ChangeLog | 10 | ||||
| -rw-r--r-- | lib/matrix/eigenvalue_decomposition.rb | 206 | ||||
| -rw-r--r-- | test/matrix/test_matrix.rb | 26 |
3 files changed, 139 insertions, 103 deletions
@@ -1,3 +1,13 @@ +Fri Oct 23 10:58:41 2015 Shugo Maeda <shugo@ruby-lang.org> + + * lib/matrix/eigenvalue_decomposition.rb (tridiagonalize): fix + indentation to avoid a warning when the command line option -w of + ruby is specified. + + * lib/matrix/eigenvalue_decomposition.rb (hessenberg_to_real_schur): + change the name of a block parameter to avoid a warning when the + command line option -w of ruby is specified. + Fri Oct 23 10:49:36 2015 Nobuyoshi Nakada <nobu@ruby-lang.org> * compile.c (iseq_compile_each): support safe navigation of simple diff --git a/lib/matrix/eigenvalue_decomposition.rb b/lib/matrix/eigenvalue_decomposition.rb index a92b49c1e0..6c3b794f33 100644 --- a/lib/matrix/eigenvalue_decomposition.rb +++ b/lib/matrix/eigenvalue_decomposition.rb @@ -119,113 +119,113 @@ class Matrix # Auto. Comp., Vol.ii-Linear Algebra, and the corresponding # Fortran subroutine in EISPACK. - @size.times do |j| - @d[j] = @v[@size-1][j] - end + @size.times do |j| + @d[j] = @v[@size-1][j] + end - # Householder reduction to tridiagonal form. + # Householder reduction to tridiagonal form. - (@size-1).downto(0+1) do |i| + (@size-1).downto(0+1) do |i| - # Scale to avoid under/overflow. + # 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 + 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. + # 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 + 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. + # 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 + 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 - @d[i] = h end + @d[i] = h + end - # Accumulate transformations. + # 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(@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| - @d[k] = @v[k][i+1] / h + g += @v[k][i+1] * @v[k][j] 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 + 0.upto(i) do |k| + @v[k][j] -= g * @d[k] 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 + 0.upto(i) do |k| + @v[k][i+1] = 0.0 end - @v[@size-1][@size-1] = 1.0 - @e[0] = 0.0 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. @@ -727,20 +727,20 @@ class Matrix return end - (nn-1).downto(0) do |n| - p = @d[n] - q = @e[n] + (nn-1).downto(0) do |k| + p = @d[k] + q = @e[k] # Real vector if (q == 0) - l = n - @h[n][n] = 1.0 - (n-1).downto(0) do |i| + l = k + @h[k][k] = 1.0 + (k-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] + l.upto(k) do |j| + r += @h[i][j] * @h[j][k] end if (@e[i] < 0.0) z = w @@ -749,9 +749,9 @@ class Matrix l = i if (@e[i] == 0.0) if (w != 0.0) - @h[i][n] = -r / w + @h[i][k] = -r / w else - @h[i][n] = -r / (eps * norm) + @h[i][k] = -r / (eps * norm) end # Solve real equations @@ -761,20 +761,20 @@ class Matrix 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 + @h[i][k] = t if (x.abs > z.abs) - @h[i+1][n] = (-r - w * t) / x + @h[i+1][k] = (-r - w * t) / x else - @h[i+1][n] = (-s - y * t) / z + @h[i+1][k] = (-s - y * t) / z end end # Overflow control - t = @h[i][n].abs + t = @h[i][k].abs if ((eps * t) * t > 1) - i.upto(n) do |j| - @h[j][n] = @h[j][n] / t + i.upto(k) do |j| + @h[j][k] = @h[j][k] / t end end end diff --git a/test/matrix/test_matrix.rb b/test/matrix/test_matrix.rb index e2530f7366..a578aeb1c8 100644 --- a/test/matrix/test_matrix.rb +++ b/test/matrix/test_matrix.rb @@ -582,4 +582,30 @@ class TestMatrix < Test::Unit::TestCase assert_equal @e2, @e2.vstack(@e2) assert_equal SubMatrix, SubMatrix.vstack(@e1).class end + + def test_eigenvalues_and_eigenvectors_symmetric + m = Matrix[ + [8, 1], + [1, 8] + ] + values = m.eigensystem.eigenvalues + assert_in_epsilon(7.0, values[0]) + assert_in_epsilon(9.0, values[1]) + vectors = m.eigensystem.eigenvectors + assert_in_epsilon(-vectors[0][0], vectors[0][1]) + assert_in_epsilon(vectors[1][0], vectors[1][1]) + end + + def test_eigenvalues_and_eigenvectors_nonsymmetric + m = Matrix[ + [8, 1], + [4, 5] + ] + values = m.eigensystem.eigenvalues + assert_in_epsilon(9.0, values[0]) + assert_in_epsilon(4.0, values[1]) + vectors = m.eigensystem.eigenvectors + assert_in_epsilon(vectors[0][0], vectors[0][1]) + assert_in_epsilon(-4 * vectors[1][0], vectors[1][1]) + end end |