/*
 * 
 * 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_PPCON_HPP
#define BOOST_NUMERIC_BINDINGS_LAPACK_PPCON_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 {

    /*
     *  DPPCON estimates the reciprocal of the condition number (in the
     *  1-norm) of a real symmetric positive definite packed matrix using
     *  the Cholesky factorization A = U**T*U or A = L*L**T computed by
     *  DPPTRF.
     *
     *  An estimate is obtained for norm(inv(A)), and the reciprocal of the
     *  condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).    
     */

    namespace detail {

      inline
      void ppcon (char uplo, const int n, const float* ap, const float anorm, 
		  float &rcond, float *work, int *iwork, int &info)
      {
	LAPACK_SPPCON(&uplo, &n, ap, &anorm, &rcond, work, iwork, &info);
      }

      inline
      void ppcon (char uplo, const int n, const double* ap, const double anorm, 
		  double &rcond, double *work, int *iwork, int &info)
      {
	LAPACK_DPPCON(&uplo, &n, ap, &anorm, &rcond, work, iwork, &info);
      }

      inline
      void ppcon (char uplo, const int n, const traits::complex_f* ap, const float anorm, 
		  float &rcond, traits::complex_f *work, float *iwork, int &info)
      {
	LAPACK_CPPCON(&uplo, &n, traits::complex_ptr (ap), &anorm, &rcond, traits::complex_ptr (work), iwork, &info);
      }

      inline
      void ppcon (char uplo, const int n, const traits::complex_d* ap, const double anorm, 
		  double &rcond, traits::complex_d *work, double *iwork, int &info)
      {
	LAPACK_ZPPCON(&uplo, &n, traits::complex_ptr (ap), &anorm, &rcond, traits::complex_ptr (work), iwork, &info);
      }

     
      template <typename SymmMatrA, typename T, typename Work>
      inline
      int ppcon (char const uplo, const SymmMatrA&a, const T anorm, T& rcond, Work& work) {
	typedef typename SymmMatrA::value_type value_type ;
	typedef typename traits::type_traits<value_type>::real_type real_type;

#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
	BOOST_STATIC_ASSERT((boost::is_same<T, real_type>::value));
#endif
	assert(uplo == 'U' || uplo == 'L');  
	
	int const n = traits::matrix_size1 (a);
	assert( n == traits::matrix_size2(a) );
	assert( traits::vector_size(work) >= 3 * n);
	// FIXME: use workspace2 to let user give the workspace?
	traits::detail::array< typename ir_workspace_traits<value_type>::type > iwork(n);
	
	int info; 
	
	detail::ppcon (uplo, n, traits::matrix_storage (a), anorm,
		       rcond,
		       traits::vector_storage(work), 
 		       traits::vector_storage(iwork),
		       info);
	return info; 
      }


      template <typename SymmMatrA, typename T, typename Work>
      inline
      T
      ppcon (const SymmMatrA& a, const T anorm, Work& work) {
	typedef typename SymmMatrA::value_type value_type;
	typedef typename traits::type_traits<value_type>::real_type real_type;

#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
	BOOST_STATIC_ASSERT((boost::is_same<
			     typename traits::matrix_traits<SymmMatrA>::matrix_structure,
			     typename traits::detail::symm_herm_pack_t<
			     typename traits::matrix_traits<SymmMatrA>::value_type
			     >::type
			    >::value));
	BOOST_STATIC_ASSERT((boost::is_same<T, real_type>::value));
#endif
	  
	char uplo = traits::matrix_uplo_tag (a);
	int info; 
	real_type rcond;
	info =  ppcon( uplo, a, anorm, rcond, work);
	return (info == 0) ? rcond : info;
      }
    } // namespace detail



    template <typename SymmMatrA, typename T>
    inline
    T
    ppcon(const SymmMatrA &a, const T anorm, optimal_workspace = optimal_workspace()) {
      typedef typename SymmMatrA::value_type value_type;
      const int n = traits::matrix_size1 (a);
      traits::detail::array<value_type> work(3*n);
      return detail::ppcon(a, anorm, work);
    }

    template <typename SymmMatrA, typename T, typename Work>
    inline
    T
    ppcon(const SymmMatrA &a, const T anorm, detail::workspace1<Work> workspace) {
      return detail::ppcon(a, anorm, workspace.w_);
    }



  }

}}}

#endif