text: /*
   * 
   * Copyright (c) Toon Knapen & Kresimir Fresl 2003
   *
   * Permission to copy, modify, use and distribute this software 
   * for any non-commercial or commercial purpose is granted provided 
   * that this license appear on all copies of the software source code.
   *
   * Authors assume no responsibility whatsoever for its use and makes 
   * no guarantees about its quality, correctness or reliability.
   *
   * KF acknowledges the support of the Faculty of Civil Engineering, 
   * University of Zagreb, Croatia.
   *
   */
  
  #ifndef BOOST_NUMERIC_BINDINGS_LAPACK_TREXC_HPP
  #define BOOST_NUMERIC_BINDINGS_LAPACK_TREXC_HPP
  
  #include <complex>
  #include <boost/numeric/bindings/traits/traits.hpp>
  #include <boost/numeric/bindings/traits/type_traits.hpp>
  #include <boost/numeric/bindings/lapack/lapack.h>
  #include <boost/numeric/bindings/traits/detail/array.hpp>
  // #include <boost/numeric/bindings/traits/std_vector.hpp>
  
  #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
  #  include <boost/static_assert.hpp>
  #  include <boost/type_traits.hpp>
  #endif 
  
  
  namespace boost { namespace numeric { namespace bindings { 
  
    namespace lapack {
  
      ///////////////////////////////////////////////////////////////////
      //
      // Reorder the Schur factorization of a matrix.
      // 
      ///////////////////////////////////////////////////////////////////
  
      /* 
       * trexc()  reorders the Schur factorization of a matrix
       * A =  Q*T*Q**T, so that the diagonal block of T with row
       * index IFST is  moved to row ILST.
       */ 
  
      namespace detail {
        inline 
        int trexc_work_size( int const n, float ) {return n;}
  
        inline 
        int trexc_work_size( int const n, double ) {return n;}
  
        inline 
        int trexc_work_size( int const n, std::complex<float> ) {return 0;}
  
        inline 
        int trexc_work_size( int const n, std::complex<double> ) {return 0;}
      }
  
      // Get the minimum size of the work array.
      template <typename MatrT>
      int trexc_work_size(const MatrT& t) {
         return detail::trexc_work_size( traits::matrix_size1(t), typename MatrT::value_type() );
      }
  
      namespace detail {
        inline 
        void trexc (char const compq, int const n,
                   float* t, int const ldt, float* q, int const ldq, int& ifst, int& ilst,
  		 float* work, int& info) 
        {
          LAPACK_STREXC (&compq, &n, t, &ldt, q, &ldq, &ifst, &ilst, work, &info);
        }
  
        inline 
        void trexc (char const compq, int const n,
                   double* t, int const ldt, double* q, int const ldq, int& ifst, int& ilst,
  		 double* work, int& info) 
        {
          LAPACK_DTREXC (&compq, &n, t, &ldt, q, &ldq, &ifst, &ilst, work, &info);
        }
  
        inline 
        void trexc (char const compq, int const n,
                   std::complex<float>* t, int const ldt, std::complex<float>* q, int const ldq, int& ifst, int& ilst,
  		 float* work, int& info) 
        {
          LAPACK_CTREXC (&compq, &n, traits::complex_ptr(t), &ldt, traits::complex_ptr(q), &ldq, &ifst, &ilst, &info);
        }
  
        inline 
        void trexc (char const compq, int const n,
                   std::complex<double>* t, int const ldt, std::complex<double>* q, int const ldq, int& ifst, int& ilst,
  		 double* work, int& info) 
        {
          LAPACK_ZTREXC (&compq, &n, traits::complex_ptr(t), &ldt, traits::complex_ptr(q), &ldq, &ifst, &ilst, &info);
        }
  
      } 
  
      // Reorder Schur factorization with Schur vectors
      template <typename MatrT, typename Q, typename Work>
      inline
      int trexc (char const compq, MatrT& t, Q& q, int& ifst, int& ilst, Work& work) {
  
  #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
        BOOST_STATIC_ASSERT((boost::is_same<
          typename traits::matrix_traits<MatrT>::matrix_structure, 
          traits::general_t
        >::value)); 
        BOOST_STATIC_ASSERT((boost::is_same<
          typename traits::matrix_traits<Q>::matrix_structure, 
          traits::general_t
        >::value)); 
  #endif 
  
        int const n = traits::matrix_size1 (t);
        assert (n == traits::matrix_size2 (t)); 
        assert (n == traits::matrix_size1 (q)); 
        assert (n == traits::matrix_size2 (q)); 
        assert (trexc_work_size(t) <= traits::vector_size (work)); 
  
        int info; 
        detail::trexc (compq, n,
                      traits::matrix_storage (t), 
                      traits::leading_dimension (t),
                      traits::matrix_storage (q),
                      traits::leading_dimension (q),
  		    ifst, ilst,
                      traits::vector_storage (work),
                      info);
        return info; 
      }
  
    }
  
  }}}
  
  #endif 
  