Squashed 'third_party/eigen/' content from commit 61d72f6
Change-Id: Iccc90fa0b55ab44037f018046d2fcffd90d9d025
git-subtree-dir: third_party/eigen
git-subtree-split: 61d72f6383cfa842868c53e30e087b0258177257
diff --git a/bench/sparse_lu.cpp b/bench/sparse_lu.cpp
new file mode 100644
index 0000000..5c75001
--- /dev/null
+++ b/bench/sparse_lu.cpp
@@ -0,0 +1,132 @@
+
+// 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;
+}
+