/* * * 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 #include #include #include #include #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK # include # include #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 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::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 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 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 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 work(3*n+1); return geqp3( a, jpvt, tau, work ); } // Computation of the QR factorization. // Workspace is taken from the array in workspace. template inline int geqp3 (A& a, IVec &jpvt, Tau& tau, detail::workspace1 workspace ) { return geqp3( a, jpvt, tau, workspace.w_ ); } } }}} #endif