#!/usr/bin/env ruby
# Ruby/GSL implementation of GSL "eigen/test.c"
require("gsl")
require("../gsl_test.rb")
include GSL::Test

def create_random_symm_matrix(size1, size2, rng, lower, upper)
  m = GSL::Matrix.alloc(size1, size2)
  for i in 0...size1 do
    for j in i...size2 do
      x = rng.uniform()*(upper-lower) + lower
      m[i,j] = x
      m[j,i] = x
    end
  end
  m
end

def create_random_herm_matrix(size1, size2, rng, lower, upper)
  m = GSL::Matrix::Complex.alloc(size1, size2)
  for i in 0...size1 do
    for j in i...size2 do
      re = rng.uniform()* (upper - lower) + lower
      if i == j
      	im = 0.0
      else
        im = rng.uniform()* (upper - lower) + lower
      end
      z = GSL::Complex.alloc(re, im)
      m[i,j] = z
      m[j,i] = z.conjugate
    end
  end
  m
end

def create_random_posdef_matrix(size1, size2, rng)
  m = GSL::Matrix.alloc(size1, size2)
  x = rng.uniform()
  for i in 0...size1 do
    for j in i...size2 do
      a = pow(x, (j - i).to_f)
      m[i,j] = a
      m[j,i] = a
    end
  end
  m
end

def create_random_complex_posdef_matrix(size1, size2, rng)
  m = GSL::Matrix::Complex.calloc(size1, size2)
  n = size1
  
  for i in 0...n do
    x = rng.uniform()
    z = GSL::Complex.alloc(x, 0.0)
    m[i,i] = z
  end
  
  work = GSL::Vector::Complex.alloc(n)
  for i in 0...n do
    for j in 0...n do
      x = 2.0*rng.uniform() - 1.0
      y = 2.0*rng.uniform() - 1.0			
      z = GSL::Complex.alloc(x, y)
      work[j] 
    end
    tau = GSL::Linalg::Complex::householder_transform(work)
    GSL::Linalg::Complex::householder_hm(tau, work, m)
    GSL::Linalg::Complex::householder_mh(tau.conjugate, work, m)		
  end
  m
end

def create_random_nonsymm_matrix(size1, size2, rng, lower, upper)
  m = GSL::Matrix.alloc(size1, size2)
  for i in 0...size1 do
    for j in 0...size2 do
      m[i,j] = rng.uniform()*(upper - lower) + lower
    end
  end
  m
end

def test_eigen_results(n, m, eval, evec, desc, desc2)
  x = GSL::Vector.alloc(n)
  y = GSL::Vector.alloc(n)
  
  # check eigenvalues
  for i in 0...n
    ei = eval[i]
    vi = evec.col(i)
    GSL::Vector.memcpy(x, vi)
    y = GSL::Blas.dgemv(GSL::Blas::NoTrans, 1.0, m, x, 0.0, y)
#    y = GSL::Blas.dgemv(GSL::Blas::NoTrans, 1.0, m, x)
    for j in 0...n
      xj = x[j]
      yj = y[j]
      desc3 = sprintf("%s, eigenvalue(%d,%d), %s", desc, i, j, desc2)
      GSL::Test::test_rel(yj, ei*xj, 1e8*DBL_EPSILON, desc3)
    end
  end
# check eigenvectors are orthonormal
  for i in 0...n
    vi = evec.col(i)
    nrm_v = GSL::Blas.dnrm2(vi)
    desc3 = sprintf("%s, normalized(%d), %s", desc, i, desc2)
    GSL::Test::test_rel(nrm_v, 1.0, n*DBL_EPSILON, desc3)
  end

  for i in 0...n
    vi = evec.col(i)
    for j in (i+1)...n
      vj = evec.col(j)
      vivj = GSL::Blas.ddot(vi, vj)
      desc3 = sprintf("%s, orthogonal(%d,%d), %s", desc, i, j, desc2)
      GSL::Test::test_abs(vivj, 0.0, n*DBL_EPSILON, desc3)
    end
  end
end

def test_eigen_complex_results(n, m, eval, evec, desc, desc2)
  x = GSL::Vector::Complex.alloc(n)
  y = GSL::Vector::Complex.alloc(n)

# check eigenvalues
  for i in 0...n
    ei = eval[i]
    vi = evec.col(i)
    GSL::Vector::Complex.memcpy(x, vi)
    c1 = GSL::Complex.alloc(1.0, 0.0)
    c0 = GSL::Complex.alloc(0.0, 0.0)
    y = GSL::Blas.zgemv(GSL::Blas::NoTrans, c1, m, x, c0, y)
    for j in 0...n
      xj = x[j]
      yj = y[j]
      desc3 = sprintf("%s, eigenvalue(%d,%d), real, %s", desc, i, j, desc2)
      GSL::Test::test_rel(yj.re, ei*xj.re, 1e8*DBL_EPSILON, desc3)
      desc3 = sprintf("%s, eigenvalue(%d,%d), imag, %s", desc, i, j, desc2)
      GSL::Test::test_rel(yj.im, ei*xj.im, 1e8*DBL_EPSILON, desc3)
    end
  end
# check eigenvectors are orthonormal
  for i in 0...n
    vi = evec.col(i)
    nrm_v = GSL::Blas.dznrm2(vi)
    desc3 = sprintf("%s, normalized(%d), %s", desc, i, desc2)
    GSL::Test::test_rel(nrm_v, 1.0, n*DBL_EPSILON, desc3)
  end

  for i in 0...n
    vi = evec.col(i)
    for j in (i+1)...n
      vj = evec.col(j)
      vivj = GSL::Blas.zdotc(vi, vj)
      desc3 = sprintf("%s, orthogonal(%d,%d), %s", desc, i, j, desc2)
      GSL::Test::test_abs(vivj.abs, 0.0, n*DBL_EPSILON, desc3)
    end
  end
end

def test_eigenvalues(n, eval, eval2, desc, desc2)
  for i in 0...n
    ei = eval[i]
    e2i = eval2[i]
    desc3 = sprintf("%s, direct eigenvalue(%d), %s", desc, i, desc2)
    GSL::Test::test_rel(ei, e2i, GSL::DBL_EPSILON, desc3)
  end
end

def chop_subnormals(x)
  # Chop any subnormal values */
  return x.abs < GSL::DBL_MIN ? 0 : x
end

def test_eigenvalues_real(eval, eval2, desc, desc2)
  n = eval.size
  emax = 0
  # check eigenvalues 
  for i in 0...n do
    ei = eval[i]
    if ei.abs > emax
      emax = ei.abs
    end
  end
  
  for i in 0...n do
    ei = eval[i]
    e2i = chop_subnormals(eval2[i])
    GSL::Test::test_abs(ei, e2i, emax * 1e8 * GSL::DBL_EPSILON,  "#{desc}, direct eigenvalue(#{i}), #{desc2}")
  end
end

def test_eigenvalues_complex(eval, eval2, desc, desc2)
  n = eval.size
  for i in 0...n do
    ei = eval[i]
    e2i = eval2[i]
    GSL::Test::test_rel(ei.real, e2i.real, 10*n*GSL::DBL_EPSILON, 
                        "#{desc}, direct eigenvalue(#{i}) real, #{desc2}")
    
    GSL::Test::test_rel(ei.imag, e2i.imag, 10*n*GSL::DBL_EPSILON, 
                        "#{desc}, direct eigenvalue(#{i}) imag, #{desc2}")		
  end
end

def test_eigen_schur(a, s, q, z, count, desc, desc2)
  n = a.size1
  
  t1 = a*z
  t2 = q*s
  for i in 0...n do
    for j in 0...n do
      x = t1[i,j]
      y = t2[i,j]
      GSL::test_abs(x, y, 1.0e8 * GSL::DBL_EPSILON,
                    "#{desc}(N=#{n},cnt=#{count}), #{desc2}, schur(#{i},#{j})")
    end
  end
end