diff options
author | 2014-01-19 14:14:05 +0000 | |
---|---|---|
committer | 2014-01-19 14:14:05 +0000 | |
commit | e897e114be4e1e7b15581fd0ff37f42bbaba51d6 (patch) | |
tree | 5c0cd9a0fbd5010806741eae98016cae85281a40 | |
parent | bcbd00db40543ed9bf3744be6c0270edd9834a40 (diff) |
Fixed direct solver bug. Using direct solver up to matrices of 4x4 now. for bigger matrices gauss-seidel is used. No wn.
-rw-r--r-- | src/emu/netlist/analog/nld_solver.c | 72 | ||||
-rw-r--r-- | src/emu/netlist/analog/nld_solver.h | 7 |
2 files changed, 53 insertions, 26 deletions
diff --git a/src/emu/netlist/analog/nld_solver.c b/src/emu/netlist/analog/nld_solver.c index 6333da2f65d..926536464cc 100644 --- a/src/emu/netlist/analog/nld_solver.c +++ b/src/emu/netlist/analog/nld_solver.c @@ -297,7 +297,7 @@ ATTR_HOT inline void netlist_matrix_solver_direct_t<m_N, _storage_N>::gauss_LE( #endif /* Singular matrix? */ double f = A[i][i]; - //if (fabs(A[i][i]) < 1e-20) printf("Singular!"); + //if (fabs(f) < 1e-20) printf("Singular!"); f = 1.0 / f; /* Eliminate column i from row j */ @@ -312,7 +312,6 @@ ATTR_HOT inline void netlist_matrix_solver_direct_t<m_N, _storage_N>::gauss_LE( RHS[j] -= RHS[i] * f1; } } - /* back substitution */ for (int j = N() - 1; j >= 0; j--) { @@ -322,12 +321,21 @@ ATTR_HOT inline void netlist_matrix_solver_direct_t<m_N, _storage_N>::gauss_LE( x[j] = (RHS[j] - tmp) / A[j][j]; } +#if 0 + for (int i = 0; i < N(); i++) + { + for (int k = 0; k < N(); k++) + printf("%f ", A[i][k]); + printf("| %f = %f \n", x[i], RHS[i]); + } + printf("\n"); +#endif } template <int m_N, int _storage_N> -ATTR_HOT inline double netlist_matrix_solver_direct_t<m_N, _storage_N>::delta_and_store( - double (* RESTRICT RHS), - double (* RESTRICT V)) +ATTR_HOT inline double netlist_matrix_solver_direct_t<m_N, _storage_N>::delta( + const double (* RESTRICT RHS), + const double (* RESTRICT V)) { double cerr = 0; double cerr2 = 0; @@ -337,13 +345,29 @@ ATTR_HOT inline double netlist_matrix_solver_direct_t<m_N, _storage_N>::delta_an double e2 = (RHS[i] - this->m_RHS[i]); cerr += e * e; cerr2 += e2 * e2; - this->m_nets[i]->m_cur.Analog = this->m_nets[i]->m_new.Analog = V[i]; - this->m_RHS[i] = RHS[i]; } return (cerr + cerr2*(100000.0 * 100000.0)) / this->N(); } template <int m_N, int _storage_N> +ATTR_HOT inline void netlist_matrix_solver_direct_t<m_N, _storage_N>::store( + const double (* RESTRICT RHS), + const double (* RESTRICT V)) +{ + for (int i = 0; i < this->N(); i++) + { + this->m_nets[i]->m_cur.Analog = this->m_nets[i]->m_new.Analog = V[i]; + } + if (RHS != NULL) + { + for (int i = 0; i < this->N(); i++) + { + this->m_RHS[i] = RHS[i]; + } + } +} + +template <int m_N, int _storage_N> ATTR_HOT int netlist_matrix_solver_direct_t<m_N, _storage_N>::solve_non_dynamic() { double A[_storage_N][_storage_N] = { { 0.0 } }; @@ -356,13 +380,17 @@ ATTR_HOT int netlist_matrix_solver_direct_t<m_N, _storage_N>::solve_non_dynamic( if (this->is_dynamic()) { - double err = delta_and_store(RHS, new_v); + double err = delta(RHS, new_v); + + store(RHS, new_v); if (err > this->m_accuracy * this->m_accuracy) { return 2; } + return 1; } + store(NULL, new_v); // ==> No need to store RHS return 1; } @@ -439,12 +467,17 @@ ATTR_HOT int netlist_matrix_solver_direct2_t::solve_non_dynamic() new_val[1] = a / (a*d - b*c) * (RHS[1] - c / a * RHS[0]); new_val[0] = (RHS[0] - b * new_val[1]) / a; - if (is_dynamic() && (delta_and_store(RHS, new_val) > m_accuracy * m_accuracy)) + if (is_dynamic()) { - return 2; + double err = delta(RHS, new_val); + store(RHS, new_val); + if (err > m_accuracy * m_accuracy) + return 2; + else + return 1; } - else - return 1; + store(NULL, new_val); + return 1; } // ---------------------------------------------------------------------------------------- @@ -729,20 +762,13 @@ ATTR_COLD void NETLIB_NAME(solver)::post_start() case 2: ms = new netlist_matrix_solver_direct2_t(); break; - /* FIXME: There is an issue (singularity during gaussian elemination - * with the direct solver - use gauss_seidel - */ case 3: - //ms = new netlist_matrix_solver_direct_t<3,3>(); - ms = new netlist_matrix_solver_gauss_seidel_t<3,3>(); + ms = new netlist_matrix_solver_direct_t<3,3>(); + //ms = new netlist_matrix_solver_gauss_seidel_t<3,3>(); break; case 4: - //ms = new netlist_matrix_solver_direct_t<4,4>(); - ms = new netlist_matrix_solver_gauss_seidel_t<4,4>(); - break; - case 5: - //ms = new netlist_matrix_solver_direct_t<5,5>(); - ms = new netlist_matrix_solver_gauss_seidel_t<5,5>(); + ms = new netlist_matrix_solver_direct_t<4,4>(); + //ms = new netlist_matrix_solver_gauss_seidel_t<4,4>(); break; default: //ms = new netlist_matrix_solver_direct_t<0,16>(); diff --git a/src/emu/netlist/analog/nld_solver.h b/src/emu/netlist/analog/nld_solver.h index 18aab7e0b9c..ad4a0b6e4db 100644 --- a/src/emu/netlist/analog/nld_solver.h +++ b/src/emu/netlist/analog/nld_solver.h @@ -86,9 +86,10 @@ protected: ATTR_HOT inline void gauss_LE(double (* RESTRICT A)[_storage_N], double (* RESTRICT RHS), double (* RESTRICT x)); - ATTR_HOT inline double delta_and_store( - double (* RESTRICT RHS), - double (* RESTRICT V)); + ATTR_HOT inline double delta( + const double (* RESTRICT RHS), + const double (* RESTRICT V)); + ATTR_HOT inline void store(const double (* RESTRICT RHS), const double (* RESTRICT V)); double m_RHS[_storage_N]; // right hand side - contains currents |