Boost logo

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