diff options
Diffstat (limited to 'src/lib/netlist/solver/nld_ms_w.h')
-rw-r--r-- | src/lib/netlist/solver/nld_ms_w.h | 585 |
1 files changed, 283 insertions, 302 deletions
diff --git a/src/lib/netlist/solver/nld_ms_w.h b/src/lib/netlist/solver/nld_ms_w.h index a01fc239938..b5b58f4ebe7 100644 --- a/src/lib/netlist/solver/nld_ms_w.h +++ b/src/lib/netlist/solver/nld_ms_w.h @@ -1,379 +1,360 @@ -// license:GPL-2.0+ +// license:BSD-3-Clause // copyright-holders:Couriersud -/* - * nld_ms_direct.h - * - * - * Woodbury Solver - * - * Computes the updated solution of A given that the change in A is - * - * A <- A + (U x transpose(V)) U,V matrices - * - * The approach is describes in "Numerical Recipes in C", Second edition, Page 75ff - * - * Whilst the book proposes to invert the matrix R=(I+transpose(V)*Z) we define - * - * w = transpose(V)*y - * a = R^-1 * w - * - * and consequently - * - * R * a = w - * - * And solve for a using Gaussian elimination. This is a lot faster. - * - * One fact omitted in the book is the fact that actually the matrix Z which contains - * in it's columns the solutions of - * - * A * zk = uk - * - * for uk being unit vectors for full rank (max(k) == n) is identical to the - * inverse of A. - * - * The approach performs relatively well for matrices up to n ~ 40 (kidniki using frontiers). - * Kidniki without frontiers has n==88. Here, the average number of Newton-Raphson - * loops increase to 20. It looks like that the approach for larger matrices - * introduces numerical instability. - */ #ifndef NLD_MS_W_H_ #define NLD_MS_W_H_ -#include "nld_matrix_solver.h" -#include "nld_solver.h" +// Names +// spell-checker: words Woodbury, Raphson, +// +// Specific technical terms +// spell-checker: words Cgso, Cgdo, Cgbo, Cjsw, Mjsw, Ucrit, Uexp, Utra, Neff, Tnom, capval, Udsat, Utst + + +/// +/// \file nld_ms_direct.h +/// +/// Woodbury Solver +/// +/// Computes the updated solution of A given that the change in A is +/// +/// A <- A + (U x transpose(V)) U,V matrices +/// +/// The approach is describes in "Numerical Recipes in C", Second edition, Page 75ff +/// +/// Whilst the book proposes to invert the matrix R=(I+transpose(V)*Z) we define +/// +/// w = transpose(V)*y +/// a = R^-1 * w +/// +/// and consequently +/// +/// R * a = w +/// +/// And solve for a using Gaussian elimination. This is a lot faster. +/// +/// One fact omitted in the book is the fact that actually the matrix Z which contains +/// in it's columns the solutions of +/// +/// A * zk = uk +/// +/// for uk being unit vectors for full rank (max(k) == n) is identical to the +/// inverse of A. +/// +/// The approach performs relatively well for matrices up to n ~ 40 (`kidniki` using frontiers). +/// `Kidniki` without frontiers has n==88. Here, the average number of Newton-Raphson +/// loops increase to 20. It looks like that the approach for larger matrices +/// introduces numerical instability. +/// + +#include "nld_matrix_solver_ext.h" #include "plib/vector_ops.h" #include <algorithm> -namespace netlist +namespace netlist::solver { - namespace devices - { - -template <typename FT, int SIZE> -class matrix_solver_w_t: public matrix_solver_t -{ - friend class matrix_solver_t; - -public: - using float_ext_type = FT; - using float_type = FT; - - // FIXME: dirty hack to make this compile - static constexpr const std::size_t storage_N = 100; - - matrix_solver_w_t(netlist_state_t &anetlist, const pstring &name, const solver_parameters_t *params, const std::size_t size); - - void vsetup(analog_net_t::list_t &nets) override; - void reset() override { matrix_solver_t::reset(); } - -protected: - unsigned vsolve_non_dynamic(const bool newton_raphson) override; - unsigned solve_non_dynamic(const bool newton_raphson); - - constexpr std::size_t size() const { return m_dim; } - - void LE_invert(); - - template <typename T> - void LE_compute_x(T & x); - - - template <typename T1, typename T2> - float_ext_type &A(const T1 &r, const T2 &c) { return m_A[r][c]; } - template <typename T1, typename T2> - float_ext_type &W(const T1 &r, const T2 &c) { return m_W[r][c]; } - - /* access to Ainv for fixed columns over row, there store transposed */ - template <typename T1, typename T2> - float_ext_type &Ainv(const T1 &r, const T2 &c) { return m_Ainv[c][r]; } - template <typename T1> - float_ext_type &RHS(const T1 &r) { return m_RHS[r]; } + template <typename FT, int SIZE> + class matrix_solver_w_t: public matrix_solver_ext_t<FT, SIZE> + { + public: + using float_ext_type = FT; + using float_type = FT; + + // FIXME: dirty hack to make this compile + static constexpr const std::size_t storage_N = 100; + + matrix_solver_w_t(devices::nld_solver &main_solver, const pstring &name, + const matrix_solver_t::net_list_t &nets, + const solver_parameters_t *params, const std::size_t size) + : matrix_solver_ext_t<FT, SIZE>(main_solver, name, nets, params, size) + , m_cnt(0) + { + this->build_mat_ptr(m_A); + } - template <typename T1, typename T2> - float_ext_type &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; } + void reset() override { matrix_solver_t::reset(); } + protected: + void upstream_solve_non_dynamic() override; -private: - static constexpr std::size_t m_pitch = ((( storage_N) + 7) / 8) * 8; - float_ext_type m_A[storage_N][m_pitch]; - float_ext_type m_Ainv[storage_N][m_pitch]; - float_ext_type m_W[storage_N][m_pitch]; - std::array<float_ext_type, storage_N> m_RHS; // right hand side - contains currents + void LE_invert(); - float_ext_type m_lA[storage_N][m_pitch]; + template <typename T> + void LE_compute_x(T & x); - /* temporary */ - float_type H[storage_N][m_pitch] ; - std::array<unsigned, storage_N> rows; - unsigned cols[storage_N][m_pitch]; - std::array<unsigned, storage_N> colcount; + template <typename T1, typename T2> + float_ext_type &A(const T1 &r, const T2 &c) { return m_A[r][c]; } + template <typename T1, typename T2> + float_ext_type &W(const T1 &r, const T2 &c) { return m_W[r][c]; } - unsigned m_cnt; + // access to the inverted matrix for fixed columns over row, values stored transposed + template <typename T1, typename T2> + float_ext_type &Ainv(const T1 &r, const T2 &c) { return m_Ainv[c][r]; } + template <typename T1> + float_ext_type &RHS(const T1 &r) { return this->m_RHS[r]; } - //float_ext_type m_RHSx[storage_N]; - const std::size_t m_dim; + template <typename T1, typename T2> + float_ext_type &lA(const T1 &r, const T2 &c) { return m_lA[r][c]; } -}; -// ---------------------------------------------------------------------------------------- -// matrix_solver_direct -// ---------------------------------------------------------------------------------------- + private: + void solve_non_dynamic(); -template <typename FT, int SIZE> -void matrix_solver_w_t<FT, SIZE>::vsetup(analog_net_t::list_t &nets) -{ - matrix_solver_t::setup_base(nets); + template <typename T, std::size_t N, std::size_t M> + using array2D = std::array<std::array<T, M>, N>; + static constexpr std::size_t m_pitch = ((( storage_N) + 7) / 8) * 8; + array2D<float_ext_type, storage_N, m_pitch> m_A; + array2D<float_ext_type, storage_N, m_pitch> m_Ainv; + array2D<float_ext_type, storage_N, m_pitch> m_W; - // FIXME: This shouldn't be necessary, recalculate on each entry ... - for (std::size_t k = 0; k < size(); k++) - state().save(*this, RHS(k), this->name(), plib::pfmt("RHS.{1}")(k)); -} + array2D<float_ext_type, storage_N, m_pitch> m_lA; + // temporary + array2D<float_ext_type, storage_N, m_pitch> H; + std::array<unsigned, storage_N> rows; + array2D<unsigned, storage_N, m_pitch> cols; + std::array<unsigned, storage_N> col_count; + unsigned m_cnt; + }; -template <typename FT, int SIZE> -void matrix_solver_w_t<FT, SIZE>::LE_invert() -{ - const std::size_t kN = size(); + // ---------------------------------------------------------------------------------------- + // matrix_solver_direct + // ---------------------------------------------------------------------------------------- - for (std::size_t i = 0; i < kN; i++) + template <typename FT, int SIZE> + void matrix_solver_w_t<FT, SIZE>::LE_invert() { - for (std::size_t j = 0; j < kN; j++) - { - W(i,j) = lA(i,j) = A(i,j); - Ainv(i,j) = 0.0; - } - Ainv(i,i) = 1.0; - } - /* down */ - for (std::size_t i = 0; i < kN; i++) - { - /* FIXME: Singular matrix? */ - const float_type f = 1.0 / W(i,i); - const auto * const p = m_terms[i]->m_nzrd.data(); - const size_t e = m_terms[i]->m_nzrd.size(); - - /* Eliminate column i from row j */ + const std::size_t kN = this->size(); - const auto * const pb = m_terms[i]->m_nzbd.data(); - const size_t eb = m_terms[i]->m_nzbd.size(); - for (std::size_t jb = 0; jb < eb; jb++) + for (std::size_t i = 0; i < kN; i++) { - const auto j = pb[jb]; - const float_type f1 = - W(j,i) * f; - if (f1 != 0.0) + for (std::size_t j = 0; j < kN; j++) { - for (std::size_t k = 0; k < e; k++) - W(j,p[k]) += W(i,p[k]) * f1; - for (std::size_t k = 0; k <= i; k ++) - Ainv(j,k) += Ainv(i,k) * f1; + W(i,j) = lA(i,j) = A(i,j); + Ainv(i,j) = plib::constants<FT>::zero(); } + Ainv(i,i) = plib::constants<FT>::one(); } - } - /* up */ - for (std::size_t i = kN; i-- > 0; ) - { - /* FIXME: Singular matrix? */ - const float_type f = 1.0 / W(i,i); - for (std::size_t j = i; j-- > 0; ) + // down + for (std::size_t i = 0; i < kN; i++) { - const float_type f1 = - W(j,i) * f; - if (f1 != 0.0) + // FIXME: Singular matrix? + const float_type f = plib::reciprocal(W(i,i)); + const auto * const p = this->m_terms[i].m_nzrd.data(); + const size_t e = this->m_terms[i].m_nzrd.size(); + + // Eliminate column i from row j + + const auto * const pb = this->m_terms[i].m_nzbd.data(); + const size_t eb = this->m_terms[i].m_nzbd.size(); + for (std::size_t jb = 0; jb < eb; jb++) { - for (std::size_t k = i; k < kN; k++) - W(j,k) += W(i,k) * f1; - for (std::size_t k = 0; k < kN; k++) - Ainv(j,k) += Ainv(i,k) * f1; + const auto j = pb[jb]; + const float_type f1 = - W(j,i) * f; + // FIXME: comparison to zero + if (f1 != plib::constants<float_type>::zero()) + { + for (std::size_t k = 0; k < e; k++) + W(j,p[k]) += W(i,p[k]) * f1; + for (std::size_t k = 0; k <= i; k ++) + Ainv(j,k) += Ainv(i,k) * f1; + } } } - for (std::size_t k = 0; k < kN; k++) + // up + for (std::size_t i = kN; i-- > 0; ) { - Ainv(i,k) *= f; + // FIXME: Singular matrix? + const float_type f = plib::reciprocal(W(i,i)); + for (std::size_t j = i; j-- > 0; ) + { + const float_type f1 = - W(j,i) * f; + // FIXME: comparison to zero + if (f1 != plib::constants<float_type>::zero()) + { + for (std::size_t k = i; k < kN; k++) + W(j,k) += W(i,k) * f1; + for (std::size_t k = 0; k < kN; k++) + Ainv(j,k) += Ainv(i,k) * f1; + } + } + for (std::size_t k = 0; k < kN; k++) + { + Ainv(i,k) *= f; + } } } -} - -template <typename FT, int SIZE> -template <typename T> -void matrix_solver_w_t<FT, SIZE>::LE_compute_x( - T & x) -{ - const std::size_t kN = size(); - - for (std::size_t i=0; i<kN; i++) - x[i] = 0.0; - for (std::size_t k=0; k<kN; k++) + template <typename FT, int SIZE> + template <typename T> + void matrix_solver_w_t<FT, SIZE>::LE_compute_x( + T & x) { - const float_type f = RHS(k); + const std::size_t kN = this->size(); for (std::size_t i=0; i<kN; i++) - x[i] += Ainv(i,k) * f; - } -} + x[i] = plib::constants<FT>::zero(); + for (std::size_t k=0; k<kN; k++) + { + const float_type f = RHS(k); -template <typename FT, int SIZE> -unsigned matrix_solver_w_t<FT, SIZE>::solve_non_dynamic(const bool newton_raphson) -{ - const auto iN = size(); + for (std::size_t i=0; i<kN; i++) + x[i] += Ainv(i,k) * f; + } + } - std::array<float_type, storage_N> new_V; // = { 0.0 }; - if ((m_cnt % 50) == 0) + template <typename FT, int SIZE> + void matrix_solver_w_t<FT, SIZE>::solve_non_dynamic() { - /* complete calculation */ - this->LE_invert(); - this->LE_compute_x(new_V); - } - else - { - /* Solve Ay = b for y */ - this->LE_compute_x(new_V); - - /* determine changed rows */ + const auto iN = this->size(); - unsigned rowcount=0; - #define VT(r,c) (A(r,c) - lA(r,c)) + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) + std::array<float_type, storage_N> t; // FIXME: convert to member + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) + std::array<float_type, storage_N> w; - for (unsigned row = 0; row < iN; row ++) + if ((m_cnt % 50) == 0) { - unsigned cc=0; - auto &nz = m_terms[row]->m_nz; - for (auto & col : nz) - { - if (A(row,col) != lA(row,col)) - cols[rowcount][cc++] = col; - } - if (cc > 0) - { - colcount[rowcount] = cc; - rows[rowcount++] = row; - } + // complete calculation + this->LE_invert(); + this->LE_compute_x(this->m_new_V); } - if (rowcount > 0) + else { - /* construct w = transform(V) * y - * dim: rowcount x iN - * */ - std::array<float_type, storage_N> w; - for (unsigned i = 0; i < rowcount; i++) + // Solve Ay = b for y + this->LE_compute_x(this->m_new_V); + + // determine changed rows + + unsigned row_count=0; + #define VT(r,c) (A(r,c) - lA(r,c)) + + for (unsigned row = 0; row < iN; row ++) { - const unsigned r = rows[i]; - double tmp = 0.0; - for (unsigned k = 0; k < iN; k++) - tmp += VT(r,k) * new_V[k]; - w[i] = tmp; + unsigned cc=0; + auto &nz = this->m_terms[row].m_nz; + for (auto & col : nz) + { + if (A(row,col) != lA(row,col)) + cols[row_count][cc++] = col; + } + if (cc > 0) + { + col_count[row_count] = cc; + rows[row_count++] = row; + } } + if (row_count > 0) + { + // construct w = transform(V) * y + // dim: row_count x iN + // + for (unsigned i = 0; i < row_count; i++) + { + const unsigned r = rows[i]; + FT tmp = plib::constants<FT>::zero(); + for (unsigned k = 0; k < iN; k++) + tmp += VT(r,k) * this->m_new_V[k]; + w[i] = tmp; + } + + for (unsigned i = 0; i < row_count; i++) + for (unsigned k=0; k< row_count; k++) + H[i][k] = plib::constants<FT>::zero(); - for (unsigned i = 0; i < rowcount; i++) - for (unsigned k=0; k< rowcount; k++) - H[i][k] = 0.0; + for (unsigned i = 0; i < row_count; i++) + H[i][i] = plib::constants<FT>::one(); + // Construct H = (I + VT*Z) + for (unsigned i = 0; i < row_count; i++) + for (unsigned k=0; k< col_count[i]; k++) + { + const unsigned col = cols[i][k]; + float_type f = VT(rows[i],col); + // FIXME: comparison to zero + if (f != plib::constants<float_type>::zero()) + for (unsigned j= 0; j < row_count; j++) + H[i][j] += f * Ainv(col,rows[j]); + } - for (unsigned i = 0; i < rowcount; i++) - H[i][i] = 1.0; - /* Construct H = (I + VT*Z) */ - for (unsigned i = 0; i < rowcount; i++) - for (unsigned k=0; k< colcount[i]; k++) + // Gaussian elimination of H + for (unsigned i = 0; i < row_count; i++) { - const unsigned col = cols[i][k]; - float_type f = VT(rows[i],col); - if (f!=0.0) - for (unsigned j= 0; j < rowcount; j++) - H[i][j] += f * Ainv(col,rows[j]); + // FIXME: comparison to zero + if (H[i][i] == plib::constants<float_type>::zero()) + plib::perrlogger("{} H singular\n", this->name()); + const float_type f = plib::reciprocal(H[i][i]); + for (unsigned j = i+1; j < row_count; j++) + { + const float_type f1 = - f * H[j][i]; + + // FIXME: comparison to zero + if (f1 != plib::constants<float_type>::zero()) + { + float_type *pj = &H[j][i+1]; + const float_type *pi = &H[i][i+1]; + for (unsigned k = 0; k < row_count-i-1; k++) + pj[k] += f1 * pi[k]; + //H[j][k] += f1 * H[i][k]; + w[j] += f1 * w[i]; + } + } } - - /* Gaussian elimination of H */ - for (unsigned i = 0; i < rowcount; i++) - { - if (H[i][i] == 0.0) - plib::perrlogger("{} H singular\n", this->name()); - const float_type f = 1.0 / H[i][i]; - for (unsigned j = i+1; j < rowcount; j++) + // Back substitution + //inv(H) w = t w = H t + for (unsigned j = row_count; j-- > 0; ) { - const float_type f1 = - f * H[j][i]; + float_type tmp = 0; + const float_type *pj = &H[j][j+1]; + const float_type *tj = &t[j+1]; + for (unsigned k = 0; k < row_count-j-1; k++) + tmp += pj[k] * tj[k]; + //tmp += H[j][k] * t[k]; + t[j] = (w[j] - tmp) / H[j][j]; + } - if (f1!=0.0) + // x = y - Zt + for (unsigned i=0; i<iN; i++) + { + float_type tmp = plib::constants<FT>::zero(); + for (unsigned j=0; j<row_count;j++) { - float_type *pj = &H[j][i+1]; - const float_type *pi = &H[i][i+1]; - for (unsigned k = 0; k < rowcount-i-1; k++) - pj[k] += f1 * pi[k]; - //H[j][k] += f1 * H[i][k]; - w[j] += f1 * w[i]; + const unsigned row = rows[j]; + tmp += Ainv(i,row) * t[j]; } + this->m_new_V[i] -= tmp; } } - /* Back substitution */ - //inv(H) w = t w = H t - std::array<float_type, storage_N> t; // FIXME: convert to member - for (unsigned j = rowcount; j-- > 0; ) - { - float_type tmp = 0; - const float_type *pj = &H[j][j+1]; - const float_type *tj = &t[j+1]; - for (unsigned k = 0; k < rowcount-j-1; k++) - tmp += pj[k] * tj[k]; - //tmp += H[j][k] * t[k]; - t[j] = (w[j] - tmp) / H[j][j]; - } + } + m_cnt++; - /* x = y - Zt */ + if (false) // NOLINT for (unsigned i=0; i<iN; i++) { - float_type tmp = 0.0; - for (unsigned j=0; j<rowcount;j++) + float_type tmp = plib::constants<FT>::zero(); + for (unsigned j=0; j<iN; j++) { - const unsigned row = rows[j]; - tmp += Ainv(i,row) * t[j]; + tmp += A(i,j) * this->m_new_V[j]; } - new_V[i] -= tmp; + if (plib::abs(tmp-RHS(i)) > static_cast<float_type>(1e-6)) + plib::perrlogger("{} failed on row {}: {} RHS: {}\n", this->name(), i, plib::abs(tmp-RHS(i)), RHS(i)); } - } } - m_cnt++; - - if (false) - for (unsigned i=0; i<iN; i++) - { - float_type tmp = 0.0; - for (unsigned j=0; j<iN; j++) - { - tmp += A(i,j) * new_V[j]; - } - if (std::abs(tmp-RHS(i)) > 1e-6) - plib::perrlogger("{} failed on row {}: {} RHS: {}\n", this->name(), i, std::abs(tmp-RHS(i)), RHS(i)); - } - const float_type err = (newton_raphson ? delta(new_V) : 0.0); - store(new_V); - return (err > this->m_params.m_accuracy) ? 2 : 1; -} + template <typename FT, int SIZE> + void matrix_solver_w_t<FT, SIZE>::upstream_solve_non_dynamic() + { + this->clear_square_mat(this->m_A); + this->fill_matrix_and_rhs(); -template <typename FT, int SIZE> -unsigned matrix_solver_w_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_raphson) -{ - this->build_LE_A(*this); - this->build_LE_RHS(*this); - - this->m_stat_calculations++; - return this->solve_non_dynamic(newton_raphson); -} - -template <typename FT, int SIZE> -matrix_solver_w_t<FT, SIZE>::matrix_solver_w_t(netlist_state_t &anetlist, const pstring &name, - const solver_parameters_t *params, const std::size_t size) - : matrix_solver_t(anetlist, name, NOSORT, params) - , m_cnt(0) - , m_dim(size) -{ -} + this->solve_non_dynamic(); + } - } //namespace devices -} // namespace netlist +} // namespace netlist::solver -#endif /* NLD_MS_DIRECT_H_ */ +#endif // NLD_MS_DIRECT_H_ |