Re: [ublas] inner_prod

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 
away.

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 
'my_matrix'.

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>
{
public:
        my_matrix () : 
                matrix<double>(5,5)
        {}
        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 
forward.

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

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