diff --git a/arbor/backends/gpu/math_cu.hpp b/arbor/backends/gpu/math_cu.hpp
index 13abd7eaf2a8bd95df9f62dec43c11d626499eef..5937ec86aeb3eb2bd67f0ed3f3f2681f92a056a9 100644
--- a/arbor/backends/gpu/math_cu.hpp
+++ b/arbor/backends/gpu/math_cu.hpp
@@ -1,4 +1,5 @@
#pragma once
+#include <cfloat>
// Implementations of mathematical operations required
// by generated CUDA mechanisms.
@@ -6,6 +7,14 @@
namespace arb {
namespace gpu {
+__device__
+inline double safeinv(double x) {
+ if (1.0+x == 1.0) {
+ return 1/DBL_EPSILON;
+ }
+ return 1/x;
+}
+
__device__
inline double exprelr(double x) {
if (1.0+x == 1.0) {
diff --git a/arbor/include/arbor/math.hpp b/arbor/include/arbor/math.hpp
index f8fba17e39868ade737002f55ffa4a8b2d0f1c98..08a5ff4e73a200c72c0608eea1c2ccaf99ae4481 100644
--- a/arbor/include/arbor/math.hpp
+++ b/arbor/include/arbor/math.hpp
@@ -116,6 +116,14 @@ C round_up(T v, U b) {
return v-m+signum(m)*impl::abs_if_signed(b, Signed{});
}
+// Returns 1/x if x != 0; 0 otherwise
+template <typename T>
+inline
+T safeinv(T x) {
+ // If abs(x) is less than epsilon return epsilon, else calculate the result directly.
+ return (T(1)==T(1)+x)? 1/std::numeric_limits<T>::epsilon(): 1/x;
+}
+
// Value of x/(exp(x)-1) with care taken to handle x=0 case
template <typename T>
inline
diff --git a/modcc/expression.cpp b/modcc/expression.cpp
index 18e5a53178496952419036ef47a81314d47fdca7..9ba6c58ea1c1a3716080163456688f8792db1981 100644
--- a/modcc/expression.cpp
+++ b/modcc/expression.cpp
@@ -996,6 +996,9 @@ void LogUnaryExpression::accept(Visitor *v) {
void AbsUnaryExpression::accept(Visitor *v) {
v->visit(this);
}
+void SafeInvUnaryExpression::accept(Visitor *v) {
+ v->visit(this);
+}
void ExprelrUnaryExpression::accept(Visitor *v) {
v->visit(this);
}
@@ -1077,6 +1080,8 @@ expression_ptr unary_expression( Location loc,
return make_expression<AbsUnaryExpression>(loc, std::move(e));
case tok::exprelr :
return make_expression<ExprelrUnaryExpression>(loc, std::move(e));
+ case tok::safeinv :
+ return make_expression<SafeInvUnaryExpression>(loc, std::move(e));
default :
std::cerr << yellow(token_string(op))
<< " is not a valid unary operator"
diff --git a/modcc/expression.hpp b/modcc/expression.hpp
index fe3eab9fb438ffab450c0292f82cd194a7a1e6b9..1da38c35247e7cb3fc0f2cf5f2c49cfa714e7baa 100644
--- a/modcc/expression.hpp
+++ b/modcc/expression.hpp
@@ -1226,6 +1226,15 @@ public:
void accept(Visitor *v) override;
};
+class SafeInvUnaryExpression : public UnaryExpression {
+public:
+ SafeInvUnaryExpression(Location loc, expression_ptr e)
+ : UnaryExpression(loc, tok::safeinv, std::move(e))
+ {}
+
+ void accept(Visitor *v) override;
+};
+
// exprel reciprocal unary expression,
// i.e. x/(exp(x)-1)=x/expm1(x) with exprelr(0)=1
class ExprelrUnaryExpression : public UnaryExpression {
diff --git a/modcc/parser.cpp b/modcc/parser.cpp
index efeb1cac10905c4a883f0b8c6b056293dc1c9ee9..bb59db9ebcd1fb3b01801cc56d96c56a076ab0bb 100644
--- a/modcc/parser.cpp
+++ b/modcc/parser.cpp
@@ -1346,6 +1346,7 @@ expression_ptr Parser::parse_unaryop() {
case tok::cos :
case tok::log :
case tok::abs :
+ case tok::safeinv:
case tok::exprelr:
get_token(); // consume operator (exp, sin, cos or log)
if(token_.type!=tok::lparen) {
diff --git a/modcc/printer/cexpr_emit.cpp b/modcc/printer/cexpr_emit.cpp
index a3be027494849905622f689e83c34b287061f619..0c88fb766e85c8714189198c2ada86b4c50e885b 100644
--- a/modcc/printer/cexpr_emit.cpp
+++ b/modcc/printer/cexpr_emit.cpp
@@ -53,7 +53,8 @@ void CExprEmitter::visit(UnaryExpression* e) {
{tok::sin, "sin"},
{tok::log, "log"},
{tok::abs, "abs"},
- {tok::exprelr, "exprelr"}
+ {tok::exprelr, "exprelr"},
+ {tok::safeinv, "safeinv"}
};
if (!unaryop_tbl.count(e->op())) {
diff --git a/modcc/printer/cprinter.cpp b/modcc/printer/cprinter.cpp
index dc0b81a52b0781746c89037f669d99dcbaa291d4..c42a16e80c5fec5804c66c5c4d61a7d445ec5700 100644
--- a/modcc/printer/cprinter.cpp
+++ b/modcc/printer/cprinter.cpp
@@ -154,6 +154,7 @@ std::string emit_cpp_source(const Module& module_, const printer_options& opt) {
"using size_type = base::size_type;\n"
"using index_type = base::index_type;\n"
"using ::arb::math::exprelr;\n"
+ "using ::arb::math::safeinv;\n"
"using ::std::abs;\n"
"using ::std::cos;\n"
"using ::std::exp;\n"
@@ -190,6 +191,11 @@ std::string emit_cpp_source(const Module& module_, const printer_options& opt) {
out <<
"using simd_value = S::simd<::arb::fvm_value_type, simd_width_, " << abi << ">;\n"
"using simd_index = S::simd<::arb::fvm_index_type, simd_width_, " << abi << ">;\n"
+ "\n"
+ "inline simd_value safeinv(simd_value x) {\n"
+ " S::where(x+1==1, x) = DBL_EPSILON;\n"
+ " return 1/x;\n"
+ "}\n"
"\n";
}
diff --git a/modcc/printer/cudaprinter.cpp b/modcc/printer/cudaprinter.cpp
index 783eb9f6573a0bb0cc8173f085b014d5706e6ac5..3aac1146c523f000d8fedd82c7ee6a80a0ec3d3e 100644
--- a/modcc/printer/cudaprinter.cpp
+++ b/modcc/printer/cudaprinter.cpp
@@ -253,6 +253,7 @@ std::string emit_cuda_cu_source(const Module& module_, const printer_options& op
out << "namespace {\n\n"; // place inside an anonymous namespace
out << "using ::arb::gpu::exprelr;\n";
+ out << "using ::arb::gpu::safeinv;\n";
out << "using ::arb::gpu::min;\n";
out << "using ::arb::gpu::max;\n\n";
diff --git a/modcc/solvers.cpp b/modcc/solvers.cpp
index 19e3700ae872df4eea6bd53a2bcbd715df869150..cd792ab092812af40ca33b2961d398611f61c302 100644
--- a/modcc/solvers.cpp
+++ b/modcc/solvers.cpp
@@ -4,7 +4,6 @@
#include <string>
#include <vector>
-#include "astmanip.hpp"
#include "expression.hpp"
#include "parser.hpp"
#include "solvers.hpp"
@@ -171,6 +170,84 @@ static expression_ptr as_expression(symge::symbol_term_diff diff) {
}
}
+std::vector<local_assignment> SystemSolver::generate_row_updates(scope_ptr scope, std::vector<symge::symbol> row_sym) {
+ std::vector<local_assignment> S_;
+ for (const auto& s: row_sym) {
+ if (primitive(s)) continue;
+
+ auto expr = as_expression(definition(s));
+ auto local_t_term = make_unique_local_assign(scope, expr.get(), "t_");
+ auto t_ = local_t_term.id->is_identifier()->spelling();
+ symtbl_.name(s, t_);
+ S_.push_back(std::move(local_t_term));
+ }
+ return S_;
+}
+
+local_assignment SystemSolver::generate_normalizing_term(scope_ptr scope, std::vector<symge::symbol> row_sym) {
+ // Get the max element in the row
+ expression_ptr max;
+ for (auto &s: row_sym) {
+ auto elem = make_expression<IdentifierExpression>(Location(), symtbl_.name(s));
+ auto abs = make_expression<AbsUnaryExpression>(elem->location(), elem->is_identifier()->clone());
+ if (!max) {
+ max = std::move(abs);
+ } else {
+ max = make_expression<MaxBinaryExpression>(elem->location(), max->clone(), std::move(abs));
+ }
+ }
+ // Safely invert max
+ auto inv = make_expression<SafeInvUnaryExpression>(max->location(), std::move(max));
+
+ // Create a symbol for inv
+ auto local_inv_term = make_unique_local_assign(scope, inv.get(), "inv_");
+
+ return local_inv_term;
+}
+
+std::vector<expression_ptr> SystemSolver::generate_normalizing_assignments(expression_ptr normalizer, std::vector<symge::symbol> row_sym) {
+ std::vector<expression_ptr> S_;
+ for (auto &s: row_sym) {
+ auto elem = make_expression<IdentifierExpression>(Location(), symtbl_.name(s));
+ auto ratio = make_expression<MulBinaryExpression>(elem->location(), elem->clone(), normalizer->clone());
+ auto assign = make_expression<AssignmentExpression>(elem->location(), elem->clone(), std::move(ratio));
+ S_.push_back(std::move(assign));
+ }
+ return S_;
+}
+
+std::vector<expression_ptr> SystemSolver::generate_solution_assignments(std::vector<std::string> lhs_vars) {
+ std::vector<expression_ptr> U_;
+
+ // State variable updates given by rhs/diagonal for reduced matrix.
+ Location loc;
+ auto nrow = A_.nrow();
+ for (unsigned i = 0; i < nrow; ++i) {
+ const symge::sym_row& row = A_[i];
+ unsigned rhs_col = A_.augcol();
+ unsigned lhs_col = -1;
+ for (unsigned r = 0; r < A_.nrow(); ++r) {
+ if (row[r]) {
+ lhs_col = r;
+ break;
+ }
+ }
+
+ if (lhs_col==unsigned(-1)) {
+ throw std::logic_error("zero row in matrix solver");
+ }
+ auto expr = make_expression<AssignmentExpression>(loc,
+ make_expression<IdentifierExpression>(loc, lhs_vars[lhs_col]),
+ make_expression<DivBinaryExpression>(loc,
+ make_expression<IdentifierExpression>(loc, symge::name(A_[i][rhs_col])),
+ make_expression<IdentifierExpression>(loc, symge::name(A_[i][lhs_col]))));
+
+ U_.push_back(std::move(expr));
+ }
+ return U_;
+}
+
+
void SparseSolverVisitor::visit(BlockExpression* e) {
// Do a first pass to extract variables comprising ODE system
// lhs; can't really trust 'STATE' block.
@@ -212,9 +289,8 @@ void SparseSolverVisitor::visit(CompartmentExpression *e) {
}
void SparseSolverVisitor::visit(AssignmentExpression *e) {
- if (A_.empty()) {
- unsigned n = dvars_.size();
- A_ = symge::sym_matrix(n, n);
+ if (system_.empty()) {
+ system_.create_square_matrix(dvars_.size());
}
auto loc = e->location();
@@ -237,7 +313,7 @@ void SparseSolverVisitor::visit(AssignmentExpression *e) {
return;
}
- if (conserve_ && !A_[deq_index_].empty()) {
+ if (conserve_ && !system_.empty_row(deq_index_)) {
deq_index_++;
return;
}
@@ -299,15 +375,14 @@ void SparseSolverVisitor::visit(AssignmentExpression *e) {
statements_.push_back(std::move(local_a_term.local_decl));
statements_.push_back(std::move(local_a_term.assignment));
- A_[deq_index_].push_back({j, symtbl_.define(a_)});
+ system_.add_entry({deq_index_, j}, a_);
}
++deq_index_;
}
void SparseSolverVisitor::visit(ConserveExpression *e) {
- if (A_.empty()) {
- unsigned n = dvars_.size();
- A_ = symge::sym_matrix(n, n);
+ if (system_.empty()) {
+ system_.create_square_matrix(dvars_.size());
}
conserve_ = true;
@@ -334,7 +409,7 @@ void SparseSolverVisitor::visit(ConserveExpression *e) {
}
// Replace that row with the conserve statement
- A_[row_idx].clear();
+ system_.clear_row(row_idx);
for (unsigned j = 0; j < dvars_.size(); ++j) {
auto state = dvars_[j];
@@ -355,7 +430,7 @@ void SparseSolverVisitor::visit(ConserveExpression *e) {
statements_.push_back(std::move(local_a_term.local_decl));
statements_.push_back(std::move(local_a_term.assignment));
- A_[row_idx].push_back({j, symtbl_.define(a_)});
+ system_.add_entry({(unsigned)row_idx, j}, a_);
}
}
@@ -378,63 +453,46 @@ void SparseSolverVisitor::finalize() {
error({"Conserve statement(s) missing in steady-state solver", {}});
}
- std::vector<symge::symbol> rhs;
+ std::vector<std::string> rhs;
for (const auto& var: dvars_) {
auto v = solve_variant_ == solverVariant::steadystate? steadystate_rhs_ : var;
- rhs.push_back(symtbl_.define(v));
+ rhs.push_back(v);
}
if (conserve_) {
for (unsigned i = 0; i < conserve_idx_.size(); ++i) {
- rhs[conserve_idx_[i]] = symtbl_.define(conserve_rhs_[i]);
+ rhs[conserve_idx_[i]] = conserve_rhs_[i];
}
}
- A_.augment(rhs);
-
- symge::gj_reduce(A_, symtbl_);
+ system_.augment(rhs);
- // Create and assign intermediate variables.
- for (unsigned i = 0; i<symtbl_.size(); ++i) {
- symge::symbol s = symtbl_[i];
+ // Reduce the system
+ auto row_symbols = system_.reduce();
- if (primitive(s)) continue;
-
- auto expr = as_expression(definition(s));
- auto local_t_term = make_unique_local_assign(block_scope_, expr.get(), "t_");
- auto t_ = local_t_term.id->is_identifier()->spelling();
- symtbl_.name(s, t_);
-
- statements_.push_back(std::move(local_t_term.local_decl));
- statements_.push_back(std::move(local_t_term.assignment));
- }
-
- // State variable updates given by rhs/diagonal for reduced matrix.
- Location loc;
- auto nrow = A_.nrow();
- for (unsigned i = 0; i<nrow; ++i) {
- const symge::sym_row& row = A_[i];
- unsigned rhs_col = A_.augcol();
- unsigned lhs_col = -1;
- for (unsigned r = 0; r<A_.nrow(); ++r) {
- if (row[r]) {
- lhs_col = r;
- break;
- }
+ // Row by row:
+ // Generate entries of the system and declare and assign as local variables
+ // Generate normalizing terms and normalize the row
+ for (auto& row: row_symbols) {
+ auto entries = system_.generate_row_updates(block_scope_, row);
+ for (auto& l: entries) {
+ statements_.push_back(std::move(l.local_decl));
+ statements_.push_back(std::move(l.assignment));
}
- if (lhs_col==unsigned(-1)) {
- throw std::logic_error("zero row in sparse solver matrix");
- }
+ // If size of system > 5 normalize the row updates
+ if (system_.size() > 5) {
+ auto norm_term = system_.generate_normalizing_term(block_scope_, row);
+ auto norm_assigns = system_.generate_normalizing_assignments(norm_term.id->clone(), row);
- auto expr =
- make_expression<AssignmentExpression>(loc,
- make_expression<IdentifierExpression>(loc, dvars_[lhs_col]),
- make_expression<DivBinaryExpression>(loc,
- make_expression<IdentifierExpression>(loc, symge::name(A_[i][rhs_col])),
- make_expression<IdentifierExpression>(loc, symge::name(A_[i][lhs_col]))));
-
- statements_.push_back(std::move(expr));
+ statements_.push_back(std::move(norm_term.local_decl));
+ statements_.push_back(std::move(norm_term.assignment));
+ std::move(std::begin(norm_assigns), std::end(norm_assigns), std::back_inserter(statements_));
+ }
}
+ // Update the state variables
+ auto updates = system_.generate_solution_assignments(dvars_);
+ std::move(std::begin(updates), std::end(updates), std::back_inserter(statements_));
+
BlockRewriterBase::finalize();
}
@@ -451,9 +509,8 @@ void LinearSolverVisitor::visit(LinearExpression *e) {
auto loc = e->location();
scope_ptr scope = e->scope();
- if (A_.empty()) {
- unsigned n = dvars_.size();
- A_ = symge::sym_matrix(n, n);
+ if (system_.empty()) {
+ system_.create_square_matrix(dvars_.size());
}
linear_test_result r = linear_test(e->lhs(), dvars_);
@@ -473,60 +530,43 @@ void LinearSolverVisitor::visit(LinearExpression *e) {
auto a_ = expr->is_identifier()->spelling();
- A_[deq_index_].push_back({j, symtbl_.define(a_)});
+ system_.add_entry({deq_index_, j}, a_);
}
- rhs_.push_back(symtbl_.define(e->rhs()->is_identifier()->spelling()));
+ rhs_.push_back(e->rhs()->is_identifier()->spelling());
++deq_index_;
}
void LinearSolverVisitor::finalize() {
- A_.augment(rhs_);
-
- symge::gj_reduce(A_, symtbl_);
-
- // Create and assign intermediate variables.
- for (unsigned i = 0; i<symtbl_.size(); ++i) {
- symge::symbol s = symtbl_[i];
+ system_.augment(rhs_);
- if (primitive(s)) continue;
-
- auto expr = as_expression(definition(s));
- auto local_t_term = make_unique_local_assign(block_scope_, expr.get(), "t_");
- auto t_ = local_t_term.id->is_identifier()->spelling();
- symtbl_.name(s, t_);
-
- statements_.push_back(std::move(local_t_term.local_decl));
- statements_.push_back(std::move(local_t_term.assignment));
- }
-
- // State variable updates given by rhs/diagonal for reduced matrix.
- Location loc;
- auto nrow = A_.nrow();
- for (unsigned i = 0; i < nrow; ++i) {
- const symge::sym_row& row = A_[i];
- unsigned rhs = A_.augcol();
- unsigned lhs = -1;
- for (unsigned r = 0; r < A_.nrow(); ++r) {
- if (row[r]) {
- lhs = r;
- break;
- }
- }
+ // Reduce the system
+ auto row_symbols = system_.reduce();
- if (lhs==unsigned(-1)) {
- throw std::logic_error("zero row in linear solver matrix");
+ // Row by row:
+ // Generate entries of the system and declare and assign as local variables
+ // Generate normalizing terms and normalize the row
+ for (auto& row: row_symbols) {
+ auto entries = system_.generate_row_updates(block_scope_, row);
+ for (auto& l: entries) {
+ statements_.push_back(std::move(l.local_decl));
+ statements_.push_back(std::move(l.assignment));
}
- auto expr =
- make_expression<AssignmentExpression>(loc,
- make_expression<IdentifierExpression>(loc, dvars_[lhs]),
- make_expression<DivBinaryExpression>(loc,
- make_expression<IdentifierExpression>(loc, symge::name(A_[i][rhs])),
- make_expression<IdentifierExpression>(loc, symge::name(A_[i][lhs]))));
+ // If size of system > 5 normalize the row updates
+ if (system_.size() > 5) {
+ auto norm_term = system_.generate_normalizing_term(block_scope_, row);
+ auto norm_assigns = system_.generate_normalizing_assignments(norm_term.id->clone(), row);
- statements_.push_back(std::move(expr));
+ statements_.push_back(std::move(norm_term.local_decl));
+ statements_.push_back(std::move(norm_term.assignment));
+ std::move(std::begin(norm_assigns), std::end(norm_assigns), std::back_inserter(statements_));
+ }
}
+ // Update the state variables
+ auto updates = system_.generate_solution_assignments(dvars_);
+ std::move(std::begin(updates), std::end(updates), std::back_inserter(statements_));
+
BlockRewriterBase::finalize();
}
@@ -591,9 +631,8 @@ void SparseNonlinearSolverVisitor::visit(CompartmentExpression *e) {
}
void SparseNonlinearSolverVisitor::visit(AssignmentExpression *e) {
- if (A_.empty()) {
- unsigned n = dvars_.size();
- A_ = symge::sym_matrix(n, n);
+ if (system_.empty()) {
+ system_.create_square_matrix(dvars_.size());
}
auto loc = e->location();
@@ -705,7 +744,7 @@ void SparseNonlinearSolverVisitor::visit(AssignmentExpression *e) {
statements_.push_back(std::move(local_j_term.local_decl));
J_.push_back(std::move(local_j_term.assignment));
- A_[deq_index_].push_back({j, symtbl_.define(j_)});
+ system_.add_entry({deq_index_, j}, j_);
}
++deq_index_;
}
@@ -713,60 +752,49 @@ void SparseNonlinearSolverVisitor::visit(AssignmentExpression *e) {
void SparseNonlinearSolverVisitor::finalize() {
// Create rhs of A_
- std::vector<symge::symbol> rhs;
+ std::vector<std::string> rhs;
for (const auto& var: F_) {
auto id = var->is_assignment()->lhs()->is_identifier()->spelling();
- rhs.push_back(symtbl_.define(id));
+ rhs.push_back(id);
}
- A_.augment(rhs);
- symge::gj_reduce(A_, symtbl_);
+ system_.augment(rhs);
- // Create and assign intermediate variables.
- // Save gaussian elimination solution statements in S_
+ // Reduce the system
+ auto row_symbols = system_.reduce();
+
+ // Row by row:
+ // Generate entries of the system and declare and assign as local variables
+ // Generate normalizing terms and normalize the row
std::vector<expression_ptr> S_;
- for (unsigned i = 0; i < symtbl_.size(); ++i) {
- symge::symbol s = symtbl_[i];
+ for (auto& row: row_symbols) {
+ auto entries = system_.generate_row_updates(block_scope_, row);
+ for (auto& l: entries) {
+ statements_.push_back(std::move(l.local_decl));
+ S_.push_back(std::move(l.assignment));
+ }
- if (primitive(s)) continue;
+ // If size of system > 5 normalize the row updates
+ if (system_.size() > 5) {
+ auto norm_term = system_.generate_normalizing_term(block_scope_, row);
+ auto norm_assigns = system_.generate_normalizing_assignments(norm_term.id->clone(), row);
- auto expr = as_expression(definition(s));
- auto local_t_term = make_unique_local_assign(block_scope_, expr.get(), "s_");
- auto t_ = local_t_term.id->is_identifier()->spelling();
- symtbl_.name(s, t_);
-
- statements_.push_back(std::move(local_t_term.local_decl));
- S_.push_back(std::move(local_t_term.assignment));
+ statements_.push_back(std::move(norm_term.local_decl));
+ S_.push_back(std::move(norm_term.assignment));
+ std::move(std::begin(norm_assigns), std::end(norm_assigns), std::back_inserter(S_));
+ }
}
+ // Update the state variables
+ auto U_ = system_.generate_solution_assignments(dvar_temp_);
+
// Create the statements that update the temporary state variables
// (dvar_temp) after a Newton's iteration and save them in U_
- std::vector<expression_ptr> U_;
-
- Location loc;
- for (unsigned i = 0; i < A_.nrow(); ++i) {
- const symge::sym_row &row = A_[i];
- unsigned rhs_col = A_.augcol();
- int lhs_col = -1;
- for (unsigned r = 0; r < A_.nrow(); r++) {
- if (row[r]) {
- lhs_col = r;
- break;
- }
- }
-
- if (lhs_col == -1) {
- error({"Something went wrong in solving the sparse matrix", loc});
- return;
- }
-
- auto var = make_expression<IdentifierExpression>(loc, dvar_temp_[lhs_col]);
- auto expr = make_expression<AssignmentExpression>(loc, var->clone(),
- make_expression<SubBinaryExpression>(loc, var->clone(),
- make_expression<DivBinaryExpression>(loc, make_expression<IdentifierExpression>(loc, symge::name(A_[i][rhs_col])),
- make_expression<IdentifierExpression>(loc, symge::name(A_[i][lhs_col])))));
-
- U_.push_back(std::move(expr));
+ for (auto& u: U_) {
+ auto lhs = u->is_assignment()->lhs();
+ auto rhs = u->is_assignment()->rhs();
+ u = make_expression<AssignmentExpression>(u->location(), lhs->clone(),
+ make_expression<SubBinaryExpression>(u->location(), lhs->clone(), rhs->clone()));
}
// Do 3 Newton iterations
@@ -786,6 +814,7 @@ void SparseNonlinearSolverVisitor::finalize() {
}
}
+ Location loc;
// Finally update the state variables, dvars_ = dvar_temp_
for (unsigned i = 0; i < dvars_.size(); ++i) {
auto expr = make_expression<AssignmentExpression>(loc,
diff --git a/modcc/solvers.hpp b/modcc/solvers.hpp
index b375edc9b88f93dfb98c50236b6c006b16783b13..718c0900ffcf3f4f51a6ad972c4b5337ceddcf7e 100644
--- a/modcc/solvers.hpp
+++ b/modcc/solvers.hpp
@@ -7,6 +7,7 @@
#include <string>
#include <vector>
+#include "astmanip.hpp"
#include "expression.hpp"
#include "symdiff.hpp"
#include "symge.hpp"
@@ -66,9 +67,75 @@ public:
virtual void visit(AssignmentExpression *e) override;
};
+class SystemSolver {
+protected:
+ // Symbolic matrix for backwards Euler step.
+ symge::sym_matrix A_;
+
+ // 'Symbol table' for initial variables.
+ symge::symbol_table symtbl_;
+
+public:
+ struct system_loc {
+ unsigned row, col;
+ };
+
+ explicit SystemSolver() {}
+
+ void reset() {
+ A_.clear();
+ symtbl_.clear();
+ }
+ unsigned size() const {
+ return A_.size();
+ }
+ bool empty() const {
+ return A_.empty();
+ }
+ bool empty_row(unsigned i) const {
+ return A_[i].empty();
+ }
+ void clear_row(unsigned i) {
+ A_[i].clear();
+ }
+ void create_square_matrix(unsigned n) {
+ A_ = symge::sym_matrix(n, n);
+ }
+ void add_entry(system_loc loc, std::string name) {
+ A_[loc.row].push_back({loc.col, symtbl_.define(name)});
+ }
+ void augment(std::vector<std::string> rhs) {
+ std::vector<symge::symbol> rhs_sym;
+ for (unsigned r = 0; r < rhs.size(); ++r) {
+ rhs_sym.push_back(symtbl_.define(rhs[r]));
+ }
+ A_.augment(rhs_sym);
+ }
+
+ // Returns a vector of rows of symbols
+ // Needed for normalization
+ std::vector<std::vector<symge::symbol>> reduce() {
+ return symge::gj_reduce(A_, symtbl_);
+ }
+
+ // Returns a vector of local assignments for row updates during system reduction
+ std::vector<local_assignment> generate_row_updates(scope_ptr scope, std::vector<symge::symbol> row_sym);
+
+ // Given a row of symbols, generates local assignment for a normalizing term
+ local_assignment generate_normalizing_term(scope_ptr scope, std::vector<symge::symbol> row_sym);
+
+ // Given a row of symbols, generates expressions normalizing row updates
+ std::vector<expression_ptr> generate_normalizing_assignments(expression_ptr normalizer, std::vector<symge::symbol> row_sym);
+
+ // Returns solution assignment of lhs_vars
+ std::vector<expression_ptr> generate_solution_assignments(std::vector<std::string> lhs_vars);
+
+};
+
class SparseSolverVisitor : public SolverVisitorBase {
protected:
solverVariant solve_variant_;
+
// 'Current' differential equation is for variable with this
// index in `dvars`.
unsigned deq_index_ = 0;
@@ -77,12 +144,6 @@ protected:
// calculations.
substitute_map local_expr_;
- // Symbolic matrix for backwards Euler step.
- symge::sym_matrix A_;
-
- // 'Symbol table' for symbolic manipulation.
- symge::symbol_table symtbl_;
-
// Flag to indicate whether conserve statements are part of the system
bool conserve_ = false;
@@ -95,6 +156,10 @@ protected:
// rhs of steadstate
std::string steadystate_rhs_;
+
+ // System Solver helper
+ SystemSolver system_;
+
public:
using SolverVisitorBase::visit;
@@ -110,13 +175,12 @@ public:
virtual void reset() override {
deq_index_ = 0;
local_expr_.clear();
- A_.clear();
- symtbl_.clear();
conserve_ = false;
scale_factor_.clear();
conserve_rhs_.clear();
conserve_idx_.clear();
steadystate_rhs_.clear();
+ system_.reset();
SolverVisitorBase::reset();
}
};
@@ -139,15 +203,12 @@ protected:
std::vector<std::string> dvar_temp_;
std::vector<std::string> dvar_init_;
- // Symbolic matrix for backwards Euler step.
- symge::sym_matrix A_;
-
- // 'Symbol table' for symbolic manipulation.
- symge::symbol_table symtbl_;
-
// State variable multiplier/divider
std::vector<expression_ptr> scale_factor_;
+ // System Solver helper
+ SystemSolver system_;
+
public:
using SolverVisitorBase::visit;
@@ -164,8 +225,8 @@ public:
local_expr_.clear();
F_.clear();
J_.clear();
- symtbl_.clear();
scale_factor_.clear();
+ system_.reset();
SolverVisitorBase::reset();
}
};
@@ -180,14 +241,11 @@ protected:
// calculations.
substitute_map local_expr_;
- // Symbolic matrix for backwards Euler step.
- symge::sym_matrix A_;
-
- // RHS
- std::vector<symge::symbol> rhs_;
+ // Stores the rhs symbols of the linear system
+ std::vector<std::string> rhs_;
- // 'Symbol table' for symbolic manipulation.
- symge::symbol_table symtbl_;
+ // System Solver helper
+ SystemSolver system_;
public:
using SolverVisitorBase::visit;
@@ -204,9 +262,8 @@ public:
virtual void reset() override {
deq_index_ = 0;
local_expr_.clear();
- A_.clear();
rhs_.clear();
- symtbl_.clear();
+ system_.reset();
SolverVisitorBase::reset();
}
};
diff --git a/modcc/symge.cpp b/modcc/symge.cpp
index 51a9008cd9b774eed6eebe1baceb2980f66b5aaf..704936c3af718a6524a1fa01a5f247aa2c1ea053 100644
--- a/modcc/symge.cpp
+++ b/modcc/symge.cpp
@@ -80,7 +80,10 @@ double estimate_cost(const sym_matrix& A, pivot p) {
// The reduction is division-free: the result will have non-zero terms
// that are symbols that are either primitive, or defined (in the symbol
// table) as products or differences of products of other symbols.
-void gj_reduce(sym_matrix& A, symbol_table& table) {
+// Returns a vector of vectors of symbols, partitioned by row of the matrix
+std::vector<std::vector<symge::symbol>> gj_reduce(sym_matrix& A, symbol_table& table) {
+ std::vector<std::vector<symge::symbol>> row_symbols;
+
if (A.nrow()>A.ncol()) throw std::runtime_error("improper matrix for reduction");
auto define_sym = [&table](symbol_term_diff t) { return table.define(t); };
@@ -123,12 +126,18 @@ void gj_reduce(sym_matrix& A, symbol_table& table) {
for (unsigned i = 0; i<A.nrow(); ++i) {
if (i==p.row || A[i].index(p.col)==msparse::row_npos) continue;
A[i] = row_reduce(p.col, A[i], A[p.row], define_sym);
+
+ std::vector<symge::symbol> row;
+ std::transform(A[i].begin(), A[i].end(), std::back_inserter(row),
+ [](auto&& entry){ return entry.value; });
+ row_symbols.emplace_back(std::move(row));
}
if (remaining_rows.empty()) {
break;
}
}
+ return row_symbols;
}
} // namespace symge
diff --git a/modcc/symge.hpp b/modcc/symge.hpp
index edd37b9031d6cad6640088661b04e544eb10692b..9343609255a8e4a6371345d312f1d4f68760b4da 100644
--- a/modcc/symge.hpp
+++ b/modcc/symge.hpp
@@ -153,6 +153,7 @@ using sym_matrix = msparse::matrix<symbol>;
// Perform Gauss-Jordan reduction on a (possibly augmented) symbolic matrix, with
// pivots taken from the diagonal elements. New symbol definitions due to fill-in
// will be added via the provided symbol table.
-void gj_reduce(sym_matrix& A, symbol_table& table);
+// Returns a vector of vectors of symbols, partitioned by row of the matrix
+std::vector<std::vector<symge::symbol>> gj_reduce(sym_matrix& A, symbol_table& table);
} // namespace symge
diff --git a/modcc/token.cpp b/modcc/token.cpp
index 6900d91cf40ecf10755388c0ba0994e6557ab057..3cdcfa4eeb966f277d39acf9a4ff984187c1b1ae 100644
--- a/modcc/token.cpp
+++ b/modcc/token.cpp
@@ -66,6 +66,7 @@ static Keyword keywords[] = {
{"log", tok::log},
{"fabs", tok::abs},
{"exprelr", tok::exprelr},
+ {"safeinv", tok::safeinv},
{"CONDUCTANCE", tok::conductance},
{nullptr, tok::reserved},
};
@@ -138,6 +139,7 @@ static TokenString token_strings[] = {
{"log", tok::log},
{"fabs", tok::abs},
{"exprelr", tok::exprelr},
+ {"safeinv", tok::safeinv},
{"cos", tok::cos},
{"sin", tok::sin},
{"cnexp", tok::cnexp},
diff --git a/modcc/token.hpp b/modcc/token.hpp
index 21467063dce1a2359f2a0f4a847b9426fd37c35f..3c5344a49b461c1c8d2fb827e83f924ac7dd99e1 100644
--- a/modcc/token.hpp
+++ b/modcc/token.hpp
@@ -69,7 +69,7 @@ enum class tok {
min, max,
// unary operators
- exp, sin, cos, log, abs,
+ exp, sin, cos, log, abs, safeinv,
exprelr, // equivalent to x/(exp(x)-1) with exprelr(0)=1
// logical keywords