/*
 * 
 * Copyright (c) Kian Ming A. Chai 2009
 *
 * 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)
 *
 * KMAC acknowledges the support DSO NL, Singapore
 *
 * Created based on geqrf.hpp by Toon Knapen, Karl Meerbergen & Kresimir Fresl
 */

#ifndef BOOST_NUMERIC_BINDINGS_LAPACK_TRCON_HPP
#define BOOST_NUMERIC_BINDINGS_LAPACK_TRCON_HPP

#include <complex>
#include <boost/numeric/bindings/traits/traits.hpp>
#include <boost/numeric/bindings/lapack/lapack.h>
#include <boost/numeric/bindings/lapack/workspace.hpp>
#include <boost/numeric/bindings/traits/detail/array.hpp>
#include <boost/numeric/bindings/traits/detail/symm_herm_traits.hpp>

#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
#  include <boost/static_assert.hpp>
#  include <boost/type_traits/is_same.hpp>
#endif 

namespace boost { namespace numeric { namespace bindings { 

  namespace lapack {

    /*
     *
     *  STRCON estimates the reciprocal of the condition number of a
     *  triangular matrix A, in either the 1-norm or the infinity-norm.
     *
     *  The norm of A is computed and an estimate is obtained for
     *  norm(inv(A)), then the reciprocal of the condition number is
     *  computed as
     *     RCOND = 1 / ( norm(A) * norm(inv(A)) ).
     *
     */

    namespace detail {

      inline 
      void trcon (char const norm, char const uplo, char const diag,
		  int const n, const float* a, int const lda, 
		  float &rcond, float *work, int *iwork, int &info) 
      {
        LAPACK_STRCON (&norm, &uplo, &diag, &n, a, &lda, &rcond, work, iwork, &info);
      }

      inline 
      void trcon (char const norm, char const uplo, char const diag,
		  int const n, const double* a, int const lda, 
		  double &rcond, double *work, int *iwork, int &info) 
      {
        LAPACK_DTRCON (&norm, &uplo, &diag, &n, a, &lda, &rcond, work, iwork, &info);
      }

      inline 
      void trcon (char const norm, char const uplo, char const diag,
		  int const n, const traits::complex_f* a, int const lda, 
		  float &rcond, traits::complex_f *work, float *iwork, int &info) 
      {
        LAPACK_CTRCON (&norm, &uplo, &diag, &n, traits::complex_ptr (a), &lda, &rcond, traits::complex_ptr (work), iwork, &info);
      }

      inline 
      void trcon (char const norm, char const uplo, char const diag,
		  int const n, const traits::complex_d* a, int const lda, 
		  double &rcond, traits::complex_d *work, double *iwork, int &info) 
      {
        LAPACK_ZTRCON (&norm, &uplo, &diag, &n, traits::complex_ptr (a), &lda, &rcond, traits::complex_ptr (work), iwork, &info);
      }

     
      template <typename TriA, typename T, typename Work>
      inline
      int trcon (char const norm, char const uplo, const TriA&a, T& rcond, Work& work) {
	typedef typename TriA::value_type value_type ;
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
	BOOST_STATIC_ASSERT((boost::is_same<T, typename traits::type_traits<value_type>::real_type>::value));
#endif
	assert(norm == '1' || norm == 'O' || norm == 'I');  
	assert(uplo == 'U' || uplo == 'L');       

	int const n = traits::matrix_size2 (a);
	int const lda = traits::leading_dimension(a);
	assert( lda >= n);
	assert( traits::vector_size(work) >= 3 * n);

	traits::detail::array< typename ir_workspace_traits<value_type>::type > iwork(n);
	
	int info; 
	
	detail::trcon (norm, uplo, 'N',
		       n, traits::matrix_storage (a), 
		       lda,
		       rcond, 
		       traits::vector_storage(work), 
		       traits::vector_storage(iwork),
		       info);
	return info; 
      }


      template <typename TriA, typename Work>
      inline
      typename traits::type_traits<typename TriA::value_type>::real_type
      trcon (char const norm, const TriA& a, Work& work) {
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
	BOOST_STATIC_ASSERT((boost::is_same<
			     typename traits::matrix_traits<TriA>::matrix_structure, 
			     traits::symmetric_t // FIXME: we shold probably have an adapter for triangular matrices
			     >::value));
#endif
	
	typedef typename TriA::value_type value_type ;
	typedef typename traits::type_traits<value_type>::real_type real_type;
	
	char uplo = traits::matrix_uplo_tag (a);
	int info; 
	real_type rcond;
	info =  trcon( norm, uplo, a, rcond, work);
	return (info == 0) ? rcond : info;
      }
    } // namespace detail



    template <typename TriA>
    inline
    typename traits::type_traits<typename TriA::value_type>::real_type
    trcon(char const norm, const TriA &a, optimal_workspace = optimal_workspace()) {
      typedef typename TriA::value_type value_type;
      const int n = traits::matrix_size2 (a);
      traits::detail::array<value_type> work(3*n);
      return detail::trcon(norm, a, work);
    }

    template <typename TriA, typename Work>
    inline
    typename traits::type_traits<typename TriA::value_type>::real_type
    trcon(char const norm, const TriA &a, detail::workspace1<Work> workspace) {
      return detail::trcon(norm, a, workspace.w_);
    }



  }

}}}

#endif