/*
 * 
 * Copyright (c) Kian Ming A. Chai 2008
 *
 * 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_GEQP3_HPP
#define BOOST_NUMERIC_BINDINGS_LAPACK_GEQP3_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>

#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 {

    ///////////////////////////////////////////////////////////////////
    //
    // QR factorization of a general m x n matrix  with pivoting 
    //                         A * P = Q * R
    // 
    ///////////////////////////////////////////////////////////////////

    /* 
     * geqp3() computes a QR factorization with column pivoting of a
     * matrix A:  A*P = Q*R  using Level 3 BLAS. 
     *
     * Q is a M x min(M,N) matrix with orthogonal
     * and normalized column (i.e. herm(Q) * Q = I) and R is a
     * min(M,N) x N upper triangular matrix.
     *
     * On return of geqp3, the elements on and above the diagonal of
     * A contain the min(M,N) x N upper trapezoidal matrix R (R is
     * upper triangular if  m  >=  n);  the elements below the diagonal,
     * with the array TAU, represent the orthogonal matrix Q as a product
     * of min(M,N) elementary reflectors.
     *
     * On entry, if jpvt(j) !=0, the j-th column of A is permuted
     * to the front of A*P (a leading column); if jpvt(j)=0,
     * the j-th column of A is a free column.
     * On exit, if jpvt(j)=k, then the j-th column of A*P was the
     * the k-th column of A.
     *
     */ 

    namespace detail {

      inline 
      void geqp3 (int const m, int const n,
                 float* a, int const lda, int *jpvt,
                 float* tau, float* work, int const lwork, int& info) 
      {
        LAPACK_SGEQP3 (&m, &n, a, &lda, jpvt, tau, work, &lwork, &info);
      }

      inline 
      void geqp3 (int const m, int const n,
                 double* a, int const lda, int *jpvt,
                 double* tau, double* work, int const lwork, int& info) 
      {
        LAPACK_DGEQP3 (&m, &n, a, &lda, jpvt, tau, work, &lwork, &info);
      }

      inline 
      void geqp3 (int const m, int const n,
                  traits::complex_f* a, int const lda, int *jpvt,
                  traits::complex_f* tau, traits::complex_f* work,
		  int const lwork, int& info) 
      {
        LAPACK_CGEQP3 (&m, &n,
                      traits::complex_ptr (a), &lda, jpvt,
                      traits::complex_ptr (tau),
                      traits::complex_ptr (work), &lwork, &info );
      }
      

      inline 
      void geqp3 (int const m, int const n,
                  traits::complex_d* a, int const lda, int *jpvt,
                  traits::complex_d* tau, traits::complex_d* work,
		  int const lwork, int& info) 
      {
        LAPACK_ZGEQP3 (&m, &n,
		       traits::complex_ptr (a), &lda, jpvt,
		       traits::complex_ptr (tau),
		       traits::complex_ptr (work), &lwork, &info );
      }
      
    } 

    template <typename A, typename IVec, typename Tau, typename Work>
    inline
    int geqp3 (A& a, IVec &jpvt, Tau& tau, Work& work) {
#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
      BOOST_STATIC_ASSERT((boost::is_same<
        typename traits::matrix_traits<A>::matrix_structure, 
        traits::general_t
      >::value)); 
#endif 

      int const m = traits::matrix_size1 (a);
      int const n = traits::matrix_size2 (a);
      int const lwork = traits::vector_size (work);

      assert (std::min(m,n) <= traits::vector_size (tau)); 
      assert (n == traits::vector_size (jpvt)); 
      assert (lwork >= 3*n + 1);

      int info; 
      detail::geqp3 (m, n,
                     traits::matrix_storage (a), 
                     traits::leading_dimension (a),
                     traits::vector_storage (jpvt),  
                     traits::vector_storage (tau),  
                     traits::vector_storage (work),
                     lwork,
                     info);
      return info; 
    }

    // Computation of the QR factorization.
    // Workspace is allocated dynamically so that the optimization of
    // blas 3 calls is optimal.
    template <typename A, typename IVec, typename Tau>
    inline
    int geqp3 (A& a, IVec &jpvt, Tau& tau, optimal_workspace = optimal_workspace()) {
       typedef typename A::value_type value_type ;
       const int n = traits::matrix_size2 (a);
       traits::detail::array<value_type> work(2*n + (n + 1)*32);
       return geqp3( a, jpvt, tau, work );
    }

    // Computation of the QR factorization.
    // Workspace is allocated dynamically to its minimum size.
    // Blas 3 calls are not optimal.
    template <typename A, typename IVec, typename Tau>
    inline
    int geqp3 (A& a, IVec &jpvt, Tau& tau, minimal_workspace ) {
       typedef typename A::value_type value_type ;
       const int n = traits::matrix_size2 (a);
       traits::detail::array<value_type> work(3*n+1);
       return geqp3( a, jpvt, tau, work );
    }

    // Computation of the QR factorization.
    // Workspace is taken from the array in workspace.
    template <typename A, typename IVec, typename Tau, typename Work>
    inline
    int geqp3 (A& a, IVec &jpvt, Tau& tau, detail::workspace1<Work> workspace ) {
       return geqp3( a, jpvt, tau, workspace.w_ );
    }

  }

}}}

#endif