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

Block decomposition: in per-block files, construct sparse Jacobians more efficiently 

Use 3-column format before calling sparse().

Ensure that the 3-column matrix is constructed in column-major order: data
locality will greatly improve performance once we implement use_dll.
issue#70
Sébastien Villemot 2020-06-19 16:09:38 +02:00
parent 93c18d514b
commit 92aff91066
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
2 changed files with 121 additions and 70 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

157
src/DynamicModel.cc
View File

@ -280,25 +280,7 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
<< " //" << endl << " //" << endl
<< " % // Simulation type " << " % // Simulation type "
<< BlockSim(simulation_type) << " //" << endl << BlockSim(simulation_type) << " //" << endl
<< " % ////////////////////////////////////////////////////////////////////////" << endl << " % ////////////////////////////////////////////////////////////////////////" << endl;
<< " if stochastic_mode" << endl
<< " g1 = spalloc(" << block_size << ", " << blocks_jacob_cols_endo[blk].size() << ", " << nze_stochastic << ");" << endl
<< " g1_x = spalloc(" << block_size << ", " << blocks_jacob_cols_exo[blk].size() << ", " << nze_exo << ");" << endl
<< " g1_xd = spalloc(" << block_size << ", " << blocks_jacob_cols_exo_det[blk].size() << ", " << nze_exo_det << ");" << endl
<< " g1_o = spalloc(" << block_size << ", " << blocks_jacob_cols_other_endo[blk].size() << ", " << nze_other_endo << ");" << endl;
if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete
|| simulation_type == BlockSimulationType::solveTwoBoundariesSimple)
output << " else" << endl
<< " g1 = spalloc(" << block_mfs_size
<< ", " << 3*block_mfs_size << ", " << nze_deterministic << ");" << endl;
else if (simulation_type == BlockSimulationType::solveBackwardSimple
|| simulation_type == BlockSimulationType::solveForwardSimple
|| simulation_type == BlockSimulationType::solveBackwardComplete
|| simulation_type == BlockSimulationType::solveForwardComplete)
output << " else" << endl
<< " g1 = spalloc(" << block_mfs_size
<< ", " << block_mfs_size << ", " << nze_deterministic << ");" << endl;
output << " end" << endl;
if (simulation_type == BlockSimulationType::solveBackwardSimple if (simulation_type == BlockSimulationType::solveBackwardSimple
|| simulation_type == BlockSimulationType::solveForwardSimple || simulation_type == BlockSimulationType::solveForwardSimple
@ -381,42 +363,96 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
write_eq_tt(blocks[blk].size); write_eq_tt(blocks[blk].size);
output << " % Jacobian" << endl output << " % Jacobian" << endl
<< " if stochastic_mode" << endl; << " if stochastic_mode" << endl
<< " g1_i = zeros(" << nze_stochastic << ", 1);" << endl
<< " g1_j = zeros(" << nze_stochastic << ", 1);" << endl
<< " g1_v = zeros(" << nze_stochastic << ", 1);" << endl
<< " g1_x_i = zeros(" << nze_exo << ", 1);" << endl
<< " g1_x_j = zeros(" << nze_exo << ", 1);" << endl
<< " g1_x_v = zeros(" << nze_exo << ", 1);" << endl
<< " g1_xd_i = zeros(" << nze_exo_det << ", 1);" << endl
<< " g1_xd_j = zeros(" << nze_exo_det << ", 1);" << endl
<< " g1_xd_v = zeros(" << nze_exo_det << ", 1);" << endl
<< " g1_o_i = zeros(" << nze_other_endo << ", 1);" << endl
<< " g1_o_j = zeros(" << nze_other_endo << ", 1);" << endl
<< " g1_o_v = zeros(" << nze_other_endo << ", 1);" << endl;
ostringstream i_output, j_output, v_output;
int line_counter = 1;
for (const auto &[indices, d] : blocks_derivatives[blk]) for (const auto &[indices, d] : blocks_derivatives[blk])
{ {
auto [eq, var, lag] = indices; auto [eq, var, lag] = indices;
int jacob_col = blocks_jacob_cols_endo[blk].at({ var, lag }); int jacob_col = blocks_jacob_cols_endo[blk].at({ var, lag });
output << " g1(" << eq+1 << ", " << jacob_col+1 << ") = "; i_output << " g1_i(" << line_counter << ")="<< eq+1 << ";" << endl;
d->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); j_output << " g1_j(" << line_counter << ")="<< jacob_col+1 << ";" << endl;
output << ";" << endl; v_output << " g1_v(" << line_counter << ")=";
d->writeOutput(v_output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
v_output << ";" << endl;
line_counter++;
} }
assert(line_counter == nze_stochastic+1);
output << i_output.str() << j_output.str() << v_output.str();
i_output.str("");
j_output.str("");
v_output.str("");
line_counter = 1;
for (const auto &[indices, d] : blocks_derivatives_exo[blk]) for (const auto &[indices, d] : blocks_derivatives_exo[blk])
{ {
auto [eq, var, lag] = indices; auto [eq, var, lag] = indices;
int jacob_col = blocks_jacob_cols_exo[blk].at({ var, lag }); int jacob_col = blocks_jacob_cols_exo[blk].at({ var, lag });
output << " g1_x(" << eq+1 << ", " << jacob_col+1 << ") = "; i_output << " g1_x_i(" << line_counter << ")="<< eq+1 << ";" << endl;
d->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); j_output << " g1_x_j(" << line_counter << ")="<< jacob_col+1 << ";" << endl;
output << ";" << endl; v_output << " g1_x_v(" << line_counter << ")=";
d->writeOutput(v_output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
v_output << ";" << endl;
line_counter++;
} }
assert(line_counter == nze_exo+1);
output << i_output.str() << j_output.str() << v_output.str();
i_output.str("");
j_output.str("");
v_output.str("");
line_counter = 1;
for (const auto &[indices, d] : blocks_derivatives_exo_det[blk]) for (const auto &[indices, d] : blocks_derivatives_exo_det[blk])
{ {
auto [eq, var, lag] = indices; auto [eq, var, lag] = indices;
int jacob_col = blocks_jacob_cols_exo_det[blk].at({ var, lag }); int jacob_col = blocks_jacob_cols_exo_det[blk].at({ var, lag });
output << " g1_xd(" << eq+1 << ", " << jacob_col+1 << ") = "; i_output << " g1_xd_i(" << line_counter << ")="<< eq+1 << ";" << endl;
d->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); j_output << " g1_xd_j(" << line_counter << ")="<< jacob_col+1 << ";" << endl;
output << ";" << endl; v_output << " g1_xd_v(" << line_counter << ")=";
d->writeOutput(v_output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
v_output << ";" << endl;
line_counter++;
} }
assert(line_counter == nze_exo_det+1);
output << i_output.str() << j_output.str() << v_output.str();
i_output.str("");
j_output.str("");
v_output.str("");
line_counter = 1;
for (const auto &[indices, d] : blocks_derivatives_other_endo[blk]) for (const auto &[indices, d] : blocks_derivatives_other_endo[blk])
{ {
auto [eq, var, lag] = indices; auto [eq, var, lag] = indices;
int jacob_col = blocks_jacob_cols_other_endo[blk].at({ var, lag }); int jacob_col = blocks_jacob_cols_other_endo[blk].at({ var, lag });
output << " g1_o(" << eq+1 << ", " << jacob_col+1 << ") = "; i_output << " g1_o_i(" << line_counter << ")="<< eq+1 << ";" << endl;
d->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); j_output << " g1_o_j(" << line_counter << ")="<< jacob_col+1 << ";" << endl;
output << ";" << endl; v_output << " g1_o_v(" << line_counter << ")=";
d->writeOutput(v_output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
v_output << ";" << endl;
line_counter++;
} }
output << " varargout{1}=g1_x;" << endl assert(line_counter == nze_other_endo+1);
<< " varargout{2}=g1_xd;" << endl output << i_output.str() << j_output.str() << v_output.str();
<< " varargout{3}=g1_o;" << endl; i_output.str("");
j_output.str("");
v_output.str("");
output << " g1=sparse(g1_i, g1_j, g1_v, " << block_size << ", " << blocks_jacob_cols_endo[blk].size() << ");" << endl
<< " varargout{1}=sparse(g1_x_i, g1_x_j, g1_x_v, " << block_size << ", " << blocks_jacob_cols_exo[blk].size() << ");" << endl
<< " varargout{2}=sparse(g1_xd_i, g1_xd_j, g1_xd_v, " << block_size << ", " << blocks_jacob_cols_exo_det[blk].size() << ");" << endl
<< " varargout{3}=sparse(g1_o_i, g1_o_j, g1_o_v, " << block_size << ", " << blocks_jacob_cols_other_endo[blk].size() << ");" << endl;
switch (simulation_type) switch (simulation_type)
{ {
@ -428,43 +464,56 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
case BlockSimulationType::solveForwardSimple: case BlockSimulationType::solveForwardSimple:
case BlockSimulationType::solveBackwardComplete: case BlockSimulationType::solveBackwardComplete:
case BlockSimulationType::solveForwardComplete: case BlockSimulationType::solveForwardComplete:
output << " else" << endl; output << " else" << endl
<< " g1_i = zeros(" << nze_deterministic << ", 1);" << endl
<< " g1_j = zeros(" << nze_deterministic << ", 1);" << endl
<< " g1_v = zeros(" << nze_deterministic << ", 1);" << endl;
line_counter = 1;
for (const auto &[indices, d] : blocks_derivatives[blk]) for (const auto &[indices, d] : blocks_derivatives[blk])
{ {
auto [eq, var, lag] = indices; auto [eq, var, lag] = indices;
if (lag == 0) if (lag == 0)
{ {
output << " g1(" << eq+1 << ", " << var+1-block_recursive_size << ") = "; i_output << " g1_i(" << line_counter << ")="<< eq+1 << ";" << endl;
d->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); j_output << " g1_j(" << line_counter << ")="<< var+1-block_recursive_size << ";" << endl;
output << ";" << endl; v_output << " g1_v(" << line_counter << ")=";
d->writeOutput(v_output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
v_output << ";" << endl;
line_counter++;
} }
} }
output << " end" << endl; assert(line_counter == nze_deterministic+1);
output << i_output.str() << j_output.str() << v_output.str()
<< " g1=sparse(g1_i, g1_j, g1_v, " << block_mfs_size
<< ", " << block_mfs_size << ");" << endl
<< " end" << endl;
break; break;
case BlockSimulationType::solveTwoBoundariesSimple: case BlockSimulationType::solveTwoBoundariesSimple:
case BlockSimulationType::solveTwoBoundariesComplete: case BlockSimulationType::solveTwoBoundariesComplete:
output << " else" << endl; output << " else" << endl
<< " g1_i = zeros(" << nze_deterministic << ", 1);" << endl
<< " g1_j = zeros(" << nze_deterministic << ", 1);" << endl
<< " g1_v = zeros(" << nze_deterministic << ", 1);" << endl;
line_counter = 1;
for (const auto &[indices, d] : blocks_derivatives[blk]) for (const auto &[indices, d] : blocks_derivatives[blk])
{ {
auto [eq, var, lag] = indices; auto [eq, var, lag] = indices;
assert(lag >= -1 && lag <= 1); assert(lag >= -1 && lag <= 1);
if (eq >= block_recursive_size && var >= block_recursive_size) if (eq >= block_recursive_size && var >= block_recursive_size)
{ {
if (lag == 0) i_output << " g1_i(" << line_counter << ")="<< eq+1-block_recursive_size << ";" << endl;
output << " g1(" << eq+1-block_recursive_size << ", " j_output << " g1_j(" << line_counter << ")="<< var+1-block_recursive_size+block_mfs_size*(lag+1) << ";" << endl;
<< var+1-block_recursive_size+block_mfs_size << ") = "; v_output << " g1_v(" << line_counter << ")=";
else if (lag == 1) d->writeOutput(v_output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
output << " g1(" << eq+1-block_recursive_size << ", " v_output << ";" << endl;
<< var+1-block_recursive_size+2*block_mfs_size << ") = "; line_counter++;
else
output << " g1(" << eq+1-block_recursive_size << ", "
<< var+1-block_recursive_size << ") = ";
d->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
output << ";" << endl;
} }
} }
assert(line_counter == nze_deterministic+1);
output << " end" << endl; output << i_output.str() << j_output.str() << v_output.str()
<< " g1=sparse(g1_i, g1_j, g1_v, " << block_mfs_size
<< ", " << 3*block_mfs_size << ");" << endl
<< " end" << endl;
break; break;
default: default:
break; break;

34
src/StaticModel.cc
View File

@ -161,11 +161,6 @@ StaticModel::writeStaticPerBlockMFiles(const string &basename) const
<< BlockSim(simulation_type) << " //" << endl << BlockSim(simulation_type) << " //" << endl
<< " % ////////////////////////////////////////////////////////////////////////" << endl; << " % ////////////////////////////////////////////////////////////////////////" << endl;
if (simulation_type != BlockSimulationType::evaluateBackward
&& simulation_type != BlockSimulationType::evaluateForward)
output << " g1=spalloc(" << blocks[blk].mfs_size << "," << blocks[blk].mfs_size
<< "," << blocks_derivatives[blk].size() << ");" << endl;
if (simulation_type != BlockSimulationType::evaluateBackward if (simulation_type != BlockSimulationType::evaluateBackward
&& simulation_type != BlockSimulationType::evaluateForward) && simulation_type != BlockSimulationType::evaluateForward)
output << " residual=zeros(" << blocks[blk].mfs_size << ",1);" << endl; output << " residual=zeros(" << blocks[blk].mfs_size << ",1);" << endl;
@ -236,25 +231,32 @@ StaticModel::writeStaticPerBlockMFiles(const string &basename) const
} }
} }
// The Jacobian if we have to solve the block // The Jacobian if we have to solve the block
switch (simulation_type) if (simulation_type == BlockSimulationType::solveBackwardSimple
|| simulation_type == BlockSimulationType::solveForwardSimple
|| simulation_type == BlockSimulationType::solveBackwardComplete
|| simulation_type == BlockSimulationType::solveForwardComplete)
{ {
case BlockSimulationType::solveBackwardSimple:
case BlockSimulationType::solveForwardSimple:
case BlockSimulationType::solveBackwardComplete:
case BlockSimulationType::solveForwardComplete:
// Write temporary terms for derivatives // Write temporary terms for derivatives
write_eq_tt(blocks[blk].size); write_eq_tt(blocks[blk].size);
output << " g1_i=zeros(" << blocks_derivatives[blk].size() << ", 1);" << endl
<< " g1_j=zeros(" << blocks_derivatives[blk].size() << ", 1);" << endl
<< " g1_v=zeros(" << blocks_derivatives[blk].size() << ", 1);" << endl;
ostringstream i_output, j_output, v_output;
int line_counter = 1;
for (const auto &[indices, d] : blocks_derivatives[blk]) for (const auto &[indices, d] : blocks_derivatives[blk])
{ {
auto [eq, var, ignore] = indices; auto [eq, var, ignore] = indices;
output << " g1(" << eq+1-block_recursive_size << ", " << var+1-block_recursive_size << ") = "; i_output << " g1_i(" << line_counter << ")=" << eq+1-block_recursive_size << ";" << endl;
d->writeOutput(output, local_output_type, temporary_terms, blocks_temporary_terms_idxs); j_output << " g1_j(" << line_counter << ")=" << var+1-block_recursive_size << ";";
output << ";" << endl; v_output << " g1_v(" << line_counter << ")=";
d->writeOutput(v_output, local_output_type, temporary_terms, blocks_temporary_terms_idxs);
v_output << ";" << endl;
line_counter++;
} }
break; output << i_output.str() << j_output.str() << v_output.str()
default: << " g1=sparse(g1_i, g1_j, g1_v, " << blocks[blk].mfs_size << "," << blocks[blk].mfs_size << ");" << endl;
break;
} }
output << "end" << endl; output << "end" << endl;
output.close(); output.close();