// license:GPL-2.0+
// 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_
#include <algorithm>
#include "solver/nld_ms_direct.h"
#include "solver/nld_solver.h"
namespace netlist
{
namespace devices
{
template <unsigned m_N, unsigned storage_N>
class matrix_solver_SOR_t: public matrix_solver_direct_t<m_N, storage_N>
{
public:
matrix_solver_SOR_t(netlist_t &anetlist, const pstring &name, const solver_parameters_t *params, int size)
: matrix_solver_direct_t<m_N, storage_N>(anetlist, name, matrix_solver_t::ASCENDING, params, size)
, m_lp_fact(*this, "m_lp_fact", 0)
{
}
virtual ~matrix_solver_SOR_t() {}
virtual void vsetup(analog_net_t::list_t &nets) override;
virtual int vsolve_non_dynamic(const bool newton_raphson) override;
private:
state_var<nl_double> m_lp_fact;
};
// ----------------------------------------------------------------------------------------
// matrix_solver - Gauss - Seidel
// ----------------------------------------------------------------------------------------
template <unsigned m_N, unsigned storage_N>
void matrix_solver_SOR_t<m_N, storage_N>::vsetup(analog_net_t::list_t &nets)
{
matrix_solver_direct_t<m_N, storage_N>::vsetup(nets);
}
template <unsigned m_N, unsigned storage_N>
int matrix_solver_SOR_t<m_N, storage_N>::vsolve_non_dynamic(const bool newton_raphson)
{
const unsigned iN = this->N();
bool resched = false;
int 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))
*/
const nl_double ws = this->m_params.m_sor;
nl_double w[storage_N];
nl_double one_m_w[storage_N];
nl_double RHS[storage_N];
nl_double new_V[storage_N];
for (unsigned k = 0; k < iN; k++)
{
nl_double gtot_t = 0.0;
nl_double gabs_t = 0.0;
nl_double RHS_t = 0.0;
const unsigned term_count = this->m_terms[k]->count();
const nl_double * const RESTRICT gt = this->m_terms[k]->gt();
const nl_double * const RESTRICT go = this->m_terms[k]->go();
const nl_double * const RESTRICT Idr = this->m_terms[k]->Idr();
const nl_double * const *other_cur_analog = this->m_terms[k]->other_curanalog();
new_V[k] = this->m_nets[k]->m_cur_Analog;
for (unsigned i = 0; i < term_count; i++)
{
gtot_t = gtot_t + gt[i];
RHS_t = RHS_t + Idr[i];
}
for (unsigned i = this->m_terms[k]->m_railstart; i < term_count; i++)
RHS_t = RHS_t + go[i] * *other_cur_analog[i];
RHS[k] = RHS_t;
if (USE_GABS)
{
for (unsigned i = 0; i < term_count; i++)
gabs_t = gabs_t + std::abs(go[i]);
gabs_t *= NL_FCONST(0.5); // derived by try and error
if (gabs_t <= gtot_t)
{
w[k] = ws / gtot_t;
one_m_w[k] = NL_FCONST(1.0) - ws;
}
else
{
w[k] = NL_FCONST(1.0) / (gtot_t + gabs_t);
one_m_w[k] = NL_FCONST(1.0) - NL_FCONST(1.0) * gtot_t / (gtot_t + gabs_t);
}
}
else
{
w[k] = ws / gtot_t;
one_m_w[k] = NL_FCONST(1.0) - ws;
}
}
const nl_double accuracy = this->m_params.m_accuracy;
/* uncommenting the line below will force dynamic updates every X iterations
* althought the system has not converged yet. This is a proof of concept,
*
*/
const bool interleaved_dynamic_updates = false;
//const bool interleaved_dynamic_updates = newton_raphson;
do {
resched = false;
nl_double err = 0;
for (unsigned k = 0; k < iN; k++)
{
const int * RESTRICT net_other = this->m_terms[k]->net_other();
const unsigned railstart = this->m_terms[k]->m_railstart;
const nl_double * RESTRICT go = this->m_terms[k]->go();
nl_double Idrive = 0.0;
for (unsigned i = 0; i < railstart; i++)
Idrive = Idrive + go[i] * new_V[net_other[i]];
const nl_double new_val = new_V[k] * one_m_w[k] + (Idrive + RHS[k]) * w[k];
err = std::max(std::abs(new_val - new_V[k]), err);
new_V[k] = new_val;
}
if (err > accuracy)
resched = true;
resched_cnt++;
//} while (resched && (resched_cnt < this->m_params.m_gs_loops));
} while (resched && ((!interleaved_dynamic_updates && resched_cnt < this->m_params.m_gs_loops) || (interleaved_dynamic_updates && resched_cnt < 5 )));
this->m_iterative_total += resched_cnt;
if (resched && !interleaved_dynamic_updates)
{
// Fallback to direct solver ...
this->m_iterative_fail++;
return matrix_solver_direct_t<m_N, storage_N>::vsolve_non_dynamic(newton_raphson);
}
this->m_stat_calculations++;
if (interleaved_dynamic_updates)
{
for (unsigned k = 0; k < iN; k++)
this->m_nets[k]->m_cur_Analog += 1.0 * (new_V[k] - this->m_nets[k]->m_cur_Analog);
}
else
{
for (unsigned k = 0; k < iN; k++)
this->m_nets[k]->m_cur_Analog = new_V[k];
}
return resched_cnt;
}
} //namespace devices
} // namespace netlist
#endif /* NLD_MS_SOR_H_ */