Re: [ublas] Returning vectors or matrices from user functions

Subject: Re: [ublas] Returning vectors or matrices from user functions
From: Gunter Winkler (guwi17_at_[hidden])
Date: 2009-09-09 17:38:25

Hello Mark,

I think I can answer some of your questions

Mark Johnson schrieb:
> Thanks Jesse for your helpful advice with this.  The Curiously
> Recurring Template Pattern is very clever!
Yes, but please keep in mind to always pass matirx/vector expressions as
(const) references. The const correctness also influences the performance.

>
> I also have two further questions -- what types should one use in user
> functions that will either return a vector or a matrix, and what types
> should one use for arguments to functions that will mutate a vector or
> matrix in place?
In this case you can always return the vector/matrix by value. If you do
not want to rely on the compiler to remove temporaries in the following code

matrix<double> A( generate_matrix_A(...) );
// or
matrix<double> A =  generate_matrix_A(...) ;

the you can explicitly call assign_temporary( rhs ), e.g. like

matrix<double> A(0,0);
A.assign_temporary( generate_matrix_A(...) );

if you want to modify the value then you simple use a non-const
reference to a matrix_expression, e.g.,

template < class E >
void fill_matrix(matrix_expression<E> & me) {
  E & matrix = me();
  // do something with matrix
}
>
> For function return values, the main issue is how to specify the
> return type to avoid needless copying, I think.  Would it be
> reasonable to return a vector_expression or a const vector_expression,
> for example?  I think this would be safe if these referred to vectors
> or matrices that are instantiated in the function call.  But perhaps
> it's also possible that a modern compiler would optimize away the
> unnecessary copy operations even if one returns ordinary vectors or
> matrices?
g++ 4 is quite smart on optimizing the construction. However you have
always the option to write a custom expression similar to the
expressions generated by trans(), prod(), operator+ ... You should not
return matrix_expressions by value. At least I had some strange issues here.
>
> For function arguments that will be modified in place, I think the
> issue is whether there is a way to define the types so that non-const
> vector expressions can be conveniently passed to the function.
this is no problem: the free function project() is very handy here:

t = norm_1( project( const_vec, range(0,5) ) );
project( vec, range(0,5) ) = scalar_vector<double>(5.1.0);
// even better
noalias( project( vec, range(0,5) ) ) = scalar_vector<double>(5.1.0);

>
> To give an artificially simple example, suppose I need to normalize
> each row in a matrix m (i.e., scale the entries so that the sum of the
> values in each row is one).  One way to do this is to write a
> "my_normalize" function that operates on a single vector, and then
> call that function on each row of the matrix m, i.e.,
> my_normalize(row(m, i)).
>
can you paste a few lines of code?

mfg
Gunter