stepan
/
preprocessor
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

preprocessor: implement third-order derivatives in USE_DLL mode 

git-svn-id: https://www.dynare.org/svn/dynare/trunk@3036 ac1d8469-bf42-47a9-8791-bf33cf982152
issue#70
sebastien 2009-10-12 16:09:16 +00:00
parent 330afb32fe
commit 66684bbb9a
2 changed files with 59 additions and 32 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

84
DynamicModel.cc
View File

@ -1167,7 +1167,7 @@ DynamicModel::writeDynamicCFile(const string &dynamic_basename) const
<< "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
<< "{" << endl << "{" << endl
<< " double *y, *x, *params;" << endl << " double *y, *x, *params;" << endl
<< " double *residual, *g1, *v2;" << endl << " double *residual, *g1, *v2, *v3;" << endl
<< " int nb_row_x, it_;" << endl << " int nb_row_x, it_;" << endl
<< endl << endl
<< " /* Create a pointer to the input matrix y. */" << endl << " /* Create a pointer to the input matrix y. */" << endl
@ -1210,12 +1210,18 @@ DynamicModel::writeDynamicCFile(const string &dynamic_basename) const
<< " /* Set the output pointer to the output matrix v2. */" << endl << " /* Set the output pointer to the output matrix v2. */" << endl
<< " plhs[2] = mxCreateDoubleMatrix(" << NNZDerivatives[1] << ", " << 3 << " plhs[2] = mxCreateDoubleMatrix(" << NNZDerivatives[1] << ", " << 3
<< ", mxREAL);" << endl << ", mxREAL);" << endl
<< " /* Create a C pointer to a copy of the output matrix g1. */" << endl
<< " v2 = mxGetPr(plhs[2]);" << endl << " v2 = mxGetPr(plhs[2]);" << endl
<< " }" << endl << " }" << endl
<< endl << endl
<< " if (nlhs >= 4)" << endl
<< " {" << endl
<< " /* Set the output pointer to the output matrix v3. */" << endl
<< " plhs[3] = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 3 << ", mxREAL);" << endl
<< " v3 = mxGetPr(plhs[3]);" << endl
<< " }" << endl
<< endl
<< " /* Call the C subroutines. */" << endl << " /* Call the C subroutines. */" << endl
<< " Dynamic(y, x, nb_row_x, params, it_, residual, g1, v2);" << endl << " Dynamic(y, x, nb_row_x, params, it_, residual, g1, v2, v3);" << endl
<< "}" << endl; << "}" << endl;
mDynamicModelFile.close(); mDynamicModelFile.close();
} }
@ -1735,16 +1741,13 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
int col_nb = id1 * dynJacobianColsNbr + id2; int col_nb = id1 * dynJacobianColsNbr + id2;
int col_nb_sym = id2 * dynJacobianColsNbr + id1; int col_nb_sym = id2 * dynJacobianColsNbr + id1;
hessian_output << "v2"; sparseHelper(2, hessian_output, k, 0, output_type);
hessianHelper(hessian_output, k, 0, output_type);
hessian_output << "=" << eq + 1 << ";" << endl; hessian_output << "=" << eq + 1 << ";" << endl;
hessian_output << "v2"; sparseHelper(2, hessian_output, k, 1, output_type);
hessianHelper(hessian_output, k, 1, output_type);
hessian_output << "=" << col_nb + 1 << ";" << endl; hessian_output << "=" << col_nb + 1 << ";" << endl;
hessian_output << "v2"; sparseHelper(2, hessian_output, k, 2, output_type);
hessianHelper(hessian_output, k, 2, output_type);
hessian_output << "="; hessian_output << "=";
d2->writeOutput(hessian_output, output_type, temporary_terms); d2->writeOutput(hessian_output, output_type, temporary_terms);
hessian_output << ";" << endl; hessian_output << ";" << endl;
@ -1754,18 +1757,15 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
// Treating symetric elements // Treating symetric elements
if (id1 != id2) if (id1 != id2)
{ {
hessian_output << "v2"; sparseHelper(2, hessian_output, k, 0, output_type);
hessianHelper(hessian_output, k, 0, output_type);
hessian_output << "=" << eq + 1 << ";" << endl; hessian_output << "=" << eq + 1 << ";" << endl;
hessian_output << "v2"; sparseHelper(2, hessian_output, k, 1, output_type);
hessianHelper(hessian_output, k, 1, output_type);
hessian_output << "=" << col_nb_sym + 1 << ";" << endl; hessian_output << "=" << col_nb_sym + 1 << ";" << endl;
hessian_output << "v2"; sparseHelper(2, hessian_output, k, 2, output_type);
hessianHelper(hessian_output, k, 2, output_type); hessian_output << "=";
hessian_output << "=v2"; sparseHelper(2, hessian_output, k-1, 2, output_type);
hessianHelper(hessian_output, k-1, 2, output_type);
hessian_output << ";" << endl; hessian_output << ";" << endl;
k++; k++;
@ -1790,10 +1790,18 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
// Reference column number for the g3 matrix // Reference column number for the g3 matrix
int ref_col = id1 * hessianColsNbr + id2 * dynJacobianColsNbr + id3; int ref_col = id1 * hessianColsNbr + id2 * dynJacobianColsNbr + id3;
third_derivatives_output << " v3(" << ++k << ",:) = [" sparseHelper(3, third_derivatives_output, k, 0, output_type);
<< eq + 1 << "," << ref_col + 1 << ","; third_derivatives_output << "=" << eq + 1 << ";" << endl;
sparseHelper(3, third_derivatives_output, k, 1, output_type);
third_derivatives_output << "=" << ref_col + 1 << ";" << endl;
sparseHelper(3, third_derivatives_output, k, 2, output_type);
third_derivatives_output << "=";
d3->writeOutput(third_derivatives_output, output_type, temporary_terms); d3->writeOutput(third_derivatives_output, output_type, temporary_terms);
third_derivatives_output << "];" << endl; third_derivatives_output << ";" << endl;
k++;
// Compute the column numbers for the 5 other permutations of (id1,id2,id3) and store them in a set (to avoid duplicates if two indexes are equal) // Compute the column numbers for the 5 other permutations of (id1,id2,id3) and store them in a set (to avoid duplicates if two indexes are equal)
set<int> cols; set<int> cols;
@ -1806,10 +1814,20 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
int k2 = 0; // Keeps the offset of the permutation relative to k int k2 = 0; // Keeps the offset of the permutation relative to k
for (set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++) for (set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++)
if (*it2 != ref_col) if (*it2 != ref_col)
third_derivatives_output << " v3(" << k + ++k2 << ",:) = [" {
<< eq + 1 << "," << *it2 + 1 << "," sparseHelper(3, third_derivatives_output, k+k2, 0, output_type);
<< "v3(" << k << ",3)];" << endl; third_derivatives_output << "=" << eq + 1 << ";" << endl;
sparseHelper(3, third_derivatives_output, k+k2, 1, output_type);
third_derivatives_output << "=" << *it2 + 1 << ";" << endl;
sparseHelper(3, third_derivatives_output, k+k2, 2, output_type);
third_derivatives_output << "=";
sparseHelper(3, third_derivatives_output, k, 2, output_type);
third_derivatives_output << ";" << endl;
k2++;
}
k += k2; k += k2;
} }
@ -1864,7 +1882,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
} }
else else
{ {
DynamicOutput << "void Dynamic(double *y, double *x, int nb_row_x, double *params, int it_, double *residual, double *g1, double *v2)" << endl DynamicOutput << "void Dynamic(double *y, double *x, int nb_row_x, double *params, int it_, double *residual, double *g1, double *v2, double *v3)" << endl
<< "{" << endl << "{" << endl
<< " double lhs, rhs;" << endl << " double lhs, rhs;" << endl
<< endl << endl
@ -1879,7 +1897,6 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
<< " }" << endl; << " }" << endl;
if (second_derivatives.size()) if (second_derivatives.size())
{
DynamicOutput << " /* Hessian for endogenous and exogenous variables */" << endl DynamicOutput << " /* Hessian for endogenous and exogenous variables */" << endl
<< " if (v2 == NULL)" << endl << " if (v2 == NULL)" << endl
<< " return;" << endl << " return;" << endl
@ -1887,7 +1904,16 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
<< " {" << endl << " {" << endl
<< hessian_output.str() << hessian_output.str()
<< " }" << endl; << " }" << endl;
}
if (third_derivatives.size())
DynamicOutput << " /* Third derivatives for endogenous and exogenous variables */" << endl
<< " if (v3 == NULL)" << endl
<< " return;" << endl
<< " else" << endl
<< " {" << endl
<< third_derivatives_output.str()
<< " }" << endl;
DynamicOutput << "}" << endl << endl; DynamicOutput << "}" << endl << endl;
} }
} }
@ -2820,13 +2846,13 @@ DynamicModel::jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOut
} }
void void
DynamicModel::hessianHelper(ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const DynamicModel::sparseHelper(int order, ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const
{ {
output << LEFT_ARRAY_SUBSCRIPT(output_type); output << "v" << order << LEFT_ARRAY_SUBSCRIPT(output_type);
if (IS_MATLAB(output_type)) if (IS_MATLAB(output_type))
output << row_nb + 1 << ", " << col_nb + 1; output << row_nb + 1 << ", " << col_nb + 1;
else else
output << row_nb + col_nb * NNZDerivatives[1]; output << row_nb + col_nb * NNZDerivatives[order-1];
output << RIGHT_ARRAY_SUBSCRIPT(output_type); output << RIGHT_ARRAY_SUBSCRIPT(output_type);
} }

7
DynamicModel.hh
View File

@ -137,9 +137,10 @@ private:
/*! Writes either (i+1,j+1) or [i+j*no_eq] */ /*! Writes either (i+1,j+1) or [i+j*no_eq] */
void jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const; void jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const;
//! Helper for writing the sparse Hessian elements in MATLAB and C //! Helper for writing the sparse Hessian or third derivatives in MATLAB and C
/*! Writes either (i+1,j+1) or [i+j*NNZDerivatives[1]] */ /*! If order=2, writes either v2(i+1,j+1) or v2[i+j*NNZDerivatives[1]]
void hessianHelper(ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const; If order=3, writes either v3(i+1,j+1) or v3[i+j*NNZDerivatives[2]] */
void sparseHelper(int order, ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const;
//! Write chain rule derivative of a recursive equation w.r. to a variable //! Write chain rule derivative of a recursive equation w.r. to a variable
void writeChainRuleDerivative(ostream &output, int eq, int var, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const; void writeChainRuleDerivative(ostream &output, int eq, int var, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;