|
|
Ublas : |
From: James N. Knight (nate_at_[hidden])
Date: 2006-10-13 01:43:46
I just realized all this has already been done. The specializations are in numeric/traits/ublas_vector2.hpp.
I should have read the docs more carefully. If you include this file you can call the lapack and blas
functions with a vector in place of a matrix. It treats the vector as a matrix of size nx1 or 1xn.
Nate
James N. Knight wrote:
> I think a more general way to handle this is to place something like
>
> // Added by nate
> template <typename T, typename A>
> inline
> typename vector_traits<ublas::vector<T,A> >::pointer matrix_storage (ublas::vector<T,A>& m) {
> return vector_traits<ublas::vector<T,A> >::storage (m);
> }
>
> template <typename T, typename A>
> inline
> int matrix_size1 (ublas::vector<T,A>& m) { return vector_traits<ublas::vector<T,A> >::size (m); }
>
> template <typename T, typename A>
> inline
> int matrix_size2 (ublas::vector<T,A>& m) { return 1; }
>
> template <typename T, typename A>
> inline
> int leading_dimension (ublas::vector<T,A>& m) {
> return vector_traits<ublas::vector<T,A> >::size (m);
> }
> //
>
> in the numeric/traits/ublas_vector.hpp file. Its should go in the boost::numeric::bindings::traits namespace.
> You might need to add a few more specializations to make everything work.
>
> Specifically, I have to compile the following with -DBOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK because I'm not sure how to specialize matrix_structure for a vector. Any ideas?
>
> #include <boost/numeric/ublas/vector.hpp>
> #include <boost/numeric/ublas/matrix.hpp>
> #include <boost/numeric/ublas/io.hpp>
>
> #include <boost/numeric/bindings/traits/ublas_matrix.hpp>
> #include <boost/numeric/bindings/traits/ublas_symmetric.hpp>
> #include <boost/numeric/bindings/traits/ublas_vector.hpp>
>
> #include <boost/numeric/bindings/lapack/gesv.hpp>
>
>
> namespace ublas = boost::numeric::ublas;
> namespace lapack = boost::numeric::bindings::lapack;
>
> int main(){
> typedef ublas::matrix<double, ublas::column_major> matrix;
>
> matrix P(2,2);
>
> P(0,0) = 2;
> P(0,1) = 0;
> P(1,0) = 1;
> P(1,1) = 1;
>
> ublas::vector<double> b(2);
>
> b(0) = 5;
> b(1) = 8;
>
> lapack::gesv(P,b);
>
> std::cout << b << std::endl;
> }
>
> Hope this helps.
>
> Nate
>
>
> Paul Thompson wrote:
>> On Friday 13 October 2006 12:29, Nico Galoppo wrote:
>>> Hi all,
>>>
>>> Often, when using the ATLAS bindings such as gesv(A,B) to solve a system A
>>> X = B, I only want to solve a system with one rhs, ie. A x = b. I find that
>>> I first have to copy my vector to a matrix and then copy the result back to
>>> a vector.
>>>
>>> Has anyone ever worked on a matrix proxy of a vector, such that this copy
>>> is not necessary anymore?
>>>
>>> Thanks!
>>>
>>> --nico
>>
>> Yes, I achieved this for the LAPACK-ATLAS Cholesky solver. I added a
>> template specialisation for ublas::vector<T,A> at the point where it becomes
>> necessary to call traits to get the assumed parameters. So then I used
>> vector traits to get the necessary terms before calling the C ATLAS function.
>>
>> It doesn't pretend to represent a vector-as-a-matrix. It does result in some
>> code duplication. Perhaps a neater solution exists??
>>
>> Also I only wrote it for boost::numeric::ublas::vector<T,A> not all the
>> vectors supported by the bindings traits.
>>
>> This type of mod would be suitable for the LAPACK-ATLAS gesv I suppose.
>>
>> The modified clapack.hpp is attached. More details below.
>>
>> Hope it helps!
>>
>> Paul T.
>>
>>
>>
>>
>>
>>
>> In more detail:
>>
>> template <typename SymmA, typename MatrB>
>> int potrs (SymmA const& a, MatrB& b)
>>
>> was specialised to:
>>
>> template <typename SymmA, typename T, typename A>
>> int potrs (SymmA const& a, boost::numeric::ublas::vector<T,A> & v)
>>
>> and:
>>
>>
>> template <typename SymmA, typename MatrB>
>> int potrs (CBLAS_UPLO const uplo, SymmA const& a, MatrB& b)
>>
>> was specialised to:
>>
>> template <typename SymmA, typename T,typename A>
>> int potrs (CBLAS_UPLO const uplo, SymmA const& a,
>> boost::numeric::ublas::vector<T,A> & v)
>>
>>
>> Users can still use the same interface, eg:
>> cholesky_substitute (SymmA const& a, MatrB& b) { return potrs (a, b); }
>> still works. The template parameter MatrB ends up as a ublas::vector
>>
>>
>>
>> ------------------------------------------------------------------------
>>
>> /*
>> *
>> * Copyright (c) Kresimir Fresl 2002
>> *
>> * 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.
>> *
>> * Author assumes no responsibility whatsoever for its use and makes
>> * no guarantees about its quality, correctness or reliability.
>> *
>> * Author acknowledges the support of the Faculty of Civil Engineering,
>> * University of Zagreb, Croatia.
>> *
>> */
>>
>> #ifndef BOOST_NUMERIC_BINDINGS_CLAPACK_HPP
>> #define BOOST_NUMERIC_BINDINGS_CLAPACK_HPP
>>
>> #include <cassert>
>> #include <new>
>>
>> #include <boost/numeric/bindings/traits/traits.hpp>
>> #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> # include <boost/numeric/bindings/traits/detail/symm_herm_traits.hpp>
>> #endif
>> #include <boost/numeric/bindings/atlas/cblas_enum.hpp>
>>
>> // see libs/numeric/bindings/atlas/doc/index.html, section 2.5.2
>> //#define BOOST_NUMERIC_BINDINGS_ATLAS_POTRF_BUG
>>
>> #include <boost/numeric/bindings/atlas/clapack_overloads.hpp>
>>
>> #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> # include <boost/static_assert.hpp>
>> # include <boost/type_traits/same_traits.hpp>
>> #endif
>>
>>
>> namespace boost { namespace numeric { namespace bindings {
>>
>> namespace atlas {
>>
>> /////////////////////////////////////////////////////////////////////
>> //
>> // general system of linear equations A * X = B
>> //
>> /////////////////////////////////////////////////////////////////////
>>
>> // gesv(): 'driver' function
>> //
>> // [comments from 'clapack_dgesv.c':]
>> /* clapack_xgesv computes the solution to a system of linear equations
>> * A * X = B,
>> * where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
>> */
>> // [but ATLAS FAQ says:]
>> /* What's the deal with the RHS in the row-major factorization/solves?
>> * Most users are confused by the row major factorization and related
>> * solves. The right-hand side vectors are probably the biggest source
>> * of confusion. The RHS array does not represent a matrix in the
>> * mathematical sense, it is instead a pasting together of the various
>> * RHS into one array for calling convenience. As such, RHS vectors are
>> * always stored contiguously, regardless of the row/col major that is
>> * chosen. This means that ldb/ldx is always independent of NRHS, and
>> * dependant on N, regardless of the row/col major setting.
>> */
>> // That is, it seems that, if B is row-major, it should be NRHS-by-N,
>> // and RHS vectors should be its rows, not columns.
>> //
>> // [comments from 'clapack_dgesv.c':]
>> /* The LU factorization used to factor A is dependent on the Order
>> * parameter, as detailed in the leading comments of clapack_dgetrf.
>> * The factored form of A is then used to solve the system of equations
>> * A * X = B.
>> * A is overwritten with the appropriate LU factorization, and B [...]
>> * is overwritten with the solution X on output.
>> */
>> // If B is row-major, solution vectors are its rows.
>> template <typename MatrA, typename MatrB, typename IVec>
>> inline
>> int gesv (MatrA& a, IVec& ipiv, MatrB& b) {
>> #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<MatrA>::matrix_structure,
>> traits::general_t
>> >::value));
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<MatrB>::matrix_structure,
>> traits::general_t
>> >::value));
>>
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<MatrA>::ordering_type,
>> typename traits::matrix_traits<MatrB>::ordering_type
>> >::value));
>> #endif
>>
>> CBLAS_ORDER const stor_ord
>> = enum_cast<CBLAS_ORDER const>
>> (storage_order<
>> #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> typename traits::matrix_traits<MatrA>::ordering_type
>> #else
>> typename MatrA::orientation_category
>> #endif
>> >::value);
>>
>> int const n = traits::matrix_size1 (a);
>> int const nrhs = stor_ord == CblasColMajor
>> ? traits::matrix_size2 (b)
>> : traits::matrix_size1 (b);
>> assert (n == traits::matrix_size2 (a));
>> assert (n == (stor_ord == CblasColMajor
>> ? traits::matrix_size1 (b)
>> : traits::matrix_size2 (b)));
>> assert (n == traits::vector_size (ipiv));
>>
>> return detail::gesv (stor_ord, n, nrhs,
>> traits::matrix_storage (a),
>> traits::leading_dimension (a),
>> traits::vector_storage (ipiv),
>> traits::matrix_storage (b),
>> traits::leading_dimension (b));
>> }
>>
>> template <typename MatrA, typename MatrB>
>> inline
>> int gesv (MatrA& a, MatrB& b) {
>> // with 'internal' pivot vector
>> #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<MatrA>::matrix_structure,
>> traits::general_t
>> >::value));
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<MatrB>::matrix_structure,
>> traits::general_t
>> >::value));
>>
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<MatrA>::ordering_type,
>> typename traits::matrix_traits<MatrB>::ordering_type
>> >::value));
>> #endif
>>
>> CBLAS_ORDER const stor_ord
>> = enum_cast<CBLAS_ORDER const>
>> (storage_order<
>> #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> typename traits::matrix_traits<MatrA>::ordering_type
>> #else
>> typename MatrA::orientation_category
>> #endif
>> >::value);
>>
>> int const n = traits::matrix_size1 (a);
>> int const nrhs = stor_ord == CblasColMajor
>> ? traits::matrix_size2 (b)
>> : traits::matrix_size1 (b);
>> assert (n == traits::matrix_size2 (a));
>> assert (n == (stor_ord == CblasColMajor
>> ? traits::matrix_size1 (b)
>> : traits::matrix_size2 (b)));
>>
>> int *ipiv = new (std::nothrow) int[n];
>> int ierr = -101;
>> // clapack_dgesv() errors:
>> // if (ierr == 0), successful
>> // if (ierr < 0), the -ierr argument had an illegal value
>> // -- we will use -101 if allocation fails
>> // if (ierr > 0), U(i-1,i-1) (or L(i-1,i-1)) is exactly zero
>>
>> if (ipiv) {
>> ierr = detail::gesv (stor_ord, n, nrhs,
>> traits::matrix_storage (a),
>> traits::leading_dimension (a),
>> ipiv,
>> traits::matrix_storage (b),
>> traits::leading_dimension (b));
>> delete[] ipiv;
>> }
>> return ierr;
>> }
>>
>> template <typename MatrA, typename MatrB>
>> inline
>> int lu_solve (MatrA& a, MatrB& b) {
>> return gesv (a, b);
>> }
>>
>>
>> // getrf(): LU factorization of A
>> // [comments from 'clapack_dgetrf.c':]
>> /* Computes one of two LU factorizations based on the setting of
>> * the Order parameter, as follows:
>> * ---------------------------------------------------------------
>> * Order == CblasColMajor
>> * Column-major factorization of form
>> * A = P * L * U
>> * where P is a row-permutation matrix, L is lower triangular with
>> * unit diagonal elements (lower trapezoidal if M > N), and U is
>> * upper triangular (upper trapezoidal if M < N).
>> * ---------------------------------------------------------------
>> * Order == CblasRowMajor
>> * Row-major factorization of form
>> * A = P * L * U
>> * where P is a column-permutation matrix, L is lower triangular
>> * (lower trapezoidal if M > N), and U is upper triangular with
>> * unit diagonals (upper trapezoidal if M < N).
>> */
>> template <typename MatrA, typename IVec>
>> inline
>> int getrf (MatrA& a, IVec& ipiv) {
>> #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<MatrA>::matrix_structure,
>> traits::general_t
>> >::value));
>> #endif
>>
>> CBLAS_ORDER const stor_ord
>> = enum_cast<CBLAS_ORDER const>
>> (storage_order<
>> #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> typename traits::matrix_traits<MatrA>::ordering_type
>> #else
>> typename MatrA::orientation_category
>> #endif
>> >::value);
>>
>> int const m = traits::matrix_size1 (a);
>> int const n = traits::matrix_size2 (a);
>> assert (traits::vector_size (ipiv) == (m < n ? m : n));
>>
>> return detail::getrf (stor_ord, m, n,
>> traits::matrix_storage (a),
>> traits::leading_dimension (a),
>> traits::vector_storage (ipiv));
>> }
>>
>> template <typename MatrA, typename IVec>
>> inline
>> int lu_factor (MatrA& a, IVec& ipiv) {
>> return getrf (a, ipiv);
>> }
>>
>>
>> // getrs(): solves a system of linear equations
>> // A * X = B or A' * X = B
>> // using the LU factorization previously computed by getrf()
>> template <typename MatrA, typename MatrB, typename IVec>
>> inline
>> int getrs (CBLAS_TRANSPOSE const Trans,
>> MatrA const& a, IVec const& ipiv, MatrB& b)
>> {
>> #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<MatrA>::matrix_structure,
>> traits::general_t
>> >::value));
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<MatrB>::matrix_structure,
>> traits::general_t
>> >::value));
>>
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<MatrA>::ordering_type,
>> typename traits::matrix_traits<MatrB>::ordering_type
>> >::value));
>> #endif
>>
>> assert (Trans == CblasNoTrans
>> || Trans == CblasTrans
>> || Trans == CblasConjTrans);
>>
>> CBLAS_ORDER const stor_ord
>> = enum_cast<CBLAS_ORDER const>
>> (storage_order<
>> #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> typename traits::matrix_traits<MatrA>::ordering_type
>> #else
>> typename MatrA::orientation_category
>> #endif
>> >::value);
>>
>> int const n = traits::matrix_size1 (a);
>> int const nrhs = stor_ord == CblasColMajor
>> ? traits::matrix_size2 (b)
>> : traits::matrix_size1 (b);
>> assert (n == traits::matrix_size2 (a));
>> assert (n == (stor_ord == CblasColMajor
>> ? traits::matrix_size1 (b)
>> : traits::matrix_size2 (b)));
>> assert (n == traits::vector_size (ipiv));
>>
>> return detail::getrs (stor_ord, Trans, n, nrhs,
>> #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
>> traits::matrix_storage (a),
>> #else
>> traits::matrix_storage_const (a),
>> #endif
>> traits::leading_dimension (a),
>> #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
>> traits::vector_storage (ipiv),
>> #else
>> traits::vector_storage_const (ipiv),
>> #endif
>> traits::matrix_storage (b),
>> traits::leading_dimension (b));
>> }
>>
>> // getrs(): solves A * X = B (after getrf())
>> template <typename MatrA, typename MatrB, typename IVec>
>> inline
>> int getrs (MatrA const& a, IVec const& ipiv, MatrB& b) {
>> return getrs (CblasNoTrans, a, ipiv, b);
>> }
>>
>> template <typename MatrA, typename MatrB, typename IVec>
>> inline
>> int lu_substitute (MatrA const& a, IVec const& ipiv, MatrB& b) {
>> return getrs (CblasNoTrans, a, ipiv, b);
>> }
>>
>>
>> // getri(): computes the inverse of a matrix A
>> // using the LU factorization previously computed by getrf()
>> template <typename MatrA, typename IVec>
>> inline
>> int getri (MatrA& a, IVec const& ipiv) {
>> #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<MatrA>::matrix_structure,
>> traits::general_t
>> >::value));
>> #endif
>>
>> CBLAS_ORDER const stor_ord
>> = enum_cast<CBLAS_ORDER const>
>> (storage_order<
>> #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> typename traits::matrix_traits<MatrA>::ordering_type
>> #else
>> typename MatrA::orientation_category
>> #endif
>> >::value);
>>
>> int const n = traits::matrix_size1 (a);
>> assert (n == traits::matrix_size2 (a));
>> assert (traits::vector_size (ipiv) == n);
>>
>> return detail::getri (stor_ord, n,
>> traits::matrix_storage (a),
>> traits::leading_dimension (a),
>> #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
>> traits::vector_storage (ipiv)
>> #else
>> traits::vector_storage_const (ipiv)
>> #endif
>> );
>> }
>>
>> template <typename MatrA, typename IVec>
>> inline
>> int lu_invert (MatrA& a, IVec& ipiv) {
>> return getri (a, ipiv);
>> }
>>
>>
>>
>> /////////////////////////////////////////////////////////////////////
>> //
>> // system of linear equations A * X = B
>> // with A symmetric or Hermitian positive definite matrix
>> //
>> /////////////////////////////////////////////////////////////////////
>>
>> #ifndef BOOST_NUMERIC_BINDINGS_ATLAS_POTRF_BUG
>> // posv(): 'driver' function
>> //
>> // [from 'dposv.f' (slightly edited):]
>> /* XPOSV computes the solution to a system of linear equations
>> * A * X = B,
>> * where A is an N-by-N symmetric/Hermitian positive definite matrix
>> * and X and B are N-by-NRHS matrices. [See also comments of gesv().]
>> *
>> * A -- On entry, the symmetric/Hermitian matrix A.
>> * If UPLO = 'U', the leading N-by-N upper triangular part of A
>> * contains the upper triangular part of the matrix A, and the
>> * strictly lower triangular part of A is not referenced.
>> * If UPLO = 'L', the leading N-by-N lower triangular part of A
>> * contains the lower triangular part of the matrix A, and the
>> * strictly upper triangular part of A is not referenced.
>> *
>> * On exit, if INFO = 0, the factor U or L from the Cholesky
>> * factorization A = U**T*U or A = L*L**T
>> * [or A = U**H*U or A = L*L**H].
>> *
>> * B -- On entry, the right hand side matrix B.
>> * On exit, if INFO = 0, the solution matrix X.
>> */
>> namespace detail {
>>
>> template <typename SymmA, typename MatrB>
>> inline
>> int posv (CBLAS_UPLO const uplo, SymmA& a, MatrB& b) {
>> #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<SymmA>::ordering_type,
>> typename traits::matrix_traits<MatrB>::ordering_type
>> >::value));
>> #endif
>>
>> CBLAS_ORDER const stor_ord
>> = enum_cast<CBLAS_ORDER const>
>> (storage_order<
>> #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> typename traits::matrix_traits<SymmA>::ordering_type
>> #else
>> typename SymmA::orientation_category
>> #endif
>> >::value);
>>
>> int const n = traits::matrix_size1 (a);
>> int const nrhs = stor_ord == CblasColMajor
>> ? traits::matrix_size2 (b)
>> : traits::matrix_size1 (b);
>> assert (n == traits::matrix_size2 (a));
>> assert (n == (stor_ord == CblasColMajor
>> ? traits::matrix_size1 (b)
>> : traits::matrix_size2 (b)));
>>
>> return posv (stor_ord, uplo, n, nrhs,
>> traits::matrix_storage (a),
>> traits::leading_dimension (a),
>> traits::matrix_storage (b),
>> traits::leading_dimension (b));
>> }
>>
>> } // detail
>>
>> template <typename SymmA, typename MatrB>
>> inline
>> int posv (CBLAS_UPLO const uplo, SymmA& a, MatrB& b) {
>> #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<SymmA>::matrix_structure,
>> traits::general_t
>> >::value));
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<MatrB>::matrix_structure,
>> traits::general_t
>> >::value));
>> #endif
>> assert (uplo == CblasUpper || uplo == CblasLower);
>> return detail::posv (uplo, a, b);
>> }
>>
>> template <typename SymmA, typename MatrB>
>> inline
>> int posv (SymmA& a, MatrB& b) {
>> #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> typedef typename traits::matrix_traits<SymmA>::value_type val_t;
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<SymmA>::matrix_structure,
>> typename traits::detail::symm_herm_t<val_t>::type
>> >::value));
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<MatrB>::matrix_structure,
>> traits::general_t
>> >::value));
>> #endif
>>
>> CBLAS_UPLO const uplo
>> = enum_cast<CBLAS_UPLO const>
>> (uplo_triang<
>> #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> typename traits::matrix_traits<SymmA>::uplo_type
>> #else
>> typename SymmA::packed_category
>> #endif
>> >::value);
>>
>> return detail::posv (uplo, a, b);
>> }
>>
>> template <typename SymmA, typename MatrB>
>> inline
>> int cholesky_solve (SymmA& a, MatrB& b) { return posv (a, b); }
>> #endif // BOOST_NUMERIC_BINDINGS_ATLAS_POTRF_BUG
>>
>>
>> // potrf(): Cholesky factorization of A
>> namespace detail {
>>
>> template <typename SymmA>
>> inline
>> int potrf (CBLAS_UPLO const uplo, SymmA& a) {
>>
>> CBLAS_ORDER const stor_ord
>> = enum_cast<CBLAS_ORDER const>
>> (storage_order<
>> #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> typename traits::matrix_traits<SymmA>::ordering_type
>> #else
>> typename SymmA::orientation_category
>> #endif
>> >::value);
>>
>> int const n = traits::matrix_size1 (a);
>> assert (n == traits::matrix_size2 (a));
>>
>> return potrf (stor_ord, uplo, n,
>> traits::matrix_storage (a),
>> traits::leading_dimension (a));
>> }
>>
>> } // detail
>>
>> template <typename SymmA>
>> inline
>> int potrf (CBLAS_UPLO const uplo, SymmA& a) {
>> #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<SymmA>::matrix_structure,
>> traits::general_t
>> >::value));
>> #endif
>> assert (uplo == CblasUpper || uplo == CblasLower);
>> return detail::potrf (uplo, a);
>> }
>>
>> template <typename SymmA>
>> inline
>> int potrf (SymmA& a) {
>> #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> typedef typename traits::matrix_traits<SymmA>::value_type val_t;
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<SymmA>::matrix_structure,
>> typename traits::detail::symm_herm_t<val_t>::type
>> >::value));
>> #endif
>>
>> CBLAS_UPLO const uplo
>> = enum_cast<CBLAS_UPLO const>
>> (uplo_triang<
>> #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> typename traits::matrix_traits<SymmA>::uplo_type
>> #else
>> typename SymmA::packed_category
>> #endif
>> >::value);
>>
>> return detail::potrf (uplo, a);
>> }
>>
>> template <typename SymmA>
>> inline
>> int cholesky_factor (SymmA& a) { return potrf (a); }
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> namespace detail
>> {
>>
>>
>>
>>
>> //! potrs specialisation for one RHS as a boost::numeric::ublas::vector
>> /*!
>> \author Paul Thompson
>> \date 06-10-06
>>
>> \todo Perhaps the input argument should be templated to accept other vector types than boost::numeric::ublas::vector<T,A> ?
>> */
>> template <typename SymmA, typename T,typename A>
>> inline
>> int potrs (CBLAS_UPLO const uplo, SymmA const& a, boost::numeric::ublas::vector<T,A> & v)
>> {
>>
>> CBLAS_ORDER const stor_ord
>> = enum_cast<CBLAS_ORDER const>
>> (storage_order<
>> #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> typename traits::matrix_traits<SymmA>::ordering_type
>> #else
>> typename SymmA::orientation_category
>> #endif
>> >::value);
>>
>> int const n = traits::matrix_size1 (a);
>> int const nrhs = 1; //! We know there's only one RHS for a vector
>> assert (n == traits::matrix_size2 (a));
>> assert (n == traits::vector_size(v));
>>
>> int ldb = traits::vector_size(v);
>>
>> #ifndef BOOST_NUMERIC_BINDINGS_ATLAS_POTRF_BUG
>> return potrs (stor_ord, uplo, n, nrhs,
>> #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
>> traits::matrix_storage (a),
>> #else
>> traits::matrix_storage_const (a),
>> #endif
>> traits::leading_dimension (a),
>> traits::vector_storage (v),
>> ldb);
>> #else // BOOST_NUMERIC_BINDINGS_ATLAS_POTRF_BUG
>> int ierr;
>> if (stor_ord == CblasColMajor)
>> ierr = potrs (stor_ord, uplo, n, nrhs,
>> #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
>> traits::matrix_storage (a),
>> #else
>> traits::matrix_storage_const (a),
>> #endif
>> traits::leading_dimension (a),
>> traits::vector_storage (v),
>> ldb);
>> else // ATLAS bug with CblasRowMajor
>> ierr = potrs_bug (stor_ord, uplo, n, nrhs,
>> #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
>> traits::matrix_storage (a),
>> #else
>> traits::matrix_storage_const (a),
>> #endif
>> traits::leading_dimension (a),
>> traits::vector_storage (v),
>> ldb);
>> return ierr;
>> #endif // BOOST_NUMERIC_BINDINGS_ATLAS_POTRF_BUG
>> }
>> //potrs detailed
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> // potrs(): solves a system of linear equations A * X = B
>> // using the Cholesky factorization computed by potrf()
>>
>>
>> template <typename SymmA, typename MatrB>
>> inline
>> int potrs (CBLAS_UPLO const uplo, SymmA const& a, MatrB& b) {
>> #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<SymmA>::ordering_type,
>> typename traits::matrix_traits<MatrB>::ordering_type
>> >::value));
>> #endif
>>
>> CBLAS_ORDER const stor_ord
>> = enum_cast<CBLAS_ORDER const>
>> (storage_order<
>> #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> typename traits::matrix_traits<SymmA>::ordering_type
>> #else
>> typename SymmA::orientation_category
>> #endif
>> >::value);
>>
>> int const n = traits::matrix_size1 (a);
>> int const nrhs = stor_ord == CblasColMajor
>> ? traits::matrix_size2 (b)
>> : traits::matrix_size1 (b);
>> assert (n == traits::matrix_size2 (a));
>> assert (n == (stor_ord == CblasColMajor
>> ? traits::matrix_size1 (b)
>> : traits::matrix_size2 (b)));
>>
>> #ifndef BOOST_NUMERIC_BINDINGS_ATLAS_POTRF_BUG
>> return potrs (stor_ord, uplo, n, nrhs,
>> #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
>> traits::matrix_storage (a),
>> #else
>> traits::matrix_storage_const (a),
>> #endif
>> traits::leading_dimension (a),
>> traits::matrix_storage (b),
>> traits::leading_dimension (b));
>> #else // BOOST_NUMERIC_BINDINGS_ATLAS_POTRF_BUG
>> int ierr;
>> if (stor_ord == CblasColMajor)
>> ierr = potrs (stor_ord, uplo, n, nrhs,
>> #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
>> traits::matrix_storage (a),
>> #else
>> traits::matrix_storage_const (a),
>> #endif
>> traits::leading_dimension (a),
>> traits::matrix_storage (b),
>> traits::leading_dimension (b));
>> else // ATLAS bug with CblasRowMajor
>> ierr = potrs_bug (stor_ord, uplo, n, nrhs,
>> #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
>> traits::matrix_storage (a),
>> #else
>> traits::matrix_storage_const (a),
>> #endif
>> traits::leading_dimension (a),
>> traits::matrix_storage (b),
>> traits::leading_dimension (b));
>> return ierr;
>> #endif // BOOST_NUMERIC_BINDINGS_ATLAS_POTRF_BUG
>> }
>>
>> } // detail
>>
>> template <typename SymmA, typename MatrB>
>> inline
>> int potrs (CBLAS_UPLO const uplo, SymmA const& a, MatrB& b) {
>> #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<SymmA>::matrix_structure,
>> traits::general_t
>> >::value));
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<MatrB>::matrix_structure,
>> traits::general_t
>> >::value));
>> #endif
>> assert (uplo == CblasUpper || uplo == CblasLower);
>> return detail::potrs (uplo, a, b);
>> }
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> template <typename SymmA, typename T, typename A>
>> inline
>> int potrs (SymmA const& a, boost::numeric::ublas::vector<T,A> & v)
>> {
>> #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> typedef typename traits::matrix_traits<SymmA>::value_type val_t;
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<SymmA>::matrix_structure,
>> typename traits::detail::symm_herm_t<val_t>::type
>> >::value));
>> #endif
>>
>> CBLAS_UPLO const uplo
>> = enum_cast<CBLAS_UPLO const>
>> (uplo_triang<
>> #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> typename traits::matrix_traits<SymmA>::uplo_type
>> #else
>> typename SymmA::packed_category
>> #endif
>> >::value);
>>
>> return detail::potrs<SymmA, T, A> (uplo, a, v);
>> }
>>
>>
>>
>>
>>
>>
>>
>>
>> template <typename SymmA, typename MatrB>
>> inline
>> int potrs (SymmA const& a, MatrB& b) {
>> #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> typedef typename traits::matrix_traits<SymmA>::value_type val_t;
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<SymmA>::matrix_structure,
>> typename traits::detail::symm_herm_t<val_t>::type
>> >::value));
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<MatrB>::matrix_structure,
>> traits::general_t
>> >::value));
>> #endif
>>
>> CBLAS_UPLO const uplo
>> = enum_cast<CBLAS_UPLO const>
>> (uplo_triang<
>> #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> typename traits::matrix_traits<SymmA>::uplo_type
>> #else
>> typename SymmA::packed_category
>> #endif
>> >::value);
>>
>> return detail::potrs (uplo, a, b);
>> }
>>
>> template <typename SymmA, typename MatrB>
>> inline
>> int cholesky_substitute (SymmA const& a, MatrB& b) { return potrs (a, b); }
>>
>>
>> #ifdef BOOST_NUMERIC_BINDINGS_ATLAS_POTRF_BUG
>> // posv(): 'driver' function
>> template <typename SymmA, typename MatrB>
>> inline
>> int posv (CBLAS_UPLO const uplo, SymmA& a, MatrB& b) {
>> int ierr = potrf (uplo, a);
>> if (ierr == 0)
>> ierr = potrs (uplo, a, b);
>> return ierr;
>> }
>>
>> template <typename SymmA, typename MatrB>
>> inline
>> int posv (SymmA& a, MatrB& b) {
>> int ierr = potrf (a);
>> if (ierr == 0)
>> ierr = potrs (a, b);
>> return ierr;
>> }
>>
>> template <typename SymmA, typename MatrB>
>> inline
>> int cholesky_solve (SymmA& a, MatrB& b) {
>> return posv (a, b);
>> }
>> #endif // BOOST_NUMERIC_BINDINGS_ATLAS_POTRF_BUG
>>
>>
>> // potri(): computes the inverse of a symmetric or Hermitian positive
>> // definite matrix A using the Cholesky factorization
>> // previously computed by potrf()
>> namespace detail {
>>
>> template <typename SymmA>
>> inline
>> int potri (CBLAS_UPLO const uplo, SymmA& a) {
>>
>> CBLAS_ORDER const stor_ord
>> = enum_cast<CBLAS_ORDER const>
>> (storage_order<
>> #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> typename traits::matrix_traits<SymmA>::ordering_type
>> #else
>> typename SymmA::orientation_category
>> #endif
>> >::value);
>>
>> int const n = traits::matrix_size1 (a);
>> assert (n == traits::matrix_size2 (a));
>>
>> return potri (stor_ord, uplo, n,
>> traits::matrix_storage (a),
>> traits::leading_dimension (a));
>> }
>>
>> } // detail
>>
>> template <typename SymmA>
>> inline
>> int potri (CBLAS_UPLO const uplo, SymmA& a) {
>> #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<SymmA>::matrix_structure,
>> traits::general_t
>> >::value));
>> #endif
>> assert (uplo == CblasUpper || uplo == CblasLower);
>> return detail::potri (uplo, a);
>> }
>>
>> template <typename SymmA>
>> inline
>> int potri (SymmA& a) {
>> #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> typedef typename traits::matrix_traits<SymmA>::value_type val_t;
>> BOOST_STATIC_ASSERT((boost::is_same<
>> typename traits::matrix_traits<SymmA>::matrix_structure,
>> typename traits::detail::symm_herm_t<val_t>::type
>> >::value));
>> #endif
>>
>> CBLAS_UPLO const uplo
>> = enum_cast<CBLAS_UPLO const>
>> (uplo_triang<
>> #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> typename traits::matrix_traits<SymmA>::uplo_type
>> #else
>> typename SymmA::packed_category
>> #endif
>> >::value);
>>
>> return detail::potri (uplo, a);
>> }
>>
>> template <typename SymmA>
>> inline
>> int cholesky_invert (SymmA& a) { return potri (a); }
>>
>>
>> } // namespace atlas
>>
>> }}}
>>
>> #endif // BOOST_NUMERIC_BINDINGS_CLAPACK_HPP
>>
>>
>> ------------------------------------------------------------------------
>>
>> 563,653d562
>> <
>> <
>> <
>> <
>> <
>> <
>> <
>> <
>> < namespace detail
>> < {
>> <
>> <
>> <
>> <
>> < //! potrs specialisation for one RHS as a boost::numeric::ublas::vector
>> < /*!
>> < \author Paul Thompson
>> < \date 06-10-06
>> <
>> < \todo Perhaps the input argument should be templated to accept other vector types than boost::numeric::ublas::vector<T,A> ?
>> < */
>> < template <typename SymmA, typename T,typename A>
>> < inline
>> < int potrs (CBLAS_UPLO const uplo, SymmA const& a, boost::numeric::ublas::vector<T,A> & v)
>> < {
>> <
>> < CBLAS_ORDER const stor_ord
>> < = enum_cast<CBLAS_ORDER const>
>> < (storage_order<
>> < #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> < typename traits::matrix_traits<SymmA>::ordering_type
>> < #else
>> < typename SymmA::orientation_category
>> < #endif
>> < >::value);
>> <
>> < int const n = traits::matrix_size1 (a);
>> < int const nrhs = 1; //! We know there's only one RHS for a vector
>> < assert (n == traits::matrix_size2 (a));
>> < assert (n == traits::vector_size(v));
>> <
>> < int ldb = traits::vector_size(v);
>> <
>> < #ifndef BOOST_NUMERIC_BINDINGS_ATLAS_POTRF_BUG
>> < return potrs (stor_ord, uplo, n, nrhs,
>> < #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
>> < traits::matrix_storage (a),
>> < #else
>> < traits::matrix_storage_const (a),
>> < #endif
>> < traits::leading_dimension (a),
>> < traits::vector_storage (v),
>> < ldb);
>> < #else // BOOST_NUMERIC_BINDINGS_ATLAS_POTRF_BUG
>> < int ierr;
>> < if (stor_ord == CblasColMajor)
>> < ierr = potrs (stor_ord, uplo, n, nrhs,
>> < #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
>> < traits::matrix_storage (a),
>> < #else
>> < traits::matrix_storage_const (a),
>> < #endif
>> < traits::leading_dimension (a),
>> < traits::vector_storage (v),
>> < ldb);
>> < else // ATLAS bug with CblasRowMajor
>> < ierr = potrs_bug (stor_ord, uplo, n, nrhs,
>> < #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
>> < traits::matrix_storage (a),
>> < #else
>> < traits::matrix_storage_const (a),
>> < #endif
>> < traits::leading_dimension (a),
>> < traits::vector_storage (v),
>> < ldb);
>> < return ierr;
>> < #endif // BOOST_NUMERIC_BINDINGS_ATLAS_POTRF_BUG
>> < }
>> < //potrs detailed
>> <
>> <
>> <
>> <
>> <
>> <
>> <
>> <
>> <
>> <
>> <
>> <
>> 656c565
>> <
>> ---
>>> namespace detail {
>> 742,788d650
>> <
>> <
>> <
>> <
>> <
>> <
>> <
>> <
>> < //! potrs specialisation for one RHS as a boost::numeric::ublas::vector
>> < /*!
>> < \author Paul Thompson
>> < \date 06-10-06
>> <
>> < \todo Perhaps the input argument should be templated to accept other vector types than boost::numeric::ublas::vector<T,A> ?
>> < */
>> < template <typename SymmA, typename T, typename A>
>> < inline
>> < int potrs (SymmA const& a, boost::numeric::ublas::vector<T,A> & v)
>> < {
>> < #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK
>> < typedef typename traits::matrix_traits<SymmA>::value_type val_t;
>> < BOOST_STATIC_ASSERT((boost::is_same<
>> < typename traits::matrix_traits<SymmA>::matrix_structure,
>> < typename traits::detail::symm_herm_t<val_t>::type
>> < >::value));
>> < #endif
>> <
>> < CBLAS_UPLO const uplo
>> < = enum_cast<CBLAS_UPLO const>
>> < (uplo_triang<
>> < #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
>> < typename traits::matrix_traits<SymmA>::uplo_type
>> < #else
>> < typename SymmA::packed_category
>> < #endif
>> < >::value);
>> <
>> < return detail::potrs<SymmA, T, A> (uplo, a, v);
>> < }
>> <
>> <
>> <
>> <
>> <
>> <
>> <
>> <
>>
>>
>> ------------------------------------------------------------------------
>>
>> _______________________________________________
>> ublas mailing list
>> ublas_at_[hidden]
>> http://lists.boost.org/mailman/listinfo.cgi/ublas
> _______________________________________________
> ublas mailing list
> ublas_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/ublas