add support for ublas

This commit is contained in:
Gael Guennebaud 2011-03-23 11:39:35 +01:00
parent ec32d2c807
commit 611fc17894
4 changed files with 70 additions and 65 deletions

View File

@ -26,7 +26,7 @@ typedef SparseMatrix<Scalar> EigenSparseMatrix;
void fillMatrix(float density, int rows, int cols, EigenSparseMatrix& dst)
{
dst.reserve(rows*cols*density);
dst.reserve(double(rows)*cols*density);
for(int j = 0; j < cols; j++)
{
for(int i = 0; i < rows; i++)
@ -122,22 +122,24 @@ void eiToCSparse(const EigenSparseMatrix& src, cs* &dst)
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/vector_of_vector.hpp>
#include <boost/numeric/ublas/operation.hpp>
// #include <boost/numeric/ublas/sparse_prod.hpp>
// using namespace boost;
// using namespace boost::numeric;
// using namespace boost::numeric::ublas;
typedef boost::numeric::ublas::compressed_matrix<Scalar,boost::numeric::ublas::column_major> UBlasSparse;
typedef boost::numeric::ublas::compressed_matrix<Scalar,boost::numeric::ublas::column_major> UblasMatrix;
void eiToUblas(const EigenSparseMatrix& src, UblasMatrix& dst)
void eiToUblas(const EigenSparseMatrix& src, UBlasSparse& dst)
{
dst.resize(src.rows(), src.cols(), false);
for (int j=0; j<src.cols(); ++j)
for (EigenSparseMatrix::InnerIterator it(src.derived(), j); it; ++it)
dst(it.index(),j) = it.value();
}
template <typename EigenType, typename UblasType>
void eiToUblasVec(const EigenType& src, UblasType& dst)
{
dst.resize(src.size());
for (int j=0; j<src.size(); ++j)
dst[j] = src.coeff(j);
}
#endif
#ifdef OSKI

View File

@ -69,4 +69,24 @@ void eiToGsl(const EigenMatrixType& src, gsl_matrix** dst)
}
#endif
#ifdef BENCH_UBLAS
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
template <typename EigenMatrixType, typename UblasMatrixType>
void eiToUblas(const EigenMatrixType& src, UblasMatrixType& dst)
{
dst.resize(src.rows(),src.cols());
for (int j=0; j<src.cols(); ++j)
for (int i=0; i<src.rows(); ++i)
dst(i,j) = src.coeff(i,j);
}
template <typename EigenType, typename UblasType>
void eiToUblasVec(const EigenType& src, UblasType& dst)
{
dst.resize(src.size());
for (int j=0; j<src.size(); ++j)
dst[j] = src.coeff(j);
}
#endif
#endif // EIGEN_BENCH_UTIL_H

View File

@ -4,7 +4,7 @@
// -DNOGMM -DNOMTL -DCSPARSE
// -I /home/gael/Coding/LinearAlgebra/CSparse/Include/ /home/gael/Coding/LinearAlgebra/CSparse/Lib/libcsparse.a
#ifndef SIZE
#define SIZE 10000
#define SIZE 650000
#endif
#ifndef DENSITY
@ -62,7 +62,8 @@ int main(int argc, char *argv[])
BenchTimer timer;
for (float density = DENSITY; density>=MINDENSITY; density*=0.5)
{
fillMatrix(density, rows, cols, sm1);
//fillMatrix(density, rows, cols, sm1);
fillMatrix2(7, rows, cols, sm1);
// dense matrices
#ifdef DENSEMATRIX
@ -76,14 +77,14 @@ int main(int argc, char *argv[])
for (int k=0; k<REPEAT; ++k)
v2 = m1 * v1;
timer.stop();
std::cout << " a * v:\t" << timer.value() << endl;
std::cout << " a * v:\t" << timer.best() << " " << double(REPEAT)/timer.best() << " * / sec " << endl;
timer.reset();
timer.start();
for (int k=0; k<REPEAT; ++k)
v2 = m1.transpose() * v1;
timer.stop();
std::cout << " a' * v:\t" << timer.value() << endl;
std::cout << " a' * v:\t" << timer.best() << endl;
}
#endif
@ -92,12 +93,12 @@ int main(int argc, char *argv[])
std::cout << "Eigen sparse\t" << sm1.nonZeros()/float(sm1.rows()*sm1.cols())*100 << "%\n";
BENCH(asm("#myc"); v2 = sm1 * v1; asm("#myd");)
std::cout << " a * v:\t" << timer.value() << endl;
std::cout << " a * v:\t" << timer.best()/REPEAT << " " << double(REPEAT)/timer.best(REAL_TIMER) << " * / sec " << endl;
BENCH( { asm("#mya"); v2 = sm1.transpose() * v1; asm("#myb"); })
std::cout << " a' * v:\t" << timer.value() << endl;
std::cout << " a' * v:\t" << timer.best()/REPEAT << endl;
}
// {
@ -130,6 +131,28 @@ int main(int argc, char *argv[])
std::cout << " a' * v:\t" << timer.value() << endl;
}
#endif
#ifndef NOUBLAS
{
std::cout << "ublas sparse\t" << density*100 << "%\n";
UBlasSparse m1(rows,cols);
eiToUblas(sm1, m1);
boost::numeric::ublas::vector<Scalar> uv1, uv2;
eiToUblasVec(v1,uv1);
eiToUblasVec(v2,uv2);
// std::vector<Scalar> gmmV1(cols), gmmV2(cols);
// Map<Matrix<Scalar,Dynamic,1> >(&gmmV1[0], cols) = v1;
// Map<Matrix<Scalar,Dynamic,1> >(&gmmV2[0], cols) = v2;
BENCH( uv2 = boost::numeric::ublas::prod(m1, uv1); )
std::cout << " a * v:\t" << timer.value() << endl;
// BENCH( boost::ublas::prod(gmm::transposed(m1), gmmV1, gmmV2); )
// std::cout << " a' * v:\t" << timer.value() << endl;
}
#endif
// MTL4
#ifndef NOMTL

View File

@ -20,6 +20,7 @@
#include <algorithm>
#include "BenchTimer.h"
#include "BenchUtil.h"
#include "BenchSparseUtil.h"
#ifndef NBTRIES
@ -228,16 +229,12 @@ int main(int argc, char *argv[])
eiToCSparse(sm1, m1);
eiToCSparse(sm2, m2);
// timer.reset();
// timer.start();
// for (int k=0; k<REPEAT; ++k)
BENCH(
{
m3 = cs_sorted_multiply(m1, m2);
if (!m3)
{
std::cerr << "cs_multiply failed\n";
// break;
}
// cs_print(m3, 0);
cs_spfree(m3);
@ -254,16 +251,11 @@ int main(int argc, char *argv[])
#ifndef NOUBLAS
{
std::cout << "ublas\t" << nnzPerCol << "%\n";
UblasMatrix m1(rows,cols), m2(rows,cols), m3(rows,cols);
UBlasSparse m1(rows,cols), m2(rows,cols), m3(rows,cols);
eiToUblas(sm1, m1);
eiToUblas(sm2, m2);
BENCH(boost::numeric::ublas::prod(m1, m2, m3););
// timer.reset();
// timer.start();
// for (int k=0; k<REPEAT; ++k)
// gmm::mult(m1, m2, gmmT3);
// timer.stop();
std::cout << " a * b:\t" << timer.value() << endl;
}
#endif
@ -277,34 +269,18 @@ int main(int argc, char *argv[])
eiToGmm(sm1, m1);
eiToGmm(sm2, m2);
timer.reset();
timer.start();
for (int k=0; k<REPEAT; ++k)
gmm::mult(m1, m2, gmmT3);
timer.stop();
BENCH(gmm::mult(m1, m2, gmmT3););
std::cout << " a * b:\t" << timer.value() << endl;
// timer.reset();
// timer.start();
// for (int k=0; k<REPEAT; ++k)
// gmm::mult(gmm::transposed(m1), m2, gmmT3);
// timer.stop();
// BENCH(gmm::mult(gmm::transposed(m1), m2, gmmT3););
// std::cout << " a' * b:\t" << timer.value() << endl;
//
// if (rows<500)
// {
// timer.reset();
// timer.start();
// for (int k=0; k<REPEAT; ++k)
// gmm::mult(gmm::transposed(m1), gmm::transposed(m2), gmmT3);
// timer.stop();
// BENCH(gmm::mult(gmm::transposed(m1), gmm::transposed(m2), gmmT3););
// std::cout << " a' * b':\t" << timer.value() << endl;
//
// timer.reset();
// timer.start();
// for (int k=0; k<REPEAT; ++k)
// gmm::mult(m1, gmm::transposed(m2), gmmT3);
// timer.stop();
// BENCH(gmm::mult(m1, gmm::transposed(m2), gmmT3););
// std::cout << " a * b':\t" << timer.value() << endl;
// }
// else
@ -323,32 +299,16 @@ int main(int argc, char *argv[])
eiToMtl(sm1, m1);
eiToMtl(sm2, m2);
timer.reset();
timer.start();
for (int k=0; k<REPEAT; ++k)
m3 = m1 * m2;
timer.stop();
BENCH(m3 = m1 * m2;);
std::cout << " a * b:\t" << timer.value() << endl;
// timer.reset();
// timer.start();
// for (int k=0; k<REPEAT; ++k)
// m3 = trans(m1) * m2;
// timer.stop();
// BENCH(m3 = trans(m1) * m2;);
// std::cout << " a' * b:\t" << timer.value() << endl;
//
// timer.reset();
// timer.start();
// for (int k=0; k<REPEAT; ++k)
// m3 = trans(m1) * trans(m2);
// timer.stop();
// BENCH(m3 = trans(m1) * trans(m2););
// std::cout << " a' * b':\t" << timer.value() << endl;
//
// timer.reset();
// timer.start();
// for (int k=0; k<REPEAT; ++k)
// m3 = m1 * trans(m2);
// timer.stop();
// BENCH(m3 = m1 * trans(m2););
// std::cout << " a * b' :\t" << timer.value() << endl;
}
#endif