diff options
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_gmres.h')
-rw-r--r-- | src/lib/netlist/solver/nld_ms_gmres.h | 152 |
1 files changed, 69 insertions, 83 deletions
diff --git a/src/lib/netlist/solver/nld_ms_gmres.h b/src/lib/netlist/solver/nld_ms_gmres.h index 2ff515ebda7..836f6e61cbb 100644 --- a/src/lib/netlist/solver/nld_ms_gmres.h +++ b/src/lib/netlist/solver/nld_ms_gmres.h @@ -1,27 +1,24 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_ms_gmres.h - * - */ #ifndef NLD_MS_GMRES_H_ #define NLD_MS_GMRES_H_ +/// +/// \file nld_ms_gmres.h +/// + +#include "nld_matrix_solver_ext.h" #include "nld_ms_direct.h" #include "nld_solver.h" #include "plib/gmres.h" -#include "plib/mat_cr.h" #include "plib/parray.h" +#include "plib/pmatrix_cr.h" #include "plib/vector_ops.h" #include <algorithm> -#include <cmath> - -namespace netlist -{ -namespace devices +namespace netlist::solver { template <typename FT, int SIZE> @@ -31,26 +28,64 @@ namespace devices using float_type = FT; - /* Sort rows in ascending order. This should minimize fill-in and thus - * maximize the efficiency of the incomplete LUT. - * This is already preconditioning. - */ - matrix_solver_GMRES_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::PREFER_BAND_MATRIX, params, size) - //, m_ops(size, 2) - , m_ops(size, 4) + // Sort rows in ascending order. This should minimize fill-in and thus + // maximize the efficiency of the incomplete LUT. + // This is already preconditioning. + + matrix_solver_GMRES_t(devices::nld_solver &main_solver, const pstring &name, + matrix_solver_t::net_list_t &nets, + const solver::solver_parameters_t *params, + const std::size_t size) + : matrix_solver_direct_t<FT, SIZE>(main_solver, name, nets, params, size) + , m_ops(this->m_arena, size, 0) , m_gmres(size) { + const std::size_t iN = this->size(); + + std::vector<std::vector<unsigned>> fill(iN); + + for (std::size_t k=0; k<iN; k++) + { + fill[k].resize(iN, decltype(m_ops.m_mat)::FILL_INFINITY); + terms_for_net_t & row = this->m_terms[k]; + for (const auto &nz_j : row.m_nz) + { + fill[k][static_cast<matrix_type>(nz_j)] = 0; + } } - void vsetup(analog_net_t::list_t &nets) override; - unsigned vsolve_non_dynamic(const bool newton_raphson) override; + m_ops.build(fill); + this->log_fill(fill, m_ops.m_mat); + + // build pointers into the compressed row format matrix for each terminal + + for (std::size_t k=0; k<iN; k++) + { + std::size_t cnt = 0; + for (std::size_t j=0; j< this->m_terms[k].rail_start();j++) + { + for (std::size_t i = m_ops.m_mat.row_idx[k]; i<m_ops.m_mat.row_idx[k+1]; i++) + if (this->m_terms[k].m_connected_net_idx[j] == static_cast<int>(m_ops.m_mat.col_idx[i])) + { + this->m_mat_ptr[k][j] = &m_ops.m_mat.A[i]; + cnt++; + break; + } + } + nl_assert(cnt == this->m_terms[k].rail_start()); + this->m_mat_ptr[k][this->m_terms[k].rail_start()] = &m_ops.m_mat.A[m_ops.m_mat.diagonal[k]]; + } + } + + void upstream_solve_non_dynamic() override; private: - using mattype = typename plib::matrix_compressed_rows_t<FT, SIZE>::index_type; + using matrix_type = typename plib::pmatrix_cr<arena_type, FT, SIZE>::index_type; - plib::mat_precondition_ILU<FT, SIZE> m_ops; + //plib::mat_precondition_none<FT, SIZE> m_ops; + plib::mat_precondition_ILU<arena_type, FT, SIZE> m_ops; + //plib::mat_precondition_diag<FT, SIZE> m_ops; plib::gmres_t<FT, SIZE> m_gmres; }; @@ -59,87 +94,38 @@ namespace devices // ---------------------------------------------------------------------------------------- template <typename FT, int SIZE> - void matrix_solver_GMRES_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets) + void matrix_solver_GMRES_t<FT, SIZE>::upstream_solve_non_dynamic() { - matrix_solver_direct_t<FT, SIZE>::vsetup(nets); - const std::size_t iN = this->size(); - std::vector<std::vector<unsigned>> fill(iN); - - for (std::size_t k=0; k<iN; k++) - { - fill[k].resize(iN, decltype(m_ops.m_mat)::FILL_INFINITY); - terms_for_net_t * row = this->m_terms[k].get(); - for (const auto &nz_j : row->m_nz) - { - fill[k][static_cast<mattype>(nz_j)] = 0; - } - } - - m_ops.build(fill); - - /* build pointers into the compressed row format matrix for each terminal */ - - for (std::size_t k=0; k<iN; k++) - { - std::size_t cnt = 0; - for (std::size_t j=0; j< this->m_terms[k]->m_railstart;j++) - { - for (std::size_t i = m_ops.m_mat.row_idx[k]; i<m_ops.m_mat.row_idx[k+1]; i++) - if (this->m_terms[k]->m_connected_net_idx[j] == static_cast<int>(m_ops.m_mat.col_idx[i])) - { - this->m_mat_ptr[k][j] = &m_ops.m_mat.A[i]; - cnt++; - break; - } - } - nl_assert(cnt == this->m_terms[k]->m_railstart); - this->m_mat_ptr[k][this->m_terms[k]->m_railstart] = &m_ops.m_mat.A[m_ops.m_mat.diag[k]]; - } - } - - template <typename FT, int SIZE> - unsigned matrix_solver_GMRES_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson) - { - const std::size_t iN = this->size(); - - plib::parray<FT, SIZE> RHS(iN); - //float_type new_V[storage_N]; - - m_ops.m_mat.set_scalar(0.0); + m_ops.m_mat.set_scalar(plib::constants<FT>::zero()); - /* populate matrix and V for first estimate */ - this->fill_matrix(iN, this->m_mat_ptr, RHS); + // populate matrix and V for first estimate + this->fill_matrix_and_rhs(); for (std::size_t k = 0; k < iN; k++) { - this->m_new_V[k] = this->m_nets[k]->Q_Analog(); + this->m_new_V[k] = static_cast<float_type>(this->m_terms[k].getV()); } - const float_type accuracy = this->m_params.m_accuracy; + const float_type accuracy(static_cast<float_type>(this->m_params.m_accuracy)); - auto iter = std::max(plib::constants<std::size_t>::one(), this->m_params.m_gs_loops); - auto gsl = m_gmres.solve(m_ops, this->m_new_V, RHS, iter, accuracy); + auto iter = std::max(plib::constants<std::size_t>::one(), this->m_params.m_gs_loops()); + auto gsl = m_gmres.solve(m_ops, this->m_new_V, this->m_RHS, iter, accuracy); this->m_iterative_total += gsl; - this->m_stat_calculations++; if (gsl > iter) { this->m_iterative_fail++; - return matrix_solver_direct_t<FT, SIZE>::vsolve_non_dynamic(newton_raphson); + 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_GMRES_H_ */ +#endif // NLD_MS_GMRES_H_ |