Re: [eigen] Re: perf issue with vector of size 2 ?

["Christophe Prud'homme" <christophe.prudhomme@xxxxxxxxxxxxxxx>]

+the code

2010/7/21 Christophe Prud'homme <christophe.prudhomme@xxxxxxxxxxxxxxx>:
> attached the new code
>
> now the results (for C1=1,2 or 3 and C2=1) are much better even for
> the t1 vs t2 vs t3 . thanks a lot!
>
> C.
> --
> Christophe Prud'homme
> Université de Grenoble      christophe.prudhomme@xxxxxxxxxxxxxxx
> LJK - Room 55                  Tel: +33476635497
> 51, rue des Mathématiques      Fax: +33476631263
> BP53 38041 Grenoble Cedex 9
>       <http://ljk.imag.fr/membres/Christophe.Prudhomme/>
>



-- 
Christophe Prud'homme
Université de Grenoble      christophe.prudhomme@xxxxxxxxxxxxxxx
LJK - Room 55                  Tel: +33476635497
51, rue des Mathématiques      Fax: +33476631263
BP53 38041 Grenoble Cedex 9
      <http://ljk.imag.fr/membres/Christophe.Prudhomme/>

#include <boost/multi_array.hpp>
#include <Eigen/Eigen>
#include <boost/timer.hpp>

const int P = 10000;
static const int N=21;
static const int Q=25;
static const int C1=3;
static const int C2=3;

typedef Eigen::Matrix<double,C1,C2> vector_type;
typedef boost::multi_array<double,4> mad4;
typedef boost::multi_array<double,2> mad2;
typedef boost::multi_array<vector_type,2> made2;

EIGEN_DONT_INLINE
void
t1(  mad4& x1, mad4 const& x2 )
{
	for(int i = 0; i < N; ++i )
		for(int q = 0; q < Q; ++q )
		{
			for(int c1 = 0; c1 < C1; ++c1 )
				for(int c2 = 0; c2 < C2; ++c2)
				{
					x1[i][q][c1][c2] = cos(x2[i][q][c1][c2])*x2[i][q][c1][c2];
				}
		}

}
EIGEN_DONT_INLINE
void
t2(  mad4& w1, mad4 const& w2 )
{
	for(int i = 0; i < N; ++i )
	{
		for(int c1 = 0; c1 < C1; ++c1 )
			for(int c2 = 0; c2 < C2; ++c2)
				for(int q = 0; q < Q; ++q )
				{
					w1[i][c1][c2][q] = cos(w2[i][c1][c2][q])*w2[i][c1][c2][q];
				}
	}
}

EIGEN_DONT_INLINE
void
t3(  made2& y1, made2 const& y2 )
{
	for(int i = 0; i < N; ++i )
		for(int q = 0; q < Q; ++q )
		{
			y1[i][q] = y2[i][q].cwise()*y2[i][q].cwise().cos();
		}
}
EIGEN_DONT_INLINE
void
t4(  mad4& x, mad4 const& x1, mad4 const& x2 )
{
	for(int i = 0; i < N; ++i )
		for(int q = 0; q < Q; ++q )
		{
			x[i][q][0][0] = 0;
			for(int c1 = 0; c1 < C1; ++c1 )
				for(int c2 = 0; c2 < C2; ++c2)
				{
					x[i][q][0][0] += x2[i][q][c1][c2]*x1[i][q][c1][c2];
				}
		}
}
EIGEN_DONT_INLINE
void
t5(  mad4& w, mad4 const& w1, mad4 const& w2 )
{
	for(int i = 0; i < N; ++i )
	{

		for(int c1 = 0; c1 < C1; ++c1 )
			for(int c2 = 0; c2 < C2; ++c2)
			{
				for(int q = 0; q < Q; ++q )
				{

					w[i][0][0][q] += w2[i][c1][c2][q]*w1[i][c1][c2][q];
				}
			}

	}
}
EIGEN_DONT_INLINE
void
t6(  mad2& y, made2 const& y1, made2 const& y2 )
{
	for(int i = 0; i < N; ++i )
		for(int q = 0; q < Q; ++q )
		{
			y[i][q] = (y2[i][q].transpose()*y1[i][q]).trace();
		}
}
int main()
{
	boost::multi_array<double,4> x( boost::extents[N][Q][1][1] );
	boost::multi_array<double,4> x1( boost::extents[N][Q][C1][C2] );
	boost::multi_array<double,4> x2( boost::extents[N][Q][C1][C2] );
	boost::multi_array<double,4> w( boost::extents[N][1][1][Q] );
	boost::multi_array<double,4> w1( boost::extents[N][C1][C2][Q] );
	boost::multi_array<double,4> w2( boost::extents[N][C1][C2][Q] );
	boost::multi_array<double,2> y( boost::extents[N][Q] );
	boost::multi_array<vector_type,2> y1( boost::extents[N][Q] );
	boost::multi_array<vector_type,2> y2( boost::extents[N][Q] );

	std::fill( x.data(), x.data()+x.num_elements(), 1. );
	std::fill( x1.data(), x1.data()+x1.num_elements(), 2. );
	std::fill( x2.data(), x2.data()+x2.num_elements(), 3. );
	std::fill( w.data(), w.data()+w.num_elements(), -1. );
	std::fill( w1.data(), w1.data()+w1.num_elements(), 10. );
	std::fill( w2.data(), w2.data()+w2.num_elements(), -2. );
	std::fill( y.data(), y.data()+y.num_elements(), 0 );
	std::fill( y1.data(), y1.data()+y1.num_elements(), vector_type::Ones()*5 );
	std::fill( y2.data(), y2.data()+y2.num_elements(), -vector_type::Ones()*2 );

	boost::timer ti;
	for(int e = 0; e < P;++e )
		t1(x1,x2);
	std::cout << "time: " << ti.elapsed() << "\n";

	ti.restart();
	for(int e = 0; e < P;++e )
		t2( w1, w2 );
	std::cout << "time: " << ti.elapsed() << "\n";

	ti.restart();
	for(int e = 0; e < P;++e )
		t3( y1, y2 );
	std::cout << "time: " << ti.elapsed() << "\n";

	std::cout << "scalar product\n";

	ti.restart();
	for(int e = 0; e < P;++e )
		t4( x, x1, x2 );
	std::cout << "time: " << ti.elapsed() << "\n";

	ti.restart();
	for(int e = 0; e < P;++e )
		t5( w, w1, w2 );
	std::cout << "time: " << ti.elapsed() << "\n";

	ti.restart();
	for(int e = 0; e < P;++e )
		t6( y, y1, y2 );
	std::cout << "time: " << ti.elapsed() << "\n";

}