module Matrix::Jacobi

Public Class Methods

J(p, q, c, s, n) click to toggle source

Returns the Jacobi rotation matrix

# File lib/extendmatrix.rb, line 1059
def self.J(p, q, c, s, n)
  j = Matrix.I(n)

  j[p,p] = c 
  j[p,q] = s
  j[q,p] = -s 
  j[q,q] = c

  j
end
max(a) click to toggle source

Returns the index pair (p, q) with 1<= p < q <= n and A[p, q] is the maximum in absolute value

# File lib/extendmatrix.rb, line 1018
def self.max(a)
  n = a.row_size
  max = 0
  p = 0
  q = 0
  n.times{|i|
    ((i+1)...n).each{|j|
      val = a[i, j].abs
      if val > max
        max = val
        p = i
        q = j
      end   
    }
  }
  return p, q
end
off(a) click to toggle source

Returns the nurm of the off-diagonal element

# File lib/extendmatrix.rb, line 1008
def self.off(a)
  n = a.row_size
  sum = 0
  n.times{|i| n.times{|j| sum += a[i, j]**2 if j != i}}
  Math.sqrt(sum)
end
sym_schur2(a, p, q) click to toggle source

Compute the cosine-sine pair (c, s) for the element A[p, q]

# File lib/extendmatrix.rb, line 1039
def self.sym_schur2(a, p, q)
  if a[p, q] != 0
    tau = Float(a[q, q] - a[p, p])/(2 * a[p, q])
    if tau >= 0
      t = 1./(tau + Math.sqrt(1 + tau ** 2))
    else
      t = -1./(-tau + Math.sqrt(1 + tau ** 2))
    end
    c = 1./Math.sqrt(1 + t ** 2)
    s = t * c
  else
    c = 1
    s = 0
  end
  return c, s
end