reduce_to_hessenberg

# File lib/matrix/eigenvalue_decomposition.rb, line 353
    def reduce_to_hessenberg
      #  This is derived from the Algol procedures orthes and ortran,
      #  by Martin and Wilkinson, Handbook for Auto. Comp.,
      #  Vol.ii-Linear Algebra, and the corresponding
      #  Fortran subroutines in EISPACK.

      low = 0
      high = @size-1

      (low+1).upto(high-1) do |m|

        # Scale column.

        scale = 0.0
        m.upto(high) do |i|
          scale = scale + @h[i][m-1].abs
        end
        if (scale != 0.0)

          # Compute Householder transformation.

          h = 0.0
          high.downto(m) do |i|
            @ort[i] = @h[i][m-1]/scale
            h += @ort[i] * @ort[i]
          end
          g = Math.sqrt(h)
          if (@ort[m] > 0)
            g = -g
          end
          h -= @ort[m] * g
          @ort[m] = @ort[m] - g

          # Apply Householder similarity transformation
          # @h = (I-u*u'/h)*@h*(I-u*u')/h)

          m.upto(@size-1) do |j|
            f = 0.0
            high.downto(m) do |i|
              f += @ort[i]*@h[i][j]
            end
            f = f/h
            m.upto(high) do |i|
              @h[i][j] -= f*@ort[i]
            end
          end

          0.upto(high) do |i|
            f = 0.0
            high.downto(m) do |j|
              f += @ort[j]*@h[i][j]
            end
            f = f/h
            m.upto(high) do |j|
              @h[i][j] -= f*@ort[j]
            end
          end
          @ort[m] = scale*@ort[m]
          @h[m][m-1] = scale*g
        end
      end

      # Accumulate transformations (Algol's ortran).

      @size.times do |i|
        @size.times do |j|
          @v[i][j] = (i == j ? 1.0 : 0.0)
        end
      end

      (high-1).downto(low+1) do |m|
        if (@h[m][m-1] != 0.0)
          (m+1).upto(high) do |i|
            @ort[i] = @h[i][m-1]
          end
          m.upto(high) do |j|
            g = 0.0
            m.upto(high) do |i|
              g += @ort[i] * @v[i][j]
            end
            # Double division avoids possible underflow
            g = (g / @ort[m]) / @h[m][m-1]
            m.upto(high) do |i|
              @v[i][j] += g * @ort[i]
            end
          end
        end
      end
    end