tridiagonalize

# File lib/matrix/eigenvalue_decomposition.rb, line 116
    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