[前][次][番号順一覧][スレッド一覧]

ruby-changes:40154

From: shugo <ko1@a...>
Date: Fri, 23 Oct 2015 11:09:10 +0900 (JST)
Subject: [ruby-changes:40154] shugo:r52235 (trunk): * lib/matrix/eigenvalue_decomposition.rb (tridiagonalize): fix

shugo	2015-10-23 11:09:03 +0900 (Fri, 23 Oct 2015)

  New Revision: 52235

  http://svn.ruby-lang.org/cgi-bin/viewvc.cgi?view=revision&revision=52235

  Log:
    * 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.

  Modified files:
    trunk/ChangeLog
    trunk/lib/matrix/eigenvalue_decomposition.rb
    trunk/test/matrix/test_matrix.rb
Index: ChangeLog
===================================================================
--- ChangeLog	(revision 52234)
+++ ChangeLog	(revision 52235)
@@ -1,3 +1,13 @@ https://github.com/ruby/ruby/blob/trunk/ChangeLog#L1
+Fri Oct 23 10:58:41 2015  Shugo Maeda  <shugo@r...>
+
+	* 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@r...>
 
 	* compile.c (iseq_compile_each): support safe navigation of simple
Index: lib/matrix/eigenvalue_decomposition.rb
===================================================================
--- lib/matrix/eigenvalue_decomposition.rb	(revision 52234)
+++ lib/matrix/eigenvalue_decomposition.rb	(revision 52235)
@@ -119,113 +119,113 @@ class Matrix https://github.com/ruby/ruby/blob/trunk/lib/matrix/eigenvalue_decomposition.rb#L119
       #  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
+
+          # Generate Householder vector.
 
-          scale = 0.0
-          h = 0.0
           i.times do |k|
-            scale = scale + @d[k].abs
+            @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
-          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.
+          # Apply similarity transformation to remaining columns.
 
-            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
+          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
-            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
+            @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 https://github.com/ruby/ruby/blob/trunk/lib/matrix/eigenvalue_decomposition.rb#L727
         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 https://github.com/ruby/ruby/blob/trunk/lib/matrix/eigenvalue_decomposition.rb#L749
               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 https://github.com/ruby/ruby/blob/trunk/lib/matrix/eigenvalue_decomposition.rb#L761
                 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
Index: test/matrix/test_matrix.rb
===================================================================
--- test/matrix/test_matrix.rb	(revision 52234)
+++ test/matrix/test_matrix.rb	(revision 52235)
@@ -582,4 +582,30 @@ class TestMatrix < Test::Unit::TestCase https://github.com/ruby/ruby/blob/trunk/test/matrix/test_matrix.rb#L582
     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

--
ML: ruby-changes@q...
Info: http://www.atdot.net/~ko1/quickml/

[前][次][番号順一覧][スレッド一覧]