Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
Loading items

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
  • chskcau/dynare-doc-fixes
27 results
Select Git revision
Loading items
Show changes
Commits on Source (7)
Showing
with 354 additions and 318 deletions
/* /*
* Copyright (C) 2010-2011 Dynare Team * Copyright (C) 2010-2012 Dynare Team
* *
* This file is part of Dynare. * This file is part of Dynare.
* *
...@@ -23,11 +23,11 @@ ...@@ -23,11 +23,11 @@
#include "DecisionRules.hh" #include "DecisionRules.hh"
DecisionRules::DecisionRules(size_t n_arg, size_t p_arg, DecisionRules::DecisionRules(ptrdiff_t n_arg, ptrdiff_t p_arg,
const std::vector<size_t> &zeta_fwrd_arg, const std::vector<ptrdiff_t> &zeta_fwrd_arg,
const std::vector<size_t> &zeta_back_arg, const std::vector<ptrdiff_t> &zeta_back_arg,
const std::vector<size_t> &zeta_mixed_arg, const std::vector<ptrdiff_t> &zeta_mixed_arg,
const std::vector<size_t> &zeta_static_arg, const std::vector<ptrdiff_t> &zeta_static_arg,
double qz_criterium) : double qz_criterium) :
n(n_arg), p(p_arg), zeta_fwrd(zeta_fwrd_arg), zeta_back(zeta_back_arg), 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), zeta_mixed(zeta_mixed_arg), zeta_static(zeta_static_arg),
...@@ -37,25 +37,24 @@ DecisionRules::DecisionRules(size_t n_arg, size_t p_arg, ...@@ -37,25 +37,24 @@ DecisionRules::DecisionRules(size_t n_arg, size_t p_arg,
n_dynamic(n-n_static), n_dynamic(n-n_static),
S(n, n_static), S(n, n_static),
A(n, n_back_mixed + n + n_fwrd_mixed), A(n, n_back_mixed + n + n_fwrd_mixed),
D(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), E(n_fwrd + n_back + 2*n_mixed, n_fwrd + n_back + 2*n_mixed),
Z_prime(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), QR(n, n_static, n_back_mixed + n + n_fwrd_mixed),
GSD(n_fwrd + n_back + 2*n_mixed, qz_criterium), GSD(n_fwrd + n_back + 2*n_mixed, qz_criterium),
LU1(n_fwrd_mixed), Solver1(n_fwrd_mixed, n_fwrd_mixed),
LU2(n_back_mixed), Solver2(n_back_mixed, n_back_mixed),
LU3(n_static), Solver3(n_static, n_static),
Solver4(n, n),
Z21(n_fwrd_mixed, n_back_mixed), Z21(n_fwrd_mixed, n_back_mixed),
g_y_back(n_back_mixed), g_y_back(n_back_mixed, n_back_mixed),
g_y_back_tmp(n_back_mixed), g_y_fwrd(n_fwrd_mixed, n_fwrd_mixed),
g_y_static(n_static, n_back_mixed), g_y_static(n_static, n_back_mixed),
A0s(n_static), A0s(n_static, n_static),
A0d(n_static, n_dynamic), A0d(n_static, n_dynamic),
g_y_dynamic(n_dynamic, n_back_mixed), 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_tmp1(n, n_back_mixed),
g_u_tmp2(n), g_u_tmp2(n, n)
LU4(n)
{ {
assert(n == n_back + n_fwrd + n_mixed + n_static); assert(n == n_back + n_fwrd + n_mixed + n_static);
...@@ -70,7 +69,7 @@ DecisionRules::DecisionRules(size_t n_arg, size_t p_arg, ...@@ -70,7 +69,7 @@ DecisionRules::DecisionRules(size_t n_arg, size_t p_arg,
back_inserter(zeta_dynamic)); back_inserter(zeta_dynamic));
// Compute beta_back and pi_back // 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]) if (find(zeta_mixed.begin(), zeta_mixed.end(), zeta_back_mixed[i])
== zeta_mixed.end()) == zeta_mixed.end())
pi_back.push_back(i); pi_back.push_back(i);
...@@ -78,7 +77,7 @@ DecisionRules::DecisionRules(size_t n_arg, size_t p_arg, ...@@ -78,7 +77,7 @@ DecisionRules::DecisionRules(size_t n_arg, size_t p_arg,
beta_back.push_back(i); beta_back.push_back(i);
// Compute beta_fwrd and pi_fwrd // 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]) if (find(zeta_mixed.begin(), zeta_mixed.end(), zeta_fwrd_mixed[i])
== zeta_mixed.end()) == zeta_mixed.end())
pi_fwrd.push_back(i); pi_fwrd.push_back(i);
...@@ -87,115 +86,92 @@ DecisionRules::DecisionRules(size_t n_arg, size_t p_arg, ...@@ -87,115 +86,92 @@ DecisionRules::DecisionRules(size_t n_arg, size_t p_arg,
} }
void 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 assert(jacobian.rows() == n
&& jacobian.getCols() == (n_back_mixed + n + n_fwrd_mixed + p)); && jacobian.cols() == (n_back_mixed + n + n_fwrd_mixed + p));
assert(g_y.getRows() == n && g_y.getCols() == n_back_mixed); assert(g_y.rows() == n && g_y.cols() == n_back_mixed);
assert(g_u.getRows() == n && g_u.getCols() == p); assert(g_u.rows() == n && g_u.cols() == p);
// Construct S, perform QR decomposition and get A = Q*jacobian // Construct S, perform QR decomposition and get A = Q*jacobian
for (size_t i = 0; i < n_static; i++) for (ptrdiff_t i = 0; i < n_static; i++)
mat::col_copy(jacobian, n_back_mixed+zeta_static[i], S, 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); QR.computeAndLeftMultByQ(S, "T", A);
// Construct matrix D // Construct matrix D
D.setAll(0.0); D.setZero();
for (size_t i = 0; i < n_mixed; i++) for (ptrdiff_t i = 0; i < n_mixed; i++)
D(n - n_static + i, beta_back[i]) = 1.0; D(n - n_static + i, beta_back[i]) = 1.0;
for (size_t j = 0; j < n_back_mixed; j++) for (ptrdiff_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.col(j).head(n-n_static) = A.col(n_back_mixed + zeta_back_mixed[j]).tail(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); 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 // Construct matrix E
E.setAll(0.0); E.setZero();
for (size_t i = 0; i < n_mixed; i++) for (ptrdiff_t i = 0; i < n_mixed; i++)
E(n - n_static + i, n_back_mixed + beta_fwrd[i]) = 1.0; 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); E.block(0, 0, n - n_static, n_back_mixed) = -A.block(n_static, 0, n - n_static, n_back_mixed);
for (size_t j = 0; j < n_fwrd; j++) for (ptrdiff_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.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);
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
// Perform the generalized Schur // Perform the generalized Schur
size_t sdim; ptrdiff_t sdim;
GSD.compute(E, D, Z_prime, sdim); GSD.compute(E, D, Z_prime, sdim);
if (n_back_mixed != sdim) if (n_back_mixed != sdim)
throw BlanchardKahnException(true, n_fwrd_mixed, n_fwrd + n_back + 2*n_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 // Compute DR for forward variables w.r. to endogenous
MatrixView Z21_prime(Z_prime, 0, n_back_mixed, n_back_mixed, n_fwrd_mixed), Solver1.compute(Z_prime.block(n_back_mixed, n_back_mixed, n_fwrd_mixed, n_fwrd_mixed).transpose());
Z22_prime(Z_prime, n_back_mixed, n_back_mixed, n_fwrd_mixed, n_fwrd_mixed);
mat::transpose(Z21, Z21_prime);
try if (!Solver1.isInvertible())
{
LU1.invMult("T", Z22_prime, Z21);
}
catch (LUSolver::LUException &e)
{
throw BlanchardKahnException(false, n_fwrd_mixed, n_fwrd + n_back + 2*n_mixed - sdim); 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++) g_y_fwrd = -Solver1.solve(Z_prime.block(0, n_back_mixed, n_back_mixed, n_fwrd_mixed).transpose());
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 // Compute DR for backward variables w.r. to endogenous
MatrixView Z11_prime(Z_prime, 0, 0, n_back_mixed, n_back_mixed), Solver2.compute(Z_prime.block(0, 0, n_back_mixed, n_back_mixed));
T11(D, 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}
S11(E, 0, 0, n_back_mixed, n_back_mixed);
mat::set_identity(g_y_back); Solver2.compute(D.block(0, 0, n_back_mixed, n_back_mixed));
g_y_back_tmp = Z11_prime; g_y_back = Z_prime.block(0, 0, n_back_mixed, n_back_mixed) * Solver2.solve(g_y_back);
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);
// TODO: avoid to copy mixed variables again, rather test it... // TODO: avoid to copy mixed variables again, rather test it...
for (size_t i = 0; i < n_back_mixed; i++) for (ptrdiff_t i = 0; i < n_back_mixed; i++)
mat::row_copy(g_y_back, i, g_y, zeta_back_mixed[i]); g_y.row(zeta_back_mixed[i]) = g_y_back.row(i);
// Compute DR for static variables w.r. to endogenous // Compute DR for static variables w.r. to endogenous
g_y_static = MatrixView(A, 0, 0, n_static, n_back_mixed); for (ptrdiff_t i = 0; i < n_dynamic; i++)
for (size_t i = 0; i < n_dynamic; i++)
{ {
mat::row_copy(g_y, zeta_dynamic[i], g_y_dynamic, i); g_y_dynamic.row(i) = g_y.row(zeta_dynamic[i]);
mat::col_copy(A, n_back_mixed + zeta_dynamic[i], 0, n_static, A0d, i, 0); 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); for (ptrdiff_t i = 0; i < n_static; i++)
blas::gemm("N", "N", 1.0, MatrixView(A, 0, n_back_mixed + n, n_static, n_fwrd_mixed), A0s.col(i) = A.col(n_back_mixed + zeta_static[i]).head(n_static);
g_y_static_tmp, 1.0, g_y_static);
for (size_t i = 0; i < n_static; i++) Solver3.compute(A0s);
mat::col_copy(A, n_back_mixed + zeta_static[i], 0, n_static, A0s, i, 0);
LU3.invMult("N", A0s, g_y_static); g_y_static = -Solver3.solve(A.block(0, n_back_mixed + n, n_static, n_fwrd_mixed)*g_y_fwrd*g_y_back
mat::negate(g_y_static); + A0d*g_y_dynamic + A.block(0, 0, n_static, n_back_mixed));
for (size_t i = 0; i < n_static; i++) for (ptrdiff_t i = 0; i < n_static; i++)
mat::row_copy(g_y_static, i, g_y, zeta_static[i]); g_y.row(zeta_static[i]) = g_y_static.row(i);
// Compute DR for all endogenous w.r. to shocks // 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_tmp1 = jacobian.block(0, n_back_mixed + n, n, n_fwrd_mixed) * g_y_fwrd;
g_u_tmp2 = MatrixConstView(jacobian, 0, n_back_mixed, n, n);
for (size_t i = 0; i < n_back_mixed; i++) g_u_tmp2 = jacobian.block(0, n_back_mixed, n, n);
{ for (ptrdiff_t i = 0; i < n_back_mixed; i++)
VectorView c1 = mat::get_col(g_u_tmp2, zeta_back_mixed[i]), g_u_tmp2.col(zeta_back_mixed[i]) += g_u_tmp1.col(i);
c2 = mat::get_col(g_u_tmp1, i); Solver4.compute(g_u_tmp2);
vec::add(c1, c2); g_u = -Solver4.solve(jacobian.block(0, n_back_mixed + n + n_fwrd_mixed, n, p));
}
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);
} }
std::ostream & std::ostream &
......
/* /*
* Copyright (C) 2010 Dynare Team * Copyright (C) 2010-2012 Dynare Team
* *
* This file is part of Dynare. * This file is part of Dynare.
* *
...@@ -20,28 +20,29 @@ ...@@ -20,28 +20,29 @@
#include <cstdlib> #include <cstdlib>
#include <vector> #include <vector>
#include "Vector.hh" #include <Eigen/Core>
#include "Matrix.hh" #include <Eigen/QR>
#include "QRDecomposition.hh" #include "QRDecomposition.hh"
#include "GeneralizedSchurDecomposition.hh" #include "GeneralizedSchurDecomposition.hh"
#include "LUSolver.hh"
using namespace Eigen;
class DecisionRules class DecisionRules
{ {
private: private:
const size_t n, p; const ptrdiff_t n, p;
const std::vector<size_t> zeta_fwrd, zeta_back, zeta_mixed, zeta_static; const std::vector<ptrdiff_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; const ptrdiff_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, std::vector<ptrdiff_t> zeta_fwrd_mixed, zeta_back_mixed, zeta_dynamic,
beta_back, beta_fwrd, pi_back, pi_fwrd; beta_back, beta_fwrd, pi_back, pi_fwrd;
Matrix S, A, D, E, Z_prime; MatrixXd S, A, D, E, Z_prime;
QRDecomposition QR; QRDecomposition QR;
GeneralizedSchurDecomposition GSD; GeneralizedSchurDecomposition GSD;
LUSolver LU1, LU2, LU3; ColPivHouseholderQR<MatrixXd> Solver1, Solver2, Solver3, Solver4;
Matrix Z21, g_y_back, g_y_back_tmp; MatrixXd Z21, g_y_back, g_y_fwrd;
Matrix g_y_static, A0s, A0d, g_y_dynamic, g_y_static_tmp; MatrixXd g_y_static, A0s, A0d, g_y_dynamic;
Matrix g_u_tmp1, g_u_tmp2; MatrixXd g_u_tmp1, g_u_tmp2;
LUSolver LU4;
public: public:
class BlanchardKahnException class BlanchardKahnException
{ {
...@@ -54,24 +55,25 @@ public: ...@@ -54,24 +55,25 @@ public:
/*! /*!
The zetas are supposed to follow C convention (first vector index is zero). 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, DecisionRules(ptrdiff_t n_arg, ptrdiff_t p_arg, const std::vector<ptrdiff_t> &zeta_fwrd_arg,
const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg, const std::vector<ptrdiff_t> &zeta_back_arg, const std::vector<ptrdiff_t> &zeta_mixed_arg,
const std::vector<size_t> &zeta_static_arg, double qz_criterium); const std::vector<ptrdiff_t> &zeta_static_arg, double qz_criterium);
virtual ~DecisionRules(){}; 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. \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); void compute(const MatrixXd &jacobian, MatrixXd &g_y, MatrixXd &g_u) throw (BlanchardKahnException, GeneralizedSchurDecomposition::GSDException);
template<class Vec1, class Vec2> template<typename Vec1, typename Vec2>
void getGeneralizedEigenvalues(Vec1 &eig_real, Vec2 &eig_cmplx); void getGeneralizedEigenvalues(PlainObjectBase<Vec1> &eig_real, PlainObjectBase<Vec2> &eig_cmplx);
}; };
std::ostream &operator<<(std::ostream &out, const DecisionRules::BlanchardKahnException &e); std::ostream &operator<<(std::ostream &out, const DecisionRules::BlanchardKahnException &e);
template<class Vec1, class Vec2> template<class Vec1, class Vec2>
void 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); GSD.getGeneralizedEigenvalues(eig_real, eig_cmplx);
} }
...@@ -30,18 +30,19 @@ ...@@ -30,18 +30,19 @@
/** /**
* compute the steady state (2nd stage), and computes first order approximation * 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, ModelSolution::ModelSolution(const std::string &dynamicDllFile, ptrdiff_t n_endo_arg,
const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg, ptrdiff_t n_exo_arg, const std::vector<ptrdiff_t> &zeta_fwrd_arg,
const std::vector<size_t> &zeta_static_arg, double INqz_criterium) : 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_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()), 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), 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), 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) llXsteadyState(n_jcols-n_exo)
{ {
Mx.setAll(0.0); Mx.setZero();
jacobian.setAll(0.0); jacobian.setZero();
set_union(zeta_fwrd_arg.begin(), zeta_fwrd_arg.end(), set_union(zeta_fwrd_arg.begin(), zeta_fwrd_arg.end(),
zeta_mixed_arg.begin(), zeta_mixed_arg.end(), zeta_mixed_arg.begin(), zeta_mixed_arg.end(),
......
...@@ -26,6 +26,10 @@ ...@@ -26,6 +26,10 @@
#if !defined(ModelSolution_5ADFF920_9C74_46f5_9FE9_88AD4D4BBF19__INCLUDED_) #if !defined(ModelSolution_5ADFF920_9C74_46f5_9FE9_88AD4D4BBF19__INCLUDED_)
#define 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 "DecisionRules.hh"
#include "dynamic_dll.hh" #include "dynamic_dll.hh"
...@@ -38,12 +42,14 @@ class ModelSolution ...@@ -38,12 +42,14 @@ class ModelSolution
{ {
public: public:
ModelSolution(const std::string &dynamicDllFile, size_t n_endo, size_t n_exo, const std::vector<size_t> &zeta_fwrd_arg, ModelSolution(const std::string &dynamicDllFile, ptrdiff_t n_endo, ptrdiff_t n_exo, const std::vector<ptrdiff_t> &zeta_fwrd_arg,
const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg, const std::vector<ptrdiff_t> &zeta_back_arg, const std::vector<ptrdiff_t> &zeta_mixed_arg,
const std::vector<size_t> &zeta_static_arg, double qz_criterium); const std::vector<ptrdiff_t> &zeta_static_arg, double qz_criterium);
virtual ~ModelSolution() {}; virtual ~ModelSolution() {};
template <class VEC> template <class Derived1, class Derived2>
void compute(VEC &steadyState, const VectorView &deepParams, Matrix &ghx, Matrix &ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException) void compute(PlainObjectBase<Derived1> &steadyState, const PlainObjectBase<Derived2> &deepParams,
MatrixXd &ghx, MatrixXd &ghu) throw (DecisionRules::BlanchardKahnException,
GeneralizedSchurDecomposition::GSDException)
{ {
// compute Steady State // compute Steady State
ComputeSteadyState(steadyState, deepParams); ComputeSteadyState(steadyState, deepParams);
...@@ -55,31 +61,31 @@ public: ...@@ -55,31 +61,31 @@ public:
} }
private: private:
const size_t n_endo; const ptrdiff_t n_endo;
const size_t n_exo; const ptrdiff_t n_exo;
const size_t n_jcols; // Num of Jacobian columns const ptrdiff_t n_jcols; // Num of Jacobian columns
std::vector<size_t> zeta_fwrd_mixed, zeta_back_mixed; std::vector<ptrdiff_t> zeta_fwrd_mixed, zeta_back_mixed;
Matrix jacobian; MatrixXd jacobian;
Vector residual; VectorXd residual;
Matrix Mx; MatrixXd Mx;
DecisionRules decisionRules; DecisionRules decisionRules;
DynamicModelDLL dynamicDLLp; DynamicModelDLL dynamicDLLp;
Vector llXsteadyState; VectorXd llXsteadyState;
//Matrix jacobian; //Matrix jacobian;
template <class VEC> template <class Derived1, class Derived2>
void ComputeModelSolution(VEC &steadyState, const VectorView &deepParams, void ComputeModelSolution(PlainObjectBase<Derived1> &steadyState, const PlainObjectBase<Derived2> &deepParams,
Matrix &ghx, Matrix &ghu) MatrixXd &ghx, MatrixXd &ghu)
throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
{ {
// set extended Steady State // 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]); 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); 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]); llXsteadyState(zeta_back_mixed.size() + n_endo + i) = steadyState(zeta_fwrd_mixed[i]);
//get jacobian //get jacobian
...@@ -88,8 +94,8 @@ private: ...@@ -88,8 +94,8 @@ private:
//compute rules //compute rules
decisionRules.compute(jacobian, ghx, ghu); decisionRules.compute(jacobian, ghx, ghu);
} }
template <class VEC> template <class Derived1, class Derived2>
void ComputeSteadyState(VEC &steadyState, const VectorView &deepParams) void ComputeSteadyState(PlainObjectBase<Derived1> &steadyState, const PlainObjectBase<Derived2> &deepParams)
{ {
// does nothig for time being. // does nothig for time being.
} }
......
/* /*
* Copyright (C) 2010-2011 Dynare Team * Copyright (C) 2010-2012 Dynare Team
* *
* This file is part of Dynare. * This file is part of Dynare.
* *
...@@ -24,7 +24,7 @@ ...@@ -24,7 +24,7 @@
double GeneralizedSchurDecomposition::criterium_static; 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) n(n_arg), criterium(criterium_arg)
{ {
alphar = new double[n]; alphar = new double[n];
......
/* /*
* Copyright (C) 2010 Dynare Team * Copyright (C) 2010-2012 Dynare Team
* *
* This file is part of Dynare. * This file is part of Dynare.
* *
...@@ -19,13 +19,14 @@ ...@@ -19,13 +19,14 @@
#include <dynlapack.h> #include <dynlapack.h>
#include "Vector.hh" #include <Eigen/Core>
#include "Matrix.hh"
using namespace Eigen;
class GeneralizedSchurDecomposition class GeneralizedSchurDecomposition
{ {
private: private:
const size_t n; const ptrdiff_t n;
const double criterium; const double criterium;
lapack_int lwork; lapack_int lwork;
double *alphar, *alphai, *beta, *vsl, *work; double *alphar, *alphai, *beta, *vsl, *work;
...@@ -40,38 +41,44 @@ public: ...@@ -40,38 +41,44 @@ public:
GSDException(lapack_int info_arg, lapack_int n_arg) : info(info_arg), n(n_arg) {}; 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 //! \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(); virtual ~GeneralizedSchurDecomposition();
//! \todo Add a lock around the modification of criterium_static for making it thread-safe //! \todo Add a lock around the modification of criterium_static for making it thread-safe
template<class Mat1, class Mat2, class Mat3> template<typename Mat1, typename Mat2, typename Mat3>
void compute(Mat1 &S, Mat2 &T, Mat3 &Z, size_t &sdim) throw (GSDException); void compute(PlainObjectBase<Mat1> &S, PlainObjectBase<Mat2> &T,
template<class Mat1, class Mat2, class Mat3, class Mat4, class Mat5> PlainObjectBase<Mat3> &Z, ptrdiff_t &sdim) throw (GSDException);
/*! /*!
\param[out] sdim Number of non-explosive generalized eigenvalues \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<typename Mat1, typename Mat2, typename Mat3, typename Mat4, typename Mat5>
template<class Vec1, class Vec2> void compute(const MatrixBase<Mat1> &D, const MatrixBase<Mat2> &E,
void getGeneralizedEigenvalues(Vec1 &eig_real, Vec2 &eig_cmplx); 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); 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 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 assert(S.rows() == n && S.cols() == n
&& T.getRows() == n && T.getCols() == n && T.rows() == n && T.cols() == n
&& Z.getRows() == n && Z.getCols() == n); && Z.rows() == n && Z.cols() == n);
lapack_int n2 = n; lapack_int n2 = n;
lapack_int info, sdim2; 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; 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() // 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, dgges("N", "V", "S", &selctg, &n2, S.data(), &lds, T.data(), &ldt,
&sdim2, alphar, alphai, beta, vsl, &n2, Z.getData(), &ldz, &sdim2, alphar, alphai, beta, vsl, &n2, Z.data(), &ldz,
work, &lwork, bwork, &info); work, &lwork, bwork, &info);
if (info != 0) if (info != 0)
...@@ -80,26 +87,27 @@ GeneralizedSchurDecomposition::compute(Mat1 &S, Mat2 &T, Mat3 &Z, size_t &sdim) ...@@ -80,26 +87,27 @@ GeneralizedSchurDecomposition::compute(Mat1 &S, Mat2 &T, Mat3 &Z, size_t &sdim)
sdim = sdim2; 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 void
GeneralizedSchurDecomposition::compute(const Mat1 &D, const Mat2 &E, GeneralizedSchurDecomposition::compute(const MatrixBase<Mat1> &D, const MatrixBase<Mat2> &E,
Mat3 &S, Mat4 &T, Mat5 &Z, size_t &sdim) throw (GSDException) PlainObjectBase<Mat3> &S, PlainObjectBase<Mat4> &T, PlainObjectBase<Mat5> &Z, ptrdiff_t &sdim) throw (GSDException)
{ {
assert(D.getRows() == n && D.getCols() == n assert(D.rows() == n && D.cols() == n
&& E.getRows() == n && E.getCols() == n); && E.rows() == n && E.cols() == n);
S = D; S = D;
T = E; T = E;
compute(S, T, Z, sdim); compute(S, T, Z, sdim);
} }
template<class Vec1, class Vec2> template<typename Vec1, typename Vec2>
void 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, 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) while (par < alphar + n)
{ {
*per = *par / *pb; *per = *par / *pb;
...@@ -110,7 +118,7 @@ GeneralizedSchurDecomposition::getGeneralizedEigenvalues(Vec1 &eig_real, Vec2 &e ...@@ -110,7 +118,7 @@ GeneralizedSchurDecomposition::getGeneralizedEigenvalues(Vec1 &eig_real, Vec2 &e
par++; par++;
pai++; pai++;
pb++; pb++;
per += eig_real.getStride(); per += eig_real.innerStride();
pei += eig_cmplx.getStride(); pei += eig_cmplx.innerStride();
} }
} }
/* /*
* Copyright (C) 2010 Dynare Team * Copyright (C) 2010-2012 Dynare Team
* *
* This file is part of Dynare. * This file is part of Dynare.
* *
...@@ -21,7 +21,7 @@ ...@@ -21,7 +21,7 @@
#include "QRDecomposition.hh" #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), rows(rows_arg), cols(cols_arg), mind(std::min(rows, cols)), cols2(cols2_arg),
lwork(rows*cols), lwork2(cols2) lwork(rows*cols), lwork2(cols2)
{ {
......
/* /*
* Copyright (C) 2010-2011 Dynare Team * Copyright (C) 2010-2012 Dynare Team
* *
* This file is part of Dynare. * This file is part of Dynare.
* *
...@@ -19,18 +19,22 @@ ...@@ -19,18 +19,22 @@
#include <algorithm> // For std::min() #include <algorithm> // For std::min()
#include <dynlapack.h> #include <Eigen/Core>
using namespace Eigen;
#include "Vector.hh" #include <dynlapack.h>
#include "Matrix.hh"
#include "BlasBindings.hh"
/* 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 class QRDecomposition
{ {
private: 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 //! Number of columns of the matrix to be left-multiplied by Q
const size_t cols2; const ptrdiff_t cols2;
lapack_int lwork, lwork2; lapack_int lwork, lwork2;
double *work, *work2, *tau; double *work, *work2, *tau;
public: public:
...@@ -40,7 +44,7 @@ public: ...@@ -40,7 +44,7 @@ public:
\param[in] cols_arg Number of columns of the matrix to decompose \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 \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(); virtual ~QRDecomposition();
//! Performs the QR decomposition of a matrix, and left-multiplies another matrix by Q //! Performs the QR decomposition of a matrix, and left-multiplies another matrix by Q
/*! /*!
...@@ -48,25 +52,27 @@ public: ...@@ -48,25 +52,27 @@ public:
\param[in] trans Specifies whether Q should be transposed before the multiplication, either "T" or "N" \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 \param[in,out] C The matrix to be left-multiplied by Q, modified in place
*/ */
template<class Mat1, class Mat2> template<typename Derived1, typename Derived2>
void computeAndLeftMultByQ(Mat1 &A, const char *trans, Mat2 &C); void computeAndLeftMultByQ(PlainObjectBase<Derived1> &A, const char *trans,
PlainObjectBase<Derived2> &C);
}; };
template<class Mat1, class Mat2> template<typename Derived1, typename Derived2>
void 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(A.rows() == rows && A.cols() == cols);
assert(C.getRows() == rows && C.getCols() == cols2); 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; 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); assert(info == 0);
n = cols2; n = cols2;
lapack_int k = mind, ldc = C.getLd(); lapack_int k = mind, ldc = C.outerStride();
dormqr("L", trans, &m, &n, &k, A.getData(), &lda, tau, C.getData(), &ldc, dormqr("L", trans, &m, &n, &k, A.data(), &lda, tau, C.data(), &ldc,
work2, &lwork2, &info); work2, &lwork2, &info);
assert(info == 0); assert(info == 0);
} }
check_PROGRAMS = test-qr test-gsd test-lu test-repmat 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_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
test_qr_CPPFLAGS = -I.. -I../../../ 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_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
test_gsd_CPPFLAGS = -I.. -I../../../ test_gsd_CPPFLAGS = -I.. -I../../../
......
/* /*
* Copyright (C) 2010 Dynare Team * Copyright (C) 2010-2012 Dynare Team
* *
* This file is part of Dynare. * This file is part of Dynare.
* *
...@@ -24,26 +24,22 @@ ...@@ -24,26 +24,22 @@
int int
main(int argc, char **argv) main(int argc, char **argv)
{ {
size_t n = 3; Matrix3d D, E;
double D_data[] = { 1, 2, 3, D << 1, 2, 3,
4, 5, 6, 4, 5, 6,
7, 8, 9 }; 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 E << 1, -3, 4,
mat::transpose(D); -7, 9, 1,
mat::transpose(E); -3, 4, 0;
std::cout << "D =" << std::endl << D << std::endl; std::cout << "D =" << std::endl << D << std::endl;
std::cout << "E =" << std::endl << E << 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); Matrix3d S, T, Z;
size_t sdim; ptrdiff_t sdim;
GSD.compute(D, E, S, T, Z, sdim); GSD.compute(D, E, S, T, Z, sdim);
...@@ -51,7 +47,7 @@ main(int argc, char **argv) ...@@ -51,7 +47,7 @@ main(int argc, char **argv)
std::cout << "T =" << std::endl << T << std::endl; std::cout << "T =" << std::endl << T << std::endl;
std::cout << "Z =" << std::endl << Z << 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); GSD.getGeneralizedEigenvalues(eig_real, eig_cmplx);
std::cout << "Real part of generalized eigenvalues: " << std::endl << eig_real << std::endl; std::cout << "Real part of generalized eigenvalues: " << std::endl << eig_real << std::endl;
......
/* /*
* Copyright (C) 2010 Dynare Team * Copyright (C) 2010-2012 Dynare Team
* *
* This file is part of Dynare. * This file is part of Dynare.
* *
...@@ -19,14 +19,13 @@ ...@@ -19,14 +19,13 @@
#include <iostream> #include <iostream>
#include "BlasBindings.hh"
#include "QRDecomposition.hh" #include "QRDecomposition.hh"
int int
main(int argc, char **argv) main(int argc, char **argv)
{ {
size_t m = 4, n = 3; 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); QRDecomposition QRD(m, n, m);
for (size_t i = 0; i < m; i++) for (size_t i = 0; i < m; i++)
...@@ -35,32 +34,27 @@ main(int argc, char **argv) ...@@ -35,32 +34,27 @@ main(int argc, char **argv)
std::cout << "Matrix to be decomposed:" << std::endl << S << std::endl; std::cout << "Matrix to be decomposed:" << std::endl << S << std::endl;
mat::set_identity(Q); Q.setIdentity();
S2 = S; S2 = S;
QRD.computeAndLeftMultByQ(S2, "N", Q); QRD.computeAndLeftMultByQ(S2, "N", Q);
std::cout << "Q =" << std::endl << Q << std::endl; 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; std::cout << "Q'*Q =" << std::endl << B << std::endl;
for (size_t j = 0; j < n; j++) S2 = S2.triangularView<Upper>();
mat::col_set(S2, j, j+1, m-j-1, 0);
std::cout << "R =" << std::endl << S2 << std::endl; 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; std::cout << "Q*R =" << std::endl << A << std::endl;
// Check values // Check values
Matrix B2(m); assert((B - MatrixXd::Identity(m,m)).lpNorm<Infinity>() < 1e-4);
mat::set_identity(B2);
mat::sub(B2, B);
assert(mat::nrminf(B2) < 1e-4);
mat::sub(A, S); assert((A-S).lpNorm<Infinity>() < 1e-4);
assert(mat::nrminf(A) < 1e-4);
} }
check_PROGRAMS = test-dr testModelSolution testInitKalman testKalman testPDF 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_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
test_dr_CPPFLAGS = -I.. -I../libmat -I../../ 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_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) $(LIBADD_DLOPEN)
testModelSolution_CPPFLAGS = -I.. -I../libmat -I../../ -I../utils testModelSolution_CPPFLAGS = -I.. -I../libmat -I../../ -I../utils
......
/* 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,:)
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
*/ */
/* /*
* Copyright (C) 2010-2011 Dynare Team * Copyright (C) 2010-2012 Dynare Team
* *
* This file is part of Dynare. * This file is part of Dynare.
* *
...@@ -21,14 +21,16 @@ ...@@ -21,14 +21,16 @@
* along with Dynare. If not, see <http://www.gnu.org/licenses/>. * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/ */
#include <iostream>
#include "DecisionRules.hh" #include "DecisionRules.hh"
int int
main(int argc, char **argv) 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 // y and c are purely forward
zeta_fwrd.push_back(0); zeta_fwrd.push_back(0);
zeta_fwrd.push_back(1); zeta_fwrd.push_back(1);
...@@ -132,10 +134,8 @@ main(int argc, char **argv) ...@@ -132,10 +134,8 @@ main(int argc, char **argv)
-1.000000000000000 -1.000000000000000
}; };
MatrixView jacob_tmp(jacob_data, 6, 14, 6); MatrixXd jacobian(6, 14), g_y(6, 3), g_u(6, 2);
jacobian = Map<Matrix<double,6,14> >(jacob_data);
Matrix jacobian(6, 14), g_y(6, 3), g_u(6, 2);
jacobian = jacob_tmp;
try try
{ {
...@@ -150,44 +150,32 @@ main(int argc, char **argv) ...@@ -150,44 +150,32 @@ main(int argc, char **argv)
std::cerr << e << std::endl; 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); dr.getGeneralizedEigenvalues(eig_real, eig_cmplx);
std::cout << "Eigenvalues (real part): " << eig_real std::cout << "Eigenvalues (real part): " << std::endl << eig_real << std::endl << std::endl
<< "Eigenvalues (complex part): " << eig_cmplx << std::endl << "Eigenvalues (complex part): " << std::endl << eig_cmplx << std::endl<< std::endl
<< "g_y = " << std::endl << g_y << std::endl << "g_y = " << std::endl << g_y << std::endl << std::endl
<< "g_u = " << std::endl << g_u; << "g_u = " << std::endl << g_u << std::endl;
// Check the results for g_y // Check the results for g_y
double real_g_y_data[] = { MatrixXd real_g_y(6, 3);
0.005358267364601, 1.836717147430803, 0.837085806295838, real_g_y << 0.005358267364601, 1.836717147430803, 0.837085806295838,
0.038541607674354, 0.424582606909411, -0.318740381721598, 0.038541607674354, 0.424582606909411, -0.318740381721598,
0.941816659690247, 1.419061793291772, 1.419061793291773, 0.941816659690247, 1.419061793291772, 1.419061793291773,
-0.000000000000000, 0.950000000000000, 0.025000000000000, -0.000000000000000, 0.950000000000000, 0.025000000000000,
-0.012546516642830, 0.341714987626857, 0.341714987626861, -0.012546516642830, 0.341714987626857, 0.341714987626861,
0.000000000000000, 0.025000000000000, 0.950000000000000 0.000000000000000, 0.025000000000000, 0.950000000000000;
};
MatrixView real_g_y_prime(real_g_y_data, 3, 6, 3); assert((real_g_y - g_y).lpNorm<Infinity>() < 1e-13);
Matrix real_g_y(6, 3);
mat::transpose(real_g_y, real_g_y_prime);
mat::sub(real_g_y, g_y);
assert(mat::nrminf(real_g_y) < 1e-13);
// Check the results for g_u // Check the results for g_u
double real_g_u_data[] = { MatrixXd real_g_u(6, 2);
1.911522267389459, 0.830839736432740, real_g_u << 1.911522267389459, 0.830839736432740,
0.456074274269694, -0.347518145871938, 0.456074274269694, -0.347518145871938,
1.455447993119765, 1.455447993119767, 1.455447993119765, 1.455447993119767,
1.000000000000000, 0, 1.000000000000000, 0,
0.350476910386520, 0.350476910386525, 0.350476910386520, 0.350476910386525,
0, 1.000000000000000 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); assert((real_g_u - g_u).lpNorm<Infinity>() < 1e-13);
} }
/* /*
* Copyright (C) 2010-2011 Dynare Team * Copyright (C) 2010-2012 Dynare Team
* *
* This file is part of Dynare. * This file is part of Dynare.
* *
...@@ -17,8 +17,9 @@ ...@@ -17,8 +17,9 @@
* along with Dynare. If not, see <http://www.gnu.org/licenses/>. * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/ */
// Test for ModelSolution // Test for ModelSolution; to be used with model1.mod
// Uses fs2000k2.mod and its ..._dynamic.mexw32
#include <iostream>
#include "ModelSolution.hh" #include "ModelSolution.hh"
...@@ -27,56 +28,50 @@ main(int argc, char **argv) ...@@ -27,56 +28,50 @@ main(int argc, char **argv)
{ {
if (argc < 2) 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); exit(EXIT_FAILURE);
} }
std::string modName = argv[1]; std::string modName = argv[1];
const int npar = 7; const int npar = 7;
const size_t n_endo = 15, n_exo = 2; const ptrdiff_t n_endo = 15, n_exo = 2;
std::vector<size_t> zeta_fwrd_arg; std::vector<ptrdiff_t> zeta_fwrd_arg;
std::vector<size_t> zeta_back_arg; std::vector<ptrdiff_t> zeta_back_arg;
std::vector<size_t> zeta_mixed_arg; std::vector<ptrdiff_t> zeta_mixed_arg;
std::vector<size_t> zeta_static_arg; std::vector<ptrdiff_t> zeta_static_arg;
double qz_criterium = 1.0+1.0e-9; double qz_criterium = 1.0+1.0e-9;
Vector deepParams(npar);
double dYSparams [] = { Matrix2d vCov;
1.0110, 2.2582, 0.4477, 1.0000, vCov <<
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.1960e-3, 0.0,
0.0, 0.0250e-3 0.0, 0.0250e-3;
};
int nVCVpar = 2;
MatrixView vCovVW(vcov, nVCVpar, nVCVpar, nVCVpar);
Matrix vCov(nVCVpar, nVCVpar);
vCov = vCovVW;
double dparams[] = { 0.3300, VectorXd deepParams(npar);
deepParams << 0.3300,
0.9900, 0.9900,
0.0030, 0.0030,
1.0110, 1.0110,
0.7000, 0.7000,
0.7870, 0.7870,
0.0200}; 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;
VectorView modParamsVW(dparams, npar, 1);
deepParams = modParamsVW;
VectorView steadyState(dYSparams, n_endo, 1);
std::cout << "Vector deepParams: " << std::endl << deepParams << std::endl; std::cout << "Vector deepParams: " << std::endl << deepParams << std::endl;
std::cout << "Matrix vCov: " << std::endl << vCov << std::endl; std::cout << "Matrix vCov: " << std::endl << vCov << std::endl;
std::cout << "Vector steadyState: " << std::endl << steadyState << std::endl; std::cout << "Vector steadyState: " << std::endl << steadyState << std::endl;
// Set zeta vectors [0:(n-1)] from Matlab indices [1:n] // Set zeta vectors [0:(n-1)] from Matlab indices [1:n]
//order_var = [ stat_var(:); pred_var(:); both_var(:); fwrd_var(:)]; //order_var = [ stat_var(:); pred_var(:); both_var(:); fwrd_var(:)];
size_t statc[] = { 4, 5, 6, 8, 9, 10, 11, 12, 14}; ptrdiff_t statc[] = { 4, 5, 6, 8, 9, 10, 11, 12, 14};
size_t back[] = {1, 7, 13}; ptrdiff_t back[] = {1, 7, 13};
size_t both[] = {2}; ptrdiff_t both[] = {2};
size_t fwd[] = { 3, 15}; ptrdiff_t fwd[] = { 3, 15};
for (int i = 0; i < 9; ++i) for (int i = 0; i < 9; ++i)
zeta_static_arg.push_back(statc[i]-1); zeta_static_arg.push_back(statc[i]-1);
for (int i = 0; i < 3; ++i) for (int i = 0; i < 3; ++i)
...@@ -86,8 +81,8 @@ main(int argc, char **argv) ...@@ -86,8 +81,8 @@ main(int argc, char **argv)
for (int i = 0; i < 2; ++i) for (int i = 0; i < 2; ++i)
zeta_fwrd_arg.push_back(fwd[i]-1); zeta_fwrd_arg.push_back(fwd[i]-1);
Matrix ghx(n_endo, zeta_back_arg.size() + zeta_mixed_arg.size()); MatrixXd ghx(n_endo, zeta_back_arg.size() + zeta_mixed_arg.size());
Matrix ghu(n_endo, n_exo); MatrixXd ghu(n_endo, n_exo);
ModelSolution modelSolution(modName, 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); zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium);
......
...@@ -19,12 +19,12 @@ ...@@ -19,12 +19,12 @@
#include "dynamic_dll.hh" #include "dynamic_dll.hh"
#include <iostream>
#include <sstream> #include <sstream>
using namespace std; using namespace std;
DynamicModelDLL::DynamicModelDLL(const std::string &dynamicDllFile, size_t n_exog_arg) throw (TSException) : DynamicModelDLL::DynamicModelDLL(const std::string &dynamicDllFile) throw (TSException)
n_exog(n_exog_arg)
{ {
std::string fName; std::string fName;
#if !defined(__CYGWIN32__) && !defined(_WIN32) #if !defined(__CYGWIN32__) && !defined(_WIN32)
......
...@@ -27,7 +27,9 @@ ...@@ -27,7 +27,9 @@
#endif #endif
#include <string> #include <string>
#include "Matrix.hh" #include <Eigen/Core>
using namespace Eigen;
#include "ts_exception.h" #include "ts_exception.h"
...@@ -44,7 +46,6 @@ class DynamicModelDLL ...@@ -44,7 +46,6 @@ class DynamicModelDLL
{ {
private: private:
DynamicFn Dynamic; // pointer to the Dynamic function in DLL DynamicFn Dynamic; // pointer to the Dynamic function in DLL
const size_t n_exog; // no of exogenous
#if defined(_WIN32) || defined(__CYGWIN32__) #if defined(_WIN32) || defined(__CYGWIN32__)
HINSTANCE dynamicHinstance; // DLL instance pointer in Windows HINSTANCE dynamicHinstance; // DLL instance pointer in Windows
#else #else
...@@ -53,17 +54,16 @@ private: ...@@ -53,17 +54,16 @@ private:
public: public:
// construct and load Dynamic model DLL // 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(); virtual ~DynamicModelDLL();
//! evaluate Dynamic model DLL //! evaluate Dynamic model DLL
void eval(const Vector &y, const Matrix &x, const Vector &params, VectorView &ySteady, template<class Derived1, class Derived2>
Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException); void eval(const VectorXd &y, const MatrixXd &x, const PlainObjectBase<Derived1> &modParams,
template<class VEC> PlainObjectBase<Derived2> &ySteady, VectorXd &residual, MatrixXd *g1, MatrixXd *g2,
void eval(const Vector &y, const Matrix &x, const VectorView &modParams, VEC &ySteady, MatrixXd *g3) throw (TSException)
Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException)
{ {
Dynamic(y.getData(), x.getData(), 1, modParams.getData(), ySteady.getData(), 0, residual.getData(), Dynamic(y.data(), x.data(), 1, modParams.data(), ySteady.data(), 0, residual.data(),
g1 == NULL ? NULL : g1->getData(), g2 == NULL ? NULL : g2->getData(), g3 == NULL ? NULL : g3->getData()); g1 == NULL ? NULL : g1->data(), g2 == NULL ? NULL : g2->data(), g3 == NULL ? NULL : g3->data());
}; };
}; };