Skip to content
This repository
Browse code

RbConfig and 1.9 stdlib cleanup.

* Add autoloaded warning for "Config" constant
* Clean up references to "Config" in files potentially touched in 1.9 mode
* Update missing/divergent stdlib files
* Remove all 'syck' code from 1.9 stdlib
* Add warning to 'ripper' lib
* Copy over missing matrix/ files
  • Loading branch information...
commit 19e94161e0e09e2b354e1f6bc3bed2909b6f89df 1 parent ac05224
Charles Oliver Nutter authored February 27, 2012

Showing 110 changed files with 1,342 additions and 10,485 deletions. Show diff stats Hide diff stats

  1. 2  bench/bench_parser.rb
  2. 4  bench/yarv/report.rb
  3. 4  bench/yarv/runc.rb
  4. 1  lib/ruby/1.9/cgi/.document
  5. 886  lib/ruby/1.9/matrix/eigenvalue_decomposition.rb
  6. 218  lib/ruby/1.9/matrix/lup_decomposition.rb
  7. 6  lib/ruby/1.9/minitest/autorun.rb
  8. 47  lib/ruby/1.9/minitest/benchmark.rb
  9. 13  lib/ruby/1.9/minitest/mock.rb
  10. 39  lib/ruby/1.9/minitest/pride.rb
  11. 68  lib/ruby/1.9/minitest/spec.rb
  12. 192  lib/ruby/1.9/minitest/unit.rb
  13. 2  lib/ruby/1.9/ostruct.rb
  14. 3  lib/ruby/1.9/rbconfig/obsolete.rb
  15. 232  lib/ruby/1.9/rdoc/README
  16. 340  lib/ruby/1.9/rdoc/diagram.rb
  17. 249  lib/ruby/1.9/rdoc/dot.rb
  18. 113  lib/ruby/1.9/rdoc/generator/chm.rb
  19. 100  lib/ruby/1.9/rdoc/generator/chm/chm.rb
  20. 445  lib/ruby/1.9/rdoc/generator/html.rb
  21. 24  lib/ruby/1.9/rdoc/generator/html/common.rb
  22. 92  lib/ruby/1.9/rdoc/generator/html/frameless.rb
  23. 150  lib/ruby/1.9/rdoc/generator/html/hefss.rb
  24. 769  lib/ruby/1.9/rdoc/generator/html/html.rb
  25. 151  lib/ruby/1.9/rdoc/generator/html/kilmer.rb
  26. 427  lib/ruby/1.9/rdoc/generator/html/kilmerfactory.rb
  27. 122  lib/ruby/1.9/rdoc/generator/html/one_page_html.rb
  28. 81  lib/ruby/1.9/rdoc/generator/texinfo.rb
  29. 44  lib/ruby/1.9/rdoc/generator/texinfo/class.texinfo.erb
  30. 6  lib/ruby/1.9/rdoc/generator/texinfo/file.texinfo.erb
  31. 6  lib/ruby/1.9/rdoc/generator/texinfo/method.texinfo.erb
  32. 28  lib/ruby/1.9/rdoc/generator/texinfo/texinfo.erb
  33. 117  lib/ruby/1.9/rdoc/generator/xml.rb
  34. 113  lib/ruby/1.9/rdoc/generator/xml/rdf.rb
  35. 123  lib/ruby/1.9/rdoc/generator/xml/xml.rb
  36. 337  lib/ruby/1.9/rdoc/markup/fragments.rb
  37. 152  lib/ruby/1.9/rdoc/markup/lines.rb
  38. 130  lib/ruby/1.9/rdoc/markup/preprocess.rb
  39. 185  lib/ruby/1.9/rdoc/markup/to_flow.rb
  40. 328  lib/ruby/1.9/rdoc/markup/to_latex.rb
  41. 69  lib/ruby/1.9/rdoc/markup/to_texinfo.rb
  42. 1,835  lib/ruby/1.9/rdoc/parser/f95.rb
  43. 165  lib/ruby/1.9/rdoc/parser/perl.rb
  44. 187  lib/ruby/1.9/rdoc/ri/cache.rb
  45. 156  lib/ruby/1.9/rdoc/ri/descriptions.rb
  46. 392  lib/ruby/1.9/rdoc/ri/display.rb
  47. 28  lib/ruby/1.9/rdoc/ri/gemdirs.rb
  48. 106  lib/ruby/1.9/rdoc/ri/reader.rb
  49. 79  lib/ruby/1.9/rdoc/ri/util.rb
  50. 68  lib/ruby/1.9/rdoc/ri/writer.rb
  51. 64  lib/ruby/1.9/rdoc/template.rb
  52. 52  lib/ruby/1.9/rdoc/tokenstream.rb
  53. 3  lib/ruby/1.9/rinda/.document
  54. 2  lib/ruby/1.9/ripper.rb
  55. 447  lib/ruby/1.9/syck.rb
  56. 242  lib/ruby/1.9/syck/baseemitter.rb
  57. 222  lib/ruby/1.9/syck/basenode.rb
  58. 45  lib/ruby/1.9/syck/constants.rb
  59. 35  lib/ruby/1.9/syck/encoding.rb
  60. 34  lib/ruby/1.9/syck/error.rb
  61. 14  lib/ruby/1.9/syck/loader.rb
  62. 466  lib/ruby/1.9/syck/rubytypes.rb
  63. 41  lib/ruby/1.9/syck/stream.rb
  64. 85  lib/ruby/1.9/syck/stringio.rb
  65. 16  lib/ruby/1.9/syck/syck.rb
  66. 95  lib/ruby/1.9/syck/tag.rb
  67. 192  lib/ruby/1.9/syck/types.rb
  68. 54  lib/ruby/1.9/syck/yamlnode.rb
  69. 54  lib/ruby/1.9/syck/ypath.rb
  70. 7  lib/ruby/1.9/uri/.document
  71. 20  lib/ruby/1.9/uri/common.rb
  72. 1  lib/ruby/1.9/xmlrpc/.document
  73. 22  lib/ruby/1.9/yaml/syck.rb
  74. 2  lib/ruby/shared/Win32API.rb
  75. 2  lib/ruby/shared/ant.rb
  76. 2  lib/ruby/shared/ffi/tools/Rakefile
  77. 48  lib/ruby/shared/mkmf.rb
  78. 2  lib/ruby/shared/rubygems/defaults/jruby.rb
  79. 2  rakelib/commands.rake
  80. 2  rakelib/installer.rake
  81. 2  spec/java_integration/utilities/jar_glob_spec.rb
  82. 4  spec/jruby.1.8.mspec
  83. 4  spec/jruby.1.9.mspec
  84. 4  spec/jruby.cext.mspec
  85. 2  spec/no-library-1.9.mspec
  86. 8  src/jruby/commands.rb
  87. 2  src/jruby/path_helper.rb
  88. 12  src/org/jruby/ext/rbconfig/RbConfigLibrary.java
  89. 4  test/externals/ruby1.8/dbm/test_dbm.rb
  90. 2  test/org/jruby/test/TestRbConfigLibrary.java
  91. 2  test/quiet.rb
  92. 2  test/rubicon/test_io.rb
  93. 6  test/rubicon/test_object_space.rb
  94. 4  test/test_backquote.rb
  95. 4  test/test_command_line_switches.rb
  96. 2  test/test_dir.rb
  97. 2  test/test_etc.rb
  98. 2  test/test_file.rb
  99. 6  test/test_helper.rb
  100. 4  test/test_higher_javasupport.rb
  101. 2  test/test_io.rb
  102. 10  test/test_jar_complete.rb
  103. 4  test/test_kernel.rb
  104. 2  test/test_line_endings.rb
  105. 2  test/test_load.rb
  106. 2  test/test_open3.rb
  107. 6  test/test_process.rb
  108. 6  test/test_rbconfig.rb
  109. 4  test/test_system_error.rb
  110. 4  tool/globals_1_9_3.rb
2  bench/bench_parser.rb
@@ -5,7 +5,7 @@
5 5
 
6 6
 # benchmark 100 parses of the RDoc rb parser
7 7
 ITER_COUNT = 25
8  
-filename = Config::CONFIG['rubylibdir'] + "/rdoc/parsers/parse_rb.rb"
  8
+filename = RbConfig::CONFIG['rubylibdir'] + "/rdoc/parsers/parse_rb.rb"
9 9
 src = File.read(filename)
10 10
 
11 11
 puts "file: " + filename
4  bench/yarv/report.rb
@@ -25,8 +25,8 @@ def exec_command type, file, w
25 25
 
26 26
 def benchmark cmd
27 27
   rubybin = ENV['RUBY'] || File.join(
28  
-    Config::CONFIG["bindir"],
29  
-    Config::CONFIG["ruby_install_name"] + Config::CONFIG["EXEEXT"])
  28
+    RbConfig::CONFIG["bindir"],
  29
+    RbConfig::CONFIG["ruby_install_name"] + RbConfig::CONFIG["EXEEXT"])
30 30
     
31 31
   IO.popen(rubybin, 'r+'){|io|
32 32
     io.write cmd
4  bench/yarv/runc.rb
@@ -6,8 +6,8 @@
6 6
 require 'rbconfig'
7 7
 
8 8
 $rubybin = ENV['RUBY'] || File.join(
9  
-  Config::CONFIG["bindir"],
10  
-  Config::CONFIG["ruby_install_name"] + Config::CONFIG["EXEEXT"])
  9
+  RbConfig::CONFIG["bindir"],
  10
+  RbConfig::CONFIG["ruby_install_name"] + RbConfig::CONFIG["EXEEXT"])
11 11
 
12 12
 def runfile file
13 13
   puts file
1  lib/ruby/1.9/cgi/.document
... ...
@@ -1 +0,0 @@
1  
-session.rb
886  lib/ruby/1.9/matrix/eigenvalue_decomposition.rb
... ...
@@ -0,0 +1,886 @@
  1
+class Matrix
  2
+  # Adapted from JAMA: http://math.nist.gov/javanumerics/jama/
  3
+
  4
+  # Eigenvalues and eigenvectors of a real matrix.
  5
+  #
  6
+  # Computes the eigenvalues and eigenvectors of a matrix A.
  7
+  #
  8
+  # If A is diagonalizable, this provides matrices V and D
  9
+  # such that A = V*D*V.inv, where D is the diagonal matrix with entries
  10
+  # equal to the eigenvalues and V is formed by the eigenvectors.
  11
+  #
  12
+  # If A is symmetric, then V is orthogonal and thus A = V*D*V.t
  13
+
  14
+  class EigenvalueDecomposition
  15
+
  16
+    # Constructs the eigenvalue decomposition for a square matrix +A+
  17
+    #
  18
+    def initialize(a)
  19
+      # @d, @e: Arrays for internal storage of eigenvalues.
  20
+      # @v: Array for internal storage of eigenvectors.
  21
+      # @h: Array for internal storage of nonsymmetric Hessenberg form.
  22
+      raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix)
  23
+      @size = a.row_size
  24
+      @d = Array.new(@size, 0)
  25
+      @e = Array.new(@size, 0)
  26
+
  27
+      if (@symmetric = a.symmetric?)
  28
+        @v = a.to_a
  29
+        tridiagonalize
  30
+        diagonalize
  31
+      else
  32
+        @v = Array.new(@size) { Array.new(@size, 0) }
  33
+        @h = a.to_a
  34
+        @ort = Array.new(@size, 0)
  35
+        reduce_to_hessenberg
  36
+        hessenberg_to_real_schur
  37
+      end
  38
+    end
  39
+
  40
+    # Returns the eigenvector matrix +V+
  41
+    #
  42
+    def eigenvector_matrix
  43
+      Matrix.send :new, build_eigenvectors.transpose
  44
+    end
  45
+    alias v eigenvector_matrix
  46
+
  47
+    # Returns the inverse of the eigenvector matrix +V+
  48
+    #
  49
+    def eigenvector_matrix_inv
  50
+      r = Matrix.send :new, build_eigenvectors
  51
+      r = r.transpose.inverse unless @symmetric
  52
+      r
  53
+    end
  54
+    alias v_inv eigenvector_matrix_inv
  55
+
  56
+    # Returns the eigenvalues in an array
  57
+    #
  58
+    def eigenvalues
  59
+      values = @d.dup
  60
+      @e.each_with_index{|imag, i| values[i] = Complex(values[i], imag) unless imag == 0}
  61
+      values
  62
+    end
  63
+
  64
+    # Returns an array of the eigenvectors
  65
+    #
  66
+    def eigenvectors
  67
+      build_eigenvectors.map{|ev| Vector.send :new, ev}
  68
+    end
  69
+
  70
+    # Returns the block diagonal eigenvalue matrix +D+
  71
+    #
  72
+    def eigenvalue_matrix
  73
+      Matrix.diagonal(*eigenvalues)
  74
+    end
  75
+    alias d eigenvalue_matrix
  76
+
  77
+    # Returns [eigenvector_matrix, eigenvalue_matrix, eigenvector_matrix_inv]
  78
+    #
  79
+    def to_ary
  80
+      [v, d, v_inv]
  81
+    end
  82
+    alias_method :to_a, :to_ary
  83
+
  84
+  private
  85
+    def build_eigenvectors
  86
+      # JAMA stores complex eigenvectors in a strange way
  87
+      # See http://cio.nist.gov/esd/emaildir/lists/jama/msg01021.html
  88
+      @e.each_with_index.map do |imag, i|
  89
+        if imag == 0
  90
+          Array.new(@size){|j| @v[j][i]}
  91
+        elsif imag > 0
  92
+          Array.new(@size){|j| Complex(@v[j][i], @v[j][i+1])}
  93
+        else
  94
+          Array.new(@size){|j| Complex(@v[j][i], -@v[j][i-1])}
  95
+        end
  96
+      end
  97
+    end
  98
+    # Complex scalar division.
  99
+
  100
+    def cdiv(xr, xi, yr, yi)
  101
+      if (yr.abs > yi.abs)
  102
+        r = yi/yr
  103
+        d = yr + r*yi
  104
+        [(xr + r*xi)/d, (xi - r*xr)/d]
  105
+      else
  106
+        r = yr/yi
  107
+        d = yi + r*yr
  108
+        [(r*xr + xi)/d, (r*xi - xr)/d]
  109
+      end
  110
+    end
  111
+
  112
+
  113
+    # Symmetric Householder reduction to tridiagonal form.
  114
+
  115
+    def tridiagonalize
  116
+
  117
+      #  This is derived from the Algol procedures tred2 by
  118
+      #  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
  119
+      #  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
  120
+      #  Fortran subroutine in EISPACK.
  121
+
  122
+        @size.times do |j|
  123
+          @d[j] = @v[@size-1][j]
  124
+        end
  125
+
  126
+        # Householder reduction to tridiagonal form.
  127
+
  128
+        (@size-1).downto(0+1) do |i|
  129
+
  130
+          # Scale to avoid under/overflow.
  131
+
  132
+          scale = 0.0
  133
+          h = 0.0
  134
+          i.times do |k|
  135
+            scale = scale + @d[k].abs
  136
+          end
  137
+          if (scale == 0.0)
  138
+            @e[i] = @d[i-1]
  139
+            i.times do |j|
  140
+              @d[j] = @v[i-1][j]
  141
+              @v[i][j] = 0.0
  142
+              @v[j][i] = 0.0
  143
+            end
  144
+          else
  145
+
  146
+            # Generate Householder vector.
  147
+
  148
+            i.times do |k|
  149
+              @d[k] /= scale
  150
+              h += @d[k] * @d[k]
  151
+            end
  152
+            f = @d[i-1]
  153
+            g = Math.sqrt(h)
  154
+            if (f > 0)
  155
+              g = -g
  156
+            end
  157
+            @e[i] = scale * g
  158
+            h -= f * g
  159
+            @d[i-1] = f - g
  160
+            i.times do |j|
  161
+              @e[j] = 0.0
  162
+            end
  163
+
  164
+            # Apply similarity transformation to remaining columns.
  165
+
  166
+            i.times do |j|
  167
+              f = @d[j]
  168
+              @v[j][i] = f
  169
+              g = @e[j] + @v[j][j] * f
  170
+              (j+1).upto(i-1) do |k|
  171
+                g += @v[k][j] * @d[k]
  172
+                @e[k] += @v[k][j] * f
  173
+              end
  174
+              @e[j] = g
  175
+            end
  176
+            f = 0.0
  177
+            i.times do |j|
  178
+              @e[j] /= h
  179
+              f += @e[j] * @d[j]
  180
+            end
  181
+            hh = f / (h + h)
  182
+            i.times do |j|
  183
+              @e[j] -= hh * @d[j]
  184
+            end
  185
+            i.times do |j|
  186
+              f = @d[j]
  187
+              g = @e[j]
  188
+              j.upto(i-1) do |k|
  189
+                @v[k][j] -= (f * @e[k] + g * @d[k])
  190
+              end
  191
+              @d[j] = @v[i-1][j]
  192
+              @v[i][j] = 0.0
  193
+            end
  194
+          end
  195
+          @d[i] = h
  196
+        end
  197
+
  198
+        # Accumulate transformations.
  199
+
  200
+        0.upto(@size-1-1) do |i|
  201
+          @v[@size-1][i] = @v[i][i]
  202
+          @v[i][i] = 1.0
  203
+          h = @d[i+1]
  204
+          if (h != 0.0)
  205
+            0.upto(i) do |k|
  206
+              @d[k] = @v[k][i+1] / h
  207
+            end
  208
+            0.upto(i) do |j|
  209
+              g = 0.0
  210
+              0.upto(i) do |k|
  211
+                g += @v[k][i+1] * @v[k][j]
  212
+              end
  213
+              0.upto(i) do |k|
  214
+                @v[k][j] -= g * @d[k]
  215
+              end
  216
+            end
  217
+          end
  218
+          0.upto(i) do |k|
  219
+            @v[k][i+1] = 0.0
  220
+          end
  221
+        end
  222
+        @size.times do |j|
  223
+          @d[j] = @v[@size-1][j]
  224
+          @v[@size-1][j] = 0.0
  225
+        end
  226
+        @v[@size-1][@size-1] = 1.0
  227
+        @e[0] = 0.0
  228
+      end
  229
+
  230
+
  231
+    # Symmetric tridiagonal QL algorithm.
  232
+
  233
+    def diagonalize
  234
+      #  This is derived from the Algol procedures tql2, by
  235
+      #  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
  236
+      #  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
  237
+      #  Fortran subroutine in EISPACK.
  238
+
  239
+      1.upto(@size-1) do |i|
  240
+        @e[i-1] = @e[i]
  241
+      end
  242
+      @e[@size-1] = 0.0
  243
+
  244
+      f = 0.0
  245
+      tst1 = 0.0
  246
+      eps = Float::EPSILON
  247
+      @size.times do |l|
  248
+
  249
+        # Find small subdiagonal element
  250
+
  251
+        tst1 = [tst1, @d[l].abs + @e[l].abs].max
  252
+        m = l
  253
+        while (m < @size) do
  254
+          if (@e[m].abs <= eps*tst1)
  255
+            break
  256
+          end
  257
+          m+=1
  258
+        end
  259
+
  260
+        # If m == l, @d[l] is an eigenvalue,
  261
+        # otherwise, iterate.
  262
+
  263
+        if (m > l)
  264
+          iter = 0
  265
+          begin
  266
+            iter = iter + 1  # (Could check iteration count here.)
  267
+
  268
+            # Compute implicit shift
  269
+
  270
+            g = @d[l]
  271
+            p = (@d[l+1] - g) / (2.0 * @e[l])
  272
+            r = Math.hypot(p, 1.0)
  273
+            if (p < 0)
  274
+              r = -r
  275
+            end
  276
+            @d[l] = @e[l] / (p + r)
  277
+            @d[l+1] = @e[l] * (p + r)
  278
+            dl1 = @d[l+1]
  279
+            h = g - @d[l]
  280
+            (l+2).upto(@size-1) do |i|
  281
+              @d[i] -= h
  282
+            end
  283
+            f += h
  284
+
  285
+            # Implicit QL transformation.
  286
+
  287
+            p = @d[m]
  288
+            c = 1.0
  289
+            c2 = c
  290
+            c3 = c
  291
+            el1 = @e[l+1]
  292
+            s = 0.0
  293
+            s2 = 0.0
  294
+            (m-1).downto(l) do |i|
  295
+              c3 = c2
  296
+              c2 = c
  297
+              s2 = s
  298
+              g = c * @e[i]
  299
+              h = c * p
  300
+              r = Math.hypot(p, @e[i])
  301
+              @e[i+1] = s * r
  302
+              s = @e[i] / r
  303
+              c = p / r
  304
+              p = c * @d[i] - s * g
  305
+              @d[i+1] = h + s * (c * g + s * @d[i])
  306
+
  307
+              # Accumulate transformation.
  308
+
  309
+              @size.times do |k|
  310
+                h = @v[k][i+1]
  311
+                @v[k][i+1] = s * @v[k][i] + c * h
  312
+                @v[k][i] = c * @v[k][i] - s * h
  313
+              end
  314
+            end
  315
+            p = -s * s2 * c3 * el1 * @e[l] / dl1
  316
+            @e[l] = s * p
  317
+            @d[l] = c * p
  318
+
  319
+            # Check for convergence.
  320
+
  321
+          end while (@e[l].abs > eps*tst1)
  322
+        end
  323
+        @d[l] = @d[l] + f
  324
+        @e[l] = 0.0
  325
+      end
  326
+
  327
+      # Sort eigenvalues and corresponding vectors.
  328
+
  329
+      0.upto(@size-2) do |i|
  330
+        k = i
  331
+        p = @d[i]
  332
+        (i+1).upto(@size-1) do |j|
  333
+          if (@d[j] < p)
  334
+            k = j
  335
+            p = @d[j]
  336
+          end
  337
+        end
  338
+        if (k != i)
  339
+          @d[k] = @d[i]
  340
+          @d[i] = p
  341
+          @size.times do |j|
  342
+            p = @v[j][i]
  343
+            @v[j][i] = @v[j][k]
  344
+            @v[j][k] = p
  345
+          end
  346
+        end
  347
+      end
  348
+    end
  349
+
  350
+    # Nonsymmetric reduction to Hessenberg form.
  351
+
  352
+    def reduce_to_hessenberg
  353
+      #  This is derived from the Algol procedures orthes and ortran,
  354
+      #  by Martin and Wilkinson, Handbook for Auto. Comp.,
  355
+      #  Vol.ii-Linear Algebra, and the corresponding
  356
+      #  Fortran subroutines in EISPACK.
  357
+
  358
+      low = 0
  359
+      high = @size-1
  360
+
  361
+      (low+1).upto(high-1) do |m|
  362
+
  363
+        # Scale column.
  364
+
  365
+        scale = 0.0
  366
+        m.upto(high) do |i|
  367
+          scale = scale + @h[i][m-1].abs
  368
+        end
  369
+        if (scale != 0.0)
  370
+
  371
+          # Compute Householder transformation.
  372
+
  373
+          h = 0.0
  374
+          high.downto(m) do |i|
  375
+            @ort[i] = @h[i][m-1]/scale
  376
+            h += @ort[i] * @ort[i]
  377
+          end
  378
+          g = Math.sqrt(h)
  379
+          if (@ort[m] > 0)
  380
+            g = -g
  381
+          end
  382
+          h -= @ort[m] * g
  383
+          @ort[m] = @ort[m] - g
  384
+
  385
+          # Apply Householder similarity transformation
  386
+          # @h = (I-u*u'/h)*@h*(I-u*u')/h)
  387
+
  388
+          m.upto(@size-1) do |j|
  389
+            f = 0.0
  390
+            high.downto(m) do |i|
  391
+              f += @ort[i]*@h[i][j]
  392
+            end
  393
+            f = f/h
  394
+            m.upto(high) do |i|
  395
+              @h[i][j] -= f*@ort[i]
  396
+            end
  397
+          end
  398
+
  399
+          0.upto(high) do |i|
  400
+            f = 0.0
  401
+            high.downto(m) do |j|
  402
+              f += @ort[j]*@h[i][j]
  403
+            end
  404
+            f = f/h
  405
+            m.upto(high) do |j|
  406
+              @h[i][j] -= f*@ort[j]
  407
+            end
  408
+          end
  409
+          @ort[m] = scale*@ort[m]
  410
+          @h[m][m-1] = scale*g
  411
+        end
  412
+      end
  413
+
  414
+      # Accumulate transformations (Algol's ortran).
  415
+
  416
+      @size.times do |i|
  417
+        @size.times do |j|
  418
+          @v[i][j] = (i == j ? 1.0 : 0.0)
  419
+        end
  420
+      end
  421
+
  422
+      (high-1).downto(low+1) do |m|
  423
+        if (@h[m][m-1] != 0.0)
  424
+          (m+1).upto(high) do |i|
  425
+            @ort[i] = @h[i][m-1]
  426
+          end
  427
+          m.upto(high) do |j|
  428
+            g = 0.0
  429
+            m.upto(high) do |i|
  430
+              g += @ort[i] * @v[i][j]
  431
+            end
  432
+            # Double division avoids possible underflow
  433
+            g = (g / @ort[m]) / @h[m][m-1]
  434
+            m.upto(high) do |i|
  435
+              @v[i][j] += g * @ort[i]
  436
+            end
  437
+          end
  438
+        end
  439
+      end
  440
+    end
  441
+
  442
+
  443
+
  444
+    # Nonsymmetric reduction from Hessenberg to real Schur form.
  445
+
  446
+    def hessenberg_to_real_schur
  447
+
  448
+      #  This is derived from the Algol procedure hqr2,
  449
+      #  by Martin and Wilkinson, Handbook for Auto. Comp.,
  450
+      #  Vol.ii-Linear Algebra, and the corresponding
  451
+      #  Fortran subroutine in EISPACK.
  452
+
  453
+      # Initialize
  454
+
  455
+      nn = @size
  456
+      n = nn-1
  457
+      low = 0
  458
+      high = nn-1
  459
+      eps = Float::EPSILON
  460
+      exshift = 0.0
  461
+      p=q=r=s=z=0
  462
+
  463
+      # Store roots isolated by balanc and compute matrix norm
  464
+
  465
+      norm = 0.0
  466
+      nn.times do |i|
  467
+        if (i < low || i > high)
  468
+          @d[i] = @h[i][i]
  469
+          @e[i] = 0.0
  470
+        end
  471
+        ([i-1, 0].max).upto(nn-1) do |j|
  472
+          norm = norm + @h[i][j].abs
  473
+        end
  474
+      end
  475
+
  476
+      # Outer loop over eigenvalue index
  477
+
  478
+      iter = 0
  479
+      while (n >= low) do
  480
+
  481
+        # Look for single small sub-diagonal element
  482
+
  483
+        l = n
  484
+        while (l > low) do
  485
+          s = @h[l-1][l-1].abs + @h[l][l].abs
  486
+          if (s == 0.0)
  487
+            s = norm
  488
+          end
  489
+          if (@h[l][l-1].abs < eps * s)
  490
+            break
  491
+          end
  492
+          l-=1
  493
+        end
  494
+
  495
+        # Check for convergence
  496
+        # One root found
  497
+
  498
+        if (l == n)
  499
+          @h[n][n] = @h[n][n] + exshift
  500
+          @d[n] = @h[n][n]
  501
+          @e[n] = 0.0
  502
+          n-=1
  503
+          iter = 0
  504
+
  505
+        # Two roots found
  506
+
  507
+        elsif (l == n-1)
  508
+          w = @h[n][n-1] * @h[n-1][n]
  509
+          p = (@h[n-1][n-1] - @h[n][n]) / 2.0
  510
+          q = p * p + w
  511
+          z = Math.sqrt(q.abs)
  512
+          @h[n][n] = @h[n][n] + exshift
  513
+          @h[n-1][n-1] = @h[n-1][n-1] + exshift
  514
+          x = @h[n][n]
  515
+
  516
+          # Real pair
  517
+
  518
+          if (q >= 0)
  519
+            if (p >= 0)
  520
+              z = p + z
  521
+            else
  522
+              z = p - z
  523
+            end
  524
+            @d[n-1] = x + z
  525
+            @d[n] = @d[n-1]
  526
+            if (z != 0.0)
  527
+              @d[n] = x - w / z
  528
+            end
  529
+            @e[n-1] = 0.0
  530
+            @e[n] = 0.0
  531
+            x = @h[n][n-1]
  532
+            s = x.abs + z.abs
  533
+            p = x / s
  534
+            q = z / s
  535
+            r = Math.sqrt(p * p+q * q)
  536
+            p /= r
  537
+            q /= r
  538
+
  539
+            # Row modification
  540
+
  541
+            (n-1).upto(nn-1) do |j|
  542
+              z = @h[n-1][j]
  543
+              @h[n-1][j] = q * z + p * @h[n][j]
  544
+              @h[n][j] = q * @h[n][j] - p * z
  545
+            end
  546
+
  547
+            # Column modification
  548
+
  549
+            0.upto(n) do |i|
  550
+              z = @h[i][n-1]
  551
+              @h[i][n-1] = q * z + p * @h[i][n]
  552
+              @h[i][n] = q * @h[i][n] - p * z
  553
+            end
  554
+
  555
+            # Accumulate transformations
  556
+
  557
+            low.upto(high) do |i|
  558
+              z = @v[i][n-1]
  559
+              @v[i][n-1] = q * z + p * @v[i][n]
  560
+              @v[i][n] = q * @v[i][n] - p * z
  561
+            end
  562
+
  563
+          # Complex pair
  564
+
  565
+          else
  566
+            @d[n-1] = x + p
  567
+            @d[n] = x + p
  568
+            @e[n-1] = z
  569
+            @e[n] = -z
  570
+          end
  571
+          n -= 2
  572
+          iter = 0
  573
+
  574
+        # No convergence yet
  575
+
  576
+        else
  577
+
  578
+          # Form shift
  579
+
  580
+          x = @h[n][n]
  581
+          y = 0.0
  582
+          w = 0.0
  583
+          if (l < n)
  584
+            y = @h[n-1][n-1]
  585
+            w = @h[n][n-1] * @h[n-1][n]
  586
+          end
  587
+
  588
+          # Wilkinson's original ad hoc shift
  589
+
  590
+          if (iter == 10)
  591
+            exshift += x
  592
+            low.upto(n) do |i|
  593
+              @h[i][i] -= x
  594
+            end
  595
+            s = @h[n][n-1].abs + @h[n-1][n-2].abs
  596
+            x = y = 0.75 * s
  597
+            w = -0.4375 * s * s
  598
+          end
  599
+
  600
+          # MATLAB's new ad hoc shift
  601
+
  602
+          if (iter == 30)
  603
+             s = (y - x) / 2.0
  604
+             s *= s + w
  605
+             if (s > 0)
  606
+                s = Math.sqrt(s)
  607
+                if (y < x)
  608
+                  s = -s
  609
+                end
  610
+                s = x - w / ((y - x) / 2.0 + s)
  611
+                low.upto(n) do |i|
  612
+                  @h[i][i] -= s
  613
+                end
  614
+                exshift += s
  615
+                x = y = w = 0.964
  616
+             end
  617
+          end
  618
+
  619
+          iter = iter + 1  # (Could check iteration count here.)
  620
+
  621
+          # Look for two consecutive small sub-diagonal elements
  622
+
  623
+          m = n-2
  624
+          while (m >= l) do
  625
+            z = @h[m][m]
  626
+            r = x - z
  627
+            s = y - z
  628
+            p = (r * s - w) / @h[m+1][m] + @h[m][m+1]
  629
+            q = @h[m+1][m+1] - z - r - s
  630
+            r = @h[m+2][m+1]
  631
+            s = p.abs + q.abs + r.abs
  632
+            p /= s
  633
+            q /= s
  634
+            r /= s
  635
+            if (m == l)
  636
+              break
  637
+            end
  638
+            if (@h[m][m-1].abs * (q.abs + r.abs) <
  639
+              eps * (p.abs * (@h[m-1][m-1].abs + z.abs +
  640
+              @h[m+1][m+1].abs)))
  641
+                break
  642
+            end
  643
+            m-=1
  644
+          end
  645
+
  646
+          (m+2).upto(n) do |i|
  647
+            @h[i][i-2] = 0.0
  648
+            if (i > m+2)
  649
+              @h[i][i-3] = 0.0
  650
+            end
  651
+          end
  652
+
  653
+          # Double QR step involving rows l:n and columns m:n
  654
+
  655
+          m.upto(n-1) do |k|
  656
+            notlast = (k != n-1)
  657
+            if (k != m)
  658
+              p = @h[k][k-1]
  659
+              q = @h[k+1][k-1]
  660
+              r = (notlast ? @h[k+2][k-1] : 0.0)
  661
+              x = p.abs + q.abs + r.abs
  662
+              if (x != 0.0)
  663
+                p /= x
  664
+                q /= x
  665
+                r /= x
  666
+              end
  667
+            end
  668
+            if (x == 0.0)
  669
+              break
  670
+            end
  671
+            s = Math.sqrt(p * p + q * q + r * r)
  672
+            if (p < 0)
  673
+              s = -s
  674
+            end
  675
+            if (s != 0)
  676
+              if (k != m)
  677
+                @h[k][k-1] = -s * x
  678
+              elsif (l != m)
  679
+                @h[k][k-1] = -@h[k][k-1]
  680
+              end
  681
+              p += s
  682
+              x = p / s
  683
+              y = q / s
  684
+              z = r / s
  685
+              q /= p
  686
+              r /= p
  687
+
  688
+              # Row modification
  689
+
  690
+              k.upto(nn-1) do |j|
  691
+                p = @h[k][j] + q * @h[k+1][j]
  692
+                if (notlast)
  693
+                  p += r * @h[k+2][j]
  694
+                  @h[k+2][j] = @h[k+2][j] - p * z
  695
+                end
  696
+                @h[k][j] = @h[k][j] - p * x
  697
+                @h[k+1][j] = @h[k+1][j] - p * y
  698
+              end
  699
+
  700
+              # Column modification
  701
+
  702
+              0.upto([n, k+3].min) do |i|
  703
+                p = x * @h[i][k] + y * @h[i][k+1]
  704
+                if (notlast)
  705
+                  p += z * @h[i][k+2]
  706
+                  @h[i][k+2] = @h[i][k+2] - p * r
  707
+                end
  708
+                @h[i][k] = @h[i][k] - p
  709
+                @h[i][k+1] = @h[i][k+1] - p * q
  710
+              end
  711
+
  712
+              # Accumulate transformations
  713
+
  714
+              low.upto(high) do |i|
  715
+                p = x * @v[i][k] + y * @v[i][k+1]
  716
+                if (notlast)
  717
+                  p += z * @v[i][k+2]
  718
+                  @v[i][k+2] = @v[i][k+2] - p * r
  719
+                end
  720
+                @v[i][k] = @v[i][k] - p
  721
+                @v[i][k+1] = @v[i][k+1] - p * q
  722
+              end
  723
+            end  # (s != 0)
  724
+          end  # k loop
  725
+        end  # check convergence
  726
+      end  # while (n >= low)
  727
+
  728
+      # Backsubstitute to find vectors of upper triangular form
  729
+
  730
+      if (norm == 0.0)
  731
+        return
  732
+      end
  733
+
  734
+      (nn-1).downto(0) do |n|
  735
+        p = @d[n]
  736
+        q = @e[n]
  737
+
  738
+        # Real vector
  739
+
  740
+        if (q == 0)
  741
+          l = n
  742
+          @h[n][n] = 1.0
  743
+          (n-1).downto(0) do |i|
  744
+            w = @h[i][i] - p
  745
+            r = 0.0
  746
+            l.upto(n) do |j|
  747
+              r += @h[i][j] * @h[j][n]
  748
+            end
  749
+            if (@e[i] < 0.0)
  750
+              z = w
  751
+              s = r
  752
+            else
  753
+              l = i
  754
+              if (@e[i] == 0.0)
  755
+                if (w != 0.0)
  756
+                  @h[i][n] = -r / w
  757
+                else
  758
+                  @h[i][n] = -r / (eps * norm)
  759
+                end
  760
+
  761
+              # Solve real equations
  762
+
  763
+              else
  764
+                x = @h[i][i+1]
  765
+                y = @h[i+1][i]
  766
+                q = (@d[i] - p) * (@d[i] - p) + @e[i] * @e[i]
  767
+                t = (x * s - z * r) / q
  768
+                @h[i][n] = t
  769
+                if (x.abs > z.abs)
  770
+                  @h[i+1][n] = (-r - w * t) / x
  771
+                else
  772
+                  @h[i+1][n] = (-s - y * t) / z
  773
+                end
  774
+              end
  775
+
  776
+              # Overflow control
  777
+
  778
+              t = @h[i][n].abs
  779
+              if ((eps * t) * t > 1)
  780
+                i.upto(n) do |j|
  781
+                  @h[j][n] = @h[j][n] / t
  782
+                end
  783
+              end
  784
+            end
  785
+          end
  786
+
  787
+        # Complex vector
  788
+
  789
+        elsif (q < 0)
  790
+          l = n-1
  791
+
  792
+          # Last vector component imaginary so matrix is triangular
  793
+
  794
+          if (@h[n][n-1].abs > @h[n-1][n].abs)
  795
+            @h[n-1][n-1] = q / @h[n][n-1]
  796
+            @h[n-1][n] = -(@h[n][n] - p) / @h[n][n-1]
  797
+          else
  798
+            cdivr, cdivi = cdiv(0.0, -@h[n-1][n], @h[n-1][n-1]-p, q)
  799
+            @h[n-1][n-1] = cdivr
  800
+            @h[n-1][n] = cdivi
  801
+          end
  802
+          @h[n][n-1] = 0.0
  803
+          @h[n][n] = 1.0
  804
+          (n-2).downto(0) do |i|
  805
+            ra = 0.0
  806
+            sa = 0.0
  807
+            l.upto(n) do |j|
  808
+              ra = ra + @h[i][j] * @h[j][n-1]
  809
+              sa = sa + @h[i][j] * @h[j][n]
  810
+            end
  811
+            w = @h[i][i] - p
  812
+
  813
+            if (@e[i] < 0.0)
  814
+              z = w
  815
+              r = ra
  816
+              s = sa
  817
+            else
  818
+              l = i
  819
+              if (@e[i] == 0)
  820
+                cdivr, cdivi = cdiv(-ra, -sa, w, q)
  821
+                @h[i][n-1] = cdivr
  822
+                @h[i][n] = cdivi
  823
+              else
  824
+
  825
+                # Solve complex equations
  826
+
  827
+                x = @h[i][i+1]
  828
+                y = @h[i+1][i]
  829
+                vr = (@d[i] - p) * (@d[i] - p) + @e[i] * @e[i] - q * q
  830
+                vi = (@d[i] - p) * 2.0 * q
  831
+                if (vr == 0.0 && vi == 0.0)
  832
+                  vr = eps * norm * (w.abs + q.abs +
  833
+                  x.abs + y.abs + z.abs)
  834
+                end
  835
+                cdivr, cdivi = cdiv(x*r-z*ra+q*sa, x*s-z*sa-q*ra, vr, vi)
  836
+                @h[i][n-1] = cdivr
  837
+                @h[i][n] = cdivi
  838
+                if (x.abs > (z.abs + q.abs))
  839
+                  @h[i+1][n-1] = (-ra - w * @h[i][n-1] + q * @h[i][n]) / x
  840
+                  @h[i+1][n] = (-sa - w * @h[i][n] - q * @h[i][n-1]) / x
  841
+                else
  842
+                  cdivr, cdivi = cdiv(-r-y*@h[i][n-1], -s-y*@h[i][n], z, q)
  843
+                  @h[i+1][n-1] = cdivr
  844
+                  @h[i+1][n] = cdivi
  845
+                end
  846
+              end
  847
+
  848
+              # Overflow control
  849
+
  850
+              t = [@h[i][n-1].abs, @h[i][n].abs].max
  851
+              if ((eps * t) * t > 1)
  852
+                i.upto(n) do |j|
  853
+                  @h[j][n-1] = @h[j][n-1] / t
  854
+                  @h[j][n] = @h[j][n] / t
  855
+                end
  856
+              end
  857
+            end
  858
+          end
  859
+        end
  860
+      end
  861
+
  862
+      # Vectors of isolated roots
  863
+
  864
+      nn.times do |i|
  865
+        if (i < low || i > high)
  866
+          i.upto(nn-1) do |j|
  867
+            @v[i][j] = @h[i][j]
  868
+          end
  869
+        end
  870
+      end
  871
+
  872
+      # Back transformation to get eigenvectors of original matrix
  873
+
  874
+      (nn-1).downto(low) do |j|
  875
+        low.upto(high) do |i|
  876
+          z = 0.0
  877
+          low.upto([j, high].min) do |k|
  878
+            z += @v[i][k] * @h[k][j]
  879
+          end
  880
+          @v[i][j] = z
  881
+        end
  882
+      end
  883
+    end
  884
+
  885
+  end
  886
+end
218  lib/ruby/1.9/matrix/lup_decomposition.rb
... ...
@@ -0,0 +1,218 @@
  1
+class Matrix
  2
+  # Adapted from JAMA: http://math.nist.gov/javanumerics/jama/
  3
+
  4
+  #
  5
+  # For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
  6
+  # unit lower triangular matrix L, an n-by-n upper triangular matrix U,
  7
+  # and a m-by-m permutation matrix P so that L*U = P*A.
  8
+  # If m < n, then L is m-by-m and U is m-by-n.
  9
+  #
  10
+  # The LUP decomposition with pivoting always exists, even if the matrix is
  11
+  # singular, so the constructor will never fail.  The primary use of the
  12
+  # LU decomposition is in the solution of square systems of simultaneous
  13
+  # linear equations.  This will fail if singular? returns true.
  14
+  #
  15
+
  16
+  class LUPDecomposition
  17
+    # Returns the lower triangular factor +L+
  18
+
  19
+    include Matrix::ConversionHelper
  20
+
  21
+    def l
  22
+      Matrix.build(@row_size, @col_size) do |i, j|
  23
+        if (i > j)
  24
+          @lu[i][j]
  25
+        elsif (i == j)
  26
+          1
  27
+        else
  28
+          0
  29
+        end
  30
+      end
  31
+    end
  32
+
  33
+    # Returns the upper triangular factor +U+
  34
+
  35
+    def u
  36
+      Matrix.build(@col_size, @col_size) do |i, j|
  37
+        if (i <= j)
  38
+          @lu[i][j]
  39
+        else
  40
+          0
  41
+        end
  42
+      end
  43
+    end
  44
+
  45
+    # Returns the permutation matrix +P+
  46
+
  47
+    def p
  48
+      rows = Array.new(@row_size){Array.new(@col_size, 0)}
  49
+      @pivots.each_with_index{|p, i| rows[i][p] = 1}
  50
+      Matrix.send :new, rows, @col_size
  51
+    end
  52
+
  53
+    # Returns +L+, +U+, +P+ in an array
  54
+
  55
+    def to_ary
  56
+      [l, u, p]
  57
+    end
  58
+    alias_method :to_a, :to_ary
  59
+
  60
+    # Returns the pivoting indices
  61
+
  62
+    attr_reader :pivots
  63
+
  64
+    # Returns +true+ if +U+, and hence +A+, is singular.
  65
+
  66
+    def singular? ()
  67
+      @col_size.times do |j|
  68
+        if (@lu[j][j] == 0)
  69
+          return true
  70
+        end
  71
+      end
  72
+      false
  73
+    end
  74
+
  75
+    # Returns the determinant of +A+, calculated efficiently
  76
+    # from the factorization.
  77
+
  78
+    def det
  79
+      if (@row_size != @col_size)
  80
+        Matrix.Raise Matrix::ErrDimensionMismatch unless square?
  81
+      end
  82
+      d = @pivot_sign
  83
+      @col_size.times do |j|
  84
+        d *= @lu[j][j]
  85
+      end
  86
+      d
  87
+    end
  88
+    alias_method :determinant, :det
  89
+
  90
+    # Returns +m+ so that <tt>A*m = b</tt>,
  91
+    # or equivalently so that <tt>L*U*m = P*b</tt>
  92
+    # +b+ can be a Matrix or a Vector
  93
+
  94
+    def solve b
  95
+      if (singular?)
  96
+        Matrix.Raise Matrix::ErrNotRegular, "Matrix is singular."
  97
+      end
  98
+      if b.is_a? Matrix
  99
+        if (b.row_size != @row_size)
  100
+          Matrix.Raise Matrix::ErrDimensionMismatch
  101
+        end
  102
+
  103
+        # Copy right hand side with pivoting
  104
+        nx = b.column_size
  105
+        m = @pivots.map{|row| b.row(row).to_a}
  106
+
  107
+        # Solve L*Y = P*b
  108
+        @col_size.times do |k|
  109
+          (k+1).upto(@col_size-1) do |i|
  110
+            nx.times do |j|
  111
+              m[i][j] -= m[k][j]*@lu[i][k]
  112
+            end
  113
+          end
  114
+        end
  115
+        # Solve U*m = Y
  116
+        (@col_size-1).downto(0) do |k|
  117
+          nx.times do |j|
  118
+            m[k][j] = m[k][j].quo(@lu[k][k])
  119
+          end
  120
+          k.times do |i|
  121
+            nx.times do |j|
  122
+              m[i][j] -= m[k][j]*@lu[i][k]
  123
+            end
  124
+          end
  125
+        end
  126
+        Matrix.send :new, m, nx
  127
+      else # same algorithm, specialized for simpler case of a vector
  128
+        b = convert_to_array(b)
  129
+        if (b.size != @row_size)
  130
+          Matrix.Raise Matrix::ErrDimensionMismatch
  131
+        end
  132
+
  133
+        # Copy right hand side with pivoting
  134
+        m = b.values_at(*@pivots)