|
|
Ublas : |
Subject: [ublas] How to implement functions that take matrix_expressions as arguments?
From: Heiko Schroeder (heiko_at_[hidden])
Date: 2009-10-11 06:12:19
Dear uBLAS developers and community!
Recently, I have decided to implement some numerics in C++ and after
comparing some linear algebra libraries I finally decided to use uBLAS.
Reasons were that I appreciate the high level of abstraction of uBLAS,
the nice design and the possibility to extend it to my own needs.
Unfortunately, the last part is what is giving me quite a hard time.
Certainly, this has to do with the fact that I'm still quite new to C++
and a bloody rookie to template metaprogramming. So, I appreciate any help.
But let me come to the concrete problem: I am trying to implement some
functions that I need to use. Those are the commutator (anticommutator),
trace (partial trace) and outer product of matrix expressions and I'd
like to implement them as general as possible in the spirit of the
matrix expression concept implemented in uBLAS. Moreover, I'd like to be
able to use them as part of more complicated expressions in the form of,
e.q., result = trace(matrix1) - trace(commute(matrix2,
matrix4-matrix5)). And now, I am already hitting a wall trying to
implement the commutator.
My first primitive approach was to define
/* \brief Compute the commutator of two matrix expressions. */
template<class M1, class M2, class M3>
BOOST_UBLAS_INLINE
M3 &
commute(M3 & m3,
const M1 & m1,
const M2 & m2)
{
return m3 = prod(m1, m2) - prod(m2, m1);
}
in the style of operations.hpp. This works, but I have to keep track of
the correct result types by the first argument (or so I understand the
need for putting M3 & there), which is against the use how I envisage it.
Then I tried to implement something along the lines
template <class E1, class E2>
BOOST_UBLAS_INLINE
???::result_type
commute(const matrix_expression<E1> & e1,
const matrix_expression<E2> & e2) {
typedef ???::expression_type expression_type;
return expression_type(e1(), e2());
}
and noticed that I have to determine the ??? type. Since the commutator
uses operations already defined in uBLAS, I thought it should be
possible to simply reuse the existing expression types and tried the
following, using templated typedefs:
template<class E1, class E2>
struct matrix_commutator_traits {
/* \brief templated typedef for matrix expression traits type of a
matrix product */
template<class A, class B>
struct prod {
typedef ublas::matrix_matrix_binary_traits<typename A::value_type, A,
typename B::value_type, B> type;
};
/* \brief templated typedef for matrix expression traits type of
element-wise matrix expression difference */
template<class A, class B>
struct diff {
typedef typename ublas::matrix_binary_traits<A,
B,
ublas::scalar_minus<typename A::value_type,
typename B::value_type> > type;
};
typedef typename diff<typename diff<E1, E2>::type::result_type,
typename diff<E2, E1>::type::result_type>::type type;
}; //struct
/* \brief commutator of two matrix expressions */
template <class E1, class E2>
BOOST_UBLAS_INLINE
typename matrix_commutator_traits<E1, E2>::type::result_type
commute(const ublas::matrix_expression<E1> & e1,
const ublas::matrix_expression<E2> & e2 )
{
typedef typename matrix_commutator_traits<E1,
E2>::type::expression_type expression_type;
return expression_type(e1 (), e2 ());
}
But it didn't work: error: no matching function for call to <complicated
template type>::matrix_binary(..., ...). (complete error message is
attached).
Now I'm wondering whether I'm doing it wrong on a fundamental level or
whether I'm just missing some detail. Any help is appreciated.
Sorry for the rather lengthy mail and best regards,
Heiko Schroeder
Using commute() in the following code snippet:
ublas::matrix<complex<double> > sigma[4];
//... defining the sigmas as 2x2 matrices
cout << commute(sigma[1], sigma[2]) << endl;
I get the following error:
ublas_extensions.cpp: In function ‘typename matrix_commutator_traits<E1, E2>::type::result_type commute(const boost::numeric::ublas::matrix_expression<E>&, const boost::numeric::ublas::matrix_expression<E2>&) [with E1 = boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >, E2 = boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >]’:
ublas_extensions.cpp:123: instantiated from here
ublas_extensions.cpp:98: error: no matching function for call to ‘boost::numeric::ublas::matrix_binary<boost::numeric::ublas::matrix_binary<boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >, boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >, boost::numeric::ublas::scalar_minus<std::complex<double>, std::complex<double> > >, boost::numeric::ublas::matrix_binary<boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >, boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >, boost::numeric::ublas::scalar_minus<std::complex<double>, std::complex<double> > >, boost::numeric::ublas::scalar_minus<std::complex<double>, std::complex<double> > >::matrix_binary(const boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >&, const boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >&)’
/usr/include/boost/numeric/ublas/matrix_expression.hpp:1710: note: candidates are: boost::numeric::ublas::matrix_binary<E1, E2, F>::matrix_binary(const E1&, const E2&) [with E1 = boost::numeric::ublas::matrix_binary<boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >, boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >, boost::numeric::ublas::scalar_minus<std::complex<double>, std::complex<double> > >, E2 = boost::numeric::ublas::matrix_binary<boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >, boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >, boost::numeric::ublas::scalar_minus<std::complex<double>, std::complex<double> > >, F = boost::numeric::ublas::scalar_minus<std::complex<double>, std::complex<double> >]
/usr/include/boost/numeric/ublas/matrix_expression.hpp:1684: note: boost::numeric::ublas::matrix_binary<boost::numeric::ublas::matrix_binary<boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >, boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >, boost::numeric::ublas::scalar_minus<std::complex<double>, std::complex<double> > >, boost::numeric::ublas::matrix_binary<boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >, boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >, boost::numeric::ublas::scalar_minus<std::complex<double>, std::complex<double> > >, boost::numeric::ublas::scalar_minus<std::complex<double>, std::complex<double> > >::matrix_binary(const boost::numeric::ublas::matrix_binary<boost::numeric::ublas::matrix_binary<boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >, boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >, boost::numeric::ublas::scalar_minus<std::complex<double>, std::complex<double> > >, boost::numeric::ublas::matrix_binary<boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >, boost::numeric::ublas::matrix<std::complex<double>, boost::numeric::ublas::basic_row_major<long unsigned int, long int>, boost::numeric::ublas::unbounded_array<std::complex<double>, std::allocator<std::complex<double> > > >, boost::numeric::ublas::scalar_minus<std::complex<double>, std::complex<double> > >, boost::numeric::ublas::scalar_minus<std::complex<double>, std::complex<double> > >&)