diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.cc b/mex/sources/k_order_perturbation/k_ord_dynare.cc
index 9735e5918ca8a0bb1f1a3e189226d197ba9b1c5e..f02ba300d94a5d0423aa6c9e83aaf092d04424ca 100644
--- a/mex/sources/k_order_perturbation/k_ord_dynare.cc
+++ b/mex/sources/k_order_perturbation/k_ord_dynare.cc
@@ -1,5 +1,5 @@
/*
- * Copyright (C) 2008-2009 Dynare Team
+ * Copyright (C) 2008-2010 Dynare Team
*
* This file is part of Dynare.
*
@@ -35,54 +35,27 @@
KordpDynare::KordpDynare(const char **endo, int num_endo,
const char **exo, int nexog, int npar,
- Vector *ysteady, TwoDMatrix *vcov, Vector *inParams, int nstat,
- int npred, int nforw, int nboth, const int jcols, const Vector *nnzd,
+ Vector &ysteady, TwoDMatrix &vcov, Vector &inParams, int nstat,
+ int npred, int nforw, int nboth, const int jcols, const Vector &nnzd,
const int nsteps, int norder,
Journal &jr, DynamicModelDLL &dynamicDLL, double sstol,
- const vector<int> *var_order, const TwoDMatrix *llincidence, double criterium) throw (TLException) :
+ const vector<int> &var_order, const TwoDMatrix &llincidence, double criterium) throw (TLException) :
nStat(nstat), nBoth(nboth), nPred(npred), nForw(nforw), nExog(nexog), nPar(npar),
nYs(npred + nboth), nYss(nboth + nforw), nY(num_endo), nJcols(jcols), NNZD(nnzd), nSteps(nsteps),
nOrder(norder), journal(jr), dynamicDLL(dynamicDLL), ySteady(ysteady), vCov(vcov), params(inParams),
- md(1), dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(sstol), varOrder(var_order),
+ md(1), dnl(*this, endo), denl(*this, exo), dsnl(*this, dnl, denl), ss_tol(sstol), varOrder(var_order),
ll_Incidence(llincidence), qz_criterium(criterium)
{
- dnl = new DynareNameList(*this, endo);
- denl = new DynareExogNameList(*this, exo);
- dsnl = new DynareStateNameList(*this, *dnl, *denl);
-
- JacobianIndices = ReorderDynareJacobianIndices(varOrder);
+ ReorderDynareJacobianIndices();
// Initialise ModelDerivativeContainer(*this, this->md, nOrder);
for (int iord = 1; iord <= nOrder; iord++)
- {
- FSSparseTensor *t = new FSSparseTensor(iord, nY+nYs+nYss+nExog, nY);
- md.insert(t);
- }
-
+ md.insert(new FSSparseTensor(iord, nY+nYs+nYss+nExog, nY));
}
KordpDynare::~KordpDynare()
{
- if (ySteady)
- delete ySteady;
- if (params)
- delete params;
- if (vCov)
- delete vCov;
- if (dnl)
- delete dnl;
- if (dsnl)
- delete dsnl;
- if (denl)
- delete denl;
- if (JacobianIndices)
- delete JacobianIndices;
- if (varOrder)
- delete varOrder;
- if (ll_Incidence)
- delete ll_Incidence;
- if (NNZD)
- delete NNZD;
+ // No need to manually delete tensors in "md", they are deleted by the TensorContainer destructor
}
void
@@ -112,69 +85,55 @@ KordpDynare::evaluateSystem(Vector &out, const Vector &yym, const Vector &yy,
* which performs actual calculation and reorders
***************************************************/
void
-KordpDynare::calcDerivatives(const Vector &yy, const Vector &xx) throw (DynareException)
+KordpDynare::calcDerivativesAtSteady() throw (DynareException)
{
- TwoDMatrix *g2 = 0; // NULL;
- TwoDMatrix *g3 = 0; // NULL;
- TwoDMatrix *g1 = new TwoDMatrix(nY, nJcols); // generate g1 for jacobian
- g1->zeros();
+ TwoDMatrix g1(nY, nJcols);
+ g1.zeros();
- if ((nJcols != g1->ncols()) || (nY != g1->nrows()))
- throw DynareException(__FILE__, __LINE__, "Error in calcDerivatives: Created wrong jacobian");
+ TwoDMatrix *g2p = NULL, *g3p = NULL;
if (nOrder > 1)
{
- // generate g2 space for sparse Hessian 3x NNZH = 3x NNZD[1]
- g2 = new TwoDMatrix((int)(*NNZD)[1], 3);
- g2->zeros();
+ // allocate space for sparse Hessian
+ g2p = new TwoDMatrix((int) NNZD[1], 3);
+ g2p->zeros();
}
+
if (nOrder > 2)
{
- g3 = new TwoDMatrix((int)(*NNZD)[2], 3);
- g3->zeros();
+ g3p = new TwoDMatrix((int) NNZD[2], 3);
+ g3p->zeros();
}
+
+ Vector xx(nexog());
+ xx.zeros();
+
Vector out(nY);
out.zeros();
- const Vector *llxYYp; // getting around the constantness
- if ((nJcols - nExog) > yy.length())
- llxYYp = (LLxSteady(yy));
- else
- llxYYp = &yy;
- const Vector &llxYY = *(llxYYp);
-
- dynamicDLL.eval(llxYY, xx, params, out, g1, g2, g3);
+ Vector llxSteady(nJcols-nExog);
+ LLxSteady(ySteady, llxSteady);
- if ((nJcols != g1->ncols()) || (nY != g1->nrows()))
- throw DynareException(__FILE__, __LINE__, "Error in calcDerivatives: dynamicDLL.eval returned wrong jacobian");
+ dynamicDLL.eval(llxSteady, xx, ¶ms, out, &g1, g2p, g3p);
populateDerivativesContainer(g1, 1, JacobianIndices);
- delete g1;
+
if (nOrder > 1)
{
- populateDerivativesContainer(g2, 2, JacobianIndices);
- delete g2;
+ populateDerivativesContainer(*g2p, 2, JacobianIndices);
+ delete g2p;
}
if (nOrder > 2)
{
- populateDerivativesContainer(g3, 3, JacobianIndices);
- delete g3;
+ populateDerivativesContainer(*g3p, 3, JacobianIndices);
+ delete g3p;
}
- delete llxYYp;
-}
-
-void
-KordpDynare::calcDerivativesAtSteady() throw (DynareException)
-{
- Vector xx(nexog());
- xx.zeros();
- calcDerivatives(*ySteady, xx);
}
/*******************************************************************************
* populateDerivatives to sparse Tensor and fit it in the Derivatives Container
*******************************************************************************/
void
-KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector<int> *vOrder)
+KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const vector<int> &vOrder)
{
// model derivatives FSSparseTensor instance
FSSparseTensor *mdTi = (new FSSparseTensor(ord, nJcols, nY));
@@ -183,15 +142,15 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector<i
if (ord == 1)
{
- for (int i = 0; i < g->ncols(); i++)
+ for (int i = 0; i < g.ncols(); i++)
{
- for (int j = 0; j < g->nrows(); j++)
+ for (int j = 0; j < g.nrows(); j++)
{
double x;
if (s[0] < nJcols-nExog)
- x = g->get(j, (*vOrder)[s[0]]);
+ x = g.get(j, vOrder[s[0]]);
else
- x = g->get(j, s[0]);
+ x = g.get(j, s[0]);
if (x != 0.0)
mdTi->insert(s, j, x);
}
@@ -203,11 +162,11 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector<i
int nJcols1 = nJcols-nExog;
vector<int> revOrder(nJcols1);
for (int i = 0; i < nJcols1; i++)
- revOrder[(*vOrder)[i]] = i;
- for (int i = 0; i < g->nrows(); i++)
+ revOrder[vOrder[i]] = i;
+ for (int i = 0; i < g.nrows(); i++)
{
- int j = (int) g->get(i, 0)-1; // hessian indices start with 1
- int i1 = (int) g->get(i, 1) -1;
+ int j = (int) g.get(i, 0)-1; // hessian indices start with 1
+ int i1 = (int) g.get(i, 1) -1;
int s0 = (int) floor(((double) i1)/((double) nJcols));
int s1 = i1- (nJcols*s0);
if (s0 < nJcols1)
@@ -220,7 +179,7 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector<i
s[1] = s1;
if (s[1] >= s[0])
{
- double x = g->get(i, 2);
+ double x = g.get(i, 2);
mdTi->insert(s, j, x);
}
}
@@ -231,11 +190,11 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector<i
int nJcols2 = nJcols*nJcols;
vector<int> revOrder(nJcols1);
for (int i = 0; i < nJcols1; i++)
- revOrder[(*vOrder)[i]] = i;
- for (int i = 0; i < g->nrows(); i++)
+ revOrder[vOrder[i]] = i;
+ for (int i = 0; i < g.nrows(); i++)
{
- int j = (int) g->get(i, 0)-1;
- int i1 = (int) g->get(i, 1) -1;
+ int j = (int) g.get(i, 0)-1;
+ int i1 = (int) g.get(i, 1) -1;
int s0 = (int) floor(((double) i1)/((double) nJcols2));
int i2 = i1 - nJcols2*s0;
int s1 = (int) floor(((double) i2)/((double) nJcols));
@@ -254,7 +213,7 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector<i
s[2] = s2;
if ((s[2] >= s[1]) && (s[1] >= s[0]))
{
- double x = g->get(i, 2);
+ double x = g.get(i, 2);
mdTi->insert(s, j, x);
}
}
@@ -263,6 +222,7 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector<i
// md container
md.remove(Symmetry(ord));
md.insert(mdTi);
+ // No need to delete mdTi, it will be deleted by TensorContainer destructor
}
/*********************************************************
@@ -270,26 +230,26 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector<i
* returns ySteady extended with leads and lags suitable for
* passing to <model>_dynamic DLL
*************************************************************/
-Vector *
-KordpDynare::LLxSteady(const Vector &yS) throw (DynareException, TLException)
+void
+KordpDynare::LLxSteady(const Vector &yS, Vector &llxSteady) throw (DynareException, TLException)
{
if ((nJcols-nExog) == yS.length())
throw DynareException(__FILE__, __LINE__, "ySteady already of right size");
// create temporary square 2D matrix size nEndo x nEndo (sparse)
// for the lag, current and lead blocks of the jacobian
- Vector *llxSteady = new Vector(nJcols-nExog);
- for (int ll_row = 0; ll_row < ll_Incidence->nrows(); ll_row++)
+ if (llxSteady.length() != nJcols-nExog)
+ throw DynareException(__FILE__, __LINE__, "llxSteady has wrong size");
+
+ for (int ll_row = 0; ll_row < ll_Incidence.nrows(); ll_row++)
{
// populate (non-sparse) vector with ysteady values
for (int i = 0; i < nY; i++)
{
- if (ll_Incidence->get(ll_row, i))
- (*llxSteady)[((int) ll_Incidence->get(ll_row, i))-1] = yS[i];
+ if (ll_Incidence.get(ll_row, i))
+ llxSteady[((int) ll_Incidence.get(ll_row, i))-1] = yS[i];
}
}
-
- return llxSteady;
}
/************************************
@@ -317,27 +277,27 @@ KordpDynare::LLxSteady(const Vector &yS) throw (DynareException, TLException)
exogen
************************************/
-vector<int> *
-KordpDynare::ReorderDynareJacobianIndices(const vector<int> *varOrder) throw (TLException)
+void
+KordpDynare::ReorderDynareJacobianIndices() throw (TLException)
{
// create temporary square 2D matrix size nEndo x nEndo (sparse)
// for the lag, current and lead blocks of the jacobian
- vector<int> *JacobianIndices = new vector<int>(nJcols);
+ JacobianIndices.resize(nJcols);
vector <int> tmp(nY);
int i, j, rjoff = nJcols-nExog-1;
- for (int ll_row = 0; ll_row < ll_Incidence->nrows(); ll_row++)
+ for (int ll_row = 0; ll_row < ll_Incidence.nrows(); ll_row++)
{
// reorder in orde-var order & populate temporary nEndo (sparse) vector with
// the lag, current and lead blocks of the jacobian respectively
for (i = 0; i < nY; i++)
- tmp[i] = ((int) ll_Incidence->get(ll_row, (*varOrder)[i]-1));
+ tmp[i] = ((int) ll_Incidence.get(ll_row, varOrder[i]-1));
// write the reordered blocks back to the jacobian
// in reverse order
for (j = nY-1; j >= 0; j--)
if (tmp[j])
{
- (*JacobianIndices)[rjoff] = tmp[j] -1;
+ JacobianIndices[rjoff] = tmp[j] -1;
rjoff--;
if (rjoff < 0)
break;
@@ -346,9 +306,7 @@ KordpDynare::ReorderDynareJacobianIndices(const vector<int> *varOrder) throw (TL
//add the indices for the nExog exogenous jacobians
for (j = nJcols-nExog; j < nJcols; j++)
- (*JacobianIndices)[j] = j;
-
- return JacobianIndices;
+ JacobianIndices[j] = j;
}
/**************************************************************************************/
@@ -375,7 +333,7 @@ DynareNameList::selectIndices(const vector<const char *> &ns) const throw (Dynar
DynareNameList::DynareNameList(const KordpDynare &dynare)
{
for (int i = 0; i < dynare.ny(); i++)
- names.push_back(dynare.dnl->getName(i));
+ names.push_back(dynare.dnl.getName(i));
}
DynareNameList::DynareNameList(const KordpDynare &dynare, const char **namesp)
{
@@ -386,7 +344,7 @@ DynareNameList::DynareNameList(const KordpDynare &dynare, const char **namesp)
DynareExogNameList::DynareExogNameList(const KordpDynare &dynare)
{
for (int i = 0; i < dynare.nexog(); i++)
- names.push_back(dynare.denl->getName(i));
+ names.push_back(dynare.denl.getName(i));
}
DynareExogNameList::DynareExogNameList(const KordpDynare &dynare, const char **namesp)
diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.hh b/mex/sources/k_order_perturbation/k_ord_dynare.hh
index 1f1ee91399bdfb1e1f4c90d95d63f2a94b7c2d59..10ad4f261ccca7993112de3d59bf6a76788a87a5 100644
--- a/mex/sources/k_order_perturbation/k_ord_dynare.hh
+++ b/mex/sources/k_order_perturbation/k_ord_dynare.hh
@@ -122,31 +122,31 @@ class KordpDynare : public DynamicModel
const int nYss; // nyss ={ nboth + nforw ; }
const int nY; // = num_endo={ nstat + npred + nboth + nforw ; }
const int nJcols; // no of jacobian columns= nExog+nEndo+nsPred+nsForw
- const Vector *NNZD; /* the total number of non-zero derivative elements
+ const Vector &NNZD; /* the total number of non-zero derivative elements
where hessian is 2nd : NZZD(order=2) */
const int nSteps;
const int nOrder;
Journal &journal;
- Vector *ySteady;
- Vector *params;
- TwoDMatrix *vCov;
+ Vector &ySteady;
+ Vector ¶ms;
+ TwoDMatrix &vCov;
TensorContainer<FSSparseTensor> md; // ModelDerivatives
- DynareNameList *dnl;
- DynareExogNameList *denl;
- DynareStateNameList *dsnl;
+ DynareNameList dnl;
+ DynareExogNameList denl;
+ DynareStateNameList dsnl;
const double ss_tol;
- const vector<int> *varOrder;
- const TwoDMatrix *ll_Incidence;
+ const vector<int> &varOrder;
+ const TwoDMatrix &ll_Incidence;
double qz_criterium;
- vector<int> *JacobianIndices;
+ vector<int> JacobianIndices;
public:
KordpDynare(const char **endo, int num_endo,
const char **exo, int num_exo, int num_par,
- Vector *ySteady, TwoDMatrix *vCov, Vector *params, int nstat, int nPred,
- int nforw, int nboth, const int nJcols, const Vector *NNZD,
+ Vector &ySteady, TwoDMatrix &vCov, Vector ¶ms, int nstat, int nPred,
+ int nforw, int nboth, const int nJcols, const Vector &NNZD,
const int nSteps, const int ord,
Journal &jr, DynamicModelDLL &dynamicDLL, double sstol,
- const vector<int> *varOrder, const TwoDMatrix *ll_Incidence,
+ const vector<int> &varOrder, const TwoDMatrix &ll_Incidence,
double qz_criterium) throw (TLException);
virtual ~KordpDynare();
@@ -203,27 +203,27 @@ public:
const NameList &
getAllEndoNames() const
{
- return *dnl;
+ return dnl;
}
const NameList &
getStateNames() const
{
- return *dsnl;
+ return dsnl;
}
const NameList &
getExogNames() const
{
- return *denl;
+ return denl;
}
const TwoDMatrix &
getVcov() const
{
- return *vCov;
+ return vCov;
}
Vector &
getParams()
{
- return *params;
+ return params;
}
const TensorContainer<FSSparseTensor> &
@@ -234,31 +234,31 @@ public:
const Vector &
getSteady() const
{
- return *ySteady;
+ return ySteady;
}
Vector &
getSteady()
{
- return *ySteady;
+ return ySteady;
}
void solveDeterministicSteady();
void evaluateSystem(Vector &out, const Vector &yy, const Vector &xx) throw (DynareException);
void evaluateSystem(Vector &out, const Vector &yym, const Vector &yy,
const Vector &yyp, const Vector &xx) throw (DynareException);
- void calcDerivatives(const Vector &yy, const Vector &xx) throw (DynareException);
void calcDerivativesAtSteady() throw (DynareException);
DynamicModelDLL &dynamicDLL;
DynamicModel *
clone() const
{
- return new KordpDynare(*this);
+ std::cerr << "KordpDynare::clone() not implemented" << std::endl;
+ exit(EXIT_FAILURE);
}
- Vector *LLxSteady(const Vector &yS) throw (DynareException, TLException); // returns ySteady extended with leads and lags
+ void LLxSteady(const Vector &yS, Vector &llxSteady) throw (DynareException, TLException); // Given the steady state in yS, returns in llxSteady the steady state extended with leads and lags
private:
- vector<int> *ReorderDynareJacobianIndices(const vector<int> *varOrder) throw (TLException);
- void populateDerivativesContainer(TwoDMatrix *g, int ord, const vector<int> *vOrder);
+ void ReorderDynareJacobianIndices() throw (TLException);
+ void populateDerivativesContainer(const TwoDMatrix &g, int ord, const vector<int> &vOrder);
};
#endif
diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc
index 5047dc2a26e350fccfc6b74dd567e2037419b01d..93a8e31e0993a40c280d2a2e5a3f540799cf6559 100644
--- a/mex/sources/k_order_perturbation/k_order_perturbation.cc
+++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc
@@ -130,17 +130,17 @@ extern "C" {
mxFldp = mxGetField(M_, 0, "params");
double *dparams = (double *) mxGetData(mxFldp);
int npar = (int) mxGetM(mxFldp);
- Vector *modParams = new Vector(dparams, npar);
+ Vector modParams(dparams, npar);
mxFldp = mxGetField(M_, 0, "Sigma_e");
dparams = (double *) mxGetData(mxFldp);
npar = (int) mxGetN(mxFldp);
- TwoDMatrix *vCov = new TwoDMatrix(npar, npar, dparams);
+ TwoDMatrix vCov(npar, npar, dparams);
mxFldp = mxGetField(dr, 0, "ys"); // and not in order of dr.order_var
dparams = (double *) mxGetData(mxFldp);
const int nSteady = (int) mxGetM(mxFldp);
- Vector *ySteady = new Vector(dparams, nSteady);
+ Vector ySteady(dparams, nSteady);
mxFldp = mxGetField(dr, 0, "nstatic");
const int nStat = (int) mxGetScalar(mxFldp);
@@ -172,9 +172,9 @@ extern "C" {
npar = (int) mxGetM(mxFldp);
if (npar != nEndo) //(nPar != npar)
mexErrMsgTxt("Incorrect number of input var_order vars.");
- vector<int> *var_order_vp = (new vector<int>(nEndo));
+ vector<int> var_order_vp(nEndo);
for (int v = 0; v < nEndo; v++)
- (*var_order_vp)[v] = (int)(*(dparams++));
+ var_order_vp[v] = (int)(*(dparams++));
// the lag, current and lead blocks of the jacobian respectively
mxFldp = mxGetField(M_, 0, "lead_lag_incidence");
@@ -182,14 +182,14 @@ extern "C" {
npar = (int) mxGetN(mxFldp);
int nrows = (int) mxGetM(mxFldp);
- TwoDMatrix *llincidence = new TwoDMatrix(nrows, npar, dparams);
+ TwoDMatrix llincidence(nrows, npar, dparams);
if (npar != nEndo)
mexErrMsgIdAndTxt("dynare:k_order_perturbation", "Incorrect length of lead lag incidences: ncol=%d != nEndo=%d.", npar, nEndo);
//get NNZH =NNZD(2) = the total number of non-zero Hessian elements
mxFldp = mxGetField(M_, 0, "NNZDerivatives");
dparams = (double *) mxGetData(mxFldp);
- Vector *NNZD = new Vector(dparams, (int) mxGetM(mxFldp));
+ Vector NNZD(dparams, (int) mxGetM(mxFldp));
const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
@@ -258,7 +258,7 @@ extern "C" {
app.getFoldDecisionRule().writeMMap(mm, string());
// get latest ysteady
- ySteady = (Vector *)(&dynare.getSteady());
+ ySteady = dynare.getSteady();
if (kOrder == 1)
{