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

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 test_eigen_symm(desc, m)
  n = m.size1
  a = GSL::Matrix.alloc(n, n)

  w1 = GSL::Eigen::Symm::Workspace.alloc(n)
  w2 = GSL::Eigen::Symmv::Workspace.alloc(n)

  GSL::Matrix.memcpy(a, m)
  eval, evec = a.eigen_symmv(w2)
  test_eigen_results(n, m, eval, evec, desc, "unsorted")

  GSL::Matrix.memcpy(a, m)
  eval2 = a.eigen_symm(w1)
  test_eigenvalues(n, eval, eval2, desc, "unsorted")

  GSL::Eigen::Symmv::sort(eval, evec, GSL::Eigen::SORT_VAL_ASC)
  test_eigen_results(n, m, eval, evec, desc, "val/asc")
  GSL::Eigen::Symmv::sort(eval, evec, GSL::Eigen::SORT_VAL_DESC)
  test_eigen_results(n, m, eval, evec, desc, "val/desc")
  GSL::Eigen::Symmv::sort(eval, evec, GSL::Eigen::SORT_ABS_ASC)
  test_eigen_results(n, m, eval, evec, desc, "abs/asc")
  GSL::Eigen::Symmv::sort(eval, evec, GSL::Eigen::SORT_ABS_DESC)
  test_eigen_results(n, m, eval, evec, desc, "abs/desc")
end

def test_eigen_herm(desc, m)
  n = m.size1
  a = GSL::Matrix::Complex.alloc(n, n)

  w1 = GSL::Eigen::Herm::Workspace.alloc(n)
  w2 = GSL::Eigen::Hermv::Workspace.alloc(n)

  GSL::Matrix::Complex.memcpy(a, m)
  eval, evec = a.eigen_hermv(w2)
  test_eigen_complex_results(n, m, eval, evec, desc, "unsorted")

  GSL::Matrix::Complex.memcpy(a, m)
  eval2 = a.eigen_herm(w1)
  test_eigenvalues(n, eval, eval2, desc, "unsorted")

  GSL::Eigen::Hermv::sort(eval, evec, GSL::Eigen::SORT_VAL_ASC)
  test_eigen_complex_results(n, m, eval, evec, desc, "val/asc")
  GSL::Eigen::Hermv::sort(eval, evec, GSL::Eigen::SORT_VAL_DESC)
  test_eigen_complex_results(n, m, eval, evec, desc, "val/desc")
  GSL::Eigen::Hermv::sort(eval, evec, GSL::Eigen::SORT_ABS_ASC)
  test_eigen_complex_results(n, m, eval, evec, desc, "abs/asc")
  GSL::Eigen::Hermv::sort(eval, evec, GSL::Eigen::SORT_ABS_DESC)
  test_eigen_complex_results(n, m, eval, evec, desc, "abs/desc")
end



r = GSL::Matrix.alloc([0, 0, -1, 0], [0, 1, 0, 1], [-1, 0, 0, 0], [0, 1, 0, 0])
test_eigen_symm("symm(4)", r)

c = r.to_complex
test_eigen_herm("herm(4)", c)

r = GSL::Matrix.alloc(4, 4)
r.set_diagonal([1, 2, 3, 4])
test_eigen_symm("symm(4) diag", r)

c = r.to_complex
test_eigen_herm("herm(4) diag", c)