Re: [ublas] expression templates,question 2

From: Michael Stevens (mail_at_[hidden])
Date: 2005-07-06 09:15:19

Hi Dima,

I think it is easiest to answer this questions first.
On Wednesday 06 July 2005 11:07, Dima Sorkin wrote:
> Hi.
>  The second question : extendability.
>  I find ublas-library VERY useful,but
> sometimes I need to write a small
> extensions for it (for example:
> multiplication of matrices of my-classes,
> concatenation,etc.). I always do it the
> simpliest way (no performance considerations
> at all).
>
>  However,I was trying (as an excersize) to add
> ublas new functionality (listed above) in a
> "ublas way".
>
>  How do I write in a "ublas way"?
>
> (with ETs or without, this really does not matter)

The uBLAS way is with ETs. uBLAS is built on 3 pillars. In order of the most accessible first.
1. Containers
2. Expressions (ETs)
3. Expression evaluation

>  My general request is :
> Could the internals of ublas be documented ?
> (the methods of work,not small things).
> I see a lack of documentation of how ublas works.

1. Is formally documented in the Concepts documentation
2.3. Are difficult. Although the interfaces (Concepts such as Expressions and Iterators) are formally documented
this does give enough information to extend uBLAS. Particularly the return values are a problem
I suspect the only way to make progress with documenting points hese internal uBLAS's mechanisms is everyone contributes a bit. We just need information to get people started.
Beware however that a lot of this is internal to the current uBLAS implementation and should not be set in stone.

> It is not simple to understand the mechanizm
> by parsing the bare code.

Personally I think a few examples would go a log long way.

If you get further with you experiments it would be good to contribute some basic info to the "Effective uBLAS" wiki.

To get you started here is an example from my own code. This defines a function prod_SPD that implements a quadratic product.
It is easy to use as you can write
        X = prod_SPD (X,S);
or indeed completely generally
        X = partial_expressionA prod_SPD (expression2, expression3) partial_expressionB;

template<class EX, class ES, class ET> inline
typename prod_expression_result<EX,ET>::E1E2T_type
 prod_SPD (const ublas::matrix_expression<EX>& X, const ublas::matrix_expression<ES>& S, ublas::matrix_expression<ET>& XStemp)
/*
 * Symmetric Positive (Semi) Definate product: X*(X*S)', XStemp = X*S
 *  Result symmetric if S is symmetric
 */
{
        return prod( X, trans(prod(X,S,XStemp())) );
}

The hard bit is that C++ requires a return type for this function. This is where all the nasty implementation details of ETs come in. It is
best to define the a help class that holdes everything together.

struct prod_expression_result
{	// Provide ET result E1E2T_type of prod(matrix_expression<E1>,trans(matrix_expression<E2>)
        typedef BOOST_UBLAS_TYPENAME ublas::matrix_unary2_traits<E2, ublas::scalar_identity<BOOST_UBLAS_TYPENAME E2::value_type> >::result_type  E2T_type;
        typedef BOOST_UBLAS_TYPENAME ublas::matrix_matrix_binary_traits<BOOST_UBLAS_TYPENAME E1::value_type, E1,
                                        BOOST_UBLAS_TYPENAME E2T_type::value_type, E2T_type>::result_type  E1E2T_type;
};

How do you work out what the return type is. You need to look at the individual return types for all the operators in the expression. These must then be nested
to provide a compound type.

Hope this gets you started.

Michael
-- 
___________________________________
Michael Stevens Systems Engineering
34128 Kassel, Germany
Phone/Fax: [phone]
Navigation Systems, Estimation  and
                 Bayesian Filtering
    http://bayesclasses.sf.net
___________________________________