123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132 |
- // 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;
- }
|