/* * * 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 #include #include #include #include #include #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK # include # include #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 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::real_type real_type; #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK BOOST_STATIC_ASSERT((boost::is_same::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::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 inline T ppcon (const SymmMatrA& a, const T anorm, Work& work) { typedef typename SymmMatrA::value_type value_type; typedef typename traits::type_traits::real_type real_type; #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK BOOST_STATIC_ASSERT((boost::is_same< typename traits::matrix_traits::matrix_structure, typename traits::detail::symm_herm_pack_t< typename traits::matrix_traits::value_type >::type >::value)); BOOST_STATIC_ASSERT((boost::is_same::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 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 work(3*n); return detail::ppcon(a, anorm, work); } template inline T ppcon(const SymmMatrA &a, const T anorm, detail::workspace1 workspace) { return detail::ppcon(a, anorm, workspace.w_); } } }}} #endif