summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--ChangeLog10
-rw-r--r--lib/matrix/eigenvalue_decomposition.rb206
-rw-r--r--test/matrix/test_matrix.rb26
3 files changed, 139 insertions, 103 deletions
diff --git a/ChangeLog b/ChangeLog
index a1dbff5b7c5..01753f80091 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -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 a92b49c1e0a..6c3b794f33e 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 e2530f73665..a578aeb1c86 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