//
//  Copyright Toon Knapen, Karl Meerbergen
//
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
//

#include "../../blas/test/random.hpp"

#include <boost/numeric/bindings/lapack/computational/stegr.hpp>
#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
#include <boost/numeric/bindings/traits/ublas_vector.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <algorithm>
#include <complex>
#include <limits>
#include <iostream>


namespace ublas = boost::numeric::ublas;
namespace lapack = boost::numeric::bindings::lapack;


template <typename T>
int do_value_type() {
   const int n = 10 ;

   typedef typename boost::numeric::bindings::traits::type_traits<T>::real_type real_type ;
   typedef std::complex< real_type >                                            complex_type ;

   typedef ublas::matrix<T, ublas::column_major> matrix_type ;

   // Set matrix
   matrix_type z( n, n );
   ublas::vector<real_type> d( n ), e ( n );
   ublas::vector<real_type> w( n );

   std::fill( d.begin(), d.end(), 2.0 ) ;
   std::fill( e.begin(), e.end(), -1.0 ) ;

   // Compute eigendecomposition.
   const char jobz = 'V';     // compute eigenvalues and eigenvectors
   const char range = 'A';
   integer_t m;
   real_type vl = real_type();
   real_type vu = real_type();
   integer_t il = 0;
   integer_t iu = 0;
   real_type abstol = real_type(); // unused
   ublas::vector<integer_t> isuppz( 2 * n );

   lapack::stegr(jobz, range, n, d, e, vl, vu, il, iu, abstol, m, w, z, isuppz);
   //std::cout << "eigenvalues:\n" << w << std::endl;
   //std::cout << "eigenvectors:\n" << z << std::endl;



   matrix_type a(n, n);
   a.clear();
   for (int i = 0; i < n; ++i)
   {
       a(i, i) = T(2.0);
       if (i < n - 1)
       {
           a(i + 1, i) = T(-1.0);
       }
       if (i > 0)
       {
           a(i - 1, i) = T(-1.0);
       }
   }
   matrix_type zaz = ublas::prod(ublas::herm(z), a);
   zaz = ublas::prod(zaz, z);
   //std::cout << zaz << std::endl;
   for (int i = 0; i < n; ++i)
   {
       for (int j = 0; j < n; ++j)
       {
           if (i == j)
           {
               if (std::abs(zaz(i, j) - w(i)) > 100.0 * std::numeric_limits<real_type>::epsilon())
               {
                   // std::cerr << "too large error: " << std::abs(zaz(i, j) - w(i)) << std::endl;
                   return 1;
               }
           }
           else
           {
               if (std::abs(zaz(i, j)) > 100.0 * std::numeric_limits<real_type>::epsilon())
               {
                   // std::cerr << "too large error: " << std::abs(zaz(i, j)) << std::endl;
                   return 1;
               }
           }
       }
   }

   return 0 ;
} // do_value_type()



int main() {
   // Run tests for different value_types
   if (do_value_type<float>()) return EXIT_FAILURE;
   if (do_value_type<double>()) return EXIT_FAILURE;

   // Compilation fails!
   // if (do_value_type<std::complex<float> >()) return 255;
   // if (do_value_type<std::complex<double> >()) return 255;

   std::cout << "Regression test succeeded\n" ;
   return EXIT_SUCCESS;
}