Eigen works with AD

include/ad.hpp

@@ -34,8 +34,10 @@ original variable.
 #include <mpi.h>
 #include <stdlib.h>
 #include "parallel.hpp"
+#ifdef EIGEN
 #include "eigen_codi.hpp"
 #include <Eigen/Dense>
+#endif
 //#include "linsolve.hpp"
 //
 #ifdef __INTEL_COMPILER
@@ -47,7 +49,9 @@ original variable.
 
 using namespace std;
 using namespace alg;
+#ifdef EIGEN
 using namespace Eigen;
+#endif
 
 typedef codi::RealForwardGen<double> t1s;
 typedef codi::RealForwardGen<t1s> t2s;
@@ -426,7 +430,11 @@ template <class T> void integrate(pVector<T> &x) {
   xold = x;
   sys->residual_beuler<T>(x, xold, y);
   J.zeros();
-
+#ifdef EIGEN
+  Map<Matrix<T, Dynamic, 1> > eigymap(y.get_datap(), y.dim());
+  Map<Matrix<T, Dynamic, Dynamic> > eigJmap(J.get_datap(), J.nrows(), J.ncols());
+#endif
+
   do {
     iteration = iteration + 1;
     sys->jac_beuler<T>(x, xold, J);
@@ -437,14 +445,13 @@ template <class T> void integrate(pVector<T> &x) {
 #endif
     yold = y;
     Jold = J;
-    Matrix<T, 14, 14> eigJ;
-    eigJ= Map<Matrix<T,14,14> >(J.get_datap());
-    Matrix<T, 14, 1> eigy;
-    eigy= Map<Matrix<T,14,1> >(y.get_datap());
-    eigy=eigJ.fullPivLu().solve(eigy);
+#ifdef EIGEN
+    eigymap=eigJmap.fullPivLu().solve(eigymap);
+#else
     adlinsolve<T>(J, y);
+#endif
     pVector<T> res(dim);
-    res = Jold * y - yold;
+    res = Jold * y + yold;
 #ifdef DBUG
     cout << "Norm(y): " << y.norm() << endl << " y: " << y << endl;
     cout << "|Ax - b| " << res.norm() << endl;