Re: [ublas] question regarding elementwise matrix product

Subject: Re: [ublas] question regarding elementwise matrix product
From: Kaveh Kohan (kaveh.kohan_at_[hidden])
Date: 2009-05-29 12:10:22

Thanks Kim,

I tried it and it is as fast as (slightly faster than) MATLAB. For the record, here is the program and the flags I used:

# include <boost/numeric/ublas/matrix_expression.hpp>
# include <sys/time.h>
# include <vector>
# include <algorithm>
#include <iostream>
//#include <valarray>

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/timer.hpp>

#include <boost/numeric/bindings/atlas/cblas3.hpp>
#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
namespace atlas = boost::numeric::bindings::atlas;

# include <sys/time.h>

namespace ublas = boost::numeric::ublas;
using namespace std;

double gettime() {
    struct timeval t;
    gettimeofday(&t, NULL);
    return (double)t.tv_sec+t.tv_usec*0.000001;
}

void initmatrix(ublas::matrix<double>& A) {
 for(size_t i = 0; i<A.size1(); i++)
  for(size_t j = 0; j<A.size2(); j++)
   A(i,j) = 1.0*rand();
}
void initmatrix(ublas::matrix<double,ublas::column_major>& A) {
 for(size_t i = 0; i<A.size1(); i++)
  for(size_t j = 0; j<A.size2(); j++)
   A(i,j) = 1.0*rand();
} 

int main(int argc, char *argv[]) {

     long int D,r,n ;

    //std::clock_t start;
    //double diff;
    double t1,t2 ;

    //D = 1000 ;
    D = atoi(argv[1]) ;
    //r = 1000 ;
    r = atoi(argv[2]) ;
    //n = 1000 ;
    n = atoi(argv[3]) ;

    {
        //ublas::matrix<double, ublas::column_major> A(D,r),B(D,r),C(D,r);
    ublas::matrix<double, ublas::column_major>  A(D,r),B(D,r),C(D,r);
         initmatrix(A);
        initmatrix(B);

        cout << "division:" << endl ;
        t1 = gettime() ;
        C=element_div(A,B);
        t2 = gettime() ;
        std::cout << "Time : " <<  t2 - t1 << std::endl;

    cout << "product:" << endl ;
        t1 = gettime() ;
        C=element_prod(A,B);
        t2 = gettime() ;
        std::cout << "Time : " <<  t2 - t1 << std::endl;

    //std::cout << "A:" << A << std::endl ;
    //std::cout << "B:" << B << std::endl ;
    //std::cout << "C:" << C << std::endl ;

       
   }  

}

compilation:
>> g++ Hadamard_test.cpp -o Hadamard_test -Wall -O3 -lblas -latlas -DNDEBUG

>> ./Hadamard_test 1000000 100 1
division:
Time : 1.20392
product:
Time : 0.990796

Thanks,
Kaveh

________________________________
From: Kim Kuen Tang <kuentang_at_[hidden]>
To: ublas mailing list <ublas_at_[hidden]>
Sent: Thursday, May 28, 2009 4:19:08 PM
Subject: Re: [ublas] question regarding elementwise matrix product

Hi Kaveh,

Kaveh Kohan schrieb:
> Hi Kim,
> 
> Thank you for your reply.
> I am sorry if my questions sounds naive. I ran into a problem compiling this program. Do I need to have a specific version of phoenix or spirit to compile this program. Since in my system, phoenix.hpp is located in  "boost/spirit/phoenix.hpp", I changed your program to the following:
This is not a good idea, because including only phoenix.hpp is not enough in your boost version.
( I use 1.39.)

>   symmetric_matrix<double,lower> m(3, 3), n(3,3), u(3,3);
>   generate(arg1, lambda[val(1.0)] )(m.data());
the plain data of your matrix is stored in the unbounded_array. Using the member data() you have access to this vector. Since you have another version of boost, i suggest you to change the above code to     std::fill(m.data().begin(), m.data().end,1.0 ); This will have the same effect.
>   generate(arg1, lambda[val(2.0)] )(n.data());
>   transform(arg1,n.data().begin(),u.data().begin(),lambda[arg1/arg2])(m.data());
Using transform you can perform an element wise division. In blas there is already an element_wise division called element_div.
Here is the complete code:

# include <boost/numeric/ublas/symmetric.hpp>
# include <boost/numeric/ublas/matrix_expression.hpp>
# include <boost/numeric/ublas/io.hpp>

# include <vector>
# include <algorithm>

using namespace boost::numeric::ublas;

int main () {

   symmetric_matrix<double,lower> m(3, 3), n(3,3), u(3,3);
   std::fill(m.data().begin(), m.data().end(),1.0 );
   std::fill(n.data().begin(), n.data().end(),2.0 );
   u=element_div(m,n);
   std::cout << u << std::endl;
}
> 
> Plus can you explain a little what it does because notation-wise, it is very confusing for me. Definitely ignorance frommy side but I would appreciate if you describe it a little bit.
> 
> Thanks,
> Kaveh
> 
> ------------------------------------------------------------------------
> **
_______________________________________________
ublas mailing list
ublas_at_[hidden]
http://lists.boost.org/mailman/listinfo.cgi/ublas