eigen/bench/benchCholesky.cpp
Gael Guennebaud 765219aa51 Big API change in Cholesky module:
* rename Cholesky to LLT
 * rename CholeskyWithoutSquareRoot to LDLT
 * rename MatrixBase::cholesky() to llt()
 * rename MatrixBase::choleskyNoSqrt() to ldlt()
 * make {LLT,LDLT}::solve() API consistent with other modules

Note that we are going to keep a source compatibility untill the next beta release.
E.g., the "old" Cholesky* classes, etc are still available for some time.
To be clear, Eigen beta2 should be (hopefully) source compatible with beta1,
and so beta2 will contain all the deprecated API of beta1. Those features marked
as deprecated will be removed in beta3 (or in the final 2.0 if there is no beta 3 !).

Also includes various updated in sparse Cholesky.
2008-10-13 15:53:27 +00:00

141 lines
3.4 KiB
C++

// g++ -DNDEBUG -O3 -I.. benchLLT.cpp -o benchLLT && ./benchLLT
// options:
// -DBENCH_GSL -lgsl /usr/lib/libcblas.so.3
// -DEIGEN_DONT_VECTORIZE
// -msse2
// -DREPEAT=100
// -DTRIES=10
// -DSCALAR=double
#include <Eigen/Array>
#include <Eigen/LLT>
#include <bench/BenchUtil.h>
using namespace Eigen;
#ifndef REPEAT
#define REPEAT 10000
#endif
#ifndef TRIES
#define TRIES 10
#endif
typedef float Scalar;
template <typename MatrixType>
__attribute__ ((noinline)) void benchLLT(const MatrixType& m)
{
int rows = m.rows();
int cols = m.cols();
int cost = 0;
for (int j=0; j<rows; ++j)
{
int r = std::max(rows - j -1,0);
cost += 2*(r*j+r+j);
}
int repeats = (REPEAT*1000)/(rows*rows);
typedef typename MatrixType::Scalar Scalar;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
MatrixType a = MatrixType::Random(rows,cols);
SquareMatrixType covMat = a * a.adjoint();
BenchTimer timerNoSqrt, timerSqrt;
Scalar acc = 0;
int r = ei_random<int>(0,covMat.rows()-1);
int c = ei_random<int>(0,covMat.cols()-1);
for (int t=0; t<TRIES; ++t)
{
timerNoSqrt.start();
for (int k=0; k<repeats; ++k)
{
LDLT<SquareMatrixType> cholnosqrt(covMat);
acc += cholnosqrt.matrixL().coeff(r,c);
}
timerNoSqrt.stop();
}
for (int t=0; t<TRIES; ++t)
{
timerSqrt.start();
for (int k=0; k<repeats; ++k)
{
LLT<SquareMatrixType> chol(covMat);
acc += chol.matrixL().coeff(r,c);
}
timerSqrt.stop();
}
if (MatrixType::RowsAtCompileTime==Dynamic)
std::cout << "dyn ";
else
std::cout << "fixed ";
std::cout << covMat.rows() << " \t"
<< (timerNoSqrt.value() * REPEAT) / repeats << "s \t"
<< (timerSqrt.value() * REPEAT) / repeats << "s "
<< "(" << 1e-6 * cost*repeats/timerSqrt.value() << " MFLOPS)\n";
#ifdef BENCH_GSL
if (MatrixType::RowsAtCompileTime==Dynamic)
{
timerSqrt.reset();
gsl_matrix* gslCovMat = gsl_matrix_alloc(covMat.rows(),covMat.cols());
gsl_matrix* gslCopy = gsl_matrix_alloc(covMat.rows(),covMat.cols());
eiToGsl(covMat, &gslCovMat);
for (int t=0; t<TRIES; ++t)
{
timerSqrt.start();
for (int k=0; k<repeats; ++k)
{
gsl_matrix_memcpy(gslCopy,gslCovMat);
gsl_linalg_cholesky_decomp(gslCopy);
acc += gsl_matrix_get(gslCopy,r,c);
}
timerSqrt.stop();
}
std::cout << " | \t"
<< timerSqrt.value() * REPEAT / repeats << "s";
gsl_matrix_free(gslCovMat);
}
#endif
std::cout << "\n";
// make sure the compiler does not optimize too much
if (acc==123)
std::cout << acc;
}
int main(int argc, char* argv[])
{
const int dynsizes[] = {/*4,6,8,12,16,24,32,49,64,67,128,129,130,131,132,*/256,257,258,259,260,512,900,0};
std::cout << "size no sqrt standard";
#ifdef BENCH_GSL
std::cout << " GSL (standard + double + ATLAS) ";
#endif
std::cout << "\n";
for (uint i=0; dynsizes[i]>0; ++i)
benchLLT(Matrix<Scalar,Dynamic,Dynamic>(dynsizes[i],dynsizes[i]));
// benchLLT(Matrix<Scalar,2,2>());
// benchLLT(Matrix<Scalar,3,3>());
// benchLLT(Matrix<Scalar,4,4>());
// benchLLT(Matrix<Scalar,5,5>());
// benchLLT(Matrix<Scalar,6,6>());
// benchLLT(Matrix<Scalar,7,7>());
// benchLLT(Matrix<Scalar,8,8>());
// benchLLT(Matrix<Scalar,12,12>());
// benchLLT(Matrix<Scalar,16,16>());
return 0;
}