/* * * 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_TRCON_HPP #define BOOST_NUMERIC_BINDINGS_LAPACK_TRCON_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 { /* * * STRCON estimates the reciprocal of the condition number of a * triangular matrix A, in either the 1-norm or the infinity-norm. * * The norm of A is computed and an estimate is obtained for * norm(inv(A)), then the reciprocal of the condition number is * computed as * RCOND = 1 / ( norm(A) * norm(inv(A)) ). * */ namespace detail { inline void trcon (char const norm, char const uplo, char const diag, int const n, const float* a, int const lda, float &rcond, float *work, int *iwork, int &info) { LAPACK_STRCON (&norm, &uplo, &diag, &n, a, &lda, &rcond, work, iwork, &info); } inline void trcon (char const norm, char const uplo, char const diag, int const n, const double* a, int const lda, double &rcond, double *work, int *iwork, int &info) { LAPACK_DTRCON (&norm, &uplo, &diag, &n, a, &lda, &rcond, work, iwork, &info); } inline void trcon (char const norm, char const uplo, char const diag, int const n, const traits::complex_f* a, int const lda, float &rcond, traits::complex_f *work, float *iwork, int &info) { LAPACK_CTRCON (&norm, &uplo, &diag, &n, traits::complex_ptr (a), &lda, &rcond, traits::complex_ptr (work), iwork, &info); } inline void trcon (char const norm, char const uplo, char const diag, int const n, const traits::complex_d* a, int const lda, double &rcond, traits::complex_d *work, double *iwork, int &info) { LAPACK_ZTRCON (&norm, &uplo, &diag, &n, traits::complex_ptr (a), &lda, &rcond, traits::complex_ptr (work), iwork, &info); } template inline int trcon (char const norm, char const uplo, const TriA&a, T& rcond, Work& work) { typedef typename TriA::value_type value_type ; #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK BOOST_STATIC_ASSERT((boost::is_same::real_type>::value)); #endif assert(norm == '1' || norm == 'O' || norm == 'I'); assert(uplo == 'U' || uplo == 'L'); int const n = traits::matrix_size2 (a); int const lda = traits::leading_dimension(a); assert( lda >= n); assert( traits::vector_size(work) >= 3 * n); traits::detail::array< typename ir_workspace_traits::type > iwork(n); int info; detail::trcon (norm, uplo, 'N', n, traits::matrix_storage (a), lda, rcond, traits::vector_storage(work), traits::vector_storage(iwork), info); return info; } template inline typename traits::type_traits::real_type trcon (char const norm, const TriA& a, Work& work) { #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK BOOST_STATIC_ASSERT((boost::is_same< typename traits::matrix_traits::matrix_structure, traits::symmetric_t // FIXME: we shold probably have an adapter for triangular matrices >::value)); #endif typedef typename TriA::value_type value_type ; typedef typename traits::type_traits::real_type real_type; char uplo = traits::matrix_uplo_tag (a); int info; real_type rcond; info = trcon( norm, uplo, a, rcond, work); return (info == 0) ? rcond : info; } } // namespace detail template inline typename traits::type_traits::real_type trcon(char const norm, const TriA &a, optimal_workspace = optimal_workspace()) { typedef typename TriA::value_type value_type; const int n = traits::matrix_size2 (a); traits::detail::array work(3*n); return detail::trcon(norm, a, work); } template inline typename traits::type_traits::real_type trcon(char const norm, const TriA &a, detail::workspace1 workspace) { return detail::trcon(norm, a, workspace.w_); } } }}} #endif