// 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_MAT_H_ #define NLD_MS_SOR_MAT_H_ #include #include "solver/nld_ms_direct.h" #include "solver/nld_matrix_solver.h" #include "solver/nld_solver.h" namespace netlist { namespace devices { template class matrix_solver_SOR_mat_t: public matrix_solver_direct_t { friend class matrix_solver_t; public: matrix_solver_SOR_mat_t(netlist_t &anetlist, const pstring &name, const solver_parameters_t *params, int size) : matrix_solver_direct_t(anetlist, name, matrix_solver_t::DESCENDING, params, size) , m_Vdelta(*this, "m_Vdelta", 0.0) , m_omega(*this, "m_omega", params->m_sor) , m_lp_fact(*this, "m_lp_fact", 0) , m_gs_fail(*this, "m_gs_fail", 0) , m_gs_total(*this, "m_gs_total", 0) { } virtual ~matrix_solver_SOR_mat_t() {} virtual void vsetup(analog_net_t::list_t &nets) override; virtual int vsolve_non_dynamic(const bool newton_raphson) override; private: state_var m_Vdelta; state_var m_omega; state_var m_lp_fact; state_var m_gs_fail; state_var m_gs_total; }; // ---------------------------------------------------------------------------------------- // matrix_solver - Gauss - Seidel // ---------------------------------------------------------------------------------------- template void matrix_solver_SOR_mat_t::vsetup(analog_net_t::list_t &nets) { matrix_solver_direct_t::vsetup(nets); } #if 0 //FIXME: move to solve_base template nl_double matrix_solver_SOR_mat_t::vsolve() { /* * enable linear prediction on first newton pass */ if (USE_LINEAR_PREDICTION) for (unsigned k = 0; k < this->N(); k++) { this->m_last_V[k] = this->m_nets[k]->m_cur_Analog; this->m_nets[k]->m_cur_Analog = this->m_nets[k]->m_cur_Analog + this->m_Vdelta[k] * this->current_timestep() * m_lp_fact; } else for (unsigned k = 0; k < this->N(); k++) { this->m_last_V[k] = this->m_nets[k]->m_cur_Analog; } this->solve_base(this); if (USE_LINEAR_PREDICTION) { nl_double sq = 0; nl_double sqo = 0; const nl_double rez_cts = 1.0 / this->current_timestep(); for (unsigned k = 0; k < this->N(); k++) { const analog_net_t *n = this->m_nets[k]; const nl_double nv = (n->Q_Analog() - this->m_last_V[k]) * rez_cts ; sq += nv * nv; sqo += this->m_Vdelta[k] * this->m_Vdelta[k]; this->m_Vdelta[k] = nv; } // FIXME: used to be 1e90, but this would not be compatible with float if (sqo > NL_FCONST(1e-20)) m_lp_fact = std::min(std::sqrt(sq/sqo), (nl_double) 2.0); else m_lp_fact = NL_FCONST(0.0); } return this->compute_next_timestep(); } #endif template int matrix_solver_SOR_mat_t::vsolve_non_dynamic(const bool newton_raphson) { /* The matrix based code looks a lot nicer but actually is 30% slower than * the optimized code which works directly on the data structures. * Need something like that for gaussian elimination as well. */ nl_double new_v[storage_N] = { 0.0 }; const unsigned iN = this->N(); matrix_solver_t::build_LE_A(); matrix_solver_t::build_LE_RHS(); bool resched = false; int resched_cnt = 0; #if 0 static int ws_cnt = 0; ws_cnt++; if (1 && ws_cnt % 200 == 0) { // update omega nl_double lambdaN = 0; nl_double lambda1 = 1e9; for (int k = 0; k < iN; k++) { #if 0 nl_double akk = std::abs(this->m_A[k][k]); if ( akk > lambdaN) lambdaN = akk; if (akk < lambda1) lambda1 = akk; #else nl_double akk = std::abs(this->m_A[k][k]); nl_double s = 0.0; for (int i=0; im_A[k][i]); akk = s / akk - 1.0; if ( akk > lambdaN) lambdaN = akk; if (akk < lambda1) lambda1 = akk; #endif } //printf("lambda: %f %f\n", lambda, 2.0 / (1.0 + 2 * sqrt(lambda)) ); //ws = 2.0 / (2.0 - lambdaN - lambda1); m_omega = 2.0 / (2.0 - lambda1); //printf("%f %f %f\n", m_omega, lambda1, lambdaN); } #endif for (unsigned k = 0; k < iN; k++) new_v[k] = this->m_nets[k]->m_cur_Analog; do { resched = false; nl_double cerr = 0.0; for (unsigned k = 0; k < iN; k++) { nl_double Idrive = 0; const auto *p = this->m_terms[k]->m_nz.data(); const unsigned e = this->m_terms[k]->m_nz.size(); for (unsigned i = 0; i < e; i++) Idrive = Idrive + this->A(k,p[i]) * new_v[p[i]]; const nl_double delta = m_omega * (this->RHS(k) - Idrive) / this->A(k,k); cerr = std::max(cerr, std::abs(delta)); new_v[k] += delta; } if (cerr > this->m_params.m_accuracy) { resched = true; } resched_cnt++; } while (resched && (resched_cnt < this->m_params.m_gs_loops)); this->m_stat_calculations++; this->m_gs_total += resched_cnt; if (resched) { //this->netlist().warning("Falling back to direct solver .. Consider increasing RESCHED_LOOPS"); this->m_gs_fail++; return matrix_solver_direct_t::solve_non_dynamic(newton_raphson); } else { this->store(new_v); return resched_cnt; } } } //namespace devices } // namespace netlist #endif /* NLD_MS_GAUSS_SEIDEL_H_ */