separating dynamic_dll handling into new dynamic_dll.cpp from k_order_perturbation.cpp containing now only mex_function for k_order_perturbationDLL and renaming k_order_perturbation.h to dynamic_dll.h
git-svn-id: https://www.dynare.org/svn/dynare/trunk@2864 ac1d8469-bf42-47a9-8791-bf33cf982152time-shift
parent
2006e1f9be
commit
6fc142406f
|
@ -0,0 +1,244 @@
|
||||||
|
/*
|
||||||
|
* Copyright (C) 2008-2009 Dynare Team
|
||||||
|
*
|
||||||
|
* This file is part of Dynare.
|
||||||
|
*
|
||||||
|
* Dynare is free software: you can redistribute it and/or modify
|
||||||
|
* it under the terms of the GNU General Public License as published by
|
||||||
|
* the Free Software Foundation, either version 3 of the License, or
|
||||||
|
* (at your option) any later version.
|
||||||
|
*
|
||||||
|
* Dynare is distributed in the hope that it will be useful,
|
||||||
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
* GNU General Public License for more details.
|
||||||
|
*
|
||||||
|
* You should have received a copy of the GNU General Public License
|
||||||
|
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
*/
|
||||||
|
|
||||||
|
|
||||||
|
#include "k_ord_dynare.h"
|
||||||
|
#include "dynamic_dll.h"
|
||||||
|
#include "math.h"
|
||||||
|
#include <cstring>
|
||||||
|
|
||||||
|
//////////////////////////////////////////////////////
|
||||||
|
// Convert Matlab Dynare endo and exo names array to C type array of string pointers
|
||||||
|
// Poblem is that Matlab mx function returns a long string concatenated by columns rather than rows
|
||||||
|
// hence a rather low level approach is needed
|
||||||
|
///////////////////////////////////////////////////////
|
||||||
|
const char **
|
||||||
|
DynareMxArrayToString(const mxArray *mxFldp, const int len, const int width)
|
||||||
|
{
|
||||||
|
char *cNamesCharStr = mxArrayToString(mxFldp);
|
||||||
|
const char **ret = DynareMxArrayToString(cNamesCharStr, len, width);
|
||||||
|
return ret;
|
||||||
|
}
|
||||||
|
|
||||||
|
const char **
|
||||||
|
DynareMxArrayToString(const char *cNamesCharStr, const int len, const int width)
|
||||||
|
{
|
||||||
|
char cNamesMX[len][width+1]; //
|
||||||
|
#ifdef DEBUG
|
||||||
|
mexPrintf("loop DynareMxArrayToString cNamesCharStr = %s \n", cNamesCharStr);
|
||||||
|
#endif
|
||||||
|
for (int i = 0; i < width; i++)
|
||||||
|
{
|
||||||
|
for (int j = 0; j < len; j++)
|
||||||
|
{
|
||||||
|
// Allow alphanumeric and underscores "_" only:
|
||||||
|
if (isalnum(cNamesCharStr[j+i*len]) || ('_' == cNamesCharStr[j+i*len]))
|
||||||
|
{
|
||||||
|
cNamesMX[j][i] = cNamesCharStr[j+i*len];
|
||||||
|
}
|
||||||
|
else cNamesMX[j][i] = '\0';
|
||||||
|
}
|
||||||
|
}
|
||||||
|
const char **ret = (const char **) mxCalloc(len, sizeof(char *));
|
||||||
|
for (int j = 0; j < len; j++)
|
||||||
|
{
|
||||||
|
cNamesMX[j][width] = '\0';
|
||||||
|
#ifdef DEBUG
|
||||||
|
// mexPrintf("String [%d]= %s \n", j, cNamesMX[j]);
|
||||||
|
#endif
|
||||||
|
char *token = (char *) mxCalloc(strlen(cNamesMX[j])+1, sizeof(char));
|
||||||
|
strcpy(token, cNamesMX[j]);
|
||||||
|
ret[j] = token;
|
||||||
|
#ifdef DEBUG
|
||||||
|
mexPrintf("ret [%d]= %s \n", j, ret[j]);
|
||||||
|
#endif
|
||||||
|
}
|
||||||
|
return ret;
|
||||||
|
}
|
||||||
|
|
||||||
|
/***********************************
|
||||||
|
* Members of DynamicModelDLL for handling loading and calling
|
||||||
|
* <model>_dynamic () function
|
||||||
|
**************************************/
|
||||||
|
DynamicModelDLL::DynamicModelDLL(const char *modName, const int y_length, const int j_cols,
|
||||||
|
const int n_max_lag, const int n_exog, const char *sExt)
|
||||||
|
: length(y_length), jcols(j_cols), nMax_lag(n_max_lag), nExog(n_exog)
|
||||||
|
{
|
||||||
|
char fName[MAX_MODEL_NAME];
|
||||||
|
strcpy(fName, modName);
|
||||||
|
using namespace std;
|
||||||
|
strcat(fName, "_dynamic");
|
||||||
|
#ifdef DEBUG
|
||||||
|
mexPrintf("MexPrintf: Call Load run DLL %s .\n", fName);
|
||||||
|
#endif
|
||||||
|
|
||||||
|
try
|
||||||
|
{
|
||||||
|
#ifdef WINDOWS
|
||||||
|
if (sExt == NULL) sExt = (".dll");
|
||||||
|
HINSTANCE dynamicHinstance;
|
||||||
|
// dynamicHinstance=::LoadLibraryEx(strcat(fNname,"_.dll"),NULL,DONT_RESOLVE_DLL_REFERENCES);//sExt); //"_.dll");
|
||||||
|
dynamicHinstance = ::LoadLibrary(strcat(fName, sExt)); //.dll); //"_.dll");
|
||||||
|
if (dynamicHinstance == NULL)
|
||||||
|
throw 1; //alt: return;
|
||||||
|
// (DynamicFn*) typedef void * (__stdcall *DynamicFn)();
|
||||||
|
# ifdef DEBUG
|
||||||
|
mexPrintf("MexPrintf: Call GetProcAddress %s .\n", fName);
|
||||||
|
# endif
|
||||||
|
Dynamic = (DynamicFn *) ::GetProcAddress(dynamicHinstance, "Dynamic");
|
||||||
|
|
||||||
|
#else // __linux__
|
||||||
|
if (sExt == NULL) sExt = (".so");
|
||||||
|
dynamicHinstance = dlopen(strcat(fName, sExt), RTLD_NOW);
|
||||||
|
if ((dynamicHinstance == NULL) || dlerror())
|
||||||
|
{
|
||||||
|
cerr << dlerror() << endl;
|
||||||
|
mexPrintf("MexPrintf:Error loading DLL: %s", dlerror);
|
||||||
|
throw 1;
|
||||||
|
}
|
||||||
|
Dynamic = (DynamicFn) dlsym(dynamicHinstance, "Dynamic");
|
||||||
|
if ((Dynamic == NULL) || dlerror())
|
||||||
|
{
|
||||||
|
cerr << dlerror() << endl;
|
||||||
|
mexPrintf("MexPrintf:Error finding DLL function: %s", dlerror);
|
||||||
|
throw 2;
|
||||||
|
}
|
||||||
|
#endif
|
||||||
|
|
||||||
|
}
|
||||||
|
catch (int i)
|
||||||
|
{
|
||||||
|
mexPrintf("MexPrintf: error in Load and run DLL %s , %d.\n", fName, i);
|
||||||
|
mexErrMsgTxt("Err: An error in Load and run DLL .\n");
|
||||||
|
return;
|
||||||
|
|
||||||
|
}
|
||||||
|
catch (...)
|
||||||
|
{
|
||||||
|
mexPrintf("MexPrintf: Unknown error in Call MATLAB function %s.\n", fName);
|
||||||
|
mexErrMsgTxt("Err: Unknown error in Load and run DLL .\n");
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// close DLL: If the referenced object was successfully closed,
|
||||||
|
// close() returns 0, non 0 otherwise
|
||||||
|
int
|
||||||
|
DynamicModelDLL::close()
|
||||||
|
{
|
||||||
|
#ifdef WINDOWS
|
||||||
|
// MS FreeLibrary returns non 0 if OK, 0 if fails.
|
||||||
|
bool rb = FreeLibrary(dynamicHinstance);
|
||||||
|
if (rb)
|
||||||
|
return 0;
|
||||||
|
else
|
||||||
|
return 1;
|
||||||
|
#else // linux
|
||||||
|
//If OK, dlclose() returns 0, non 0 otherwise
|
||||||
|
return dlclose(dynamicHinstance);
|
||||||
|
#endif
|
||||||
|
};
|
||||||
|
|
||||||
|
void
|
||||||
|
DynamicModelDLL::eval(const Vector &y, const TwoDMatrix &x, const Vector *modParams,
|
||||||
|
int it_, Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2)
|
||||||
|
{
|
||||||
|
|
||||||
|
double *dresidual, *dg1 = NULL, *dg2 = NULL;
|
||||||
|
//int length=y.length(); // not!
|
||||||
|
if ((jcols-nExog) != y.length())
|
||||||
|
{
|
||||||
|
// throw DLL Error
|
||||||
|
mexPrintf(" DLL Error: (jcols-nExog)!=ys.length() \n");
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
if (residual.length() < length) // dummy or insufficient
|
||||||
|
{
|
||||||
|
Vector *tempv = new Vector(length);
|
||||||
|
residual = *tempv;
|
||||||
|
delete tempv;
|
||||||
|
residual.zeros();
|
||||||
|
}
|
||||||
|
if (g1 != NULL)
|
||||||
|
{
|
||||||
|
if (g1->nrows() != length) // dummy
|
||||||
|
{
|
||||||
|
delete g1;
|
||||||
|
g1 = new TwoDMatrix(length, jcols); // and get a new one
|
||||||
|
g1->zeros();
|
||||||
|
}
|
||||||
|
dg1 = const_cast<double *>(g1->base());
|
||||||
|
}
|
||||||
|
if (g2 != NULL)
|
||||||
|
{
|
||||||
|
dg2 = const_cast<double *>(g2->base());
|
||||||
|
}
|
||||||
|
dresidual = const_cast<double *>(residual.base());
|
||||||
|
double *dy = const_cast<double *>(y.base());
|
||||||
|
double *dx = const_cast<double *>(x.base());
|
||||||
|
double *dbParams = const_cast<double *>(modParams->base());
|
||||||
|
#ifdef DEBUG
|
||||||
|
mexPrintf(" try eval Dynamic with ne g1: cols=%d , rows=%d\n",
|
||||||
|
g1->ncols(), g1->nrows());
|
||||||
|
for (int i = 0; i < modParams->length(); i++)
|
||||||
|
{
|
||||||
|
mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]);
|
||||||
|
}
|
||||||
|
for (int i = 0; i < jcols-nExog; i++)
|
||||||
|
{
|
||||||
|
mexPrintf("k_ord_perturbation: Ys[%d]= %g.\n", i, dy[i]);
|
||||||
|
}
|
||||||
|
mexPrintf("k_order_perturbation: call <model> Dynamic dParams= %g , , dy = %g dx = %f .\n",
|
||||||
|
dbParams[0], dy[0], dx[0]);
|
||||||
|
|
||||||
|
#endif
|
||||||
|
try
|
||||||
|
{
|
||||||
|
Dynamic(dy, dx, nExog, dbParams, it_, dresidual, dg1, dg2);
|
||||||
|
}
|
||||||
|
catch (...)
|
||||||
|
{
|
||||||
|
mexPrintf("MexPrintf: error in run Dynamic DLL \n");
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
void
|
||||||
|
DynamicModelDLL::eval(const Vector &y, const TwoDMatrix &x, const Vector *modParams,
|
||||||
|
Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2)
|
||||||
|
{
|
||||||
|
|
||||||
|
eval(y, x, modParams, nMax_lag, residual, g1, g2);
|
||||||
|
};
|
||||||
|
|
||||||
|
void
|
||||||
|
DynamicModelDLL::eval(const Vector &y, const Vector &x, const Vector *modParams,
|
||||||
|
Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2)
|
||||||
|
{
|
||||||
|
|
||||||
|
/** ignore given exogens and create new 2D x matrix since
|
||||||
|
* when calling <model>_dynamic(z,x,params,it_) x must be equal to
|
||||||
|
* zeros(M_.maximum_lag+1,M_.exo_nbr)
|
||||||
|
**/
|
||||||
|
TwoDMatrix &mx = *(new TwoDMatrix(nMax_lag+1, nExog));
|
||||||
|
mx.zeros(); // initialise shocks to 0s
|
||||||
|
|
||||||
|
eval(y, mx, modParams, nMax_lag, residual, g1, g2);
|
||||||
|
};
|
||||||
|
|
||||||
|
|
|
@ -22,6 +22,7 @@
|
||||||
#include <vector>
|
#include <vector>
|
||||||
#include "first_order.h"
|
#include "first_order.h"
|
||||||
#include "k_ord_dynare.h"
|
#include "k_ord_dynare.h"
|
||||||
|
#include "dynamic_dll.h"
|
||||||
|
|
||||||
#include "mex.h"
|
#include "mex.h"
|
||||||
|
|
||||||
|
@ -40,12 +41,12 @@ class KordpJacobian;
|
||||||
KordpDynare::KordpDynare(const char **endo, int num_endo,
|
KordpDynare::KordpDynare(const char **endo, int num_endo,
|
||||||
const char **exo, int nexog, int npar, //const char** par,
|
const char **exo, int nexog, int npar, //const char** par,
|
||||||
Vector *ysteady, TwoDMatrix *vcov, Vector *inParams, int nstat,
|
Vector *ysteady, TwoDMatrix *vcov, Vector *inParams, int nstat,
|
||||||
int npred, int nforw, int nboth, const int jcols, const Vector *nnzd,
|
int npred, int nforw, int nboth, const int jcols, const Vector *nnzd,
|
||||||
const int nsteps, int norder, //const char* modName,
|
const int nsteps, int norder, //const char* modName,
|
||||||
Journal &jr, DynamicModelDLL &dynamicDLL, double sstol,
|
Journal &jr, DynamicModelDLL &dynamicDLL, double sstol,
|
||||||
const vector<int> *var_order, const TwoDMatrix *llincidence, double criterium)
|
const vector<int> *var_order, const TwoDMatrix *llincidence, double criterium)
|
||||||
: nStat(nstat), nBoth(nboth), nPred(npred), nForw(nforw), nExog(nexog), nPar(npar),
|
: 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),
|
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),
|
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(NULL), denl(NULL), dsnl(NULL), ss_tol(sstol), varOrder(var_order),
|
||||||
ll_Incidence(llincidence), qz_criterium(criterium)
|
ll_Incidence(llincidence), qz_criterium(criterium)
|
||||||
|
@ -94,8 +95,8 @@ KordpDynare::KordpDynare(const char **endo, int num_endo,
|
||||||
KordpDynare::KordpDynare(const KordpDynare &dynare)
|
KordpDynare::KordpDynare(const KordpDynare &dynare)
|
||||||
: nStat(dynare.nStat), nBoth(dynare.nBoth), nPred(dynare.nPred),
|
: nStat(dynare.nStat), nBoth(dynare.nBoth), nPred(dynare.nPred),
|
||||||
nForw(dynare.nForw), nExog(dynare.nExog), nPar(dynare.nPar),
|
nForw(dynare.nForw), nExog(dynare.nExog), nPar(dynare.nPar),
|
||||||
nYs(dynare.nYs), nYss(dynare.nYss), nY(dynare.nY), nJcols(dynare.nJcols),
|
nYs(dynare.nYs), nYss(dynare.nYss), nY(dynare.nY), nJcols(dynare.nJcols),
|
||||||
NNZD(dynare.NNZD), nSteps(dynare.nSteps), nOrder(dynare.nOrder),
|
NNZD(dynare.NNZD), nSteps(dynare.nSteps), nOrder(dynare.nOrder),
|
||||||
journal(dynare.journal), dynamicDLL(dynare.dynamicDLL), //modName(dynare.modName),
|
journal(dynare.journal), dynamicDLL(dynare.dynamicDLL), //modName(dynare.modName),
|
||||||
ySteady(NULL), params(NULL), vCov(NULL), md(dynare.md),
|
ySteady(NULL), params(NULL), vCov(NULL), md(dynare.md),
|
||||||
dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(dynare.ss_tol),
|
dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(dynare.ss_tol),
|
||||||
|
@ -203,13 +204,13 @@ KordpDynare::calcDerivatives(const Vector &yy, const Vector &xx)
|
||||||
// Hessian TwoDMatrix *g2;
|
// Hessian TwoDMatrix *g2;
|
||||||
if (nOrder > 1)
|
if (nOrder > 1)
|
||||||
{
|
{
|
||||||
//g2 = new TwoDMatrix(nY, nJcols*nJcols); // generate g2 for Hessian
|
//g2 = new TwoDMatrix(nY, nJcols*nJcols); // generate g2 for Hessian
|
||||||
// generate g2 space for sparse Hessian 3x NNZH = 3x NNZD[1]
|
// generate g2 space for sparse Hessian 3x NNZH = 3x NNZD[1]
|
||||||
g2 = new TwoDMatrix((int) (*NNZD)[1],3);
|
g2 = new TwoDMatrix((int) (*NNZD)[1],3);
|
||||||
g2->zeros();
|
g2->zeros();
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
mexPrintf(" g2 cols %d rows %d \n", g2->numCols(), g2->numRows());
|
mexPrintf(" g2 cols %d rows %d \n", g2->numCols(), g2->numRows());
|
||||||
//g2->print();
|
//g2->print();
|
||||||
#endif
|
#endif
|
||||||
}
|
}
|
||||||
Vector &out = *(new Vector(nY));
|
Vector &out = *(new Vector(nY));
|
||||||
|
@ -249,12 +250,12 @@ KordpDynare::calcDerivatives(const Vector &yy, const Vector &xx)
|
||||||
populateDerivativesContainer(g1, 1, JacobianIndices);
|
populateDerivativesContainer(g1, 1, JacobianIndices);
|
||||||
if (nOrder > 1)
|
if (nOrder > 1)
|
||||||
{
|
{
|
||||||
// ReorderBlocks(g2,JacobianIndices);
|
// ReorderBlocks(g2,JacobianIndices);
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
mexPrintf(" post dll g2 cols %d rows %d \n", g2->numCols(), g2->numRows());
|
mexPrintf(" post dll g2 cols %d rows %d \n", g2->numCols(), g2->numRows());
|
||||||
for (int ii=0;ii<g2->numRows(); ii++)
|
for (int ii=0;ii<g2->numRows(); ii++)
|
||||||
mexPrintf(" g2[%d]: %d %d %f \n", ii, (int)g2->get(ii,0),(int)g2->get(ii,1),g2->get(ii,2));
|
mexPrintf(" g2[%d]: %d %d %f \n", ii, (int)g2->get(ii,0),(int)g2->get(ii,1),g2->get(ii,2));
|
||||||
//g2->print();
|
//g2->print();
|
||||||
#endif
|
#endif
|
||||||
populateDerivativesContainer(g2, 2, JacobianIndices);
|
populateDerivativesContainer(g2, 2, JacobianIndices);
|
||||||
}
|
}
|
||||||
|
@ -336,7 +337,7 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector<i
|
||||||
}
|
}
|
||||||
else if (ord == 2)
|
else if (ord == 2)
|
||||||
{
|
{
|
||||||
int nJcols1 = nJcols-nExog;
|
int nJcols1 = nJcols-nExog;
|
||||||
vector<int> revOrder(nJcols1);
|
vector<int> revOrder(nJcols1);
|
||||||
for (int i = 0; i < nJcols1; i++)
|
for (int i = 0; i < nJcols1; i++)
|
||||||
revOrder[(*vOrder)[i]] = i;
|
revOrder[(*vOrder)[i]] = i;
|
||||||
|
@ -344,7 +345,7 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector<i
|
||||||
{
|
{
|
||||||
int j = (int)g->get(i,0)-1; // hessian indices start with 1
|
int j = (int)g->get(i,0)-1; // hessian indices start with 1
|
||||||
int i1 = (int)g->get(i,1) -1;
|
int i1 = (int)g->get(i,1) -1;
|
||||||
int s0 = (int)floor(i1/nJcols);
|
int s0 = (int)floor(i1/nJcols);
|
||||||
int s1 = i1- (nJcols*s0);
|
int s1 = i1- (nJcols*s0);
|
||||||
if (s0 < nJcols1)
|
if (s0 < nJcols1)
|
||||||
s[0] = revOrder[s0];
|
s[0] = revOrder[s0];
|
||||||
|
@ -355,20 +356,20 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector<i
|
||||||
else
|
else
|
||||||
s[1] = s1;
|
s[1] = s1;
|
||||||
if (s[1] >= s[0])
|
if (s[1] >= s[0])
|
||||||
{
|
{
|
||||||
double x = g->get(i,2);
|
double x = g->get(i,2);
|
||||||
mdTi->insert(s, j, x);
|
mdTi->insert(s, j, x);
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
mexPrintf(" s[0]=%d s[1]=%d j=%d x=%f \n", s[0],s[1],j,x);
|
mexPrintf(" s[0]=%d s[1]=%d j=%d x=%f \n", s[0],s[1],j,x);
|
||||||
s.print();
|
s.print();
|
||||||
std::cout << s[0] << " " << s[1] << "\n";
|
std::cout << s[0] << " " << s[1] << "\n";
|
||||||
#endif
|
#endif
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
mexPrintf("k_ord_dynare.cpp: END populate FSSparseTensor in calcDerivatives: ord=%d \n",ord);
|
mexPrintf("k_ord_dynare.cpp: END populate FSSparseTensor in calcDerivatives: ord=%d \n",ord);
|
||||||
mdTi->print();
|
mdTi->print();
|
||||||
#endif
|
#endif
|
||||||
// md container
|
// md container
|
||||||
//md.clear();// this is to be used only for 1st order!!
|
//md.clear();// this is to be used only for 1st order!!
|
||||||
|
|
|
@ -37,7 +37,6 @@
|
||||||
#include "kord_exception.h"
|
#include "kord_exception.h"
|
||||||
#include "nlsolve.h"
|
#include "nlsolve.h"
|
||||||
#include "approximation.h"
|
#include "approximation.h"
|
||||||
#include "k_order_perturbation.h"
|
|
||||||
|
|
||||||
class KordpDynare;
|
class KordpDynare;
|
||||||
|
|
||||||
|
@ -123,8 +122,8 @@ class KordpDynare : public DynamicModel
|
||||||
const int nYs; // ={npred + nboth ; }
|
const int nYs; // ={npred + nboth ; }
|
||||||
const int nYss; // nyss ={ nboth + nforw ; }
|
const int nYss; // nyss ={ nboth + nforw ; }
|
||||||
const int nY; // = num_endo={ nstat + npred + nboth + nforw ; }
|
const int nY; // = num_endo={ nstat + npred + nboth + nforw ; }
|
||||||
const int nJcols; // no of jacobian columns= nExog+nEndo+nsPred+nsForw
|
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)
|
// where hessian is 2nd : NZZD(order=2)
|
||||||
const int nSteps;
|
const int nSteps;
|
||||||
const int nOrder;
|
const int nOrder;
|
||||||
|
@ -147,10 +146,10 @@ public:
|
||||||
KordpDynare(const char **endo, int num_endo,
|
KordpDynare(const char **endo, int num_endo,
|
||||||
const char **exo, int num_exo, int num_par, //const char** par,
|
const char **exo, int num_exo, int num_par, //const char** par,
|
||||||
Vector *ySteady, TwoDMatrix *vCov, Vector *params, int nstat, int nPred,
|
Vector *ySteady, TwoDMatrix *vCov, Vector *params, int nstat, int nPred,
|
||||||
int nforw, int nboth, const int nJcols, const Vector *NNZD,
|
int nforw, int nboth, const int nJcols, const Vector *NNZD,
|
||||||
const int nSteps, const int ord, //const char* modName,
|
const int nSteps, const int ord, //const char* modName,
|
||||||
Journal &jr, DynamicModelDLL &dynamicDLL, double sstol,
|
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);
|
double qz_criterium);
|
||||||
|
|
||||||
/** Makes a deep copy of the object. */
|
/** Makes a deep copy of the object. */
|
||||||
|
|
|
@ -39,6 +39,7 @@
|
||||||
**********************************************************/
|
**********************************************************/
|
||||||
|
|
||||||
#include "k_ord_dynare.h"
|
#include "k_ord_dynare.h"
|
||||||
|
#include "dynamic_dll.h"
|
||||||
#include "math.h"
|
#include "math.h"
|
||||||
#include <cstring>
|
#include <cstring>
|
||||||
|
|
||||||
|
@ -82,6 +83,7 @@ CK_order_perturbation::CK_order_perturbation()
|
||||||
|
|
||||||
#endif // _MSC_VER && WINDOWS
|
#endif // _MSC_VER && WINDOWS
|
||||||
|
|
||||||
|
#ifdef MATLAB_MEX_FILE // exclude mexFunction for other applications
|
||||||
extern "C" {
|
extern "C" {
|
||||||
|
|
||||||
// mexFunction: Matlab Inerface point and the main application driver
|
// mexFunction: Matlab Inerface point and the main application driver
|
||||||
|
@ -176,8 +178,8 @@ extern "C" {
|
||||||
const int nPar = (int) mxGetScalar(mxFldp);
|
const int nPar = (int) mxGetScalar(mxFldp);
|
||||||
// it_ should be set to M_.maximum_lag
|
// it_ should be set to M_.maximum_lag
|
||||||
mxFldp = mxGetField(M_, 0, "maximum_lag");
|
mxFldp = mxGetField(M_, 0, "maximum_lag");
|
||||||
const int nMax_lag = (int) mxGetScalar(mxFldp);
|
const int nMax_lag = (int) mxGetScalar(mxFldp);
|
||||||
|
|
||||||
nPred -= nBoth; // correct nPred for nBoth.
|
nPred -= nBoth; // correct nPred for nBoth.
|
||||||
|
|
||||||
mxFldp = mxGetField(dr, 0, "order_var");
|
mxFldp = mxGetField(dr, 0, "order_var");
|
||||||
|
@ -221,15 +223,15 @@ extern "C" {
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
//get NNZH =NNZD(2) = the total number of non-zero Hessian elements
|
//get NNZH =NNZD(2) = the total number of non-zero Hessian elements
|
||||||
mxFldp = mxGetField(M_, 0, "NNZDerivatives");
|
mxFldp = mxGetField(M_, 0, "NNZDerivatives");
|
||||||
dparams = (double *) mxGetData(mxFldp);
|
dparams = (double *) mxGetData(mxFldp);
|
||||||
Vector *NNZD = new Vector (dparams, (int) mxGetM(mxFldp));
|
Vector *NNZD = new Vector (dparams, (int) mxGetM(mxFldp));
|
||||||
#ifdef DEBUG
|
#ifdef DEBUG
|
||||||
mexPrintf("NNZH=%d, \n", (int) (*NNZD)[1]);
|
mexPrintf("NNZH=%d, \n", (int) (*NNZD)[1]);
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
|
const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
|
||||||
mexPrintf("k_order_perturbation: jcols= %d .\n", jcols);
|
mexPrintf("k_order_perturbation: jcols= %d .\n", jcols);
|
||||||
|
|
||||||
|
@ -308,7 +310,7 @@ extern "C" {
|
||||||
// make KordpDynare object
|
// make KordpDynare object
|
||||||
KordpDynare dynare(endoNamesMX, nEndo, exoNamesMX, nExog, nPar, // paramNames,
|
KordpDynare dynare(endoNamesMX, nEndo, exoNamesMX, nExog, nPar, // paramNames,
|
||||||
ySteady, vCov, modParams, nStat, nPred, nForw, nBoth,
|
ySteady, vCov, modParams, nStat, nPred, nForw, nBoth,
|
||||||
jcols, NNZD, nSteps, kOrder, journal, dynamicDLL,
|
jcols, NNZD, nSteps, kOrder, journal, dynamicDLL,
|
||||||
sstol, var_order_vp, llincidence, qz_criterium);
|
sstol, var_order_vp, llincidence, qz_criterium);
|
||||||
|
|
||||||
// construct main K-order approximation class
|
// construct main K-order approximation class
|
||||||
|
@ -469,222 +471,10 @@ extern "C" {
|
||||||
} //catch
|
} //catch
|
||||||
}; // end of mexFunction()
|
}; // end of mexFunction()
|
||||||
}; // end of extern C
|
}; // end of extern C
|
||||||
|
#endif // ifdef MATLAB_MEX_FILE to exclude mexFunction for other applications
|
||||||
|
|
||||||
//////////////////////////////////////////////////////
|
|
||||||
// Convert Matlab Dynare endo and exo names array to C type array of string pointers
|
|
||||||
// Poblem is that Matlab mx function returns a long string concatenated by columns rather than rows
|
|
||||||
// hence a rather low level approach is needed
|
|
||||||
///////////////////////////////////////////////////////
|
|
||||||
const char **
|
|
||||||
DynareMxArrayToString(const mxArray *mxFldp, const int len, const int width)
|
|
||||||
{
|
|
||||||
char *cNamesCharStr = mxArrayToString(mxFldp);
|
|
||||||
const char **ret = DynareMxArrayToString(cNamesCharStr, len, width);
|
|
||||||
return ret;
|
|
||||||
}
|
|
||||||
|
|
||||||
const char **
|
|
||||||
DynareMxArrayToString(const char *cNamesCharStr, const int len, const int width)
|
|
||||||
{
|
|
||||||
char cNamesMX[len][width+1]; //
|
|
||||||
#ifdef DEBUG
|
|
||||||
mexPrintf("loop DynareMxArrayToString cNamesCharStr = %s \n", cNamesCharStr);
|
|
||||||
#endif
|
|
||||||
for (int i = 0; i < width; i++)
|
|
||||||
{
|
|
||||||
for (int j = 0; j < len; j++)
|
|
||||||
{
|
|
||||||
// Allow alphanumeric and underscores "_" only:
|
|
||||||
if (isalnum(cNamesCharStr[j+i*len]) || ('_' == cNamesCharStr[j+i*len]))
|
|
||||||
{
|
|
||||||
cNamesMX[j][i] = cNamesCharStr[j+i*len];
|
|
||||||
}
|
|
||||||
else cNamesMX[j][i] = '\0';
|
|
||||||
}
|
|
||||||
}
|
|
||||||
const char **ret = (const char **) mxCalloc(len, sizeof(char *));
|
|
||||||
for (int j = 0; j < len; j++)
|
|
||||||
{
|
|
||||||
cNamesMX[j][width] = '\0';
|
|
||||||
#ifdef DEBUG
|
|
||||||
// mexPrintf("String [%d]= %s \n", j, cNamesMX[j]);
|
|
||||||
#endif
|
|
||||||
char *token = (char *) mxCalloc(strlen(cNamesMX[j])+1, sizeof(char));
|
|
||||||
strcpy(token, cNamesMX[j]);
|
|
||||||
ret[j] = token;
|
|
||||||
#ifdef DEBUG
|
|
||||||
mexPrintf("ret [%d]= %s \n", j, ret[j]);
|
|
||||||
#endif
|
|
||||||
}
|
|
||||||
return ret;
|
|
||||||
}
|
|
||||||
|
|
||||||
/***********************************
|
|
||||||
* Members of DynamicModelDLL for handling loading and calling
|
|
||||||
* <model>_dynamic () function
|
|
||||||
**************************************/
|
|
||||||
DynamicModelDLL::DynamicModelDLL(const char *modName, const int y_length, const int j_cols,
|
|
||||||
const int n_max_lag, const int n_exog, const char *sExt)
|
|
||||||
: length(y_length), jcols(j_cols), nMax_lag(n_max_lag), nExog(n_exog)
|
|
||||||
{
|
|
||||||
char fName[MAX_MODEL_NAME];
|
|
||||||
strcpy(fName, modName);
|
|
||||||
using namespace std;
|
|
||||||
strcat(fName, "_dynamic");
|
|
||||||
#ifdef DEBUG
|
|
||||||
mexPrintf("MexPrintf: Call Load run DLL %s .\n", fName);
|
|
||||||
#endif
|
|
||||||
|
|
||||||
try
|
|
||||||
{
|
|
||||||
#ifdef WINDOWS
|
|
||||||
if (sExt == NULL) sExt = (".dll");
|
|
||||||
HINSTANCE dynamicHinstance;
|
|
||||||
// dynamicHinstance=::LoadLibraryEx(strcat(fNname,"_.dll"),NULL,DONT_RESOLVE_DLL_REFERENCES);//sExt); //"_.dll");
|
|
||||||
dynamicHinstance = ::LoadLibrary(strcat(fName, sExt)); //.dll); //"_.dll");
|
|
||||||
if (dynamicHinstance == NULL)
|
|
||||||
throw 1; //alt: return;
|
|
||||||
// (DynamicFn*) typedef void * (__stdcall *DynamicFn)();
|
|
||||||
# ifdef DEBUG
|
|
||||||
mexPrintf("MexPrintf: Call GetProcAddress %s .\n", fName);
|
|
||||||
# endif
|
|
||||||
Dynamic = (DynamicFn *) ::GetProcAddress(dynamicHinstance, "Dynamic");
|
|
||||||
|
|
||||||
#else // __linux__
|
|
||||||
if (sExt == NULL) sExt = (".so");
|
|
||||||
dynamicHinstance = dlopen(strcat(fName, sExt), RTLD_NOW);
|
|
||||||
if ((dynamicHinstance == NULL) || dlerror())
|
|
||||||
{
|
|
||||||
cerr << dlerror() << endl;
|
|
||||||
mexPrintf("MexPrintf:Error loading DLL: %s", dlerror);
|
|
||||||
throw 1;
|
|
||||||
}
|
|
||||||
Dynamic = (DynamicFn) dlsym(dynamicHinstance, "Dynamic");
|
|
||||||
if ((Dynamic == NULL) || dlerror())
|
|
||||||
{
|
|
||||||
cerr << dlerror() << endl;
|
|
||||||
mexPrintf("MexPrintf:Error finding DLL function: %s", dlerror);
|
|
||||||
throw 2;
|
|
||||||
}
|
|
||||||
#endif
|
|
||||||
|
|
||||||
}
|
|
||||||
catch (int i)
|
|
||||||
{
|
|
||||||
mexPrintf("MexPrintf: error in Load and run DLL %s , %d.\n", fName, i);
|
|
||||||
mexErrMsgTxt("Err: An error in Load and run DLL .\n");
|
|
||||||
return;
|
|
||||||
|
|
||||||
}
|
|
||||||
catch (...)
|
|
||||||
{
|
|
||||||
mexPrintf("MexPrintf: Unknown error in Call MATLAB function %s.\n", fName);
|
|
||||||
mexErrMsgTxt("Err: Unknown error in Load and run DLL .\n");
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// close DLL: If the referenced object was successfully closed,
|
|
||||||
// close() returns 0, non 0 otherwise
|
|
||||||
int
|
|
||||||
DynamicModelDLL::close()
|
|
||||||
{
|
|
||||||
#ifdef WINDOWS
|
|
||||||
// MS FreeLibrary returns non 0 if OK, 0 if fails.
|
|
||||||
bool rb = FreeLibrary(dynamicHinstance);
|
|
||||||
if (rb)
|
|
||||||
return 0;
|
|
||||||
else
|
|
||||||
return 1;
|
|
||||||
#else // linux
|
|
||||||
//If OK, dlclose() returns 0, non 0 otherwise
|
|
||||||
return dlclose(dynamicHinstance);
|
|
||||||
#endif
|
|
||||||
};
|
|
||||||
|
|
||||||
void
|
|
||||||
DynamicModelDLL::eval(const Vector &y, const TwoDMatrix &x, const Vector *modParams,
|
|
||||||
int it_, Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2)
|
|
||||||
{
|
|
||||||
|
|
||||||
double *dresidual, *dg1 = NULL, *dg2 = NULL;
|
|
||||||
//int length=y.length(); // not!
|
|
||||||
if ((jcols-nExog) != y.length())
|
|
||||||
{
|
|
||||||
// throw DLL Error
|
|
||||||
mexPrintf(" DLL Error: (jcols-nExog)!=ys.length() \n");
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
if (residual.length() < length) // dummy or insufficient
|
|
||||||
{
|
|
||||||
Vector *tempv = new Vector(length);
|
|
||||||
residual = *tempv;
|
|
||||||
delete tempv;
|
|
||||||
residual.zeros();
|
|
||||||
}
|
|
||||||
if (g1 != NULL)
|
|
||||||
{
|
|
||||||
if (g1->nrows() != length) // dummy
|
|
||||||
{
|
|
||||||
delete g1;
|
|
||||||
g1 = new TwoDMatrix(length, jcols); // and get a new one
|
|
||||||
g1->zeros();
|
|
||||||
}
|
|
||||||
dg1 = const_cast<double *>(g1->base());
|
|
||||||
}
|
|
||||||
if (g2 != NULL)
|
|
||||||
{
|
|
||||||
dg2 = const_cast<double *>(g2->base());
|
|
||||||
}
|
|
||||||
dresidual = const_cast<double *>(residual.base());
|
|
||||||
double *dy = const_cast<double *>(y.base());
|
|
||||||
double *dx = const_cast<double *>(x.base());
|
|
||||||
double *dbParams = const_cast<double *>(modParams->base());
|
|
||||||
#ifdef DEBUG
|
|
||||||
mexPrintf(" try eval Dynamic with ne g1: cols=%d , rows=%d\n",
|
|
||||||
g1->ncols(), g1->nrows());
|
|
||||||
for (int i = 0; i < modParams->length(); i++)
|
|
||||||
{
|
|
||||||
mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]);
|
|
||||||
}
|
|
||||||
for (int i = 0; i < jcols-nExog; i++)
|
|
||||||
{
|
|
||||||
mexPrintf("k_ord_perturbation: Ys[%d]= %g.\n", i, dy[i]);
|
|
||||||
}
|
|
||||||
mexPrintf("k_order_perturbation: call <model> Dynamic dParams= %g , , dy = %g dx = %f .\n",
|
|
||||||
dbParams[0], dy[0], dx[0]);
|
|
||||||
|
|
||||||
#endif
|
|
||||||
try
|
|
||||||
{
|
|
||||||
Dynamic(dy, dx, nExog, dbParams, it_, dresidual, dg1, dg2);
|
|
||||||
}
|
|
||||||
catch (...)
|
|
||||||
{
|
|
||||||
mexPrintf("MexPrintf: error in run Dynamic DLL \n");
|
|
||||||
}
|
|
||||||
};
|
|
||||||
|
|
||||||
void
|
|
||||||
DynamicModelDLL::eval(const Vector &y, const TwoDMatrix &x, const Vector *modParams,
|
|
||||||
Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2)
|
|
||||||
{
|
|
||||||
|
|
||||||
eval(y, x, modParams, nMax_lag, residual, g1, g2);
|
|
||||||
};
|
|
||||||
|
|
||||||
void
|
|
||||||
DynamicModelDLL::eval(const Vector &y, const Vector &x, const Vector *modParams,
|
|
||||||
Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2)
|
|
||||||
{
|
|
||||||
|
|
||||||
/** ignore given exogens and create new 2D x matrix since
|
|
||||||
* when calling <model>_dynamic(z,x,params,it_) x must be equal to
|
|
||||||
* zeros(M_.maximum_lag+1,M_.exo_nbr)
|
|
||||||
**/
|
|
||||||
TwoDMatrix &mx = *(new TwoDMatrix(nMax_lag+1, nExog));
|
|
||||||
mx.zeros(); // initialise shocks to 0s
|
|
||||||
|
|
||||||
eval(y, mx, modParams, nMax_lag, residual, g1, g2);
|
|
||||||
};
|
|
||||||
|
|
||||||
|
|
|
@ -26,6 +26,7 @@
|
||||||
|
|
||||||
//#include "stdafx.h"
|
//#include "stdafx.h"
|
||||||
#include "k_ord_dynare.h"
|
#include "k_ord_dynare.h"
|
||||||
|
#include "dynamic_dll.h"
|
||||||
|
|
||||||
int
|
int
|
||||||
main(int argc, char *argv[])
|
main(int argc, char *argv[])
|
||||||
|
@ -79,9 +80,9 @@ main(int argc, char *argv[])
|
||||||
};
|
};
|
||||||
const int nSteady = 16; //27 //31;//29, 16 (int)mxGetM(mxFldp);
|
const int nSteady = 16; //27 //31;//29, 16 (int)mxGetM(mxFldp);
|
||||||
Vector *ySteady = new Vector(dYSparams, nSteady);
|
Vector *ySteady = new Vector(dYSparams, nSteady);
|
||||||
|
|
||||||
double nnzd[3]={ 77,217,0};
|
double nnzd[3]={ 77,217,0};
|
||||||
const Vector *NNZD = new Vector(nnzd, 3);
|
const Vector *NNZD = new Vector(nnzd, 3);
|
||||||
|
|
||||||
//mxFldp = mxGetField(dr, 0,"nstatic" );
|
//mxFldp = mxGetField(dr, 0,"nstatic" );
|
||||||
const int nStat = 7; //(int)mxGetScalar(mxFldp);
|
const int nStat = 7; //(int)mxGetScalar(mxFldp);
|
||||||
|
|
|
@ -1 +1 @@
|
||||||
gcc -DMATLAB_MEX_FILE -mthreads -shared -DWINDOWS -DPOSIX_THREADS -I"c:/Program Files"/MATLAB_SV71/extern/include -I Dynare_pp/src -I Dynare_pp/kord -I Dynare_pp/tl/cc -I Dynare_pp/utils/cc/ -I Dynare_pp/sylv/cc/ -I"f:/Pthreads/Pre-built.2/include" -I"f:/mingw/include" -o k_order_perturbation.dll k_ord_dynare.o k_order_perturbation.o nlsolve.o "dynare_pp/extern/matlab"/dynarelib.a -Wl,-L"c:/Program Files"/MATLAB_SV71/extern/lib/win32/microsoft/ -Wl,-llibmex -Wl,-llibmx -Wl,-llibmwlapack -Wl,-llibdflapack -lg2c -lmingw32 -lstdc++ -L"f:/Pthreads/Pre-built.2/lib" -lpthreadGC2
|
gcc -DMATLAB_MEX_FILE -mthreads -shared -DWINDOWS -DPOSIX_THREADS -I"c:/Program Files"/MATLAB_SV71/extern/include -I Dynare_pp/src -I Dynare_pp/kord -I Dynare_pp/tl/cc -I Dynare_pp/utils/cc/ -I Dynare_pp/sylv/cc/ -I"f:/Pthreads/Pre-built.2/include" -I"f:/mingw/include" -o k_order_perturbation.dll k_ord_dynare.o k_order_perturbation.o dynamic_dll.o nlsolve.o "dynare_pp/extern/matlab"/dynarelib.a -Wl,-L"c:/Program Files"/MATLAB_SV71/extern/lib/win32/microsoft/ -Wl,-llibmex -Wl,-llibmx -Wl,-llibmwlapack -Wl,-llibdflapack -lg2c -lmingw32 -lstdc++ -L"f:/Pthreads/Pre-built.2/lib" -lpthreadGC2
|
||||||
|
|
|
@ -1 +1 @@
|
||||||
gcc -DMATLAB_MEX_FILE -mthreads -g -DWINDOWS -DPOSIX_THREADS -I"c:/Program Files"/MATLAB_SV71/extern/include -I Dynare_pp/src -I Dynare_pp/kord -I Dynare_pp/tl/cc -I Dynare_pp/utils/cc/ -I Dynare_pp/sylv/cc/ -I"f:/Pthreads/Pre-built.2/include" -I"f:/mingw/include" -o k_orddbgtest.exe k_order_test_main.o k_ord_dynare.o k_order_perturbation.o nlsolve.o "dynare_pp/extern/matlab"/dynarelib.a -Wl,-L"c:/Program Files"/MATLAB_SV71/extern/lib/win32/microsoft/ -Wl,-llibmex -Wl,-llibmx -Wl,-llibmwlapack -Wl,-llibdflapack -lg2c -lmingw32 -lstdc++ -L"f:/Pthreads/Pre-built.2/lib" -lpthreadGC2
|
gcc -DMATLAB_MEX_FILE -mthreads -g -DWINDOWS -DPOSIX_THREADS -I"c:/Program Files"/MATLAB_SV71/extern/include -I Dynare_pp/src -I Dynare_pp/kord -I Dynare_pp/tl/cc -I Dynare_pp/utils/cc/ -I Dynare_pp/sylv/cc/ -I"f:/Pthreads/Pre-built.2/include" -I"f:/mingw/include" -o k_orddbgtest.exe k_order_test_main.o k_ord_dynare.o k_order_perturbation.o nlsolve.o dynamic_dll.o "dynare_pp/extern/matlab"/dynarelib.a -Wl,-L"c:/Program Files"/MATLAB_SV71/extern/lib/win32/microsoft/ -Wl,-llibmex -Wl,-llibmx -Wl,-llibmwlapack -Wl,-llibdflapack -lg2c -lmingw32 -lstdc++ -L"f:/Pthreads/Pre-built.2/lib" -lpthreadGC2
|
||||||
|
|
Loading…
Reference in New Issue