blob: f074e55de8920d6eb106a926109b7462487b6f79 [file] [log] [blame]
// g++ -I.. sparse_lu.cpp -O3 -g0 -I /usr/include/superlu/ -lsuperlu -lgfortran -DSIZE=1000 -DDENSITY=.05 && ./a.out
#define EIGEN_SUPERLU_SUPPORT
#define EIGEN_UMFPACK_SUPPORT
#include <Eigen/Sparse>
#define NOGMM
#define NOMTL
#ifndef SIZE
#define SIZE 10
#endif
#ifndef DENSITY
#define DENSITY 0.01
#endif
#ifndef REPEAT
#define REPEAT 1
#endif
#include "BenchSparseUtil.h"
#ifndef MINDENSITY
#define MINDENSITY 0.0004
#endif
#ifndef NBTRIES
#define NBTRIES 10
#endif
#define BENCH(X) \
timer.reset(); \
for (int _j = 0; _j < NBTRIES; ++_j) { \
timer.start(); \
for (int _k = 0; _k < REPEAT; ++_k) { \
X \
} \
timer.stop(); \
}
typedef Matrix<Scalar, Dynamic, 1> VectorX;
#include <Eigen/LU>
template <int Backend>
void doEigen(const char* name, const EigenSparseMatrix& sm1, const VectorX& b, VectorX& x, int flags = 0) {
std::cout << name << "..." << std::flush;
BenchTimer timer;
timer.start();
SparseLU<EigenSparseMatrix, Backend> lu(sm1, flags);
timer.stop();
if (lu.succeeded())
std::cout << ":\t" << timer.value() << endl;
else {
std::cout << ":\t FAILED" << endl;
return;
}
bool ok;
timer.reset();
timer.start();
ok = lu.solve(b, &x);
timer.stop();
if (ok)
std::cout << " solve:\t" << timer.value() << endl;
else
std::cout << " solve:\t"
<< " FAILED" << endl;
// std::cout << x.transpose() << "\n";
}
int main(int argc, char* argv[]) {
int rows = SIZE;
int cols = SIZE;
float density = DENSITY;
BenchTimer timer;
VectorX b = VectorX::Random(cols);
VectorX x = VectorX::Random(cols);
bool densedone = false;
// for (float density = DENSITY; density>=MINDENSITY; density*=0.5)
// float density = 0.5;
{
EigenSparseMatrix sm1(rows, cols);
fillMatrix(density, rows, cols, sm1);
// dense matrices
#ifdef DENSEMATRIX
if (!densedone) {
densedone = true;
std::cout << "Eigen Dense\t" << density * 100 << "%\n";
DenseMatrix m1(rows, cols);
eiToDense(sm1, m1);
BenchTimer timer;
timer.start();
FullPivLU<DenseMatrix> lu(m1);
timer.stop();
std::cout << "Eigen/dense:\t" << timer.value() << endl;
timer.reset();
timer.start();
lu.solve(b, &x);
timer.stop();
std::cout << " solve:\t" << timer.value() << endl;
// std::cout << b.transpose() << "\n";
// std::cout << x.transpose() << "\n";
}
#endif
#ifdef EIGEN_UMFPACK_SUPPORT
x.setZero();
doEigen<Eigen::UmfPack>("Eigen/UmfPack (auto)", sm1, b, x, 0);
#endif
#ifdef EIGEN_SUPERLU_SUPPORT
x.setZero();
doEigen<Eigen::SuperLU>("Eigen/SuperLU (nat)", sm1, b, x, Eigen::NaturalOrdering);
// doEigen<Eigen::SuperLU>("Eigen/SuperLU (MD AT+A)", sm1, b, x, Eigen::MinimumDegree_AT_PLUS_A);
// doEigen<Eigen::SuperLU>("Eigen/SuperLU (MD ATA)", sm1, b, x, Eigen::MinimumDegree_ATA);
doEigen<Eigen::SuperLU>("Eigen/SuperLU (COLAMD)", sm1, b, x, Eigen::ColApproxMinimumDegree);
#endif
}
return 0;
}