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

def test_eigen_gen_results(a, b, alpha, beta, evec, count, desc, desc2)
  n = a.size1
  ma = a.to_complex
  mb = b.to_complex
  
  for i in 0...n do
    vi = evec.column(i)
    ai = alpha[i]
    bi = beta[i]
    
    x = mb*vi
    x *= ai
    
    y = ma*vi
    y *= bi
    
    for j in 0...n do
      GSL::Test::test_abs(y[j].real, x[j].real, 1e8*GSL::DBL_EPSILON, 
                          "gen(N=#{n},cnt=#{count}), #{desc}, eigenvalue(#{i},#{j}), real, #{desc2}")
      GSL::Test::test_abs(y[j].imag, x[j].imag, 1e8*GSL::DBL_EPSILON, 
                          "gen(N=#{n},cnt=#{count}), #{desc}, eigenvalue(#{i},#{j}), imag, #{desc2}")
    end
  end
end

def test_eigen_gen_pencil(a, b, count, desc, test_schur, w, wv)
  n = a.size1
  
  aa = a.clone
  bb = b.clone
  
  if test_schur == 1
    alphav, betav, evec, q, z = GSL::Eigen::genv_QZ(aa, bb, wv)
    test_eigen_schur(a, aa, q, z, count, "genv/A", desc)
    test_eigen_schur(b, bb, q, z, count, "genv/B", desc)
  else
    alphav, betav, evec = GSL::Eigen::genv(aa, bb, wv)
  end
  
  test_eigen_gen_results(a, b, alphav, betav, evec, count, desc, "unsorted")
  
  aa = a.clone
  bb = b.clone
  if test_schur == 1
    GSL::Eigen::gen_params(1, 1, 0, w)
    alpha, beta, q, z = GSL::Eigen::gen_QZ(aa, bb, w)
    test_eigen_schur(a, aa, q, z, count, "gen/A", desc)
    test_eigen_schur(b, bb, q, z, count, "gen/B", desc)
  else
    GSL::Eigen::gen_params(0, 0, 0, w)
    alpha, beta = GSL::Eigen::gen(aa, bb, w)
  end
  
  eval = GSL::Vector::Complex.alloc(n)
  evalv = GSL::Vector::Complex.alloc(n)	
  for i in 0...n do
    ai = alpha[i]
    bi = beta[i]
    z = GSL::Complex.alloc(ai.real/bi, ai.imag/bi)
    eval[i] = z
    
    ai = alphav[i]
    bi = betav[i]
    z = GSL::Complex.alloc(ai.real/bi, ai.imag/bi)
    evalv[i] = z
  end
  
  GSL::Eigen::nonsymmv_sort(eval, nil, GSL::EIGEN_SORT_ABS_ASC)
  GSL::Eigen::nonsymmv_sort(evalv, nil, GSL::EIGEN_SORT_ABS_ASC)
  test_eigenvalues_complex(evalv, eval, "gen", desc)
  
  GSL::Eigen::genv_sort(alphav, betav, evec, GSL::EIGEN_SORT_ABS_ASC)
  test_eigen_gen_results(a, b, alphav, betav, evec, count, desc, "abs/asc")
  GSL::Eigen::genv_sort(alphav, betav, evec, GSL::EIGEN_SORT_ABS_DESC)
  test_eigen_gen_results(a, b, alphav, betav, evec, count, desc, "abs/desc")
end

def test_eigen_gen()
  n_max = 20
  rng = GSL::Rng.alloc()
  for n in 1..n_max do
    w = GSL::Eigen::Gen.alloc(n)
    wv = GSL::Eigen::Genv.alloc(n)		
    for i in 0...5 do
      a = create_random_nonsymm_matrix(n, n, rng, -10, 10)
      b = create_random_nonsymm_matrix(n, n, rng, -10, 10)
      test_eigen_gen_pencil(a, b, i, "random", 0, w, wv);
      test_eigen_gen_pencil(a, b, i, "random", 1, w, wv);
    end
  end
  ma = GSL::Matrix.alloc([1, 1, 0, 0, 0, -1, 1, 0, 0], 3, 3)
  mb = GSL::Matrix.alloc([-1, 0, -1, 0, -1, 0, 0, 0, -1], 3, 3)
  w = GSL::Eigen::Gen.alloc(3)
  wv = GSL::Eigen::Genv.alloc(3)	    
  test_eigen_gen_pencil(ma, mb, 0, "integer", 0, w, wv);
  test_eigen_gen_pencil(ma, mb, 0, "integer", 1, w, wv);
end

test_eigen_gen()