From 2482ad6ca470a905928f31de901d9accdd26acaa Mon Sep 17 00:00:00 2001
From: Michel Schanen <mschanen@anl.gov>
Date: Tue, 17 Apr 2018 13:36:19 -0500
Subject: [PATCH] analyze sparse pattern once

---
 include/ad.hpp | 14 +++++++++++++-
 1 file changed, 13 insertions(+), 1 deletion(-)

diff --git a/include/ad.hpp b/include/ad.hpp
index 7065581..b5c9516 100644
--- a/include/ad.hpp
+++ b/include/ad.hpp
@@ -448,6 +448,19 @@ template <class T> void integrate(pVector<T> &x) {
   std::vector<Eigen::Triplet<T> > tripletList;
   tripletList.reserve(y.dim()*y.dim());
   SparseLU<SparseMatrix<T> >  solver;
+#ifdef EIGEN_SPARSE
+  sys->jac_beuler<T>(x, xold, J);
+  for(size_t i = 0; i < J.nrows(); ++i) {
+    for (size_t j = 0; j < J.ncols(); ++j) {
+      if(J[i][j] != 0.0) {
+        tripletList.push_back(Eigen::Triplet<T>(i,j, J[i][j]));
+      }
+    }
+  }
+  sJ.setFromTriplets(tripletList.begin(), tripletList.end());
+  tripletList.clear();
+  solver.analyzePattern(sJ); 
+#endif
 #endif
   do {
     iteration = iteration + 1;
@@ -469,7 +482,6 @@ template <class T> void integrate(pVector<T> &x) {
     }
     sJ.setFromTriplets(tripletList.begin(), tripletList.end());
     tripletList.clear();
-    solver.analyzePattern(sJ); 
     solver.compute(sJ);
     eigymap = solver.solve(eigymap);
 #elif EIGEN