summaryrefslogtreecommitdiffstatshomepage
path: root/src/lib/netlist/solver/nld_ms_sor.h
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_sor.h')
-rw-r--r--src/lib/netlist/solver/nld_ms_sor.h249
1 files changed, 119 insertions, 130 deletions
diff --git a/src/lib/netlist/solver/nld_ms_sor.h b/src/lib/netlist/solver/nld_ms_sor.h
index c31aaa6d46a..ca5a4839ce3 100644
--- a/src/lib/netlist/solver/nld_ms_sor.h
+++ b/src/lib/netlist/solver/nld_ms_sor.h
@@ -1,175 +1,164 @@
-// license:GPL-2.0+
+// license:BSD-3-Clause
// copyright-holders:Couriersud
-/*
- * nld_ms_sor.h
- *
- * Generic successive over relaxation solver.
- *
- * Fow w==1 we will do the classic Gauss-Seidel approach
- *
- */
#ifndef NLD_MS_SOR_H_
#define NLD_MS_SOR_H_
+// Names
+// spell-checker: words Seidel
+
+///
+/// \file nld_ms_sor.h
+///
+/// Generic successive over relaxation solver.
+///
+/// Fow w==1 we will do the classic Gauss-Seidel approach.
+///
+
+#include "nld_matrix_solver_ext.h"
#include "nld_ms_direct.h"
-#include "nld_solver.h"
#include <algorithm>
-namespace netlist
-{
- namespace devices
-{
-
-template <typename FT, int SIZE>
-class matrix_solver_SOR_t: public matrix_solver_direct_t<FT, SIZE>
+namespace netlist::solver
{
-public:
- using float_type = FT;
-
- matrix_solver_SOR_t(netlist_state_t &anetlist, const pstring &name, const solver_parameters_t *params, const std::size_t size)
- : matrix_solver_direct_t<FT, SIZE>(anetlist, name, matrix_solver_t::ASCENDING, params, size)
- , m_lp_fact(*this, "m_lp_fact", 0)
- , w(size, 0.0)
- , one_m_w(size, 0.0)
- , RHS(size, 0.0)
- //, new_V(size, 0.0)
- {
- }
+ template <typename FT, int SIZE>
+ class matrix_solver_SOR_t: public matrix_solver_direct_t<FT, SIZE>
+ {
+ public:
- void vsetup(analog_net_t::list_t &nets) override;
- unsigned vsolve_non_dynamic(const bool newton_raphson) override;
+ using float_type = FT;
-private:
- state_var<float_type> m_lp_fact;
- std::vector<float_type> w;
- std::vector<float_type> one_m_w;
- std::vector<float_type> RHS;
- //std::vector<float_type> new_V;
-};
+ matrix_solver_SOR_t(devices::nld_solver &main_solver, const pstring &name,
+ matrix_solver_t::net_list_t &nets,
+ const solver_parameters_t *params, const std::size_t size)
+ : matrix_solver_direct_t<FT, SIZE>(main_solver, name, nets, params, size)
+ , m_lp_fact(*this, "m_lp_fact", 0)
+ , w(size, plib::constants<FT>::zero())
+ , one_m_w(size, plib::constants<FT>::zero())
+ {
+ }
-// ----------------------------------------------------------------------------------------
-// matrix_solver - Gauss - Seidel
-// ----------------------------------------------------------------------------------------
+ void upstream_solve_non_dynamic() override;
+ private:
+ state_var<float_type> m_lp_fact;
+ std::vector<float_type> w;
+ std::vector<float_type> one_m_w;
+ };
-template <typename FT, int SIZE>
-void matrix_solver_SOR_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets)
-{
- matrix_solver_direct_t<FT, SIZE>::vsetup(nets);
-}
+ ///
+ /// \brief Gauss - Seidel matrix_solver
+ ///
+ ///
+ template <typename FT, int SIZE>
+ void matrix_solver_SOR_t<FT, SIZE>::upstream_solve_non_dynamic()
+ {
+ const std::size_t iN = this->size();
+ bool resched = false;
+ unsigned resched_cnt = 0;
-template <typename FT, int SIZE>
-unsigned matrix_solver_SOR_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson)
-{
- const std::size_t iN = this->size();
- bool resched = false;
- unsigned resched_cnt = 0;
+ // ideally, we could get an estimate for the spectral radius of
+ // Inv(D - L) * U
+ //
+ // and estimate using
+ //
+ // omega = 2.0 / (1.0 + std::sqrt(1-rho))
+ //
- /* ideally, we could get an estimate for the spectral radius of
- * Inv(D - L) * U
- *
- * and estimate using
- *
- * omega = 2.0 / (1.0 + std::sqrt(1-rho))
- */
+ const auto ws(static_cast<float_type>(this->m_params.m_gs_sor));
- const float_type ws = this->m_params.m_gs_sor;
+ for (std::size_t k = 0; k < iN; k++)
+ {
- for (std::size_t k = 0; k < iN; k++)
- {
- float_type gtot_t = 0.0;
- float_type gabs_t = 0.0;
- float_type RHS_t = 0.0;
+ const std::size_t term_count = this->m_terms[k].count();
+ const auto * const gt = this->m_gtn[k];
+ const auto * const go = this->m_gonn[k];
+ const auto * const Idr = this->m_Idrn[k];
+ auto other_cur_analog = this->m_connected_net_Vn[k];
- const std::size_t term_count = this->m_terms[k]->count();
- const float_type * const gt = this->m_gtn[k];
- const float_type * const go = this->m_gonn[k];
- const float_type * const Idr = this->m_Idrn[k];
- auto other_cur_analog = this->m_connected_net_Vn[k];
+ using fpaggtype = std::remove_reference_t<std::remove_cv_t<decltype(this->m_gtn[0][0])>>;
- this->m_new_V[k] = this->m_nets[k]->Q_Analog();
+ fpaggtype gtot_t = nlconst_base<fpaggtype>::zero();
+ fpaggtype gabs_t = nlconst_base<fpaggtype>::zero();
+ fpaggtype RHS_t = nlconst_base<fpaggtype>::zero();
- for (std::size_t i = 0; i < term_count; i++)
- {
- gtot_t = gtot_t + gt[i];
- RHS_t = RHS_t + Idr[i];
- }
+ this->m_new_V[k] = static_cast<float_type>(this->m_terms[k].getV());
- for (std::size_t i = this->m_terms[k]->m_railstart; i < term_count; i++)
- RHS_t = RHS_t - go[i] * *other_cur_analog[i];
+ for (std::size_t i = 0; i < term_count; i++)
+ {
+ gtot_t = gtot_t + gt[i];
+ RHS_t = RHS_t + Idr[i];
+ }
- RHS[k] = RHS_t;
+ for (std::size_t i = this->m_terms[k].rail_start(); i < term_count; i++)
+ RHS_t = RHS_t - go[i] * *other_cur_analog[i];
- if (this->m_params.m_use_gabs)
- {
- for (std::size_t i = 0; i < term_count; i++)
- gabs_t = gabs_t + std::abs(go[i]);
+ this->m_RHS[k] = static_cast<float_type>(RHS_t);
- gabs_t *= plib::constants<nl_double>::cast(0.5); // derived by try and error
- if (gabs_t <= gtot_t)
+ if (this->m_params.m_use_gabs)
{
- w[k] = ws / gtot_t;
- one_m_w[k] = plib::constants<FT>::one() - ws;
+ for (std::size_t i = 0; i < term_count; i++)
+ gabs_t = gabs_t + plib::abs(go[i]);
+
+ gabs_t *= nlconst::half(); // derived by try and error
+ if (gabs_t <= gtot_t)
+ {
+ w[k] = ws / static_cast<float_type>(gtot_t);
+ one_m_w[k] = plib::constants<FT>::one() - ws;
+ }
+ else
+ {
+ w[k] = plib::reciprocal(static_cast<float_type>(gtot_t + gabs_t));
+ one_m_w[k] = plib::constants<FT>::one() - plib::constants<FT>::one() * static_cast<FT>(gtot_t / (gtot_t + gabs_t));
+ }
}
else
{
- w[k] = plib::constants<FT>::one() / (gtot_t + gabs_t);
- one_m_w[k] = plib::constants<FT>::one() - plib::constants<FT>::one() * gtot_t / (gtot_t + gabs_t);
+ w[k] = ws / static_cast<float_type>(gtot_t);
+ one_m_w[k] = plib::constants<FT>::one() - ws;
}
}
- else
- {
- w[k] = ws / gtot_t;
- one_m_w[k] = plib::constants<FT>::one() - ws;
- }
- }
- const float_type accuracy = this->m_params.m_accuracy;
+ const auto accuracy(static_cast<float_type>(this->m_params.m_accuracy));
- do {
- resched = false;
- float_type err = 0;
- for (std::size_t k = 0; k < iN; k++)
- {
- const int * net_other = this->m_terms[k]->m_connected_net_idx.data();
- const std::size_t railstart = this->m_terms[k]->m_railstart;
- const float_type * go = this->m_gonn[k];
+ do {
+ resched = false;
+ float_type err = 0;
+ for (std::size_t k = 0; k < iN; k++)
+ {
+ const int * net_other = this->m_terms[k].m_connected_net_idx.data();
+ const std::size_t rail_start = this->m_terms[k].rail_start();
+ const auto * go = this->m_gonn[k];
- float_type Idrive = 0.0;
- for (std::size_t i = 0; i < railstart; i++)
- Idrive = Idrive - go[i] * this->m_new_V[static_cast<std::size_t>(net_other[i])];
+ float_type Idrive = plib::constants<float_type>::zero();
+ for (std::size_t i = 0; i < rail_start; i++)
+ Idrive = Idrive - static_cast<float_type>(go[i]) * this->m_new_V[static_cast<std::size_t>(net_other[i])];
- const float_type new_val = this->m_new_V[k] * one_m_w[k] + (Idrive + RHS[k]) * w[k];
+ const float_type new_val = this->m_new_V[k] * one_m_w[k] + (Idrive + this->m_RHS[k]) * w[k];
- err = std::max(std::abs(new_val - this->m_new_V[k]), err);
- this->m_new_V[k] = new_val;
- }
+ err = std::max(plib::abs(new_val - this->m_new_V[k]), err);
+ this->m_new_V[k] = new_val;
+ }
- if (err > accuracy)
- resched = true;
+ if (err > accuracy)
+ resched = true;
- resched_cnt++;
- } while (resched && (resched_cnt < this->m_params.m_gs_loops));
+ resched_cnt++;
+ } while (resched && (resched_cnt < this->m_params.m_gs_loops));
- this->m_iterative_total += resched_cnt;
- this->m_stat_calculations++;
+ this->m_iterative_total += resched_cnt;
- if (resched)
- {
- // Fallback to direct solver ...
- this->m_iterative_fail++;
- return matrix_solver_direct_t<FT, SIZE>::vsolve_non_dynamic(newton_raphson);
- }
+ if (resched)
+ {
+ // Fallback to direct solver ...
+ this->m_iterative_fail++;
+ matrix_solver_direct_t<FT, SIZE>::upstream_solve_non_dynamic();
+ }
- const float_type err = (newton_raphson ? this->delta(this->m_new_V) : 0.0);
- this->store(this->m_new_V);
- return (err > this->m_params.m_accuracy) ? 2 : 1;
-}
+ }
- } //namespace devices
-} // namespace netlist
+} // namespace netlist::solver
-#endif /* NLD_MS_SOR_H_ */
+#endif // NLD_MS_SOR_H_