|
|
Ublas : |
From: Ian Fellows (ifellows_at_[hidden])
Date: 2007-05-30 20:23:39
I have also needed to do element wise arithmetic in addition to matrix
arithmetic. ublas provides element wise multiplication, division, and
addition, as well as scalar multiplication and division. I don't know of any
built in function to multiply/add the rows/columns of a matrix by a vector,
or to add a scalar to a to a vector/matrix.
thank you very much for your useful dialog. I have incorporated "noalias"
into a number of crude functions that I had written, however I ran into
trouble in two places. Note that I am a beginner in C++ so don't wince to
much my work. I'm planning on making the functions templates after the
numerics have been settled. When doing a row/column wise product assign with
noalias, I get the error :
error C2676: binary '*=' : 'boost::numeric::ublas::noalias_proxy<C>' does
not define this operator or a conversion to a type acceptable to the
predefined operator
here is the function :
void row_product_assign(matrix<float> &mat,vector<float> &vec)
{
#ifndef NDEBUG
if(vec.size()!=mat.size2())
std::cout<<"Error: In row_product_assign. vector length not equal to
matrix.\n";
#endif
for(unsigned i = 0; i < mat.size2 (); ++ i)
{
//I get errors here with noalias : error C2676: binary '*=' :
'boost::numeric::ublas::noalias_proxy<C>' does not define this operator or a
conversion to a type acceptable to the predefined operator
(column(mat,i))*=vec(i);
}
}
Also, I was wondering if anyone had any suggestions on how to add a scalar
to a vector/matrix efficiently? My crude attempt is:
vector<float> element_sum(const vector<float> &vec, const float &scalar)
{
vector<float> target_vec(vec.size());
for(unsigned i = 0; i < vec.size (); ++ i)
{
//error with noalias
target_vec(i) = vec(i)+scalar;
}
return (target_vec);
}
Any other suggestions/comments about how to improve the speed of these
operations would be greatly appreciated. Would this functionality be useful
to add to ublas?
Thanks
Ian
p.s. Below is the full code for the row/column wise arithmetic and scalar
addition
void row_product_assign_alias(matrix<float> &mat,vector<float> &vec)
{
#ifndef NDEBUG
if(vec.size()!=mat.size2())
std::cout<<"Error: In row_product_assign. vector length not equal to
matrix.\n";
#endif
for(unsigned i = 0; i < mat.size2 (); ++ i)
{
column(mat,i)*=vec(i);
}
}
void row_product_assign(matrix<float> &mat,vector<float> &vec)
{
#ifndef NDEBUG
if(vec.size()!=mat.size2())
std::cout<<"Error: In row_product_assign. vector length not equal to
matrix.\n";
#endif
for(unsigned i = 0; i < mat.size2 (); ++ i)
{
//I get errors here with noalias : error C2676: binary '*=' :
'boost::numeric::ublas::noalias_proxy<C>' does not define this operator or a
conversion to a type acceptable to the predefined operator
(column(mat,i))*=vec(i);
}
}
void col_product_assign_alias(matrix<float> &mat,vector<float> &vec)
{
#ifndef NDEBUG
if(vec.size()!=mat.size1())
std::cout<<"Error: In row_product_assign. vector length not equal to
matrix.\n";
#endif
for(unsigned i = 0; i < mat.size1 (); ++ i)
{
row(mat,i) *= vec(i);
}
}
void col_product_assign(matrix<float> &mat,vector<float> &vec)
{
#ifndef NDEBUG
if(vec.size()!=mat.size1())
std::cout<<"Error: In row_product_assign. vector length not equal to
matrix.\n";
#endif
for(unsigned i = 0; i < mat.size1 (); ++ i)
{
//Error here with noalias error C2676: binary '*=' :
'boost::numeric::ublas::noalias_proxy<C>' does not define this operator or a
conversion to a type acceptable to the predefined operator
(row(mat,i)) *= vec(i);
}
}
void row_product_alias(matrix<float> &mat,vector<float> &vec,matrix<float>
&target_matrix)
{
#ifndef NDEBUG
if(vec.size()!=mat.size2())
std::cout<<"Error: In row_product_assign. vector length not equal to
matrix.\n";
#endif
for(unsigned i = 0; i < mat.size2 (); ++ i)
{
column(target_matrix,i)=column(mat,i)*vec(i);
}
}
void row_product(matrix<float> &mat,vector<float> &vec,matrix<float>
&target_matrix)
{
#ifndef NDEBUG
if(vec.size()!=mat.size2())
std::cout<<"Error: In row_product_assign. vector length not equal to
matrix.\n";
#endif
for(unsigned i = 0; i < mat.size2 (); ++ i)
{
noalias(column(target_matrix,i))=column(mat,i)*vec(i);
}
}
void col_product_alias(matrix<float> &mat, vector<float> &vec, matrix<float>
&target_matrix)
{
#ifndef NDEBUG
if(vec.size()!=mat.size1())
std::cout<<"Error: In row_product_assign. vector length not equal to
matrix.\n";
#endif
for(unsigned i = 0; i < mat.size1 (); ++ i)
{
row(target_matrix,i) = row(mat,i)*vec(i);
}
}
void col_product(matrix<float> &mat, vector<float> &vec, matrix<float>
&target_matrix)
{
#ifndef NDEBUG
if(vec.size()!=mat.size1())
std::cout<<"Error: In row_product_assign. vector length not equal to
matrix.\n";
#endif
for(unsigned i = 0; i < mat.size1 (); ++ i)
{
noalias(row(target_matrix,i)) = row(mat,i)*vec(i);
}
}
matrix<float> row_product_alias(const matrix<float> &mat,const
vector<float> &vec)
{
matrix<float> target_matrix(mat.size1(),mat.size2());
#ifndef NDEBUG
if(vec.size()!=mat.size2())
std::cout<<"Error: In row_product_assign. vector length not equal to
matrix.\n";
#endif
for(unsigned i = 0; i < mat.size2 (); ++ i)
{
column(target_matrix,i)=column(mat,i)*vec(i);
}
return(target_matrix);
}
matrix<float> row_product(const matrix<float> &mat,const vector<float> &vec)
{
matrix<float> target_matrix(mat.size1(),mat.size2());
#ifndef NDEBUG
if(vec.size()!=mat.size2())
std::cout<<"Error: In row_product_assign. vector length not equal to
matrix.\n";
#endif
for(unsigned i = 0; i < mat.size2 (); ++ i)
{
noalias(column(target_matrix,i))=column(mat,i)*vec(i);
}
return(target_matrix);
}
matrix<float> col_product_alias(const matrix<float> &mat, const
vector<float> &vec)
{
matrix<float> target_matrix(mat.size1(),mat.size2());
#ifndef NDEBUG
if(vec.size()!=mat.size1())
std::cout<<"Error: In row_product_assign. vector length not equal to
matrix.\n";
#endif
for(unsigned i = 0; i < mat.size1 (); ++ i)
{
row(target_matrix,i) = row(mat,i)*vec(i);
}
return (target_matrix);
}
matrix<float> col_product(const matrix<float> &mat, const vector<float>
&vec)
{
matrix<float> target_matrix(mat.size1(),mat.size2());
#ifndef NDEBUG
if(vec.size()!=mat.size1())
std::cout<<"Error: In row_product_assign. vector length not equal to
matrix.\n";
#endif
for(unsigned i = 0; i < mat.size1 (); ++ i)
{
noalias(row(target_matrix,i)) = row(mat,i)*vec(i);
}
return (target_matrix);
}
matrix<float> element_row_sum_alias(const matrix<float> &mat,const
vector<float> &vec)
{
matrix<float> target_matrix(mat.size1(),mat.size2());
#ifndef NDEBUG
if(vec.size()!=mat.size1())
std::cout<<"Error: In element_row_sum. vector length not equal to
matrix.\n";
#endif
for(unsigned i = 0; i < mat.size1 (); ++ i)
{
row(target_matrix,i)=row(mat,i)+vec;
}
return(target_matrix);
}
matrix<float> element_row_sum(const matrix<float> &mat,const vector<float>
&vec)
{
matrix<float> target_matrix(mat.size1(),mat.size2());
#ifndef NDEBUG
if(vec.size()!=mat.size1())
std::cout<<"Error: In element_row_sum. vector length not equal to
matrix.\n";
#endif
for(unsigned i = 0; i < mat.size1 (); ++ i)
{
noalias(row(target_matrix,i))=row(mat,i)+vec;
}
return(target_matrix);
}
matrix<float> element_col_sum_alias(const matrix<float> &mat, const
vector<float> &vec)
{
matrix<float> target_matrix(mat.size1(),mat.size2());
#ifndef NDEBUG
if(vec.size()!=mat.size2())
std::cout<<"Error: In element_col_sum. vector length not equal to
matrix.\n";
#endif
for(unsigned i = 0; i < mat.size2 (); ++ i)
{
column(target_matrix,i) = column(mat,i)+vec;
}
return (target_matrix);
}
matrix<float> element_col_sum(const matrix<float> &mat, const vector<float>
&vec)
{
matrix<float> target_matrix(mat.size1(),mat.size2());
#ifndef NDEBUG
if(vec.size()!=mat.size2())
std::cout<<"Error: In element_col_sum. vector length not equal to
matrix.\n";
#endif
for(unsigned i = 0; i < mat.size2 (); ++ i)
{
noalias(column(target_matrix,i)) = column(mat,i)+vec;
}
return (target_matrix);
}
matrix<float> element_sum(const matrix<float> &mat, const float &scalar)
{
matrix<float> target_matrix(mat.size1(),mat.size2());
for(unsigned i = 0; i < mat.size1 (); ++ i)
{
for(unsigned j = 0; j < mat.size2(); ++j)
{
//error with noalias
target_matrix(i,j) = mat(i,j)+scalar;
}
}
return (target_matrix);
}
vector<float> element_sum(const vector<float> &vec, const float &scalar)
{
vector<float> target_vec(vec.size());
for(unsigned i = 0; i < vec.size (); ++ i)
{
//error with noalias
target_vec(i) = vec(i)+scalar;
}
return (target_vec);
}
-----Original Message-----
From: ublas-bounces_at_[hidden]
[mailto:ublas-bounces_at_[hidden]]On Behalf Of Gunter Winkler
Sent: Wednesday, May 30, 2007 3:13 PM
To: ublas mailing list
Subject: Re: [ublas] Matrix-vector and vector-vector operations
Am Mittwoch, 30. Mai 2007 21:37 schrieb Albert Strasheim:
> Do you have any other suggestions for how to accomplish this without
> needing extra memory?
Another note. The noalias speeds up both version:
* without noalias() (g++ -O2 -DNDEBUG)
for loop: 1.4 s
outer_prod: failed
* with noalias()
for loop: 1.1 s
outer_prod: 1.1 s
(I see _no_ time difference)
#include <boost/timer.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector.hpp>
using namespace boost::numeric::ublas;
int main(int argc, char* argv[])
{
matrix<double> A(scalar_matrix<double>(196608, 1024, 1.0));
vector<double> b(scalar_vector<double>(A.size2(),1.0));
scalar_vector<double> e(A.size1(), 1.0);
boost::timer t;
t.restart();
#if 1
noalias( A ) += outer_prod(e, b);
#else
for (std::size_t i = 0; i < A.size1(); ++i)
{
noalias( row(A, i) ) += b;
}
#endif
std::cout << "time " << t.elapsed() << " s\n";
return 0;
}
mfg
Gunter