kotovalexarian-likes-github/ruby--ruby
1
0
Fork 0
mirror of https://github.com/ruby/ruby.git synced 2022-11-09 12:17:21 -05:00
Code Releases Activity

* lib/matrix/eigenvalue_decomposition.rb (tridiagonalize): fix

Browse source 
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.

git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@52235 b2dd03c8-39d4-4d8f-98ff-823fe69b080e
This commit is contained in:
shugo 2015-10-23 02:09:03 +00:00
parent 0b7d473734
commit d07ce082a6
3 changed files with 158 additions and 122 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

10
ChangeLog
View file

@ -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> Fri Oct 23 10:49:36 2015 Nobuyoshi Nakada <nobu@ruby-lang.org>
* compile.c (iseq_compile_each): support safe navigation of simple * compile.c (iseq_compile_each): support safe navigation of simple

244
lib/matrix/eigenvalue_decomposition.rb
View file

@ -119,114 +119,114 @@ class Matrix
# Auto. Comp., Vol.ii-Linear Algebra, and the corresponding # Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
# Fortran subroutine in EISPACK. # Fortran subroutine in EISPACK.
@size.times do |j| @size.times do |j|
@d[j] = @v[@size-1][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 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. # Symmetric tridiagonal QL algorithm.
@ -727,20 +727,20 @@ class Matrix
return return
end end
(nn-1).downto(0) do |n| (nn-1).downto(0) do |k|
p = @d[n] p = @d[k]
q = @e[n] q = @e[k]
# Real vector # Real vector
if (q == 0) if (q == 0)
l = n l = k
@h[n][n] = 1.0 @h[k][k] = 1.0
(n-1).downto(0) do |i| (k-1).downto(0) do |i|
w = @h[i][i] - p w = @h[i][i] - p
r = 0.0 r = 0.0
l.upto(n) do |j| l.upto(k) do |j|
r += @h[i][j] * @h[j][n] r += @h[i][j] * @h[j][k]
end end
if (@e[i] < 0.0) if (@e[i] < 0.0)
z = w z = w
@ -749,9 +749,9 @@ class Matrix
l = i l = i
if (@e[i] == 0.0) if (@e[i] == 0.0)
if (w != 0.0) if (w != 0.0)
@h[i][n] = -r / w @h[i][k] = -r / w
else else
@h[i][n] = -r / (eps * norm) @h[i][k] = -r / (eps * norm)
end end
# Solve real equations # Solve real equations
@ -761,20 +761,20 @@ class Matrix
y = @h[i+1][i] y = @h[i+1][i]
q = (@d[i] - p) * (@d[i] - p) + @e[i] * @e[i] q = (@d[i] - p) * (@d[i] - p) + @e[i] * @e[i]
t = (x * s - z * r) / q t = (x * s - z * r) / q
@h[i][n] = t @h[i][k] = t
if (x.abs > z.abs) if (x.abs > z.abs)
@h[i+1][n] = (-r - w * t) / x @h[i+1][k] = (-r - w * t) / x
else else
@h[i+1][n] = (-s - y * t) / z @h[i+1][k] = (-s - y * t) / z
end end
end end
# Overflow control # Overflow control
t = @h[i][n].abs t = @h[i][k].abs
if ((eps * t) * t > 1) if ((eps * t) * t > 1)
i.upto(n) do |j| i.upto(k) do |j|
@h[j][n] = @h[j][n] / t @h[j][k] = @h[j][k] / t
end end
end end
end end

26
test/matrix/test_matrix.rb
View file

@ -582,4 +582,30 @@ class TestMatrix < Test::Unit::TestCase
assert_equal @e2, @e2.vstack(@e2) assert_equal @e2, @e2.vstack(@e2)
assert_equal SubMatrix, SubMatrix.vstack(@e1).class assert_equal SubMatrix, SubMatrix.vstack(@e1).class
end 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 end