mex/sources/estimation/DecisionRules.cc

 /*
- * Copyright (C) 2010-2011 Dynare Team
+ * Copyright (C) 2010-2012 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -23,11 +23,11 @@
 
 #include "DecisionRules.hh"
 
-DecisionRules::DecisionRules(size_t n_arg, size_t p_arg,
-                             const std::vector<size_t> &zeta_fwrd_arg,
-                             const std::vector<size_t> &zeta_back_arg,
-                             const std::vector<size_t> &zeta_mixed_arg,
-                             const std::vector<size_t> &zeta_static_arg,
+DecisionRules::DecisionRules(ptrdiff_t n_arg, ptrdiff_t p_arg,
+                             const std::vector<ptrdiff_t> &zeta_fwrd_arg,
+                             const std::vector<ptrdiff_t> &zeta_back_arg,
+                             const std::vector<ptrdiff_t> &zeta_mixed_arg,
+                             const std::vector<ptrdiff_t> &zeta_static_arg,
                              double qz_criterium) :
   n(n_arg), p(p_arg), zeta_fwrd(zeta_fwrd_arg), zeta_back(zeta_back_arg),
   zeta_mixed(zeta_mixed_arg), zeta_static(zeta_static_arg),
@@ -37,25 +37,24 @@ DecisionRules::DecisionRules(size_t n_arg, size_t p_arg,
   n_dynamic(n-n_static),
   S(n, n_static),
   A(n, n_back_mixed + n + n_fwrd_mixed),
-  D(n_fwrd + n_back + 2*n_mixed),
-  E(n_fwrd + n_back + 2*n_mixed),
-  Z_prime(n_fwrd + n_back + 2*n_mixed),
+  D(n_fwrd + n_back + 2*n_mixed, n_fwrd + n_back + 2*n_mixed),
+  E(n_fwrd + n_back + 2*n_mixed, n_fwrd + n_back + 2*n_mixed),
+  Z_prime(n_fwrd + n_back + 2*n_mixed, n_fwrd + n_back + 2*n_mixed),
   QR(n, n_static, n_back_mixed + n + n_fwrd_mixed),
   GSD(n_fwrd + n_back + 2*n_mixed, qz_criterium),
-  LU1(n_fwrd_mixed),
-  LU2(n_back_mixed),
-  LU3(n_static),
+  Solver1(n_fwrd_mixed, n_fwrd_mixed),
+  Solver2(n_back_mixed, n_back_mixed),
+  Solver3(n_static, n_static),
+  Solver4(n, n),
   Z21(n_fwrd_mixed, n_back_mixed),
-  g_y_back(n_back_mixed),
-  g_y_back_tmp(n_back_mixed),
+  g_y_back(n_back_mixed, n_back_mixed),
+  g_y_fwrd(n_fwrd_mixed, n_fwrd_mixed),
   g_y_static(n_static, n_back_mixed),
-  A0s(n_static),
+  A0s(n_static, n_static),
   A0d(n_static, n_dynamic),
   g_y_dynamic(n_dynamic, n_back_mixed),
-  g_y_static_tmp(n_fwrd_mixed, n_back_mixed),
   g_u_tmp1(n, n_back_mixed),
-  g_u_tmp2(n),
-  LU4(n)
+  g_u_tmp2(n, n)
 {
   assert(n == n_back + n_fwrd + n_mixed + n_static);
 
@@ -70,7 +69,7 @@ DecisionRules::DecisionRules(size_t n_arg, size_t p_arg,
             back_inserter(zeta_dynamic));
 
   // Compute beta_back and pi_back
-  for (size_t i = 0; i < n_back_mixed; i++)
+  for (ptrdiff_t i = 0; i < n_back_mixed; i++)
     if (find(zeta_mixed.begin(), zeta_mixed.end(), zeta_back_mixed[i])
         == zeta_mixed.end())
       pi_back.push_back(i);
@@ -78,7 +77,7 @@ DecisionRules::DecisionRules(size_t n_arg, size_t p_arg,
       beta_back.push_back(i);
 
   // Compute beta_fwrd and pi_fwrd
-  for (size_t i = 0; i < n_fwrd_mixed; i++)
+  for (ptrdiff_t i = 0; i < n_fwrd_mixed; i++)
     if (find(zeta_mixed.begin(), zeta_mixed.end(), zeta_fwrd_mixed[i])
         == zeta_mixed.end())
       pi_fwrd.push_back(i);
@@ -87,115 +86,92 @@ DecisionRules::DecisionRules(size_t n_arg, size_t p_arg,
 }
 
 void
-DecisionRules::compute(const Matrix &jacobian, Matrix &g_y, Matrix &g_u) throw (BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
+DecisionRules::compute(const MatrixXd &jacobian, MatrixXd &g_y, MatrixXd &g_u) throw (BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
 {
-  assert(jacobian.getRows() == n
-         && jacobian.getCols() == (n_back_mixed + n + n_fwrd_mixed + p));
-  assert(g_y.getRows() == n && g_y.getCols() == n_back_mixed);
-  assert(g_u.getRows() == n && g_u.getCols() == p);
+  assert(jacobian.rows() == n
+         && jacobian.cols() == (n_back_mixed + n + n_fwrd_mixed + p));
+  assert(g_y.rows() == n && g_y.cols() == n_back_mixed);
+  assert(g_u.rows() == n && g_u.cols() == p);
 
   // Construct S, perform QR decomposition and get A = Q*jacobian
-  for (size_t i = 0; i < n_static; i++)
-    mat::col_copy(jacobian, n_back_mixed+zeta_static[i], S, i);
+  for (ptrdiff_t i = 0; i < n_static; i++)
+    S.col(i) = jacobian.col(n_back_mixed+zeta_static[i]);
 
-  A = MatrixConstView(jacobian, 0, 0, n, n_back_mixed + n + n_fwrd_mixed);
+  A = jacobian.block(0, 0, n, n_back_mixed + n + n_fwrd_mixed);
   QR.computeAndLeftMultByQ(S, "T", A);
 
   // Construct matrix D
-  D.setAll(0.0);
-  for (size_t i = 0; i < n_mixed; i++)
+  D.setZero();
+  for (ptrdiff_t i = 0; i < n_mixed; i++)
     D(n - n_static + i, beta_back[i]) = 1.0;
-  for (size_t j = 0; j < n_back_mixed; j++)
-    mat::col_copy(A, n_back_mixed + zeta_back_mixed[j], n_static, n - n_static,
-                  D, j, 0);
-  MatrixView(D, 0, n_back_mixed, n - n_static, n_fwrd_mixed) = MatrixView(A, n_static, n_back_mixed + n, n - n_static, n_fwrd_mixed);
+  for (ptrdiff_t j = 0; j < n_back_mixed; j++)
+    D.col(j).head(n-n_static) = A.col(n_back_mixed + zeta_back_mixed[j]).tail(n-n_static);
+  
+  D.block(0, n_back_mixed, n - n_static, n_fwrd_mixed) = A.block(n_static, n_back_mixed + n, n - n_static, n_fwrd_mixed);
 
   // Construct matrix E
-  E.setAll(0.0);
-  for (size_t i = 0; i < n_mixed; i++)
+  E.setZero();
+  for (ptrdiff_t i = 0; i < n_mixed; i++)
     E(n - n_static + i, n_back_mixed + beta_fwrd[i]) = 1.0;
-  MatrixView(E, 0, 0, n - n_static, n_back_mixed) = MatrixView(A, n_static, 0, n - n_static, n_back_mixed);
-  for (size_t j = 0; j < n_fwrd; j++)
-    mat::col_copy(A, n_back_mixed + zeta_fwrd_mixed[pi_fwrd[j]], n_static, n - n_static,
-                  E, n_back_mixed + pi_fwrd[j], 0);
-  MatrixView E_tmp(E, 0, 0, n - n_static, n_fwrd + n_back + 2*n_mixed);
-  mat::negate(E_tmp); // Here we take the opposite of some of the zeros initialized in the constructor, but it is not a problem
+  E.block(0, 0, n - n_static, n_back_mixed) = -A.block(n_static, 0, n - n_static, n_back_mixed);
+  for (ptrdiff_t j = 0; j < n_fwrd; j++)
+    E.col(n_back_mixed + pi_fwrd[j]).head(n-n_static) = -A.col(n_back_mixed + zeta_fwrd_mixed[pi_fwrd[j]]).tail(n - n_static);
 
   // Perform the generalized Schur
-  size_t sdim;
+  ptrdiff_t sdim;
   GSD.compute(E, D, Z_prime, sdim);
 
   if (n_back_mixed != sdim)
     throw BlanchardKahnException(true, n_fwrd_mixed, n_fwrd + n_back + 2*n_mixed - sdim);
 
   // Compute DR for forward variables w.r. to endogenous
-  MatrixView Z21_prime(Z_prime, 0, n_back_mixed, n_back_mixed, n_fwrd_mixed),
-    Z22_prime(Z_prime, n_back_mixed, n_back_mixed, n_fwrd_mixed, n_fwrd_mixed);
+  Solver1.compute(Z_prime.block(n_back_mixed, n_back_mixed, n_fwrd_mixed, n_fwrd_mixed).transpose());
+  
+  if (!Solver1.isInvertible())
+    throw BlanchardKahnException(false, n_fwrd_mixed, n_fwrd + n_back + 2*n_mixed - sdim);
 
-  mat::transpose(Z21, Z21_prime);
+  g_y_fwrd = -Solver1.solve(Z_prime.block(0, n_back_mixed, n_back_mixed, n_fwrd_mixed).transpose());
 
-  try
-    {
-      LU1.invMult("T", Z22_prime, Z21);
-    }
-  catch (LUSolver::LUException &e)
-    {
-      throw BlanchardKahnException(false, n_fwrd_mixed, n_fwrd + n_back + 2*n_mixed - sdim);
-    }
-  mat::negate(Z21);
-  const Matrix &g_y_fwrd = Z21;
-
-  for (size_t i = 0; i < n_fwrd_mixed; i++)
-    mat::row_copy(g_y_fwrd, i, g_y, zeta_fwrd_mixed[i]);
+  for (ptrdiff_t i = 0; i < n_fwrd_mixed; i++)
+    g_y.row(zeta_fwrd_mixed[i]) = g_y_fwrd.row(i);
 
   // Compute DR for backward variables w.r. to endogenous
-  MatrixView Z11_prime(Z_prime, 0, 0, n_back_mixed, n_back_mixed),
-    T11(D, 0, 0, n_back_mixed, n_back_mixed),
-    S11(E, 0, 0, n_back_mixed, n_back_mixed);
-  mat::set_identity(g_y_back);
-  g_y_back_tmp = Z11_prime;
-  LU2.invMult("N", g_y_back_tmp, g_y_back);
-  g_y_back_tmp = g_y_back;
-  blas::gemm("N", "N", 1.0, S11, g_y_back_tmp, 0.0, g_y_back);
-  LU2.invMult("N", T11, g_y_back);
-  g_y_back_tmp = g_y_back;
-  blas::gemm("N", "N", 1.0, Z11_prime, g_y_back_tmp, 0.0, g_y_back);
-
+  Solver2.compute(Z_prime.block(0, 0, n_back_mixed, n_back_mixed));
+  g_y_back.noalias() = E.block(0, 0, n_back_mixed, n_back_mixed) * Solver2.solve(MatrixXd::Identity(n_back_mixed, n_back_mixed)); // S_{11} * Z'_{11}^{-1}
+  
+  Solver2.compute(D.block(0, 0, n_back_mixed, n_back_mixed));
+  g_y_back = Z_prime.block(0, 0, n_back_mixed, n_back_mixed) * Solver2.solve(g_y_back);
+  
   // TODO: avoid to copy mixed variables again, rather test it...
-  for (size_t i = 0; i < n_back_mixed; i++)
-    mat::row_copy(g_y_back, i, g_y, zeta_back_mixed[i]);
+  for (ptrdiff_t i = 0; i < n_back_mixed; i++)
+    g_y.row(zeta_back_mixed[i]) = g_y_back.row(i);
 
   // Compute DR for static variables w.r. to endogenous
-  g_y_static = MatrixView(A, 0, 0, n_static, n_back_mixed);
-  for (size_t i = 0; i < n_dynamic; i++)
+  for (ptrdiff_t i = 0; i < n_dynamic; i++)
     {
-      mat::row_copy(g_y, zeta_dynamic[i], g_y_dynamic, i);
-      mat::col_copy(A, n_back_mixed + zeta_dynamic[i], 0, n_static, A0d, i, 0);
+      g_y_dynamic.row(i) = g_y.row(zeta_dynamic[i]);
+      A0d.col(i) = A.col(n_back_mixed + zeta_dynamic[i]).head(n_static);
     }
-  blas::gemm("N", "N", 1.0, A0d, g_y_dynamic, 1.0, g_y_static);
-  blas::gemm("N", "N", 1.0, g_y_fwrd, g_y_back, 0.0, g_y_static_tmp);
-  blas::gemm("N", "N", 1.0, MatrixView(A, 0, n_back_mixed + n, n_static, n_fwrd_mixed),
-             g_y_static_tmp, 1.0, g_y_static);
-  for (size_t i = 0; i < n_static; i++)
-    mat::col_copy(A, n_back_mixed + zeta_static[i], 0, n_static, A0s, i, 0);
-  LU3.invMult("N", A0s, g_y_static);
-  mat::negate(g_y_static);
-
-  for (size_t i = 0; i < n_static; i++)
-    mat::row_copy(g_y_static, i, g_y, zeta_static[i]);
 
+  for (ptrdiff_t i = 0; i < n_static; i++)
+    A0s.col(i) = A.col(n_back_mixed + zeta_static[i]).head(n_static);
+  
+  Solver3.compute(A0s);
+  
+  g_y_static = -Solver3.solve(A.block(0, n_back_mixed + n, n_static, n_fwrd_mixed)*g_y_fwrd*g_y_back
+                              + A0d*g_y_dynamic + A.block(0, 0, n_static, n_back_mixed));
+
+  for (ptrdiff_t i = 0; i < n_static; i++)
+    g_y.row(zeta_static[i]) = g_y_static.row(i);
+  
   // Compute DR for all endogenous w.r. to shocks
-  blas::gemm("N", "N", 1.0, MatrixConstView(jacobian, 0, n_back_mixed + n, n, n_fwrd_mixed), g_y_fwrd, 0.0, g_u_tmp1);
-  g_u_tmp2 = MatrixConstView(jacobian, 0, n_back_mixed, n, n);
-  for (size_t i = 0; i < n_back_mixed; i++)
-    {
-      VectorView c1 = mat::get_col(g_u_tmp2, zeta_back_mixed[i]),
-        c2 = mat::get_col(g_u_tmp1, i);
-      vec::add(c1, c2);
-    }
-  g_u = MatrixConstView(jacobian, 0, n_back_mixed + n + n_fwrd_mixed, n, p);
-  LU4.invMult("N", g_u_tmp2, g_u);
-  mat::negate(g_u);
+  g_u_tmp1 = jacobian.block(0, n_back_mixed + n, n, n_fwrd_mixed) * g_y_fwrd;
+  
+  g_u_tmp2 = jacobian.block(0, n_back_mixed, n, n);
+  for (ptrdiff_t i = 0; i < n_back_mixed; i++)
+    g_u_tmp2.col(zeta_back_mixed[i]) += g_u_tmp1.col(i);
+  Solver4.compute(g_u_tmp2);
+  g_u = -Solver4.solve(jacobian.block(0, n_back_mixed + n + n_fwrd_mixed, n, p));
 }
 
 std::ostream &

mex/sources/estimation/DecisionRules.hh

 /*
- * Copyright (C) 2010 Dynare Team
+ * Copyright (C) 2010-2012 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -20,28 +20,29 @@
 #include <cstdlib>
 #include <vector>
 
-#include "Vector.hh"
-#include "Matrix.hh"
+#include <Eigen/Core>
+#include <Eigen/QR>
+
 #include "QRDecomposition.hh"
 #include "GeneralizedSchurDecomposition.hh"
-#include "LUSolver.hh"
+
+using namespace Eigen;
 
 class DecisionRules
 {
 private:
-  const size_t n, p;
-  const std::vector<size_t> zeta_fwrd, zeta_back, zeta_mixed, zeta_static;
-  const size_t n_fwrd, n_back, n_mixed, n_static, n_back_mixed, n_fwrd_mixed, n_dynamic;
-  std::vector<size_t> zeta_fwrd_mixed, zeta_back_mixed, zeta_dynamic,
+  const ptrdiff_t n, p;
+  const std::vector<ptrdiff_t> zeta_fwrd, zeta_back, zeta_mixed, zeta_static;
+  const ptrdiff_t n_fwrd, n_back, n_mixed, n_static, n_back_mixed, n_fwrd_mixed, n_dynamic;
+  std::vector<ptrdiff_t> zeta_fwrd_mixed, zeta_back_mixed, zeta_dynamic,
     beta_back, beta_fwrd, pi_back, pi_fwrd;
-  Matrix S, A, D, E, Z_prime;
+  MatrixXd S, A, D, E, Z_prime;
   QRDecomposition QR;
   GeneralizedSchurDecomposition GSD;
-  LUSolver LU1, LU2, LU3;
-  Matrix Z21, g_y_back, g_y_back_tmp;
-  Matrix g_y_static, A0s, A0d, g_y_dynamic, g_y_static_tmp;
-  Matrix g_u_tmp1, g_u_tmp2;
-  LUSolver LU4;
+  ColPivHouseholderQR<MatrixXd> Solver1, Solver2, Solver3, Solver4;
+  MatrixXd Z21, g_y_back, g_y_fwrd;
+  MatrixXd g_y_static, A0s, A0d, g_y_dynamic;
+  MatrixXd g_u_tmp1, g_u_tmp2;
 public:
   class BlanchardKahnException
   {
@@ -54,24 +55,25 @@ public:
   /*!
     The zetas are supposed to follow C convention (first vector index is zero).
   */
-  DecisionRules(size_t n_arg, size_t p_arg, const std::vector<size_t> &zeta_fwrd_arg,
-                const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg,
-                const std::vector<size_t> &zeta_static_arg, double qz_criterium);
+  DecisionRules(ptrdiff_t n_arg, ptrdiff_t p_arg, const std::vector<ptrdiff_t> &zeta_fwrd_arg,
+                const std::vector<ptrdiff_t> &zeta_back_arg, const std::vector<ptrdiff_t> &zeta_mixed_arg,
+                const std::vector<ptrdiff_t> &zeta_static_arg, double qz_criterium);
   virtual ~DecisionRules(){};
 
   /*!
     \param jacobian First columns are backetermined vars at t-1 (in the order of zeta_back_mixed), then all vars at t (in the orig order), then forward vars at t+1 (in the order of zeta_fwrd_mixed), then exogenous vars.
   */
-  void compute(const Matrix &jacobian, Matrix &g_y, Matrix &g_u) throw (BlanchardKahnException, GeneralizedSchurDecomposition::GSDException);
-  template<class Vec1, class Vec2>
-  void getGeneralizedEigenvalues(Vec1 &eig_real, Vec2 &eig_cmplx);
+  void compute(const MatrixXd &jacobian, MatrixXd &g_y, MatrixXd &g_u) throw (BlanchardKahnException, GeneralizedSchurDecomposition::GSDException);
+  template<typename Vec1, typename Vec2>
+  void getGeneralizedEigenvalues(PlainObjectBase<Vec1> &eig_real, PlainObjectBase<Vec2> &eig_cmplx);
 };
 
 std::ostream &operator<<(std::ostream &out, const DecisionRules::BlanchardKahnException &e);
 
 template<class Vec1, class Vec2>
 void
-DecisionRules::getGeneralizedEigenvalues(Vec1 &eig_real, Vec2 &eig_cmplx)
+DecisionRules::getGeneralizedEigenvalues(PlainObjectBase<Vec1> &eig_real,
+                                         PlainObjectBase<Vec2> &eig_cmplx)
 {
   GSD.getGeneralizedEigenvalues(eig_real, eig_cmplx);
 }

mex/sources/estimation/ModelSolution.cc

@@ -30,18 +30,19 @@
 /**
  * compute the steady state (2nd stage), and computes first order approximation
  */
-ModelSolution::ModelSolution(const std::string &dynamicDllFile,  size_t n_endo_arg, size_t n_exo_arg, const std::vector<size_t> &zeta_fwrd_arg,
-                             const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg,
-                             const std::vector<size_t> &zeta_static_arg, double INqz_criterium) :
+ModelSolution::ModelSolution(const std::string &dynamicDllFile, ptrdiff_t n_endo_arg,
+                             ptrdiff_t n_exo_arg, const std::vector<ptrdiff_t> &zeta_fwrd_arg,
+                             const std::vector<ptrdiff_t> &zeta_back_arg, const std::vector<ptrdiff_t> &zeta_mixed_arg,
+                             const std::vector<ptrdiff_t> &zeta_static_arg, double INqz_criterium) :
   n_endo(n_endo_arg), n_exo(n_exo_arg),  // n_jcols = Num of Jacobian columns = nStat+2*nPred+3*nBoth+2*nForw+nExog
   n_jcols(n_exo+n_endo+ zeta_back_arg.size() /*nsPred*/ + zeta_fwrd_arg.size() /*nsForw*/ +2*zeta_mixed_arg.size()),
   jacobian(n_endo, n_jcols), residual(n_endo), Mx(1, n_exo),
   decisionRules(n_endo_arg, n_exo_arg, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, INqz_criterium),
-  dynamicDLLp(dynamicDllFile, n_exo),
+  dynamicDLLp(dynamicDllFile),
   llXsteadyState(n_jcols-n_exo)
 {
-  Mx.setAll(0.0);
-  jacobian.setAll(0.0);
+  Mx.setZero();
+  jacobian.setZero();
 
   set_union(zeta_fwrd_arg.begin(), zeta_fwrd_arg.end(),
             zeta_mixed_arg.begin(), zeta_mixed_arg.end(),

mex/sources/estimation/ModelSolution.hh

@@ -26,6 +26,10 @@
 #if !defined(ModelSolution_5ADFF920_9C74_46f5_9FE9_88AD4D4BBF19__INCLUDED_)
 #define ModelSolution_5ADFF920_9C74_46f5_9FE9_88AD4D4BBF19__INCLUDED_
 
+#include <Eigen/Core>
+
+using namespace Eigen;
+
 #include "DecisionRules.hh"
 #include "dynamic_dll.hh"
 
@@ -38,12 +42,14 @@ class ModelSolution
 {
 
 public:
-  ModelSolution(const std::string &dynamicDllFile,  size_t n_endo, size_t n_exo, const std::vector<size_t> &zeta_fwrd_arg,
-                const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg,
-                const std::vector<size_t> &zeta_static_arg, double qz_criterium);
+  ModelSolution(const std::string &dynamicDllFile, ptrdiff_t n_endo, ptrdiff_t n_exo, const std::vector<ptrdiff_t> &zeta_fwrd_arg,
+                const std::vector<ptrdiff_t> &zeta_back_arg, const std::vector<ptrdiff_t> &zeta_mixed_arg,
+                const std::vector<ptrdiff_t> &zeta_static_arg, double qz_criterium);
   virtual ~ModelSolution() {};
-  template <class VEC>
-  void compute(VEC &steadyState, const VectorView &deepParams,      Matrix &ghx, Matrix &ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
+  template <class Derived1, class Derived2>
+  void compute(PlainObjectBase<Derived1> &steadyState, const PlainObjectBase<Derived2> &deepParams,
+               MatrixXd &ghx, MatrixXd &ghu) throw (DecisionRules::BlanchardKahnException,
+                                                    GeneralizedSchurDecomposition::GSDException)
   {
     // compute Steady State
     ComputeSteadyState(steadyState, deepParams);
@@ -55,31 +61,31 @@ public:
   }
 
 private:
-  const size_t n_endo;
-  const size_t n_exo;
-  const size_t n_jcols; // Num of Jacobian columns
-  std::vector<size_t> zeta_fwrd_mixed, zeta_back_mixed;
-  Matrix jacobian;
-  Vector residual;
-  Matrix Mx;
+  const ptrdiff_t n_endo;
+  const ptrdiff_t n_exo;
+  const ptrdiff_t n_jcols; // Num of Jacobian columns
+  std::vector<ptrdiff_t> zeta_fwrd_mixed, zeta_back_mixed;
+  MatrixXd jacobian;
+  VectorXd residual;
+  MatrixXd Mx;
   DecisionRules decisionRules;
   DynamicModelDLL dynamicDLLp;
-  Vector llXsteadyState;
+  VectorXd llXsteadyState;
   //Matrix jacobian;
-  template <class VEC>
-  void ComputeModelSolution(VEC &steadyState, const VectorView &deepParams,         
-			    Matrix &ghx, Matrix &ghu) 
+  template <class Derived1, class Derived2>
+  void ComputeModelSolution(PlainObjectBase<Derived1> &steadyState, const PlainObjectBase<Derived2> &deepParams,         
+			    MatrixXd &ghx, MatrixXd &ghu) 
     throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
   {
     // set extended Steady State
 
-    for (size_t i = 0; i < zeta_back_mixed.size(); i++)
+    for (ptrdiff_t i = 0; i < (ptrdiff_t) zeta_back_mixed.size(); i++)
       llXsteadyState(i) = steadyState(zeta_back_mixed[i]);
 
-    for (size_t i = 0; i < n_endo; i++)
+    for (ptrdiff_t i = 0; i < n_endo; i++)
       llXsteadyState(zeta_back_mixed.size() + i) = steadyState(i);
 
-    for (size_t i = 0; i < zeta_fwrd_mixed.size(); i++)
+    for (ptrdiff_t i = 0; i < (ptrdiff_t) zeta_fwrd_mixed.size(); i++)
       llXsteadyState(zeta_back_mixed.size() + n_endo + i) = steadyState(zeta_fwrd_mixed[i]);
 
     //get jacobian
@@ -88,8 +94,8 @@ private:
     //compute rules
     decisionRules.compute(jacobian, ghx, ghu);
   }
-  template <class VEC>
-  void ComputeSteadyState(VEC &steadyState, const VectorView &deepParams)
+  template <class Derived1, class Derived2>
+  void ComputeSteadyState(PlainObjectBase<Derived1> &steadyState, const PlainObjectBase<Derived2> &deepParams)
   {
     // does nothig for time being.
   }

mex/sources/estimation/libmat/GeneralizedSchurDecomposition.cc

 /*
- * Copyright (C) 2010-2011 Dynare Team
+ * Copyright (C) 2010-2012 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -24,7 +24,7 @@
 
 double GeneralizedSchurDecomposition::criterium_static;
 
-GeneralizedSchurDecomposition::GeneralizedSchurDecomposition(size_t n_arg, double criterium_arg) :
+GeneralizedSchurDecomposition::GeneralizedSchurDecomposition(ptrdiff_t n_arg, double criterium_arg) :
   n(n_arg), criterium(criterium_arg)
 {
   alphar = new double[n];

mex/sources/estimation/libmat/GeneralizedSchurDecomposition.hh

 /*
- * Copyright (C) 2010 Dynare Team
+ * Copyright (C) 2010-2012 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -19,13 +19,14 @@
 
 #include <dynlapack.h>
 
-#include "Vector.hh"
-#include "Matrix.hh"
+#include <Eigen/Core>
+
+using namespace Eigen;
 
 class GeneralizedSchurDecomposition
 {
 private:
-  const size_t n;
+  const ptrdiff_t n;
   const double criterium;
   lapack_int lwork;
   double *alphar, *alphai, *beta, *vsl, *work;
@@ -40,38 +41,44 @@ public:
     GSDException(lapack_int info_arg, lapack_int n_arg) : info(info_arg), n(n_arg) {};
   };
   //! \todo Replace heuristic choice for workspace size by a query to determine the optimal size
-  GeneralizedSchurDecomposition(size_t n_arg, double criterium_arg);
+  GeneralizedSchurDecomposition(ptrdiff_t n_arg, double criterium_arg);
   virtual ~GeneralizedSchurDecomposition();
   //! \todo Add a lock around the modification of criterium_static for making it thread-safe
-  template<class Mat1, class Mat2, class Mat3>
-  void compute(Mat1 &S, Mat2 &T, Mat3 &Z, size_t &sdim) throw (GSDException);
-  template<class Mat1, class Mat2, class Mat3, class Mat4, class Mat5>
+  template<typename Mat1, typename Mat2, typename Mat3>
+  void compute(PlainObjectBase<Mat1> &S, PlainObjectBase<Mat2> &T,
+               PlainObjectBase<Mat3> &Z, ptrdiff_t &sdim) throw (GSDException);
+
   /*!
     \param[out] sdim Number of non-explosive generalized eigenvalues
   */
-  void compute(const Mat1 &D, const Mat2 &E, Mat3 &S, Mat4 &T, Mat5 &Z, size_t &sdim) throw (GSDException);
-  template<class Vec1, class Vec2>
-  void getGeneralizedEigenvalues(Vec1 &eig_real, Vec2 &eig_cmplx);
+  template<typename Mat1, typename Mat2, typename Mat3, typename Mat4, typename Mat5>
+  void compute(const MatrixBase<Mat1> &D, const MatrixBase<Mat2> &E,
+               PlainObjectBase<Mat3> &S, PlainObjectBase<Mat4> &T,
+               PlainObjectBase<Mat5> &Z, ptrdiff_t &sdim) throw (GSDException);
+ 
+  template<typename Vec1, typename Vec2>
+  void getGeneralizedEigenvalues(PlainObjectBase<Vec1> &eig_real,
+                                 PlainObjectBase<Vec2> &eig_cmplx);
 };
 
 std::ostream &operator<<(std::ostream &out, const GeneralizedSchurDecomposition::GSDException &e);
 
-template<class Mat1, class Mat2, class Mat3>
+template<typename Mat1, typename Mat2, typename Mat3>
 void
-GeneralizedSchurDecomposition::compute(Mat1 &S, Mat2 &T, Mat3 &Z, size_t &sdim) throw (GSDException)
+GeneralizedSchurDecomposition::compute(PlainObjectBase<Mat1> &S, PlainObjectBase<Mat2> &T, PlainObjectBase<Mat3> &Z, ptrdiff_t &sdim) throw (GSDException)
 {
-  assert(S.getRows() == n && S.getCols() == n
-         && T.getRows() == n && T.getCols() == n
-         && Z.getRows() == n && Z.getCols() == n);
+  assert(S.rows() == n && S.cols() == n
+         && T.rows() == n && T.cols() == n
+         && Z.rows() == n && Z.cols() == n);
 
   lapack_int n2 = n;
   lapack_int info, sdim2;
-  lapack_int lds = S.getLd(), ldt = T.getLd(), ldz = Z.getLd();
+  lapack_int lds = S.outerStride(), ldt = T.outerStride(), ldz = Z.outerStride();
 
   criterium_static = criterium;
   // Here we are forced to give space for left Schur vectors, even if we don't use them, because of a bug in dgges()
-  dgges("N", "V", "S", &selctg, &n2, S.getData(), &lds, T.getData(), &ldt,
-        &sdim2, alphar, alphai, beta, vsl, &n2, Z.getData(), &ldz,
+  dgges("N", "V", "S", &selctg, &n2, S.data(), &lds, T.data(), &ldt,
+        &sdim2, alphar, alphai, beta, vsl, &n2, Z.data(), &ldz,
         work, &lwork, bwork, &info);
 
   if (info != 0)
@@ -80,26 +87,27 @@ GeneralizedSchurDecomposition::compute(Mat1 &S, Mat2 &T, Mat3 &Z, size_t &sdim)
   sdim = sdim2;
 }
 
-template<class Mat1, class Mat2, class Mat3, class Mat4, class Mat5>
+template<typename Mat1, typename Mat2, typename Mat3, typename Mat4, typename Mat5>
 void
-GeneralizedSchurDecomposition::compute(const Mat1 &D, const Mat2 &E,
-                                       Mat3 &S, Mat4 &T, Mat5 &Z, size_t &sdim) throw (GSDException)
+GeneralizedSchurDecomposition::compute(const MatrixBase<Mat1> &D, const MatrixBase<Mat2> &E,
+                                       PlainObjectBase<Mat3> &S, PlainObjectBase<Mat4> &T, PlainObjectBase<Mat5> &Z, ptrdiff_t &sdim) throw (GSDException)
 {
-  assert(D.getRows() == n && D.getCols() == n
-         && E.getRows() == n && E.getCols() == n);
+  assert(D.rows() == n && D.cols() == n
+         && E.rows() == n && E.cols() == n);
   S = D;
   T = E;
   compute(S, T, Z, sdim);
 }
 
-template<class Vec1, class Vec2>
+template<typename Vec1, typename Vec2>
 void
-GeneralizedSchurDecomposition::getGeneralizedEigenvalues(Vec1 &eig_real, Vec2 &eig_cmplx)
+GeneralizedSchurDecomposition::getGeneralizedEigenvalues(PlainObjectBase<Vec1> &eig_real,
+                                                         PlainObjectBase<Vec2> &eig_cmplx)
 {
-  assert(eig_real.getSize() == n && eig_cmplx.getSize() == n);
+  assert(eig_real.size() == n && eig_cmplx.size() == n);
 
   double *par = alphar, *pai = alphai, *pb = beta,
-    *per = eig_real.getData(), *pei = eig_cmplx.getData();
+    *per = eig_real.data(), *pei = eig_cmplx.data();
   while (par < alphar + n)
     {
       *per = *par / *pb;
@@ -110,7 +118,7 @@ GeneralizedSchurDecomposition::getGeneralizedEigenvalues(Vec1 &eig_real, Vec2 &e
       par++;
       pai++;
       pb++;
-      per += eig_real.getStride();
-      pei += eig_cmplx.getStride();
+      per += eig_real.innerStride();
+      pei += eig_cmplx.innerStride();
     }
 }

mex/sources/estimation/libmat/QRDecomposition.cc

 /*
- * Copyright (C) 2010 Dynare Team
+ * Copyright (C) 2010-2012 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -21,7 +21,7 @@
 
 #include "QRDecomposition.hh"
 
-QRDecomposition::QRDecomposition(size_t rows_arg, size_t cols_arg, size_t cols2_arg) :
+QRDecomposition::QRDecomposition(ptrdiff_t rows_arg, ptrdiff_t cols_arg, ptrdiff_t cols2_arg) :
   rows(rows_arg), cols(cols_arg), mind(std::min(rows, cols)), cols2(cols2_arg),
   lwork(rows*cols), lwork2(cols2)
 {

mex/sources/estimation/libmat/QRDecomposition.hh

 /*
- * Copyright (C) 2010-2011 Dynare Team
+ * Copyright (C) 2010-2012 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -19,18 +19,22 @@
 
 #include <algorithm> // For std::min()
 
-#include <dynlapack.h>
+#include <Eigen/Core>
+
+using namespace Eigen;
 
-#include "Vector.hh"
-#include "Matrix.hh"
-#include "BlasBindings.hh"
+#include <dynlapack.h>
 
+/* Note that Eigen provides built-in QR decomposers (classes *HouseHolderQR).
+   However, we need to have our own because we need the ability to
+   left-multiply arbitrary matrices by Q, and this is not offered in a
+   convenient way by the Eigen classes. */
 class QRDecomposition
 {
 private:
-  const size_t rows, cols, mind;
+  const ptrdiff_t rows, cols, mind;
   //! Number of columns of the matrix to be left-multiplied by Q
-  const size_t cols2;
+  const ptrdiff_t cols2;
   lapack_int lwork, lwork2;
   double *work, *work2, *tau;
 public:
@@ -40,7 +44,7 @@ public:
     \param[in] cols_arg Number of columns of the matrix to decompose
     \param[in] cols2_arg Number of columns of the matrix to be multiplied by Q
   */
-  QRDecomposition(size_t rows_arg, size_t cols_arg, size_t cols2_arg);
+  QRDecomposition(ptrdiff_t rows_arg, ptrdiff_t cols_arg, ptrdiff_t cols2_arg);
   virtual ~QRDecomposition();
   //! Performs the QR decomposition of a matrix, and left-multiplies another matrix by Q
   /*!
@@ -48,25 +52,27 @@ public:
     \param[in] trans Specifies whether Q should be transposed before the multiplication, either "T" or "N"
     \param[in,out] C The matrix to be left-multiplied by Q, modified in place
   */
-  template<class Mat1, class Mat2>
-  void computeAndLeftMultByQ(Mat1 &A, const char *trans, Mat2 &C);
+  template<typename Derived1, typename Derived2>
+  void computeAndLeftMultByQ(PlainObjectBase<Derived1> &A, const char *trans,
+                             PlainObjectBase<Derived2> &C);
 };
 
-template<class Mat1, class Mat2>
+template<typename Derived1, typename Derived2>
 void
-QRDecomposition::computeAndLeftMultByQ(Mat1 &A, const char *trans, Mat2 &C)
+QRDecomposition::computeAndLeftMultByQ(PlainObjectBase<Derived1> &A, const char *trans,
+                                       PlainObjectBase<Derived2> &C)
 {
-  assert(A.getRows() == rows && A.getCols() == cols);
-  assert(C.getRows() == rows && C.getCols() == cols2);
+  assert(A.rows() == rows && A.cols() == cols);
+  assert(C.rows() == rows && C.cols() == cols2);
 
-  lapack_int m = rows, n = cols, lda = A.getLd();
+  lapack_int m = rows, n = cols, lda = A.outerStride();
   lapack_int info;
-  dgeqrf(&m, &n, A.getData(), &lda, tau, work, &lwork, &info);
+  dgeqrf(&m, &n, A.data(), &lda, tau, work, &lwork, &info);
   assert(info == 0);
 
   n = cols2;
-  lapack_int k = mind, ldc = C.getLd();
-  dormqr("L", trans, &m, &n, &k, A.getData(), &lda, tau, C.getData(), &ldc,
+  lapack_int k = mind, ldc = C.outerStride();
+  dormqr("L", trans, &m, &n, &k, A.data(), &lda, tau, C.data(), &ldc,
          work2, &lwork2, &info);
   assert(info == 0);
 }

mex/sources/estimation/libmat/tests/Makefile.am

 check_PROGRAMS = test-qr test-gsd test-lu test-repmat
 
-test_qr_SOURCES = ../Matrix.cc ../Vector.cc ../QRDecomposition.cc test-qr.cc
+test_qr_SOURCES = ../QRDecomposition.cc test-qr.cc
 test_qr_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
 test_qr_CPPFLAGS = -I.. -I../../../
 
-test_gsd_SOURCES = ../Matrix.cc ../Vector.cc ../GeneralizedSchurDecomposition.cc test-gsd.cc
+test_gsd_SOURCES = ../GeneralizedSchurDecomposition.cc test-gsd.cc
 test_gsd_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
 test_gsd_CPPFLAGS = -I.. -I../../../
 

mex/sources/estimation/libmat/tests/test-gsd.cc

 /*
- * Copyright (C) 2010 Dynare Team
+ * Copyright (C) 2010-2012 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -24,26 +24,22 @@
 int
 main(int argc, char **argv)
 {
-  size_t n = 3;
-  double D_data[] = { 1, 2, 3,
-                      4, 5, 6,
-                      7, 8, 9 };
-  double E_data[] = { 1, -3, 4,
-                      -7, 9, 1,
-                      -3, 4, 0 };
-  MatrixView D(D_data, n, n, n), E(E_data, n, n, n);
-
-  // Need to transpose because internally matrices are in column-major order
-  mat::transpose(D);
-  mat::transpose(E);
+  Matrix3d D, E;
+  D << 1, 2, 3,
+    4, 5, 6,
+    7, 8, 9;
+  
+  E << 1, -3, 4,
+    -7, 9, 1,
+    -3, 4, 0;
 
   std::cout << "D =" << std::endl << D << std::endl;
   std::cout << "E =" << std::endl << E << std::endl;
 
-  GeneralizedSchurDecomposition GSD(n, 1.00001);
+  GeneralizedSchurDecomposition GSD(3, 1.00001);
 
-  Matrix S(n), T(n), Z(n);
-  size_t sdim;
+  Matrix3d S, T, Z;
+  ptrdiff_t sdim;
 
   GSD.compute(D, E, S, T, Z, sdim);
 
@@ -51,7 +47,7 @@ main(int argc, char **argv)
   std::cout << "T =" << std::endl << T << std::endl;
   std::cout << "Z =" << std::endl << Z << std::endl;
 
-  Vector eig_real(n), eig_cmplx(n);
+  Vector3d eig_real, eig_cmplx;
   GSD.getGeneralizedEigenvalues(eig_real, eig_cmplx);
 
   std::cout << "Real part of generalized eigenvalues: " << std::endl << eig_real << std::endl;

mex/sources/estimation/libmat/tests/test-qr.cc

 /*
- * Copyright (C) 2010 Dynare Team
+ * Copyright (C) 2010-2012 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -19,14 +19,13 @@
 
 #include <iostream>
 
-#include "BlasBindings.hh"
 #include "QRDecomposition.hh"
 
 int
 main(int argc, char **argv)
 {
   size_t m = 4, n = 3;
-  Matrix S(m, n), Q(m), A(m, n), B(m), S2(m, n);
+  MatrixXd S(m, n), Q(m, m), A(m, n), S2(m, n), B;
   QRDecomposition QRD(m, n, m);
 
   for (size_t i = 0; i < m; i++)
@@ -35,32 +34,27 @@ main(int argc, char **argv)
 
   std::cout << "Matrix to be decomposed:" << std::endl << S << std::endl;
 
-  mat::set_identity(Q);
+  Q.setIdentity();
 
   S2 = S;
   QRD.computeAndLeftMultByQ(S2, "N", Q);
 
   std::cout << "Q =" << std::endl << Q << std::endl;
 
-  blas::gemm("T", "N", 1.0, Q, Q, 0.0, B);
-
+  B = Q.transpose()*Q;
+  
   std::cout << "Q'*Q =" << std::endl << B << std::endl;
 
-  for (size_t j = 0; j < n; j++)
-    mat::col_set(S2, j, j+1, m-j-1, 0);
+  S2 = S2.triangularView<Upper>();
 
   std::cout << "R =" << std::endl << S2 << std::endl;
 
-  blas::gemm("N", "N", 1.0, Q, S2, 0.0, A);
+  A = Q*S2;
 
   std::cout << "Q*R =" << std::endl << A << std::endl;
 
   // Check values
-  Matrix B2(m);
-  mat::set_identity(B2);
-  mat::sub(B2, B);
-  assert(mat::nrminf(B2) < 1e-4);
+  assert((B - MatrixXd::Identity(m,m)).lpNorm<Infinity>() < 1e-4);
 
-  mat::sub(A, S);
-  assert(mat::nrminf(A) < 1e-4);
+  assert((A-S).lpNorm<Infinity>() < 1e-4);
 }

mex/sources/estimation/tests/Makefile.am

 check_PROGRAMS = test-dr testModelSolution testInitKalman testKalman testPDF
 
-test_dr_SOURCES = ../libmat/Matrix.cc ../libmat/Vector.cc ../libmat/QRDecomposition.cc ../libmat/GeneralizedSchurDecomposition.cc ../libmat/LUSolver.cc ../DecisionRules.cc test-dr.cc
+test_dr_SOURCES = ../libmat/GeneralizedSchurDecomposition.cc ../libmat/QRDecomposition.cc ../DecisionRules.cc test-dr.cc
 test_dr_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
 test_dr_CPPFLAGS = -I.. -I../libmat -I../../
 
-testModelSolution_SOURCES = ../libmat/Matrix.cc ../libmat/Vector.cc ../libmat/QRDecomposition.cc ../libmat/GeneralizedSchurDecomposition.cc ../libmat/LUSolver.cc ../utils/dynamic_dll.cc ../DecisionRules.cc ../ModelSolution.cc testModelSolution.cc
+testModelSolution_SOURCES = ../libmat/QRDecomposition.cc ../libmat/GeneralizedSchurDecomposition.cc ../utils/dynamic_dll.cc ../DecisionRules.cc ../ModelSolution.cc testModelSolution.cc
 testModelSolution_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) $(LIBADD_DLOPEN)
 testModelSolution_CPPFLAGS = -I.. -I../libmat -I../../ -I../utils
 

mex/sources/estimation/tests/model1.mod

+/* To be used with testModelSolution; derived from fs2000.mod */
+
+var m P c e W R k d n l gy_obs gp_obs y dA ;
+varexo e_a e_m;
+
+parameters alp bet gam mst rho psi del;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+
+model (use_dll);
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
+P*c = m;
+m-1+d = l;
+e = exp(e_a);
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+initval;
+m = mst;
+P = 2.25;
+c = 0.45;
+e = 1;
+W = 4;
+R = 1.02;
+k = 6;
+d = 0.85;
+n = 0.19;
+l = 0.86;
+y = 0.6;
+gy_obs = exp(gam);
+gp_obs = exp(-gam); 
+dA = exp(gam);
+end;
+
+shocks;
+var e_a; stderr 0.196e-3;
+var e_m; stderr 0.0250e-3;
+end;
+
+steady;
+
+stoch_simul(nograph,noprint,order=1);
+
+// Output ghx and ghu using the ordering used by the DLL
+
+[s,i]=sort(oo_.dr.order_var(oo_.dr.nstatic+1:oo_.dr.npred+oo_.dr.nstatic));
+oo_.dr.ghx(oo_.dr.inv_order_var,i)
+
+oo_.dr.ghu(oo_.dr.inv_order_var,:)

mex/sources/estimation/tests/test-dr.cc

@@ -3,7 +3,7 @@
  */
 
 /*
- * Copyright (C) 2010-2011 Dynare Team
+ * Copyright (C) 2010-2012 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -21,14 +21,16 @@
  * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include <iostream>
+
 #include "DecisionRules.hh"
 
 int
 main(int argc, char **argv)
 {
-  size_t endo_nbr = 6, exo_nbr = 2;
+  ptrdiff_t endo_nbr = 6, exo_nbr = 2;
 
-  std::vector<size_t> zeta_fwrd, zeta_back, zeta_mixed, zeta_static;
+  std::vector<ptrdiff_t> zeta_fwrd, zeta_back, zeta_mixed, zeta_static;
   // y and c are purely forward
   zeta_fwrd.push_back(0);
   zeta_fwrd.push_back(1);
@@ -132,10 +134,8 @@ main(int argc, char **argv)
     -1.000000000000000
   };
 
-  MatrixView jacob_tmp(jacob_data, 6, 14, 6);
-
-  Matrix jacobian(6, 14), g_y(6, 3), g_u(6, 2);
-  jacobian = jacob_tmp;
+  MatrixXd jacobian(6, 14), g_y(6, 3), g_u(6, 2);
+  jacobian = Map<Matrix<double,6,14> >(jacob_data);
 
   try
     {
@@ -150,44 +150,32 @@ main(int argc, char **argv)
       std::cerr << e << std::endl;
     }
 
-  Vector eig_real(6), eig_cmplx(6);
+  VectorXd eig_real(6), eig_cmplx(6);
   dr.getGeneralizedEigenvalues(eig_real, eig_cmplx);
-  std::cout << "Eigenvalues (real part): " << eig_real
-            << "Eigenvalues (complex part): " << eig_cmplx << std::endl
-            << "g_y = " << std::endl << g_y << std::endl
-            << "g_u = " << std::endl << g_u;
+  std::cout << "Eigenvalues (real part): " << std::endl << eig_real << std::endl << std::endl
+            << "Eigenvalues (complex part): " << std::endl << eig_cmplx << std::endl<< std::endl
+            << "g_y = " << std::endl << g_y << std::endl << std::endl
+            << "g_u = " << std::endl << g_u << std::endl;
 
   // Check the results for g_y
-  double real_g_y_data[] = {
-    0.005358267364601, 1.836717147430803, 0.837085806295838,
+  MatrixXd real_g_y(6, 3);
+  real_g_y << 0.005358267364601, 1.836717147430803, 0.837085806295838,
     0.038541607674354, 0.424582606909411, -0.318740381721598,
     0.941816659690247, 1.419061793291772, 1.419061793291773,
     -0.000000000000000, 0.950000000000000, 0.025000000000000,
     -0.012546516642830, 0.341714987626857, 0.341714987626861,
-    0.000000000000000, 0.025000000000000, 0.950000000000000
-  };
-
-  MatrixView real_g_y_prime(real_g_y_data, 3, 6, 3);
-  Matrix real_g_y(6, 3);
-  mat::transpose(real_g_y, real_g_y_prime);
-  mat::sub(real_g_y, g_y);
+    0.000000000000000, 0.025000000000000, 0.950000000000000;
 
-  assert(mat::nrminf(real_g_y) < 1e-13);
+  assert((real_g_y - g_y).lpNorm<Infinity>() < 1e-13);
 
   // Check the results for g_u
-  double real_g_u_data[] = {
-    1.911522267389459, 0.830839736432740,
+  MatrixXd real_g_u(6, 2);
+  real_g_u << 1.911522267389459, 0.830839736432740,
     0.456074274269694, -0.347518145871938,
     1.455447993119765, 1.455447993119767,
     1.000000000000000, 0,
     0.350476910386520, 0.350476910386525,
-    0, 1.000000000000000
-  };
-
-  MatrixView real_g_u_prime(real_g_u_data, 2, 6, 2);
-  Matrix real_g_u(6, 2);
-  mat::transpose(real_g_u, real_g_u_prime);
-  mat::sub(real_g_u, g_u);
-
-  assert(mat::nrminf(real_g_u) < 1e-13);
+    0, 1.000000000000000;
+  
+  assert((real_g_u - g_u).lpNorm<Infinity>() < 1e-13);
 }

mex/sources/estimation/tests/testModelSolution.cc

 /*
- * Copyright (C) 2010-2011 Dynare Team
+ * Copyright (C) 2010-2012 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -17,8 +17,9 @@
  * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-// Test  for ModelSolution
-// Uses fs2000k2.mod and its ..._dynamic.mexw32
+// Test for ModelSolution; to be used with model1.mod
+
+#include <iostream>
 
 #include "ModelSolution.hh"
 
@@ -27,56 +28,50 @@ main(int argc, char **argv)
 {
   if (argc < 2)
     {
-      std::cerr << argv[0] << ": please provide as argument the name of the dynamic DLL generated from fs2000k2.mod (typically fs2000k2_dynamic.mex*)" << std::endl;
+      std::cerr << argv[0] << ": please provide as argument the name of the dynamic DLL generated from model1.mod (typically model1_dynamic.mex*)" << std::endl;
       exit(EXIT_FAILURE);
     }
 
   std::string modName = argv[1];
   const int npar = 7;
-  const size_t n_endo = 15, n_exo = 2;
-  std::vector<size_t> zeta_fwrd_arg;
-  std::vector<size_t> zeta_back_arg;
-  std::vector<size_t> zeta_mixed_arg;
-  std::vector<size_t> zeta_static_arg;
+  const ptrdiff_t n_endo = 15, n_exo = 2;
+  std::vector<ptrdiff_t> zeta_fwrd_arg;
+  std::vector<ptrdiff_t> zeta_back_arg;
+  std::vector<ptrdiff_t> zeta_mixed_arg;
+  std::vector<ptrdiff_t> zeta_static_arg;
   double qz_criterium = 1.0+1.0e-9;
-  Vector deepParams(npar);
 
-  double dYSparams [] = {
+  Matrix2d vCov;
+  vCov <<
+    0.1960e-3, 0.0,
+    0.0, 0.0250e-3;
+  
+  VectorXd deepParams(npar);
+  deepParams << 0.3300,
+    0.9900,
+    0.0030,
+    1.0110,
+    0.7000,
+    0.7870,
+    0.0200;
+
+  VectorXd steadyState(n_endo);
+  steadyState <<
     1.0110,  2.2582,  0.4477,  1.0000,
     4.5959,  1.0212,  5.8012,  0.8494,
     0.1872,  0.8604,  1.0030,  1.0080,
-    0.5808,  1.0030,  2.2093
-  };
-  double vcov[] = {
-    0.1960e-3, 0.0,
-    0.0, 0.0250e-3
-  };
-  int nVCVpar = 2;
-  MatrixView vCovVW(vcov, nVCVpar, nVCVpar, nVCVpar);
-  Matrix vCov(nVCVpar, nVCVpar);
-  vCov = vCovVW;
-
-  double dparams[] = { 0.3300,
-                       0.9900,
-                       0.0030,
-                       1.0110,
-                       0.7000,
-                       0.7870,
-                       0.0200};
+    0.5808,  1.0030,  2.2093;
 
-  VectorView modParamsVW(dparams, npar, 1);
-  deepParams = modParamsVW;
-  VectorView steadyState(dYSparams, n_endo, 1);
   std::cout << "Vector deepParams: " << std::endl << deepParams << std::endl;
   std::cout << "Matrix vCov: " << std::endl << vCov << std::endl;
   std::cout << "Vector steadyState: " << std::endl << steadyState << std::endl;
 
   // Set zeta vectors [0:(n-1)] from Matlab indices [1:n]
   //order_var = [ stat_var(:); pred_var(:); both_var(:); fwrd_var(:)];
-  size_t statc[] = { 4, 5, 6, 8, 9, 10, 11, 12, 14};
-  size_t back[] = {1, 7, 13};
-  size_t both[] = {2};
-  size_t fwd[] = { 3, 15};
+  ptrdiff_t statc[] = { 4, 5, 6, 8, 9, 10, 11, 12, 14};
+  ptrdiff_t back[] = {1, 7, 13};
+  ptrdiff_t both[] = {2};
+  ptrdiff_t fwd[] = { 3, 15};
   for (int i = 0; i < 9; ++i)
     zeta_static_arg.push_back(statc[i]-1);
   for (int i = 0; i < 3; ++i)
@@ -86,8 +81,8 @@ main(int argc, char **argv)
   for (int i = 0; i < 2; ++i)
     zeta_fwrd_arg.push_back(fwd[i]-1);
 
-  Matrix ghx(n_endo, zeta_back_arg.size() + zeta_mixed_arg.size());
-  Matrix ghu(n_endo, n_exo);
+  MatrixXd ghx(n_endo, zeta_back_arg.size() + zeta_mixed_arg.size());
+  MatrixXd ghu(n_endo, n_exo);
 
   ModelSolution modelSolution(modName, n_endo, n_exo,
                               zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium);

mex/sources/estimation/utils/dynamic_dll.cc

@@ -19,12 +19,12 @@
 
 #include "dynamic_dll.hh"
 
+#include <iostream>
 #include <sstream>
 
 using namespace std;
 
-DynamicModelDLL::DynamicModelDLL(const std::string &dynamicDllFile, size_t n_exog_arg) throw (TSException) :
-  n_exog(n_exog_arg)
+DynamicModelDLL::DynamicModelDLL(const std::string &dynamicDllFile) throw (TSException)
 {
   std::string fName;
 #if !defined(__CYGWIN32__) && !defined(_WIN32)

mex/sources/estimation/utils/dynamic_dll.hh

@@ -27,7 +27,9 @@
 #endif
 
 #include <string>
-#include "Matrix.hh"
+#include <Eigen/Core>
+
+using namespace Eigen;
 
 #include "ts_exception.h"
 
@@ -44,7 +46,6 @@ class DynamicModelDLL
 {
 private:
   DynamicFn Dynamic; // pointer to the Dynamic function in DLL
-  const size_t n_exog; // no of exogenous
 #if defined(_WIN32) || defined(__CYGWIN32__)
   HINSTANCE dynamicHinstance;  // DLL instance pointer in Windows
 #else
@@ -53,17 +54,16 @@ private:
 
 public:
   // construct and load Dynamic model DLL
-  DynamicModelDLL(const std::string &dynamicDllFile, size_t n_exog_arg) throw (TSException);
+  DynamicModelDLL(const std::string &dynamicDllFile) throw (TSException);
   virtual ~DynamicModelDLL();
 
   //! evaluate Dynamic model DLL
-  void eval(const Vector &y, const Matrix &x, const Vector &params, VectorView &ySteady,
-            Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException);
-  template<class VEC>
-  void eval(const Vector &y, const Matrix &x, const VectorView &modParams, VEC &ySteady,
-                      Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException)
+  template<class Derived1, class Derived2>
+  void eval(const VectorXd &y, const MatrixXd &x, const PlainObjectBase<Derived1> &modParams,
+            PlainObjectBase<Derived2> &ySteady, VectorXd &residual, MatrixXd *g1, MatrixXd *g2,
+            MatrixXd *g3) throw (TSException)
   {
-    Dynamic(y.getData(), x.getData(), 1, modParams.getData(), ySteady.getData(), 0, residual.getData(),
-	    g1 == NULL ? NULL : g1->getData(), g2 == NULL ? NULL : g2->getData(), g3 == NULL ? NULL : g3->getData());
+    Dynamic(y.data(), x.data(), 1, modParams.data(), ySteady.data(), 0, residual.data(),
+	    g1 == NULL ? NULL : g1->data(), g2 == NULL ? NULL : g2->data(), g3 == NULL ? NULL : g3->data());
   };
 };