Re: [ublas] matrix proxy of vector

From: James N. Knight (nate_at_[hidden])
Date: 2006-10-13 01:30:15

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