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

expand extended preprocessor + first implementation of Petsc interface

issue#70
Michel Juillard 2015-08-27 10:00:27 +02:00
parent 73a1ca92c0
commit a538e29f1e
5 changed files with 112 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

102
DynamicModel.cc
View File

@ -4467,7 +4467,7 @@ DynamicModel::writeResidualsC(const string &basename, bool cuda) const
deriv_node_temp_terms_t tef_terms;
ostringstream model_output; // Used for storing model equations
writeModelEquations(model_output, oCDynamicModel);
writeModelEquations(model_output, oCDynamic2Model);
mDynamicModelFile << " double lhs, rhs;" << endl
<< endl
@ -4540,6 +4540,106 @@ DynamicModel::writeFirstDerivativesC(const string &basename, bool cuda) const
// using compressed sparse row format (CSR)
void
DynamicModel::writeFirstDerivativesC_csr(const string &basename, bool cuda) const
{
string filename = basename + "_first_derivatives.c";
ofstream mDynamicModelFile, mDynamicMexFile;
mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
if (!mDynamicModelFile.is_open())
{
cerr << "Error: Can't open file " << filename << " for writing" << endl;
exit(EXIT_FAILURE);
}
mDynamicModelFile << "/*" << endl
<< " * " << filename << " : Computes first order derivatives of the model for Dynare" << endl
<< " *" << endl
<< " * Warning : this file is generated automatically by Dynare" << endl
<< " * from model " << basename << "(.mod)" << endl
<< " */" << endl
<< "#include <math.h>" << endl;
mDynamicModelFile << "#include <stdlib.h>" << endl;
mDynamicModelFile << "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl
<< "#define min(a, b) (((a) > (b)) ? (b) : (a))" << endl;
// Write function definition if oPowerDeriv is used
writePowerDerivCHeader(mDynamicModelFile);
mDynamicModelFile << "void FirstDerivatives(const double *y, double *x, int nb_row_x, double *params, double *steady_state, int it_, double *residual, int *row_ptr, int *col_ptr, double *value)" << endl
<< "{" << endl;
int cols_nbr = 3*symbol_table.endo_nbr() + symbol_table.exo_nbr() + symbol_table.exo_det_nbr();
// this is always empty here, but needed by d1->writeOutput
deriv_node_temp_terms_t tef_terms;
// Indexing derivatives in column order
vector<derivative> D;
for (first_derivatives_t::const_iterator it = first_derivatives.begin();
it != first_derivatives.end(); it++)
{
int eq = it->first.first;
int dynvar = it->first.second;
int lag = getLagByDerivID(dynvar);
int symb = getSymbIDByDerivID(dynvar);
SymbolType type = getTypeByDerivID(dynvar);
int col_id;
switch(type)
{
case eEndogenous:
col_id = symb+(lag+1)*symbol_table.endo_nbr();
break;
case eExogenous:
col_id = symb+3*symbol_table.endo_nbr();
break;
case eExogenousDet:
col_id = symb+3*symbol_table.endo_nbr()+symbol_table.exo_nbr();
break;
default:
std::cerr << "This case shouldn't happen" << std::endl;
exit(1);
}
derivative deriv(col_id + eq*cols_nbr,col_id,eq,it->second);
D.push_back(deriv);
}
sort(D.begin(), D.end(), derivative_less_than() );
// writing sparse Jacobian
vector<int> row_ptr(equations.size());
fill(row_ptr.begin(),row_ptr.end(),0.0);
int k = 0;
for(vector<derivative>::const_iterator it = D.begin(); it != D.end(); ++it)
{
row_ptr[it->row_nbr]++;
mDynamicModelFile << "col_ptr[" << k << "] "
<< "=" << it->col_nbr << ";" << endl;
mDynamicModelFile << "value[" << k << "] = ";
// oCstaticModel makes reference to the static variables
it->value->writeOutput(mDynamicModelFile, oCDynamic2Model, temporary_terms, tef_terms);
mDynamicModelFile << ";" << endl;
k++;
}
// row_ptr must point to the relative address of the first element of the row
int cumsum = 0;
mDynamicModelFile << "int row_ptr_data[" << row_ptr.size() + 1 << "] = { 0";
for (vector<int>::iterator it=row_ptr.begin(); it != row_ptr.end(); ++it)
{
cumsum += *it;
mDynamicModelFile << ", " << cumsum;
}
mDynamicModelFile << "};" << endl
<< "int i;" << endl
<< "for (i=0; i < " << row_ptr.size() + 1 << "; i++) row_ptr[i] = row_ptr_data[i];" << endl;
mDynamicModelFile << "}" << endl;
// writePowerDeriv(mDynamicModelFile, true);
mDynamicModelFile.close();
}
void
DynamicModel::writeSecondDerivativesC_csr(const string &basename, bool cuda) const
{

2
DynamicModel.hh
View File

@ -480,6 +480,8 @@ public:
void writeResidualsC(const string &basename, bool cuda) const;
//! Writes C file containing first order derivatives of model evaluated at steady state
void writeFirstDerivativesC(const string &basename, bool cuda) const;
//! Writes C file containing first order derivatives of model evaluated at steady state (conpressed sparse column)
void writeFirstDerivativesC_csr(const string &basename, bool cuda) const;
//! Writes C file containing second order derivatives of model evaluated at steady state (compressed sparse column)
void writeSecondDerivativesC_csr(const string &basename, bool cuda) const;
//! Writes C file containing third order derivatives of model evaluated at steady state (compressed sparse column)

6
ExprNode.cc
View File

@ -659,6 +659,10 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
i = datatree.getDynJacobianCol(datatree.getDerivID(symb_id, lag)) + ARRAY_SUBSCRIPT_OFFSET(output_type);
output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
break;
case oCDynamic2Model:
i = symb_id + (lag+1)*datatree.symbol_table.endo_nbr() + ARRAY_SUBSCRIPT_OFFSET(output_type);
output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
break;
case oCStaticModel:
case oMatlabStaticModel:
case oMatlabStaticModelSparse:
@ -710,6 +714,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
output << "x(it_, " << i << ")";
break;
case oCDynamicModel:
case oCDynamic2Model:
if (lag == 0)
output << "x[it_+" << i << "*nb_row_x]";
else if (lag > 0)
@ -755,6 +760,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
output << "x(it_, " << i << ")";
break;
case oCDynamicModel:
case oCDynamic2Model:
if (lag == 0)
output << "x[it_+" << i << "*nb_row_x]";
else if (lag > 0)

2
ExprNode.hh
View File

@ -65,6 +65,7 @@ enum ExprNodeOutputType
oMatlabStaticModelSparse, //!< Matlab code, static block decomposed model
oMatlabDynamicModelSparse, //!< Matlab code, dynamic block decomposed model
oCDynamicModel, //!< C code, dynamic model
oCDynamic2Model, //!< C code, dynamic model, alternative numbering of endogenous variables
oCStaticModel, //!< C code, static model
oMatlabOutsideModel, //!< Matlab code, outside model block (for example in initval)
oLatexStaticModel, //!< LaTeX code, static model
@ -87,6 +88,7 @@ enum ExprNodeOutputType
|| (output_type) == oSteadyStateFile)
#define IS_C(output_type) ((output_type) == oCDynamicModel \
|| (output_type) == oCDynamic2Model \
|| (output_type) == oCStaticModel \
|| (output_type) == oCDynamicSteadyStateOperator \
|| (output_type) == oCSteadyStateFile)

2
ModFile.cc
View File

@ -972,7 +972,7 @@ ModFile::writeExternalFilesCC(const string &basename, FileOutputType output) con
// dynamic_model.writeResidualsC(basename, cuda);
// dynamic_model.writeParamsDerivativesFileC(basename, cuda);
dynamic_model.writeResidualsC(basename, cuda);
dynamic_model.writeFirstDerivativesC(basename, cuda);
dynamic_model.writeFirstDerivativesC_csr(basename, cuda);
if (output == second)
dynamic_model.writeSecondDerivativesC_csr(basename, cuda);