summaryrefslogtreecommitdiffstatshomepage
path: root/src/lib/netlist/solver/nld_ms_direct_lu.h
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_direct_lu.h')
-rw-r--r--src/lib/netlist/solver/nld_ms_direct_lu.h170
1 files changed, 86 insertions, 84 deletions
diff --git a/src/lib/netlist/solver/nld_ms_direct_lu.h b/src/lib/netlist/solver/nld_ms_direct_lu.h
index e7cedc1dd29..c51ccac4da5 100644
--- a/src/lib/netlist/solver/nld_ms_direct_lu.h
+++ b/src/lib/netlist/solver/nld_ms_direct_lu.h
@@ -1,13 +1,18 @@
-// license:GPL-2.0+
+// license:BSD-3-Clause
// copyright-holders:Couriersud
-/*
- * nld_ms_direct.h
- *
- */
+
+///
+/// \file nld_ms_direct.h
+///
+///
#if 0
#ifndef NLD_MS_DIRECT_H_
#define NLD_MS_DIRECT_H_
+// Names
+// spell-checker: words Seidel,Crout
+
+
#include "solver/nld_solver.h"
#include "solver/nld_matrix_solver.h"
@@ -38,12 +43,12 @@ public:
unsigned N() const { if (m_N == 0) return m_dim; else return m_N; }
- int vsolve_non_dynamic(const bool newton_raphson);
+ int upstream_solve_non_dynamic(bool newton_raphson);
protected:
virtual void add_term(int net_idx, terminal_t *term) override;
- int solve_non_dynamic(const bool newton_raphson);
+ int solve_non_dynamic(bool newton_raphson);
void build_LE_A();
void build_LE_RHS(nl_double * RESTRICT rhs);
@@ -132,11 +137,11 @@ protected:
nl_double delta(const nl_double * RESTRICT V);
void store(const nl_double * RESTRICT V);
- /* bring the whole system to the current time
- * Don't schedule a new calculation time. The recalculation has to be
- * triggered by the caller after the netlist element was changed.
- */
- nl_double compute_next_timestep();
+ // bring the whole system to the current time
+ // Don't schedule a new calculation time. The recalculation has to be
+ // triggered by the caller after the netlist element was changed.
+
+ nl_double compute_next_time_step();
template <typename T1, typename T2>
nl_ext_double &A(const T1 r, const T2 c) { return m_A[r][c]; }
@@ -164,65 +169,64 @@ matrix_solver_direct_t<m_N, storage_N>::~matrix_solver_direct_t()
}
template <unsigned m_N, unsigned storage_N>
-nl_double matrix_solver_direct_t<m_N, storage_N>::compute_next_timestep()
+nl_double matrix_solver_direct_t<m_N, storage_N>::compute_next_time_step()
{
- nl_double new_solver_timestep = m_params.m_max_timestep;
+ nl_double new_solver_time_step = m_params.m_max_time_step;
if (m_params.m_dynamic_ts)
{
- /*
- * FIXME: We should extend the logic to use either all nets or
- * only output nets.
- */
+ //
+ // FIXME: We should extend the logic to use either all nets or
+ // only output nets.
for (unsigned k = 0, iN=N(); k < iN; k++)
{
analog_net_t *n = m_nets[k];
const nl_double DD_n = (n->Q_Analog() - m_last_V[k]);
- const nl_double hn = current_timestep();
+ const nl_double hn = current_time_step();
nl_double DD2 = (DD_n / hn - n->m_DD_n_m_1 / n->m_h_n_m_1) / (hn + n->m_h_n_m_1);
- nl_double new_net_timestep;
+ nl_double new_net_time_step;
n->m_h_n_m_1 = hn;
n->m_DD_n_m_1 = DD_n;
- if (std::abs(DD2) > NL_FCONST(1e-30)) // avoid div-by-zero
- new_net_timestep = std::sqrt(m_params.m_dynamic_lte / std::abs(NL_FCONST(0.5)*DD2));
+ if (plib::abs(DD2) > NL_FCONST(1e-30)) // avoid div-by-zero
+ new_net_time_step = std::sqrt(m_params.m_dynamic_lte / plib::abs(NL_FCONST(0.5)*DD2));
else
- new_net_timestep = m_params.m_max_timestep;
+ new_net_time_step = m_params.m_max_time_step;
- if (new_net_timestep < new_solver_timestep)
- new_solver_timestep = new_net_timestep;
+ if (new_net_time_step < new_solver_time_step)
+ new_solver_time_step = new_net_time_step;
}
- if (new_solver_timestep < m_params.m_min_timestep)
- new_solver_timestep = m_params.m_min_timestep;
- if (new_solver_timestep > m_params.m_max_timestep)
- new_solver_timestep = m_params.m_max_timestep;
+ if (new_solver_time_step < m_params.m_min_time_step)
+ new_solver_time_step = m_params.m_min_time_step;
+ if (new_solver_time_step > m_params.m_max_time_step)
+ new_solver_time_step = m_params.m_max_time_step;
}
- //if (new_solver_timestep > 10.0 * hn)
- // new_solver_timestep = 10.0 * hn;
- return new_solver_timestep;
+ //#if (new_solver_time_step > 10.0 * hn)
+ //# new_solver_time_step = 10.0 * hn;
+ return new_solver_time_step;
}
template <unsigned m_N, unsigned storage_N>
void matrix_solver_direct_t<m_N, storage_N>::add_term(int k, terminal_t *term)
{
- if (term->m_otherterm->net().isRailNet())
+ if (term->m_other_terminal->net().isRailNet())
{
m_rails_temp[k].add(term, -1, false);
}
else
{
- int ot = get_net_idx(&term->m_otherterm->net());
+ int ot = get_net_idx(&term->m_other_terminal->net());
if (ot>=0)
{
m_terms[k]->add(term, ot, true);
}
- /* Should this be allowed ? */
+ // Should this be allowed ?
else // if (ot<0)
{
m_rails_temp[k].add(term, ot, true);
- netlist().error("found term with missing othernet {1}\n", term->name());
+ netlist().error("found term with missing other net {1}\n", term->name());
}
}
}
@@ -244,7 +248,7 @@ void matrix_solver_direct_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
for (unsigned k = 0; k < N(); k++)
{
- m_terms[k]->m_railstart = m_terms[k]->count();
+ m_terms[k]->m_rail_start = m_terms[k]->count();
for (unsigned i = 0; i < m_rails_temp[k].count(); i++)
this->m_terms[k]->add(m_rails_temp[k].terms()[i], m_rails_temp[k].connected_net_idx()[i], false);
@@ -254,32 +258,32 @@ void matrix_solver_direct_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
#if 1
- /* Sort in descending order by number of connected matrix voltages.
- * The idea is, that for Gauss-Seidel algo the first voltage computed
- * depends on the greatest number of previous voltages thus taking into
- * account the maximum amout of information.
- *
- * This actually improves performance on popeye slightly. Average
- * GS computations reduce from 2.509 to 2.370
- *
- * Smallest to largest : 2.613
- * Unsorted : 2.509
- * Largest to smallest : 2.370
- *
- * Sorting as a general matrix pre-conditioning is mentioned in
- * literature but I have found no articles about Gauss Seidel.
- *
- * For Gaussian Elimination however increasing order is better suited.
- * FIXME: Even better would be to sort on elements right of the matrix diagonal.
- *
- */
+ // Sort in descending order by number of connected matrix voltages.
+ // The idea is, that for Gauss-Seidel algo the first voltage computed
+ // depends on the greatest number of previous voltages thus taking into
+ // account the maximum amount of information.
+ //
+ // This actually improves performance on popeye slightly. Average
+ // GS computations reduce from 2.509 to 2.370
+ //
+ // Smallest to largest : 2.613
+ // Unsorted : 2.509
+ // Largest to smallest : 2.370
+ //
+ // Sorting as a general matrix pre-conditioning is mentioned in
+ // literature but I have found no articles about Gauss Seidel.
+ //
+ // For Gaussian Elimination however increasing order is better suited.
+ // FIXME: Even better would be to sort on elements right of the matrix diagonal.
+ //
+ //
int sort_order = (type() == GAUSS_SEIDEL ? 1 : -1);
for (unsigned k = 0; k < N() / 2; k++)
for (unsigned i = 0; i < N() - 1; i++)
{
- if ((m_terms[i]->m_railstart - m_terms[i+1]->m_railstart) * sort_order < 0)
+ if ((m_terms[i]->m_rail_start - m_terms[i+1]->m_rail_start) * sort_order < 0)
{
std::swap(m_terms[i],m_terms[i+1]);
std::swap(m_nets[i], m_nets[i+1]);
@@ -291,19 +295,18 @@ void matrix_solver_direct_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
int *other = m_terms[k]->connected_net_idx();
for (unsigned i = 0; i < m_terms[k]->count(); i++)
if (other[i] != -1)
- other[i] = get_net_idx(&m_terms[k]->terms()[i]->m_otherterm->net());
+ other[i] = get_net_idx(&m_terms[k]->terms()[i]->m_other_terminal->net());
}
#endif
- /* create a list of non zero elements right of the diagonal
- * These list anticipate the population of array elements by
- * Gaussian elimination.
- */
+ // create a list of non zero elements right of the diagonal
+ // These list anticipate the population of array elements by
+ // Gaussian elimination.
+
for (unsigned k = 0; k < N(); k++)
{
terms_for_net_t * t = m_terms[k];
- /* pretty brutal */
int *other = t->connected_net_idx();
t->m_nz.clear();
@@ -325,7 +328,7 @@ void matrix_solver_direct_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
for (unsigned j = 0; j < N(); j++)
{
- for (unsigned i = 0; i < t->m_railstart; i++)
+ for (unsigned i = 0; i < t->m_rail_start; i++)
{
if (!t->m_nzrd.contains(other[i]) && other[i] >= (int) (k + 1))
t->m_nzrd.add(other[i]);
@@ -349,9 +352,10 @@ void matrix_solver_direct_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
log("\n");
}
- /*
- * save states
- */
+ //
+ // save states
+ //
+
save(NLNAME(m_RHS));
save(NLNAME(m_last_V));
@@ -378,7 +382,7 @@ void matrix_solver_direct_t<m_N, storage_N>::build_LE_A()
nl_double akk = 0.0;
const unsigned terms_count = m_terms[k]->count();
- const unsigned railstart = m_terms[k]->m_railstart;
+ const unsigned rail_start = m_terms[k]->m_rail_start;
const nl_double * RESTRICT gt = m_terms[k]->gt();
const nl_double * RESTRICT go = m_terms[k]->go();
const int * RESTRICT net_other = m_terms[k]->connected_net_idx();
@@ -388,7 +392,7 @@ void matrix_solver_direct_t<m_N, storage_N>::build_LE_A()
A(k,k) += akk;
- for (unsigned i = 0; i < railstart; i++)
+ for (unsigned i = 0; i < rail_start; i++)
A(k, net_other[i]) -= go[i];
}
}
@@ -410,8 +414,8 @@ void matrix_solver_direct_t<m_N, storage_N>::build_LE_RHS(nl_double * RESTRICT r
for (int i = 0; i < terms_count; i++)
rhsk_a = rhsk_a + Idr[i];
- for (int i = m_terms[k]->m_railstart; i < terms_count; i++)
- //rhsk = rhsk + go[i] * terms[i]->m_otherterm->net().as_analog().Q_Analog();
+ for (int i = m_terms[k]->m_rail_start; i < terms_count; i++)
+ //#rhsk = rhsk + go[i] * terms[i]->m_other_terminal->net().as_analog().Q_Analog();
rhsk_b = rhsk_b + go[i] * *other_cur_analog[i];
rhs[k] = rhsk_a + rhsk_b;
@@ -426,8 +430,8 @@ void matrix_solver_direct_t<m_N, storage_N>::LE_solve()
{
const unsigned kN = N();
- ATTR_UNUSED int imax;
- ATTR_UNUSED double big,temp;
+ [[maybe_unused]] int imax;
+ [[maybe_unused]] double big,temp;
#if 0
double vv[storage_N];
@@ -438,7 +442,7 @@ void matrix_solver_direct_t<m_N, storage_N>::LE_solve()
for (j=0;j<kN;j++)
if ((temp=fabs(m_A[i][j])) > big)
big=temp;
- //if (big == 0.0) nrerror("Singular matrix in routine LUDCMP");
+ //#if (big == 0.0) nrerror("Singular matrix in routine LUDCMP");
vv[i]=1.0/big;
}
#endif
@@ -489,7 +493,6 @@ void matrix_solver_direct_t<m_N, storage_N>::LE_solve()
m_A[imax][k]=m_A[j][k];
m_A[j][k]=dum;
}
- //*d = -(*d);
vv[imax]=vv[j];
}
indx[j]=imax;
@@ -509,7 +512,7 @@ void matrix_solver_direct_t<m_N, storage_N>::LE_back_subst(
{
const unsigned kN = N();
- /* back substitution */
+ // back substitution
// int ip;
// ii=-1
@@ -541,15 +544,14 @@ template <unsigned m_N, unsigned storage_N>
nl_double matrix_solver_direct_t<m_N, storage_N>::delta(
const nl_double * RESTRICT V)
{
- /* FIXME: Ideally we should also include currents (RHS) here. This would
- * need a revaluation of the right hand side after voltages have been updated
- * and thus belong into a different calculation. This applies to all solvers.
- */
+ // FIXME: Ideally we should also include currents (RHS) here. This would
+ // need a revaluation of the right hand side after voltages have been updated
+ // and thus belong into a different calculation. This applies to all solvers.
const unsigned iN = this->N();
nl_double cerr = 0;
for (unsigned i = 0; i < iN; i++)
- cerr = std::fmax(cerr, std::abs(V[i] - this->m_nets[i]->m_cur_Analog));
+ cerr = std::fmax(cerr, plib::abs(V[i] - this->m_nets[i]->m_cur_Analog));
return cerr;
}
@@ -565,7 +567,7 @@ void matrix_solver_direct_t<m_N, storage_N>::store(
template <unsigned m_N, unsigned storage_N>
-unsigned matrix_solver_direct_t<m_N, storage_N>::solve_non_dynamic(const bool newton_raphson)
+unsigned matrix_solver_direct_t<m_N, storage_N>::solve_non_dynamic(bool newton_raphson)
{
nl_double new_V[storage_N]; // = { 0.0 };
@@ -587,7 +589,7 @@ unsigned matrix_solver_direct_t<m_N, storage_N>::solve_non_dynamic(const bool ne
}
template <unsigned m_N, unsigned storage_N>
-int matrix_solver_direct_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson)
+int matrix_solver_direct_t<m_N, storage_N>::upstream_solve_non_dynamic(bool newton_raphson)
{
this->build_LE_A();
this->build_LE_RHS(m_RHS);
@@ -621,5 +623,5 @@ matrix_solver_direct_t<m_N, storage_N>::matrix_solver_direct_t(const eSolverType
} //namespace devices
} // namespace netlist
-#endif /* NLD_MS_DIRECT_H_ */
+#endif // NLD_MS_DIRECT_H_
#endif