text: /*
   * 
   * Copyright (c) Toon Knapen, Karl Meerbergen & 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_ORMQR_HPP
  #define BOOST_NUMERIC_BINDINGS_LAPACK_ORMQR_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/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 {
  
      ///////////////////////////////////////////////////////////////////
      //
      // Apply the orthogonal transformation abtained by geqrf() to
      // a general matrix.
      // 
      ///////////////////////////////////////////////////////////////////
  
      /* 
       * ormqr() overwrites the general M by N matrix C with
       *
       *                         SIDE = 'L'     SIDE = 'R'
       *  TRANS = 'N':             Q * C          C * Q
       *  TRANS = 'C' || 'T':      Q**H * C       C * Q**H
       *
       *  where Q is a complex unitary matrix defined as the product of k
       *  elementary reflectors
       *
       *        Q = H(1) H(2) . . . H(k)
       *
       *  as returned by geqrf(). Q is of order M if SIDE = 'L' and of order N
       *  if SIDE = 'R'.
       *
       *  Workspace is organized following the arguments in the calling sequence.
       *   optimal_workspace() : for optimizing use of blas 3 kernels
       *   minimal_workspace() : minimum size of workarrays, but does not allow for optimization
       *                         of blas 3 kernels
       *   workspace( work ) where is an array that is used as auxiliary memory with the same value_type
       *                     as C.
       *                     We must have that vector_size( work ) >= matrix_size2( c )
       *                     if SIDE=='L' otherwise  vector_size( work ) >= matrix_size1( c )
       */ 
  
      namespace detail {
  
        inline 
        void ormqr (char const side, char const trans, int const m, int const n,
  		 int const k, const float* a, int const lda,
  		 const float* tau, float* c,
  		 int const ldc, float* work, int const lwork,
                   int& info) 
        {
          LAPACK_SORMQR (&side, &trans, &m, &n, &k,
  		      a, &lda,
  		      tau,
  		      c, &ldc,
  		      work, &lwork,
  		      &info);
        }
  
        inline 
        void ormqr (char const side, char const trans, int const m, int const n,
  		 int const k, const double* a, int const lda,
  		 const double* tau, double* c,
  		 int const ldc, double* work, int const lwork,
                   int& info) 
        {
          LAPACK_DORMQR (&side, &trans, &m, &n, &k,
  		      a, &lda,
  		      tau,
  		      c, &ldc,
  		      work, &lwork,
  		      &info);
        }
  
        inline 
        void ormqr (char const side, char const trans, int const m, int const n,
  		 int const k, const std::complex<float>* a, int const lda,
  		 const std::complex<float>* tau, std::complex<float>* c,
  		 int const ldc, std::complex<float>* work, int const lwork,
                   int& info) 
        {
          LAPACK_CUNMQR (&side, &trans, &m, &n, &k,
  		      reinterpret_cast<const fcomplex_t*>(a), &lda,
  		      reinterpret_cast<const fcomplex_t*>(tau),
  		      reinterpret_cast<fcomplex_t*>(c), &ldc,
  		      reinterpret_cast<fcomplex_t*>(work), &lwork,
  		      &info);
        }
  
        inline 
        void ormqr (char const side, char const trans, int const m, int const n,
  		 int const k, const std::complex<double>* a, int const lda,
  		 const std::complex<double>* tau, std::complex<double>* c,
  		 int const ldc, std::complex<double>* work, int const lwork,
                   int& info) 
        {
          LAPACK_ZUNMQR (&side, &trans, &m, &n, &k,
  		      reinterpret_cast<const dcomplex_t*>(a), &lda,
  		      reinterpret_cast<const dcomplex_t*>(tau),
  		      reinterpret_cast<dcomplex_t*>(c), &ldc,
  		      reinterpret_cast<dcomplex_t*>(work), &lwork,
  		      &info);
        }
  
  
        template <typename A, typename Tau, typename C, typename Work>
        inline
        int ormqr (char side, char trans, const A& a, const Tau& tau, C& c,
                   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)); 
          BOOST_STATIC_ASSERT((boost::is_same<
            typename traits::matrix_traits<C>::matrix_structure, 
            traits::general_t
          >::value)); 
  #endif 
  
          int const m = traits::matrix_size1 (c);
          int const n = traits::matrix_size2 (c);
          int const k = traits::vector_size (tau);
          int const lwork = traits::vector_size (work);
  
          assert ( side=='L' || side=='R' );
          assert ( trans=='N' || trans=='C' || trans=='T' );
          assert ( (side=='L' ?  m >= k : n >= k ) );
  
          assert ( (side=='L' ?
                    m == traits::matrix_size1 (a) :
                    n == traits::matrix_size1 (a) ) );
          assert (traits::matrix_size2 (a)==k); 
  
          assert ( (side=='L' ?
                    lwork >= n : lwork >= m ) );
  
          int info; 
          ormqr (side, trans, m, n, k,
                         traits::matrix_storage (a), 
                         traits::leading_dimension (a),
                         traits::vector_storage (tau),  
                         traits::matrix_storage (c), 
                         traits::leading_dimension (c),
                         traits::vector_storage (work),
                         lwork,
                         info);
          return info; 
        }
      } // namespace detail
  
  
      // Function that allocates temporary arrays for optimal execution time.
      template <typename A, typename Tau, typename C>
      inline
      int ormqr (char side, char trans, const A& a, const Tau& tau, C& c, optimal_workspace ) {
         typedef typename A::value_type                              value_type ;
  
         int const n_w = (side=='L' ? traits::matrix_size2 (c)
                                    : traits::matrix_size1 (c) );
  
         traits::detail::array<value_type> work( std::max(1,n_w*32) );
  
         return detail::ormqr( side, trans, a, tau, c, work );
      }
  
  
      // Function that allocates temporary arrays with minimal size
      template <typename A, typename Tau, typename C>
      inline
      int ormqr (char side, char trans, const A& a, const Tau& tau, C& c, minimal_workspace ) {
         typedef typename A::value_type                              value_type ;
  
         int const n_w = (side=='L' ? traits::matrix_size2 (c)
                                    : traits::matrix_size1 (c) );
  
         traits::detail::array<value_type> work( std::max(1,n_w) );
  
         return detail::ormqr( side, trans, a, tau, c, work );
      }
  
  
      // Function that uses auxiliary array in workspace
      // The calling sequence is ormqr(side, trans, a, tau, c, workspace( work_array ) )
      // where work_array is an array with the same value_type as a and c .
      template <typename A, typename Tau, typename C, typename Work>
      inline
      int ormqr (char side, char trans, const A& a, const Tau& tau, C& c, detail::workspace1<Work> workspace ) {
         typedef typename A::value_type                              value_type ;
  
         return detail::ormqr( side, trans, a, tau, c, workspace.w_ );
      }
  
    }
  
  }}}
  
  #endif 
  