1 // license:GPL-2.0+ 2 // copyright-holders:Couriersud 3 4 5 #include "nl_factory.h" 6 #include "core/setup.h" 7 #include "nl_errstr.h" 8 #include "nl_setup.h" // FIXME: only needed for splitter code 9 #include "nld_matrix_solver.h" 10 #include "nld_ms_direct.h" 11 #include "nld_ms_direct1.h" 12 #include "nld_ms_direct2.h" 13 #include "nld_ms_gcr.h" 14 #include "nld_ms_gmres.h" 15 #include "nld_ms_sm.h" 16 #include "nld_ms_sor.h" 17 #include "nld_ms_sor_mat.h" 18 #include "nld_ms_w.h" 19 #include "nld_solver.h" 20 #include "plib/pomp.h" 21 #include "plib/ptimed_queue.h" 22 23 #include <algorithm> 24 #include <type_traits> 25 26 namespace netlist 27 { 28 namespace devices 29 { 30 31 // ---------------------------------------------------------------------------------------- 32 // solver 33 // ---------------------------------------------------------------------------------------- 34 NETLIB_RESET(solver)35 NETLIB_RESET(solver) 36 { 37 if (exec().stats_enabled()) 38 m_fb_step.set_delegate(NETLIB_DELEGATE(fb_step<true>)); 39 for (auto &s : m_mat_solvers) 40 s->reset(); 41 for (auto &s : m_mat_solvers) 42 m_queue.push<false>({netlist_time_ext::zero(), s.get()}); 43 } 44 NETLIB_NAME(solver)45 void NETLIB_NAME(solver)::stop() 46 { 47 for (auto &s : m_mat_solvers) 48 s->log_stats(); 49 } 50 51 #if 1 52 53 template<bool KEEP_STATS> NETLIB_HANDLER(solver,fb_step)54 NETLIB_HANDLER(solver, fb_step) 55 { 56 const netlist_time_ext now(exec().time()); 57 const std::size_t nthreads = m_params.m_parallel() < 2 ? 1 : std::min(static_cast<std::size_t>(m_params.m_parallel()), plib::omp::get_max_threads()); 58 const netlist_time_ext sched(now + (nthreads <= 1 ? netlist_time_ext::zero() : netlist_time_ext::from_nsec(100))); 59 plib::uninitialised_array<solver::matrix_solver_t *, config::MAX_SOLVER_QUEUE_SIZE::value> tmp; //NOLINT 60 plib::uninitialised_array<netlist_time, config::MAX_SOLVER_QUEUE_SIZE::value> nt; //NOLINT 61 std::size_t p=0; 62 63 while (!m_queue.empty()) 64 { 65 const auto t = m_queue.top().exec_time(); 66 auto *o = m_queue.top().object(); 67 if (t != now) 68 if (t > sched) 69 break; 70 tmp[p++] = o; 71 m_queue.pop(); 72 } 73 74 // FIXME: Disabled for now since parallel processing will decrease performance 75 // for tested applications. More testing required here 76 if (true || nthreads < 2) 77 { 78 if (!KEEP_STATS) 79 { 80 for (std::size_t i = 0; i < p; i++) 81 nt[i] = tmp[i]->solve(now, "no-parallel"); 82 } 83 else 84 { 85 stats()->m_stat_total_time.stop(); 86 for (std::size_t i = 0; i < p; i++) 87 { 88 tmp[i]->stats()->m_stat_call_count.inc(); 89 auto g(tmp[i]->stats()->m_stat_total_time.guard()); 90 nt[i] = tmp[i]->solve(now, "no-parallel"); 91 } 92 stats()->m_stat_total_time.start(); 93 } 94 95 for (std::size_t i = 0; i < p; i++) 96 { 97 if (nt[i] != netlist_time::zero()) 98 m_queue.push<false>({now + nt[i], tmp[i]}); 99 tmp[i]->update_inputs(); 100 } 101 } 102 else 103 { 104 plib::omp::set_num_threads(nthreads); 105 plib::omp::for_static(static_cast<std::size_t>(0), p, [&tmp, &nt,now](std::size_t i) 106 { 107 nt[i] = tmp[i]->solve(now, "parallel"); 108 }); 109 for (std::size_t i = 0; i < p; i++) 110 { 111 if (nt[i] != netlist_time::zero()) 112 m_queue.push<false>({now + nt[i], tmp[i]}); 113 tmp[i]->update_inputs(); 114 } 115 } 116 if (!m_queue.empty()) 117 m_Q_step.net().toggle_and_push_to_queue(static_cast<netlist_time>(m_queue.top().exec_time() - now)); 118 } 119 NETLIB_NAME(solver)120 void NETLIB_NAME(solver) :: reschedule(solver::matrix_solver_t *solv, netlist_time ts) 121 { 122 const netlist_time_ext now(exec().time()); 123 const netlist_time_ext sched(now + ts); 124 m_queue.remove<false>(solv); 125 m_queue.push<false>({sched, solv}); 126 127 if (m_Q_step.net().is_queued()) 128 { 129 if (m_Q_step.net().next_scheduled_time() > sched) 130 m_Q_step.net().toggle_and_push_to_queue(ts); 131 } 132 else 133 m_Q_step.net().toggle_and_push_to_queue(ts); 134 } 135 #else NETLIB_HANDLER(solver,fb_step)136 NETLIB_HANDLER(solver, fb_step) 137 { 138 if (m_params.m_dynamic_ts) 139 return; 140 141 netlist_time_ext now(exec().time()); 142 // force solving during start up if there are no time-step devices 143 // FIXME: Needs a more elegant solution 144 bool force_solve = (now < netlist_time_ext::from_fp<decltype(m_params.m_max_timestep)>(2 * m_params.m_max_timestep)); 145 146 std::size_t nthreads = std::min(static_cast<std::size_t>(m_params.m_parallel()), plib::omp::get_max_threads()); 147 148 std::vector<solver_entry *> &solvers = (force_solve ? m_mat_solvers_all : m_mat_solvers_timestepping); 149 150 if (nthreads > 1 && solvers.size() > 1) 151 { 152 plib::omp::set_num_threads(nthreads); 153 plib::omp::for_static(static_cast<std::size_t>(0), solvers.size(), [&solvers, now](std::size_t i) 154 { 155 const netlist_time ts = solvers[i]->ptr->solve(now); 156 plib::unused_var(ts); 157 }); 158 } 159 else 160 for (auto & solver : solvers) 161 { 162 const netlist_time ts = solver->ptr->solve(now); 163 plib::unused_var(ts); 164 } 165 166 for (auto & solver : solvers) 167 solver->ptr->update_inputs(); 168 169 // step circuit 170 if (!m_Q_step.net().is_queued()) 171 { 172 m_Q_step.net().toggle_and_push_to_queue(netlist_time::from_fp(m_params.m_max_timestep)); 173 } 174 } 175 #endif 176 177 // FIXME: should be created in device space 178 template <class C> NETLIB_NAME(solver)179 NETLIB_NAME(solver)::solver_ptr create_it(NETLIB_NAME(solver) &main_solver, pstring name, 180 NETLIB_NAME(solver)::net_list_t &nets, 181 const solver::solver_parameters_t *params, std::size_t size) 182 { 183 return plib::make_unique<C, device_arena>(main_solver, name, nets, params, size); 184 } 185 186 template <typename FT, int SIZE> NETLIB_NAME(solver)187 NETLIB_NAME(solver)::solver_ptr NETLIB_NAME(solver)::create_solver(std::size_t size, 188 const pstring &solvername, const solver::solver_parameters_t *params, 189 NETLIB_NAME(solver)::net_list_t &nets) 190 { 191 switch (params->m_method()) 192 { 193 case solver::matrix_type_e::MAT_CR: 194 return create_it<solver::matrix_solver_GCR_t<FT, SIZE>>(*this, solvername, nets, params, size); 195 case solver::matrix_type_e::MAT: 196 return create_it<solver::matrix_solver_direct_t<FT, SIZE>>(*this, solvername, nets, params, size); 197 case solver::matrix_type_e::GMRES: 198 return create_it<solver::matrix_solver_GMRES_t<FT, SIZE>>(*this, solvername, nets, params, size); 199 #if (NL_USE_ACADEMIC_SOLVERS) 200 case solver::matrix_type_e::SOR: 201 return create_it<solver::matrix_solver_SOR_t<FT, SIZE>>(*this, solvername, nets, params, size); 202 case solver::matrix_type_e::SOR_MAT: 203 return create_it<solver::matrix_solver_SOR_mat_t<FT, SIZE>>(*this, solvername, nets, params, size); 204 case solver::matrix_type_e::SM: 205 // Sherman-Morrison Formula 206 return create_it<solver::matrix_solver_sm_t<FT, SIZE>>(*this, solvername, nets, params, size); 207 case solver::matrix_type_e::W: 208 // Woodbury Formula 209 return create_it<solver::matrix_solver_w_t<FT, SIZE>>(*this, solvername, nets, params, size); 210 #else 211 //case solver::matrix_type_e::GMRES: 212 case solver::matrix_type_e::SOR: 213 case solver::matrix_type_e::SOR_MAT: 214 case solver::matrix_type_e::SM: 215 case solver::matrix_type_e::W: 216 state().log().warning(MW_SOLVER_METHOD_NOT_SUPPORTED(params->m_method().name(), "MAT_CR")); 217 return create_it<solver::matrix_solver_GCR_t<FT, SIZE>>(*this, solvername, nets, params, size); 218 #endif 219 } 220 return solver_ptr(); 221 } 222 223 template <typename FT> NETLIB_NAME(solver)224 NETLIB_NAME(solver)::solver_ptr NETLIB_NAME(solver)::create_solvers( 225 const pstring &sname, const solver::solver_parameters_t *params, 226 net_list_t &nets) 227 { 228 std::size_t net_count = nets.size(); 229 switch (net_count) 230 { 231 #if !defined(__EMSCRIPTEN__) 232 case 1: 233 return plib::make_unique<solver::matrix_solver_direct1_t<FT>, device_arena>(*this, sname, nets, params); 234 case 2: 235 return plib::make_unique<solver::matrix_solver_direct2_t<FT>, device_arena>(*this, sname, nets, params); 236 case 3: 237 return create_solver<FT, 3>(3, sname, params, nets); 238 case 4: 239 return create_solver<FT, 4>(4, sname, params, nets); 240 case 5: 241 return create_solver<FT, 5>(5, sname, params, nets); 242 case 6: 243 return create_solver<FT, 6>(6, sname, params, nets); 244 case 7: 245 return create_solver<FT, 7>(7, sname, params, nets); 246 case 8: 247 return create_solver<FT, 8>(8, sname, params, nets); 248 #endif 249 default: 250 log().info(MI_NO_SPECIFIC_SOLVER(net_count)); 251 if (net_count <= 16) 252 { 253 return create_solver<FT, -16>(net_count, sname, params, nets); 254 } 255 if (net_count <= 32) 256 { 257 return create_solver<FT, -32>(net_count, sname, params, nets); 258 } 259 if (net_count <= 64) 260 { 261 return create_solver<FT, -64>(net_count, sname, params, nets); 262 } 263 if (net_count <= 128) 264 { 265 return create_solver<FT, -128>(net_count, sname, params, nets); 266 } 267 if (net_count <= 256) 268 { 269 return create_solver<FT, -256>(net_count, sname, params, nets); 270 } 271 if (net_count <= 512) 272 { 273 return create_solver<FT, -512>(net_count, sname, params, nets); 274 } 275 return create_solver<FT, 0>(net_count, sname, params, nets); 276 } 277 } 278 279 struct net_splitter 280 { runnetlist::devices::net_splitter281 void run(netlist_state_t &nlstate) 282 { 283 for (auto & net : nlstate.nets()) 284 { 285 nlstate.log().verbose("processing {1}", net->name()); 286 if (!net->is_rail_net() && !nlstate.core_terms(*net).empty()) 287 { 288 nlstate.log().verbose(" ==> not a rail net"); 289 // Must be an analog net 290 auto &n = dynamic_cast<analog_net_t &>(*net); 291 if (!already_processed(n)) 292 { 293 groupspre.emplace_back(NETLIB_NAME(solver)::net_list_t()); 294 process_net(nlstate, n); 295 } 296 } 297 } 298 for (auto &g : groupspre) 299 if (!g.empty()) 300 groups.push_back(g); 301 } 302 303 std::vector<NETLIB_NAME(solver)::net_list_t> groups; 304 305 private: 306 already_processednetlist::devices::net_splitter307 bool already_processed(const analog_net_t &n) const 308 { 309 // no need to process rail nets - these are known variables 310 if (n.is_rail_net()) 311 return true; 312 // if it's already processed - no need to continue 313 for (const auto & grp : groups) 314 if (plib::container::contains(grp, &n)) 315 return true; 316 return false; 317 } 318 check_if_processed_and_joinnetlist::devices::net_splitter319 bool check_if_processed_and_join(const analog_net_t &n) 320 { 321 // no need to process rail nets - these are known variables 322 if (n.is_rail_net()) 323 return true; 324 // First check if it is in a previous group. 325 // In this case we need to merge this group into the current group 326 if (groupspre.size() > 1) 327 { 328 for (std::size_t i = 0; i<groupspre.size() - 1; i++) 329 if (plib::container::contains(groupspre[i], &n)) 330 { 331 // copy all nets 332 for (auto & cn : groupspre[i]) 333 if (!plib::container::contains(groupspre.back(), cn)) 334 groupspre.back().push_back(cn); 335 // clear 336 groupspre[i].clear(); 337 return true; 338 } 339 } 340 // if it's already processed - no need to continue 341 if (!groupspre.empty() && plib::container::contains(groupspre.back(), &n)) 342 return true; 343 return false; 344 } 345 346 // NOLINTNEXTLINE(misc-no-recursion) process_netnetlist::devices::net_splitter347 void process_net(netlist_state_t &nlstate, analog_net_t &n) 348 { 349 // ignore empty nets. FIXME: print a warning message 350 nlstate.log().verbose("Net {}", n.name()); 351 if (!nlstate.core_terms(n).empty()) 352 { 353 // add the net 354 groupspre.back().push_back(&n); 355 // process all terminals connected to this net 356 for (auto &term : nlstate.core_terms(n)) 357 { 358 nlstate.log().verbose("Term {} {}", term->name(), static_cast<int>(term->type())); 359 // only process analog terminals 360 if (term->is_type(detail::terminal_type::TERMINAL)) 361 { 362 auto &pt = dynamic_cast<terminal_t &>(*term); 363 // check the connected terminal 364 const auto *const connected_terminals = nlstate.setup().get_connected_terminals(pt); 365 for (auto ct = connected_terminals->begin(); *ct != nullptr; ct++) 366 { 367 analog_net_t &connected_net = (*ct)->net(); 368 nlstate.log().verbose(" Connected net {}", connected_net.name()); 369 if (!check_if_processed_and_join(connected_net)) 370 process_net(nlstate, connected_net); 371 } 372 } 373 } 374 } 375 } 376 377 std::vector<NETLIB_NAME(solver)::net_list_t> groupspre; 378 }; 379 NETLIB_NAME(solver)380 void NETLIB_NAME(solver)::post_start() 381 { 382 log().verbose("Scanning net groups ..."); 383 // determine net groups 384 385 net_splitter splitter; 386 387 splitter.run(state()); 388 log().verbose("Found {1} net groups in {2} nets\n", splitter.groups.size(), state().nets().size()); 389 390 int num_errors = 0; 391 392 log().verbose("checking net consistency ..."); 393 for (const auto &grp : splitter.groups) 394 { 395 int railterms = 0; 396 pstring nets_in_grp; 397 for (const auto &n : grp) 398 { 399 nets_in_grp += (n->name() + " "); 400 if (!n->is_analog()) 401 { 402 state().log().error(ME_SOLVER_CONSISTENCY_NOT_ANALOG_NET(n->name())); 403 num_errors++; 404 } 405 if (n->is_rail_net()) 406 { 407 state().log().error(ME_SOLVER_CONSISTENCY_RAIL_NET(n->name())); 408 num_errors++; 409 } 410 for (const auto &t : state().core_terms(*n)) 411 { 412 if (!t->has_net()) 413 { 414 state().log().error(ME_SOLVER_TERMINAL_NO_NET(t->name())); 415 num_errors++; 416 } 417 else 418 { 419 auto *otherterm = dynamic_cast<terminal_t *>(t); 420 if (otherterm != nullptr) 421 if (state().setup().get_connected_terminal(*otherterm)->net().is_rail_net()) 422 railterms++; 423 } 424 } 425 } 426 if (railterms == 0) 427 { 428 state().log().error(ME_SOLVER_NO_RAIL_TERMINAL(nets_in_grp)); 429 num_errors++; 430 } 431 } 432 if (num_errors > 0) 433 throw nl_exception(MF_SOLVER_CONSISTENCY_ERRORS(num_errors)); 434 435 436 // setup the solvers 437 for (auto & grp : splitter.groups) 438 { 439 solver_ptr ms; 440 pstring sname = plib::pfmt("Solver_{1}")(m_mat_solvers.size()); 441 params_uptr params = plib::make_unique<solver::solver_parameters_t, solver_arena>(*this, sname + ".", m_params); 442 443 switch (params->m_fp_type()) 444 { 445 case solver::matrix_fp_type_e::FLOAT: 446 if (!config::use_float_matrix()) 447 log().info("FPTYPE {1} not supported. Using DOUBLE", params->m_fp_type().name()); 448 ms = create_solvers<std::conditional_t<config::use_float_matrix::value, float, double>>(sname, params.get(), grp); 449 break; 450 case solver::matrix_fp_type_e::DOUBLE: 451 ms = create_solvers<double>(sname, params.get(), grp); 452 break; 453 case solver::matrix_fp_type_e::LONGDOUBLE: 454 if (!config::use_long_double_matrix()) 455 log().info("FPTYPE {1} not supported. Using DOUBLE", params->m_fp_type().name()); 456 ms = create_solvers<std::conditional_t<config::use_long_double_matrix::value, long double, double>>(sname, params.get(), grp); 457 break; 458 case solver::matrix_fp_type_e::FLOATQ128: 459 #if (NL_USE_FLOAT128) 460 ms = create_solvers<FLOAT128>(sname, params.get(), grp); 461 #else 462 log().info("FPTYPE {1} not supported. Using DOUBLE", params->m_fp_type().name()); 463 ms = create_solvers<double>(sname, params.get(), grp); 464 #endif 465 break; 466 } 467 468 log().verbose("Solver {1}", ms->name()); 469 log().verbose(" ==> {1} nets", grp.size()); 470 log().verbose(" has {1} dynamic elements", ms->dynamic_device_count()); 471 log().verbose(" has {1} timestep elements", ms->timestep_device_count()); 472 for (auto &n : grp) 473 { 474 log().verbose("Net {1}", n->name()); 475 for (const auto &t : state().core_terms(*n)) 476 { 477 log().verbose(" {1}", t->name()); 478 } 479 } 480 481 m_mat_params.push_back(std::move(params)); 482 m_mat_solvers.push_back(std::move(ms)); 483 } 484 485 } 486 NETLIB_NAME(solver)487 solver::static_compile_container NETLIB_NAME(solver)::create_solver_code(solver::static_compile_target target) 488 { 489 solver::static_compile_container mp; 490 for (auto & s : m_mat_solvers) 491 { 492 auto r = s->create_solver_code(target); 493 if (!r.first.empty()) // ignore solvers not supporting static compile 494 mp.push_back(r); 495 } 496 return mp; 497 } 498 NETLIB_NAME(solver)499 std::size_t NETLIB_NAME(solver)::get_solver_id(const solver::matrix_solver_t *net) const 500 { 501 for (std::size_t i=0; i < m_mat_solvers.size(); i++) 502 if (m_mat_solvers[i].get() == net) 503 return i; 504 return std::numeric_limits<std::size_t>::max(); 505 } 506 NETLIB_NAME(solver)507 solver::matrix_solver_t * NETLIB_NAME(solver)::solver_by_id(std::size_t id) const 508 { 509 return m_mat_solvers[id].get(); 510 } 511 512 513 NETLIB_DEVICE_IMPL(solver, "SOLVER", "FREQ") 514 515 } // namespace devices 516 } // namespace netlist 517