|
|
Ublas : |
From: Michael Stevens (mail_at_[hidden])
Date: 2005-10-19 11:39:44
On Wednesday, 19. October 2005 17:44, Shawn D. Pautz wrote:
> Karl,
>
> Good point. I had tried to distill the code down; perhaps I went too
> far. Here's a less distilled version that I can't compile:
>
> template <class Matrix>
> typename Matrix::value_type determinant(const Matrix& A)
> {
> assert(A.size1() == A.size2());
>
> typedef typename Matrix::value_type value_type;
> typedef typename Matrix::size_type size_type;
>
> Matrix LU = A; // mutable type?
> ublas::permutation_matrix<size_type> pivots(LU.size1());
> if (lu_factorize(LU, pivots) != 0) // this line won't compile in
> some situations
> return 0;
>SNIP
> return det;
> }
>
SNIP
> determinant(ublas::prod(A,B)); // won't compile
> determinant(ublas::matrix<double>(ublas::prod(A,B))); // compiles
> determinant(static_cast<ublas::matrix<double>
>
> >(ublas::prod(A,B))); // compiles
>
> ublas::matrix<double> tempMatrix(ublas::prod(A,B)); // compiles
> determinant(tempMatrix); // compiles
> }
>
> I note that only one of the lines above will not compile now; thanks for
> catching the const issue. But what do I need to do to get the most
> compact line above to compile?
> Is the problem in my local template
> code?
Yes! Hopefully I can explain what goes wrong here.
First thing to do is to make the name of the function template parameter more meaningful...
template <class CouldbeAnything>
typename Matrix::value_type determinant(const CouldbeAnything>& A)
{
This is important, calling it 'Matrix' is very missleading! The actual type is often something very different. In the case of
the 'prod' function its return value is an instance of an expression template. The expression template is a type that represents the operation.
The instance will hold generalised references to the operands (A,B).
Since the result of 'prod' is a Matrix Expression the actual type will be
matrix_expression<something_horribly_complex>
You can do several with this. Principly you can iterate over it's elements, computing the values as you go along.
The matrix_expression instance is however not mutable, just as 'a+b' is not mutable.
The determinate function then fails later because the line
CouldbeAnything LU = A; // mutable type?
copy constructs an imputable type.
uBLAS has a much better way to write this kind of function. You should use...
template <class E>
typename matrix_expression<E>::value_type determinant (const matrix_expression<E>& A)
{
This will match any Matrix Expression. In this case you will need to create a temporary
Matrix LU = A;
If you want to avoid the inefficiency of this copy for Matrix types you can specialises using ..
template <class C>
typename matrix_container<C>::value_type determinant (matrix_container<C>& A)
{
> Is there a problem in the ublas conversions? More generally,
> should I be able to string together various matrix operations, not just
> the basic mathematical ones like operator+ but algorithms with a
> matrix-like return type?
The issue of return values is difficult in C++. For 'determinate' it is OK as it is simply the value_type.
If you want to create new operations that return Matrix or Vector Expression then you have to go through several hoops.
The easist way for me to explain this with my own code...
/*
* prod_SPD - uBLAS Expression templates for X*X'
*/
template <class E1, class E2>
struct prod_expression_result
{ // Provide ET result E1E2T_type of prod(matrix_expression<E1>,trans(matrix_expression<E2>)
typedef typename ublas::matrix_unary2_traits<E2, ublas::scalar_identity<typename E2::value_type> >::result_type E2T_type;
typedef typename ublas::matrix_matrix_binary_traits<typename E1::value_type, E1,
typename E2T_type::value_type, E2T_type>::result_type E1E2T_type;
};
template <class E>
typename prod_expression_result<E,E>::E1E2T_type
prod_SPD (const ublas::matrix_expression<E>& X)
/*
* Symmetric Positive (Semi) Definate product: X*X'
*/
{
return prod( X, trans(X));
}
As you can see 90% of the code is in deducing the return type. This cannot be done automatically in the current C++ standard.
Hope this makes sense,
Michael
--
___________________________________
Michael Stevens Systems Engineering
34128 Kassel, Germany
Phone/Fax: +49 561 5218038
Navigation Systems, Estimation and
Bayesian Filtering
http://bayesclasses.sf.net
___________________________________