Boost logo

Ublas :

From: Michael Stevens (mail_at_[hidden])
Date: 2005-08-26 09:02:59

On Saturday 20 August 2005 12:49, Gunter Winkler wrote:
> Am Freitag, 19. August 2005 22:00 schrieb Shawn D. Pautz:
> > Is there a way to use inner_prod (or something similar) for vectors
> > of more complicated types?
> > What's the best way to perform the above operation without resorting
> > to writing loops?
> Currently this is not directly supported. The inner_prod uses the
> operator * (x, y) function to compute the products which should work
> for scalar times Matrix. However, the return type of the expression can
> not be calculated automatically. You have to create a specialization of
> promote_traits.
> typedef double SCALAR;
> typedef matrix<SCALAR> MATRIX;
> template<>
> struct promote_traits< SCALAR, MATRIX > {
> typedef MATRIX base_type;
> typedef MATRIX promote_type;
> };
> Unfortunatly this only fixes the return type deduction. The second
> problem is that the return type has to be constructible from (int) 0.
> This can be achieved by using a type derived from MATRIX with an
> appropriate constructor.
> class my_matrix : public MATRIX {
> public: my_matrix(const size_t zero) :
> MATRIX( zero_matrix<double>(5,5) ) { }
> }
> I hope you got the idea.

It is worth noting that the two problems have different root causes. The need
for promote_traits is a restriction of the current return type deduction.
Some time in the future (when C++ supports this) this restriction will go

The zero construction problem is fundamental and is caused by the fundamental
properties of matrices. uBLAS provides container types based on the objects
properties : dense/sparse, vector/matrix etc. It does not associate a type
with order (or dimension) of these objects. Therefore it is not possible to
construct a conformant object based on type.

The current solution is to require that the element type defines a 'zero'
based on calling the constructor with (0). This follows from the properties
of C++ builtin types and is easy to achieve for scalar types.

In the case of using vectors as the elements involved in an inner_prod we need
to enforce this and effectively encode the order of the matrix in the type

Based on Gunter's remarks I have put together the following example code which
is complete.

 * Test nested LA objects for inner_prod
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>

#include <iostream>

using namespace boost::numeric::ublas;

template <class T>
class my_matrix : public matrix<T>
        my_matrix () :
        my_matrix (const int s) :
                matrix<double>( scalar_matrix<T>(5,5) )

namespace boost { namespace numeric { namespace ublas {
// Nasty intrusion into the ublas namespace!
        template <class SCALAR>
        struct promote_traits<SCALAR, my_matrix<SCALAR> > {
            typedef my_matrix<SCALAR> base_type;
            typedef base_type promote_type;

int main()
        vector<double> vec;
        vector<my_matrix<double> > vecOfMatrices;
        my_matrix<double> result;
        result = inner_prod (vec, vecOfMatrices);

> mfg
> Gunter
> PS: How far is the work for type_traits<MATRIX>::zero ?

I can't say I like the idea of type_traits<>:zero . In fact I dislike the
whole type_traits mechanism! I think it is: error prone, hard to use, hard to
document and lacks 'encapsulation'.

So I don't think I am like to put forward a patch to implement it! Of course,
it would be good if someone who does feel this is the way to go put something

If type_traits are the way to go, then something like
  type_traits<T>::clone_zero (const &T);
may be better then just zero.
This would be much more useful in the context of implementing algorithms such
as inner_prod.

Alternatively it may be possible to solve the problem by implementing uBLAS so
the zero elements are not constructed in any algorithms. Using the operation
        Z = 0* A where 0* multiples elements of A by the scalar 0
we can generate Z such that its assignment and additively conformant with A.

Many algorithms could use this or simply special case the first element in a
loop and use copy construction to generate a conformant temporary. Of course
it would remain to be seen if compilers can optimise so there is no overhead
in the case of scalar elements. This could be assisted with special cases for
bultin types.


Michael Stevens Systems Engineering
34128 Kassel, Germany
Phone/Fax: +49 561 5218038
Navigation Systems, Estimation  and
                 Bayesian Filtering