tridiagonalize

# File lib/matrix/eigenvalue_decomposition.rb, line 115
    def tridiagonalize

      #  This is derived from the Algol procedures tred2 by
      #  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
      #  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
      #  Fortran subroutine in EISPACK.

        @size.times do |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