namespace boost {
namespace numeric {
namespace bindings {
namespace lapack {

//$DESCRIPTION

// overloaded functions to call lapack
namespace detail {
    inline void pteqr( char const compz, integer_t const n, float* d,
            float* e, float* z, integer_t const ldz, float* work,
            integer_t& info ) {
        LAPACK_SPTEQR( &compz, &n, d, e, z, &ldz, work, &info );
    }
    inline void pteqr( char const compz, integer_t const n, double* d,
            double* e, double* z, integer_t const ldz, double* work,
            integer_t& info ) {
        LAPACK_DPTEQR( &compz, &n, d, e, z, &ldz, work, &info );
    }
    inline void pteqr( char const compz, integer_t const n, float* d,
            float* e, traits::complex_f* z, integer_t const ldz, float* work,
            integer_t& info ) {
        LAPACK_CPTEQR( &compz, &n, d, e, traits::complex_ptr(z), &ldz, work,
                &info );
    }
    inline void pteqr( char const compz, integer_t const n, double* d,
            double* e, traits::complex_d* z, integer_t const ldz,
            double* work, integer_t& info ) {
        LAPACK_ZPTEQR( &compz, &n, d, e, traits::complex_ptr(z), &ldz, work,
                &info );
    }
}

// value-type based template
template< typename ValueType, typename Enable = void >
struct pteqr_impl{};

// real specialization
template< typename ValueType >
struct pteqr_impl< ValueType, typename boost::enable_if< traits::is_real<ValueType> >::type > {

    typedef ValueType value_type;
    typedef typename traits::type_traits<ValueType>::real_type real_type;

    // user-defined workspace specialization
    template< typename VectorD, typename VectorE, typename MatrixZ,
            typename WORK >
    static void invoke( char const compz, integer_t const n, VectorD& d,
            VectorE& e, MatrixZ& z, integer_t& info, detail::workspace1<
            WORK > work ) {
        BOOST_STATIC_ASSERT( (boost::is_same< typename traits::vector_traits<
                VectorD >::value_type, typename traits::vector_traits<
                VectorE >::value_type >::value) );
        BOOST_STATIC_ASSERT( (boost::is_same< typename traits::vector_traits<
                VectorD >::value_type, typename traits::matrix_traits<
                MatrixZ >::value_type >::value) );
        BOOST_ASSERT( compz == 'N' || compz == 'V' || compz == 'I' );
        BOOST_ASSERT( n >= 0 );
        BOOST_ASSERT( traits::vector_size(d) >= n );
        BOOST_ASSERT( traits::vector_size(e) >= n-1 );
        BOOST_ASSERT( traits::vector_size(work.select(real_type())) >=
                min_size_work( n ));
        detail::pteqr( compz, n, traits::vector_storage(d),
                traits::vector_storage(e), traits::matrix_storage(z),
                traits::leading_dimension(z),
                traits::vector_storage(work.select(real_type())), info );
    }

    // minimal workspace specialization
    template< typename VectorD, typename VectorE, typename MatrixZ >
    static void invoke( char const compz, integer_t const n, VectorD& d,
            VectorE& e, MatrixZ& z, integer_t& info, minimal_workspace work ) {
        traits::detail::array< real_type > tmp_work( min_size_work( n ) );
        invoke( compz, n, d, e, z, info, workspace( tmp_work ) );
    }

    // optimal workspace specialization
    template< typename VectorD, typename VectorE, typename MatrixZ >
    static void invoke( char const compz, integer_t const n, VectorD& d,
            VectorE& e, MatrixZ& z, integer_t& info, optimal_workspace work ) {
        invoke( compz, n, d, e, z, info, minimal_workspace() );
    }

    static integer_t min_size_work( integer_t const n ) {
        return 4*n;
    }
};

// complex specialization
template< typename ValueType >
struct pteqr_impl< ValueType, typename boost::enable_if< traits::is_complex<ValueType> >::type > {

    typedef ValueType value_type;
    typedef typename traits::type_traits<ValueType>::real_type real_type;

    // user-defined workspace specialization
    template< typename VectorD, typename VectorE, typename MatrixZ,
            typename WORK >
    static void invoke( char const compz, integer_t const n, VectorD& d,
            VectorE& e, MatrixZ& z, integer_t& info, detail::workspace1<
            WORK > work ) {
        BOOST_STATIC_ASSERT( (boost::is_same< typename traits::vector_traits<
                VectorD >::value_type, typename traits::vector_traits<
                VectorE >::value_type >::value) );
        BOOST_ASSERT( compz == 'N' || compz == 'V' || compz == 'I' );
        BOOST_ASSERT( n >= 0 );
        BOOST_ASSERT( traits::vector_size(d) >= n );
        BOOST_ASSERT( traits::vector_size(e) >= n-1 );
        BOOST_ASSERT( traits::vector_size(work.select(real_type())) >=
                min_size_work( n ));
        detail::pteqr( compz, n, traits::vector_storage(d),
                traits::vector_storage(e), traits::matrix_storage(z),
                traits::leading_dimension(z),
                traits::vector_storage(work.select(real_type())), info );
    }

    // minimal workspace specialization
    template< typename VectorD, typename VectorE, typename MatrixZ >
    static void invoke( char const compz, integer_t const n, VectorD& d,
            VectorE& e, MatrixZ& z, integer_t& info, minimal_workspace work ) {
        traits::detail::array< real_type > tmp_work( min_size_work( n ) );
        invoke( compz, n, d, e, z, info, workspace( tmp_work ) );
    }

    // optimal workspace specialization
    template< typename VectorD, typename VectorE, typename MatrixZ >
    static void invoke( char const compz, integer_t const n, VectorD& d,
            VectorE& e, MatrixZ& z, integer_t& info, optimal_workspace work ) {
        invoke( compz, n, d, e, z, info, minimal_workspace() );
    }

    static integer_t min_size_work( integer_t const n ) {
        return 4*n;
    }
};


// template function to call pteqr
template< typename VectorD, typename VectorE, typename MatrixZ,
        typename Workspace >
inline integer_t pteqr( char const compz, integer_t const n, VectorD& d,
        VectorE& e, MatrixZ& z, Workspace work ) {
    typedef typename traits::matrix_traits< MatrixZ >::value_type value_type;
    integer_t info(0);
    pteqr_impl< value_type >::invoke( compz, n, d, e, z, info, work );
    return info;
}

// template function to call pteqr, default workspace type
template< typename VectorD, typename VectorE, typename MatrixZ >
inline integer_t pteqr( char const compz, integer_t const n, VectorD& d,
        VectorE& e, MatrixZ& z ) {
    typedef typename traits::matrix_traits< MatrixZ >::value_type value_type;
    integer_t info(0);
    pteqr_impl< value_type >::invoke( compz, n, d, e, z, info,
            optimal_workspace() );
    return info;
}

}}}} // namespace boost::numeric::bindings::lapack

#endif
