From 5f05a14b175a56c9f00bb6c893ab8091ce82b572 Mon Sep 17 00:00:00 2001 From: michel Date: Sun, 29 Nov 2009 20:50:39 +0000 Subject: [PATCH] v4.1: fixing k_order_perturbation DLL for order == 3 git-svn-id: https://www.dynare.org/svn/dynare/trunk@3170 ac1d8469-bf42-47a9-8791-bf33cf982152 --- .../k_order_perturbation/dynamic_dll.cpp | 95 +++++----------- .../k_order_perturbation/dynamic_dll.h | 15 ++- .../k_order_perturbation/k_ord_dynare.cpp | 106 +++++++++++------- .../k_order_perturbation.cpp | 105 ++++++++++++++--- 4 files changed, 188 insertions(+), 133 deletions(-) diff --git a/mex/sources/k_order_perturbation/dynamic_dll.cpp b/mex/sources/k_order_perturbation/dynamic_dll.cpp index 977206d27..fc87df3b8 100644 --- a/mex/sources/k_order_perturbation/dynamic_dll.cpp +++ b/mex/sources/k_order_perturbation/dynamic_dll.cpp @@ -23,60 +23,6 @@ #include "math.h" #include -////////////////////////////////////////////////////// -// 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; - cNamesMX = (char **) mxCalloc(len, sizeof(char *)); - for(int i = 0; i < len; i++) - cNamesMX[i] = (char *) mxCalloc(width+1, sizeof(char)); - -#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 - } - mxFree(cNamesMX); - return ret; -} - /*********************************** * Members of DynamicModelDLL for handling loading and calling * _dynamic () function @@ -86,7 +32,12 @@ DynamicModelDLL::DynamicModelDLL(const char *modName, const int y_length, const : length(y_length), jcols(j_cols), nMax_lag(n_max_lag), nExog(n_exog) { char fName[MAX_MODEL_NAME]; +#if defined _WIN32 strcpy(fName, modName); +#else + strcpy(fName, "./"); + strcat(fName, modName); +#endif using namespace std; strcat(fName, "_dynamic"); #ifdef DEBUG @@ -115,14 +66,14 @@ DynamicModelDLL::DynamicModelDLL(const char *modName, const int y_length, const if ((dynamicHinstance == NULL) || dlerror()) { cerr << dlerror() << endl; - mexPrintf("MexPrintf:Error loading DLL"); + // mexPrintf("MexPrintf:Error loading DLL"); throw 1; } Dynamic = (DynamicFn) dlsym(dynamicHinstance, "Dynamic"); if ((Dynamic == NULL) || dlerror()) { cerr << dlerror() << endl; - mexPrintf("MexPrintf:Error finding DLL function"); + // mexPrintf("MexPrintf:Error finding DLL function"); throw 2; } #endif @@ -130,15 +81,15 @@ DynamicModelDLL::DynamicModelDLL(const char *modName, const int y_length, const } 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"); + // 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"); + // mexPrintf("MexPrintf: Unknown error in Call MATLAB function %s.\n", fName); + // mexErrMsgTxt("Err: Unknown error in Load and run DLL .\n"); return; } } @@ -163,15 +114,16 @@ DynamicModelDLL::close() void DynamicModelDLL::eval(const Vector &y, const TwoDMatrix &x, const Vector *modParams, - int it_, Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2) + int it_, Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) { - double *dresidual, *dg1 = NULL, *dg2 = NULL; + double *dresidual, *dg1 = NULL, *dg2 = NULL, *dg3 = NULL; //int length=y.length(); // not! + //mexPrintf("y.length() = %d\n",y.length()); if ((jcols-nExog) != y.length()) { // throw DLL Error - mexPrintf(" DLL Error: (jcols-nExog)!=ys.length() \n"); + // mexPrintf(" DLL Error: (jcols-nExog)!=ys.length() \n"); return; } if (residual.length() < length) // dummy or insufficient @@ -196,6 +148,11 @@ DynamicModelDLL::eval(const Vector &y, const TwoDMatrix &x, const Vector *modPa dg2 = const_cast(g2->base()); } dresidual = const_cast(residual.base()); + if (g3 != NULL) + { + dg3 = const_cast(g3->base()); + } + dresidual = const_cast(residual.base()); double *dy = const_cast(y.base()); double *dx = const_cast(x.base()); double *dbParams = const_cast(modParams->base()); @@ -216,25 +173,25 @@ DynamicModelDLL::eval(const Vector &y, const TwoDMatrix &x, const Vector *modPa #endif try { - Dynamic(dy, dx, nExog, dbParams, it_, dresidual, dg1, dg2); + Dynamic(dy, dx, nExog, dbParams, it_, dresidual, dg1, dg2, dg3); } catch (...) { - mexPrintf("MexPrintf: error in run Dynamic DLL \n"); + // 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) + Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) { - eval(y, x, modParams, nMax_lag, residual, g1, g2); + eval(y, x, modParams, nMax_lag, residual, g1, g2, g3); }; void DynamicModelDLL::eval(const Vector &y, const Vector &x, const Vector *modParams, - Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2) + Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) { /** ignore given exogens and create new 2D x matrix since @@ -244,7 +201,7 @@ DynamicModelDLL::eval(const Vector &y, const Vector &x, const Vector *modParams, TwoDMatrix &mx = *(new TwoDMatrix(nMax_lag+1, nExog)); mx.zeros(); // initialise shocks to 0s - eval(y, mx, modParams, nMax_lag, residual, g1, g2); + eval(y, mx, modParams, nMax_lag, residual, g1, g2, g3); }; diff --git a/mex/sources/k_order_perturbation/dynamic_dll.h b/mex/sources/k_order_perturbation/dynamic_dll.h index dcefe332f..9cf3e9204 100644 --- a/mex/sources/k_order_perturbation/dynamic_dll.h +++ b/mex/sources/k_order_perturbation/dynamic_dll.h @@ -43,7 +43,7 @@ typedef void *(DynamicFn) typedef void (*DynamicFn) #endif (double *y, double *x, int nb_row_x, double *params, - int it_, double *residual, double *g1, double *g2); + int it_, double *residual, double *g1, double *g2, double *g3); //DynamicFn Dynamic; @@ -86,18 +86,17 @@ public: // evaluate Dynamic model DLL void eval(double *y, double *x, int nb_row_x, double *params, - int it_, double *residual, double *g1, double *g2) + int it_, double *residual, double *g1, double *g2, double *g3) { - Dynamic(y, x, nb_row_x, params, it_, residual, g1, g2); + Dynamic(y, x, nb_row_x, params, it_, residual, g1, g2, g3); }; void eval(const Vector &y, const Vector &x, const Vector *params, - Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2); + Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3); void eval(const Vector &y, const TwoDMatrix &x, const Vector *params, - int it_, Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2); + int it_, Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3); void eval(const Vector &y, const TwoDMatrix &x, const Vector *params, - Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2); - // void eval(const Vector&y, const TwoDMatrix&x, const Vector* params, - // Vector& residual, double *g1, double *g2); + Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3); + // close DLL: If the referenced object was successfully closed, // close() returns 0, non 0 otherwise int close(); diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.cpp b/mex/sources/k_order_perturbation/k_ord_dynare.cpp index 5c7920a06..8c358cc78 100644 --- a/mex/sources/k_order_perturbation/k_ord_dynare.cpp +++ b/mex/sources/k_order_perturbation/k_ord_dynare.cpp @@ -84,7 +84,7 @@ KordpDynare::KordpDynare(const char **endo, int num_endo, } catch (...) { - mexPrintf("k_ord_dynare: dynare constructor, error in StateNamelist construction.\n"); + // mexPrintf("k_ord_dynare: dynare constructor, error in StateNamelist construction.\n"); throw DynareException(__FILE__, __LINE__, string("Could not construct Name Lists. \n")); } } @@ -159,7 +159,7 @@ KordpDynare::evaluateSystem(Vector &out, const Vector &yy, const Vector &xx) { dynamicDLL.eval(yy, xx, //int nb_row_x, params, //int it_, - out, NULL, NULL); + out, NULL, NULL, NULL); } @@ -176,7 +176,7 @@ KordpDynare::evaluateSystem(Vector &out, const Vector &yym, const Vector &yy, #endif dynamicDLL.eval(yy, xx, //int nb_row_x, params, //int it_, - out, NULL, NULL); + out, NULL, NULL, NULL); } /************************************************ * this is main derivative calculation functin that indirectly calls dynamic.dll @@ -188,13 +188,14 @@ KordpDynare::calcDerivatives(const Vector &yy, const Vector &xx) { // Hessian TwoDMatrix *g2; TwoDMatrix *g2 = NULL; + TwoDMatrix *g3 = NULL; //Jacobian TwoDMatrix *g1 = new TwoDMatrix(nY, nJcols); // generate g1 for jacobian g1->zeros(); if ((nJcols != g1->ncols()) && (nY != g1->nrows())) { - mexPrintf("k_ord_dynare.cpp: Error in calcDerivatives: Created wrong jacobian"); + // mexPrintf("k_ord_dynare.cpp: Error in calcDerivatives: Created wrong jacobian"); return; } @@ -210,6 +211,11 @@ KordpDynare::calcDerivatives(const Vector &yy, const Vector &xx) //g2->print(); #endif } + if (nOrder > 2) + { + g3 = new TwoDMatrix((int) (*NNZD)[2],3); + g3->zeros(); + } Vector &out = *(new Vector(nY)); out.zeros(); const Vector *llxYYp; // getting around the constantness @@ -228,34 +234,25 @@ KordpDynare::calcDerivatives(const Vector &yy, const Vector &xx) #endif try { - dynamicDLL.eval(llxYY, xx, //int nb_row_x, - params, //int it_, - out, g1, g2); - // } + dynamicDLL.eval(llxYY, xx, params, out, g1, g2, g3); } catch (...) { - mexPrintf("k_ord_dynare.cpp: Error in dynamicDLL.eval in calcDerivatives"); + // mexPrintf("k_ord_dynare.cpp: Error in dynamicDLL.eval in calcDerivatives"); return; } if ((nJcols != g1->ncols()) && (nY != g1->nrows())) { - mexPrintf("k_ord_dynare.cpp: Error in calcDerivatives: dynamicDLL.eval returned wrong jacobian"); + // mexPrintf("k_ord_dynare.cpp: Error in calcDerivatives: dynamicDLL.eval returned wrong jacobian"); return; } - // ReorderCols(g1, JacobianIndices); and populate container + populateDerivativesContainer(g1, 1, JacobianIndices); if (nOrder > 1) - { - // ReorderBlocks(g2,JacobianIndices); -#ifdef DEBUG - mexPrintf(" post dll g2 cols %d rows %d \n", g2->numCols(), g2->numRows()); - for (int ii=0;iinumRows(); ii++) - mexPrintf(" g2[%d]: %d %d %f \n", ii, (int)g2->get(ii,0),(int)g2->get(ii,1),g2->get(ii,2)); - //g2->print(); -#endif populateDerivativesContainer(g2, 2, JacobianIndices); - } + if (nOrder > 2) + populateDerivativesContainer(g3, 3, JacobianIndices); + } /* This version is not currently in use */ @@ -271,7 +268,7 @@ KordpDynare::calcDerivatives(const Vector &yy, ogu::Jacobian &jacob) xx.zeros(); dynamicDLL.eval(yy, xx, //int nb_row_x, params, //int it_, - out, jj, NULL); + out, jj, NULL, NULL); // model derivatives FSSparseTensor instance FSSparseTensor &mdTi = *(new FSSparseTensor(1, jj->ncols(), jj->nrows())); for (int i = 0; i < jj->ncols(); i++) @@ -309,7 +306,7 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vectornrows())); + FSSparseTensor *mdTi = (new FSSparseTensor(ord, nJcols, nY)); IntSequence s(ord, 0); //s[0] = 0; @@ -364,14 +361,47 @@ KordpDynare::populateDerivativesContainer(TwoDMatrix *g, int ord, const vector revOrder(nJcols1); + for (int i = 0; i < nJcols1; 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 s0 = (int)floor(i1/nJcols2); + int i2 = i1 - nJcols2*s0; + int s1 = (int)floor(i2/nJcols); + int s2 = i2 - nJcols*s1; + if (s0 < nJcols1) + s[0] = revOrder[s0]; + else + s[0] = s0; + if (s1 < nJcols1) + s[1] = revOrder[s1]; + else + s[1] = s1; + if (s2 < nJcols1) + s[2] = revOrder[s2]; + else + s[2] = s2; + if ((s[2] >= s[1]) && (s[1] >= s[0])) + { + double x = g->get(i,2); + mdTi->insert(s, j, x); + } + } + } #ifdef DEBUG mexPrintf("k_ord_dynare.cpp: END populate FSSparseTensor in calcDerivatives: ord=%d \n",ord); - mdTi->print(); #endif // md container //md.clear();// this is to be used only for 1st order!! md.remove(Symmetry(ord)); - md.insert(mdTi); //(&mdTi); + md.insert(mdTi); } void @@ -406,7 +436,7 @@ KordpDynare::LLxSteady(const Vector &yS) { if ((nJcols-nExog) == yS.length()) { - mexPrintf("k_ord_dynare.cpp: Warning in LLxSteady: ySteady already. right size"); + // mexPrintf("k_ord_dynare.cpp: Warning in LLxSteady: ySteady already. right size"); return NULL; } // create temporary square 2D matrix size nEndo x nEndo (sparse) @@ -426,13 +456,13 @@ KordpDynare::LLxSteady(const Vector &yS) } catch (const TLException &e) { - mexPrintf("Caugth TL exception in LLxSteady: "); + // mexPrintf("Caugth TL exception in LLxSteady: "); e.print(); return NULL; // 255; } catch (...) { - mexPrintf(" Error in LLxSteady - wrong index?"); + // mexPrintf(" Error in LLxSteady - wrong index?"); } #ifdef DEBUG @@ -514,13 +544,13 @@ KordpDynare::ReorderDynareJacobianIndices(const vector *varOrder) } catch (const TLException &e) { - mexPrintf("Caugth TL exception in ReorderIndices: "); + // mexPrintf("Caugth TL exception in ReorderIndices: "); e.print(); return NULL; // 255; } catch (...) { - mexPrintf(" Error in ReorderIndices - wrong index?"); + // mexPrintf(" Error in ReorderIndices - wrong index?"); } //add the indices for the nExog exogenous jacobians for (j = nJcols-nExog; j < nJcols; j++) @@ -547,7 +577,7 @@ KordpDynare::ReorderCols(TwoDMatrix *tdx, const vector *vOrder) if (tdx->ncols() > vOrder->size()) { - mexPrintf(" Error in ReorderColumns - size of order var is too small"); + // mexPrintf(" Error in ReorderColumns - size of order var is too small"); return; } TwoDMatrix tmp(*tdx); // temporary 2D matrix @@ -561,13 +591,13 @@ KordpDynare::ReorderCols(TwoDMatrix *tdx, const vector *vOrder) } catch (const TLException &e) { - printf("Caugth TL exception in ReorderColumns: "); + //printf("Caugth TL exception in ReorderColumns: "); e.print(); return; // 255; } catch (...) { - mexPrintf(" Error in ReorderColumns - wrong index?"); + // mexPrintf(" Error in ReorderColumns - wrong index?"); } } void @@ -585,13 +615,13 @@ KordpDynare::ReorderCols(TwoDMatrix *tdx, const int *vOrder) } catch (const TLException &e) { - printf("Caugth TL exception in ReorderColumns: "); + //printf("Caugth TL exception in ReorderColumns: "); e.print(); return; // 255; } catch (...) { - mexPrintf(" Error in ReorderColumns - wrong index?"); + // mexPrintf(" Error in ReorderColumns - wrong index?"); } } @@ -610,7 +640,7 @@ KordpDynare::ReorderBlocks(TwoDMatrix *tdx, const vector *vOrder) int ibOrder = (int) dbOrder; if ((double) ibOrder != dbOrder || ibOrder > nOrder) { - mexPrintf(" Error in ReorderBlocks - wrong order %d", dbOrder); + // mexPrintf(" Error in ReorderBlocks - wrong order %d", dbOrder); return; } TwoDMatrix tmp(*tdx); // temporary 2D matrix @@ -633,7 +663,7 @@ KordpDynare::ReorderBlocks(TwoDMatrix *tdx, const vector *vOrder) //ReorderColumns(TwoDMatrix * subtdx, const vector * vOrder) if (tdx->ncols() > vOrder->size()) { - mexPrintf(" Error in ReorderColumns - size of order var is too small"); + // mexPrintf(" Error in ReorderColumns - size of order var is too small"); return; } // reorder the columns @@ -644,13 +674,13 @@ KordpDynare::ReorderBlocks(TwoDMatrix *tdx, const vector *vOrder) } catch (const TLException &e) { - printf("Caugth TL exception in ReorderColumns: "); + // printf("Caugth TL exception in ReorderColumns: "); e.print(); return; // 255; } catch (...) { - mexPrintf(" Error in ReorderColumns - wrong index?"); + // mexPrintf(" Error in ReorderColumns - wrong index?"); } } } diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cpp b/mex/sources/k_order_perturbation/k_order_perturbation.cpp index 6dfe27c04..4ae559d2b 100644 --- a/mex/sources/k_order_perturbation/k_order_perturbation.cpp +++ b/mex/sources/k_order_perturbation/k_order_perturbation.cpp @@ -86,10 +86,10 @@ extern "C" { mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs < 5) - mexErrMsgTxt("Must have at least 5 input parameters.\n"); - if (nlhs == 0) - mexErrMsgTxt("Must have at least 1 output parameter.\n"); + // if (nrhs < 5) + // mexErrMsgTxt("Must have at least 5 input parameters.\n"); + //if (nlhs == 0) + // mexErrMsgTxt("Must have at least 1 output parameter.\n"); const mxArray *dr = prhs[0]; const int check_flag = (int) mxGetScalar(prhs[1]); @@ -119,7 +119,14 @@ extern "C" { kOrder = (int) mxGetScalar(mxFldp); else kOrder = 1; + +// if (kOrder == 1 && nlhs != 1) +// mexErrMsgTxt("k_order_perturbation at order 1 requires exactly 1 argument in output\'n"); +// else if (kOrder > 1 && nlhs != kOrder+1) +// mexErrMsgTxt("k_order_perturbation at order > 1 requires exactly order + 1 argument in output\'n"); + + double qz_criterium = 1+1e-6; mxFldp = mxGetField(options_, 0, "qz_criterium"); if (mxIsNumeric(mxFldp)) @@ -146,9 +153,11 @@ extern "C" { // mxFldp = mxGetField(oo_, 0,"steady_state" ); // use in order of declaration mxFldp = mxGetField(dr, 0, "ys"); // and not in order of dr.order_var + //mexPrintf("mxFldp %x",mxFldp); // mxFldp = mxGetField(oo_, 0,"dyn_ys" ); // and NOT extended ys dparams = (double *) mxGetData(mxFldp); const int nSteady = (int) mxGetM(mxFldp); + // mexPrintf("nSteady %d\n",nSteady); Vector *ySteady = new Vector(dparams, nSteady); mxFldp = mxGetField(dr, 0, "nstatic"); @@ -198,6 +207,8 @@ extern "C" { dparams = (double *) mxGetData(mxFldp); npar = (int) mxGetN(mxFldp); int nrows = (int) mxGetM(mxFldp); + + #ifdef DEBUG mexPrintf("ll_Incidence nrows=%d, ncols = %d .\n", nrows, npar); #endif @@ -227,10 +238,9 @@ extern "C" { #endif const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns - mexPrintf("k_order_perturbation: jcols= %d .\n", jcols); mxFldp = mxGetField(M_, 0, "var_order_endo_names"); - mexPrintf("k_order_perturbation: Get nendo .\n"); + // mexPrintf("k_order_perturbation: Get nendo .\n"); const int nendo = (int) mxGetM(mxFldp); const int widthEndo = (int) mxGetN(mxFldp); const char **endoNamesMX = DynareMxArrayToString(mxFldp, nendo, widthEndo); @@ -301,6 +311,7 @@ extern "C" { #ifdef DEBUG mexPrintf("k_order_perturbation: Calling dynare constructor .\n"); #endif + // make KordpDynare object KordpDynare dynare(endoNamesMX, nEndo, exoNamesMX, nExog, nPar, // paramNames, ySteady, vCov, modParams, nStat, nPred, nForw, nBoth, @@ -311,6 +322,7 @@ extern "C" { #ifdef DEBUG mexPrintf("k_order_perturbation: Call Approximation constructor with qz_criterium=%f \n", qz_criterium); #endif + Approximation app(dynare, journal, nSteps, false, qz_criterium); // run stochastic steady #ifdef DEBUG @@ -367,29 +379,31 @@ extern "C" { int nb_row_x; ysteady = NULL; - if (nlhs >= 1) + if (kOrder == 1) { /* Set the output pointer to the output matrix ysteady. */ - plhs[0] = mxCreateDoubleMatrix(nEndo, 1, mxREAL); - /* Create a C pointer to a copy of the output ysteady. */ - TwoDMatrix tmp_ss(nEndo, 1, mxGetPr(plhs[0])); - tmp_ss = (const TwoDMatrix &)ss; + map::const_iterator cit = mm.begin(); + ++cit; + plhs[0] = mxCreateDoubleMatrix((*cit).second.numRows(), (*cit).second.numCols(), mxREAL); + TwoDMatrix g((*cit).second.numRows(), (*cit).second.numCols(), mxGetPr(plhs[0])); + g = (const TwoDMatrix &)(*cit).second; + // /* Create a C pointer to a copy of the output ysteady. */ + // TwoDMatrix tmp_ss(nEndo, 1, mxGetPr(plhs[0])); + // tmp_ss = (const TwoDMatrix &)ss; #ifdef DEBUG // tmp_ss.print(); !! This print Crashes??? #endif } - if (nlhs >= 2) + if (kOrder >= 2) { - /* Set the output pointer to the combined output matrix gyu = [gy gu]. */ - int ii = 1; + int ii = 0; for (map::const_iterator cit = mm.begin(); ((cit != mm.end()) && (ii < nlhs)); ++cit) { - //if ((*cit).first!="g_0") { plhs[ii] = mxCreateDoubleMatrix((*cit).second.numRows(), (*cit).second.numCols(), mxREAL); - TwoDMatrix dgyu((*cit).second.numRows(), (*cit).second.numCols(), mxGetPr(plhs[ii])); - dgyu = (const TwoDMatrix &)(*cit).second; + TwoDMatrix g_ii((*cit).second.numRows(), (*cit).second.numCols(), mxGetPr(plhs[ii])); + g_ii = (const TwoDMatrix &)(*cit).second; #ifdef DEBUG mexPrintf("k_order_perturbation: cit %d print: \n", ii); (*cit).second.print(); @@ -462,9 +476,64 @@ extern "C" { e.print(errfd); fclose(errfd); return; // 255; - } //catch + } //catch }; // end of mexFunction() }; // end of extern C + +////////////////////////////////////////////////////// +// 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; + cNamesMX = (char **) calloc(len, sizeof(char *)); + for(int i = 0; i < len; i++) + cNamesMX[i] = (char *) calloc(width+1, sizeof(char)); + +#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 **) calloc(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 *) calloc(strlen(cNamesMX[j])+1, sizeof(char)); + strcpy(token, cNamesMX[j]); + ret[j] = token; +#ifdef DEBUG + mexPrintf("ret [%d]= %s \n", j, ret[j]); +#endif + } + mxFree(cNamesMX); + return ret; +} + #endif // ifdef MATLAB_MEX_FILE to exclude mexFunction for other applications