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;
+}
+