Previous

199869 Revisions

Next

r30665 Sunday 25th May, 2014 at 21:45:51 UTC by Couriersud
Further simplification and clean up around the linear system solvers.

[src/emu/netlist]	nl_base.c nl_base.h plists.h
[src/emu/netlist/analog]	nld_solver.c nld_solver.h

trunk/src/emu/netlist/nl_base.c

r30664	r30665
771	771	, m_go(NETLIST_GMIN_DEFAULT)
772	772	, m_gt(NETLIST_GMIN_DEFAULT)
773	773	, m_otherterm(NULL)
774		, m_otherterm_ptr(NULL)
	774	, m_new_analog_ptr(NULL)
775	775	{
776	776	}
777	777

trunk/src/emu/netlist/nl_base.h

r30664	r30665
429	429	NETLIST_PREVENT_COPYING(netlist_terminal_t)
430	430	public:
431	431
432		typedef plinearlist_t<netlist_terminal_t *> list_t;
	432	typedef plinearlist_t<netlist_terminal_t * RESTRICT> list_t;
433	433
434	434	ATTR_COLD netlist_terminal_t();
435	435
r30664	r30665
462	462
463	463	// used by gauss-seidel-solver
464	464
465		double * RESTRICT m_otherterm_ptr;
	465	double * RESTRICT m_new_analog_ptr;
466	466
467	467	protected:
468	468	ATTR_COLD virtual void save_register();
r30664	r30665
732	732
733	733	//FIXME: needed by current solver code
734	734	netlist_matrix_solver_t *m_solver;
735		netlist_terminal_t::list_t m_terms;
736		netlist_terminal_t::list_t m_rails;
	735	//    netlist_terminal_t::list_t m_terms;
	736	//    netlist_terminal_t::list_t m_rails;
737	737	};
738	738
739	739	// ----------------------------------------------------------------------------------------

trunk/src/emu/netlist/analog/nld_solver.h

r30664	r30665
46	46	ATTR_COLD netlist_matrix_solver_t();
47	47	ATTR_COLD virtual ~netlist_matrix_solver_t();
48	48
49		ATTR_COLD virtual void vsetup(netlist_analog_net_t::list_t &nets, NETLIB_NAME(solver) &owner);
	49	ATTR_COLD virtual void vsetup(netlist_analog_net_t::list_t &nets,
	50	NETLIB_NAME(solver) &owner) = 0;
50	51
51	52	ATTR_HOT double solve();
52	53
r30664	r30665
64	65
65	66	netlist_solver_parameters_t m_params;
66	67
67		ATTR_COLD void log_stats();
	68	ATTR_COLD virtual void log_stats() {};
68	69
69	70	protected:
70	71
71		netlist_analog_net_t::list_t m_nets;
	72	class net_entry
	73	{
	74	NETLIST_PREVENT_COPYING(net_entry)
72	75
	76	public:
	77	net_entry(netlist_analog_net_t *net) : m_net(net) {}
	78	net_entry() : m_net(NULL) {}
73	79
	80	netlist_analog_net_t * RESTRICT m_net;
	81	netlist_terminal_t::list_t m_terms;
	82	netlist_terminal_t::list_t m_rails;
	83	};
	84
	85	ATTR_COLD virtual void setup(netlist_analog_net_t::list_t &nets,
	86	NETLIB_NAME(solver) &owner, net_entry *list);
	87
74	88	NETLIB_NAME(solver) *m_owner;
75	89
76	90	// return true if a reschedule is needed ...
77	91	ATTR_HOT virtual int vsolve_non_dynamic() = 0;
78	92
79	93	int m_calculations;
80		int m_gs_fail;
81		int m_gs_total;
82	94
83	95	private:
84	96
r30664	r30665
110	122
111	123	netlist_matrix_solver_direct_t()
112	124	: netlist_matrix_solver_t()
113		, m_term_num(0)
	125	, m_dim(0)
114	126	, m_rail_start(0)
115	127	{}
116	128
r30664	r30665
119	131	ATTR_COLD virtual void vsetup(netlist_analog_net_t::list_t &nets, NETLIB_NAME(solver) &owner);
120	132	ATTR_COLD virtual void reset() { netlist_matrix_solver_t::reset(); }
121	133
122		ATTR_HOT inline const int N() const { if (m_N == 0) return m_nets.count(); else return m_N; }
	134	ATTR_HOT inline const int N() const { if (m_N == 0) return m_dim; else return m_N; }
123	135
124	136	protected:
125	137	ATTR_HOT virtual int vsolve_non_dynamic();
126		ATTR_HOT int solve_non_dynamic(double (* RESTRICT A)[_storage_N], double (* RESTRICT RHS));
127		ATTR_HOT inline void build_LE(double (* RESTRICT A)[_storage_N], double (* RESTRICT RHS));
128		ATTR_HOT inline void gauss_LE(double (* RESTRICT A)[_storage_N],
129		double (* RESTRICT RHS),
130		double (* RESTRICT x));
	138	ATTR_HOT int solve_non_dynamic();
	139	ATTR_HOT inline void build_LE();
	140	ATTR_HOT inline void gauss_LE(double (* RESTRICT x));
131	141	ATTR_HOT inline double delta(
132		const double (* RESTRICT RHS),
133	142	const double (* RESTRICT V));
134		ATTR_HOT inline void store(const double (* RESTRICT RHS), const double (* RESTRICT V));
	143	ATTR_HOT inline void store(const double (* RESTRICT V), bool store_RHS);
135	144
136		double m_last_RHS[_storage_N]; // right hand side - contains currents
	145	net_entry m_nets[_storage_N];
137	146
	147	double m_A[_storage_N][_storage_N];
	148	double m_RHS[_storage_N];
	149	double m_last_RHS[_storage_N]; // right hand side - contains currents
	150
138	151	private:
139
140	152	ATTR_COLD int get_net_idx(netlist_net_t *net);
141	153
142	154	struct terms_t{
143		int net_this;
144		int net_other;
145		netlist_terminal_t *term;
	155
	156	terms_t(netlist_terminal_t *term, int net_this, int net_other)
	157	: m_term(term), m_net_this(net_this), m_net_other(net_other)
	158	{}
	159	terms_t()
	160	: m_term(NULL), m_net_this(-1), m_net_other(-1)
	161	{}
	162
	163	netlist_terminal_t *m_term;
	164	int m_net_this;
	165	int m_net_other;
146	166	};
147		int m_term_num;
	167
	168	int m_dim;
148	169	int m_rail_start;
149		terms_t m_terms[_storage_N * _storage_N];
	170	plinearlist_t<terms_t> m_terms;
150	171	};
151	172
152	173	template <int m_N, int _storage_N>
r30664	r30665
154	175	{
155	176	public:
156	177
157		netlist_matrix_solver_gauss_seidel_t() : netlist_matrix_solver_direct_t<m_N, _storage_N>() {}
	178	netlist_matrix_solver_gauss_seidel_t()
	179	: netlist_matrix_solver_direct_t<m_N, _storage_N>()
	180	, m_gs_fail(0)
	181	, m_gs_total(0)
	182	{}
158	183
159	184	virtual ~netlist_matrix_solver_gauss_seidel_t() {}
160	185
	186	ATTR_COLD virtual void log_stats();
	187
161	188	protected:
162	189	ATTR_HOT int vsolve_non_dynamic();
163	190
	191	private:
	192	int m_gs_fail;
	193	int m_gs_total;
	194
164	195	};
165	196
166	197	class netlist_matrix_solver_direct1_t: public netlist_matrix_solver_direct_t<1,1>

trunk/src/emu/netlist/analog/nld_solver.c

r30664	r30665
25	25	ATTR_COLD netlist_matrix_solver_t::netlist_matrix_solver_t()
26	26	: m_owner(NULL)
27	27	, m_calculations(0)
28		, m_gs_fail(0)
29		, m_gs_total(0)
30	28	{
31	29	}
32	30
r30664	r30665
36	34	delete m_inps[i];
37	35	}
38	36
39		ATTR_COLD void netlist_matrix_solver_t::vsetup(netlist_analog_net_t::list_t &nets, NETLIB_NAME(solver) &aowner)
	37	ATTR_COLD void netlist_matrix_solver_t::setup(netlist_analog_net_t::list_t &nets, NETLIB_NAME(solver) &aowner, net_entry *list)
40	38	{
41	39	m_owner = &aowner;
42	40
43	41	NL_VERBOSE_OUT(("New solver setup\n"));
44	42
45		for (netlist_analog_net_t * const * pn = nets.first(); pn != NULL; pn = nets.next(pn))
	43	for (int k = 0; k < nets.count(); k++)
	44	list[k].m_net = nets[k];
	45
	46	for (int k = 0; k < nets.count(); k++)
46	47	{
47	48	NL_VERBOSE_OUT(("setting up net\n"));
48	49
49		m_nets.add(*pn);
	50	netlist_analog_net_t *net = list[k].m_net;
50	51
51		(*pn)->m_solver = this;
	52	net->m_solver = this;
52	53
53		for (int i = 0; i < (*pn)->m_core_terms.count(); i++)
	54	for (int i = 0; i < net->m_core_terms.count(); i++)
54	55	{
55		netlist_core_terminal_t *p = (*pn)->m_core_terms[i];
56		NL_VERBOSE_OUT(("%s %s %d\n", p->name().cstr(), (*pn)->name().cstr(), (int) (*pn)->isRailNet()));
	56	netlist_core_terminal_t *p = net->m_core_terms[i];
	57	NL_VERBOSE_OUT(("%s %s %d\n", p->name().cstr(), net->name().cstr(), (int) net->isRailNet()));
57	58	switch (p->type())
58	59	{
59	60	case netlist_terminal_t::TERMINAL:
r30664	r30665
77	78	{
78	79	netlist_terminal_t *pterm = dynamic_cast<netlist_terminal_t *>(p);
79	80	// for gauss seidel
80		pterm->m_otherterm_ptr = &pterm->m_otherterm->net().as_analog().m_new_Analog;
	81	pterm->m_new_analog_ptr = &pterm->m_otherterm->net().as_analog().m_new_Analog;
81	82
82	83	if (pterm->m_otherterm->net().isRailNet())
83		(*pn)->m_rails.add(pterm);
	84	list[k].m_rails.add(pterm);
84	85	else
85		(*pn)->m_terms.add(pterm);
	86	list[k].m_terms.add(pterm);
86	87	}
87	88	NL_VERBOSE_OUT(("Added terminal\n"));
88	89	break;
r30664	r30665
114	115	break;
115	116	}
116	117	}
117		NL_VERBOSE_OUT(("added net with %d populated connections (%d railnets)\n", (*pn)->m_terms.count(), (*pn)->m_rails.count()));
	118	NL_VERBOSE_OUT(("added net with %d populated connections (%d railnets)\n", net->m_terms.count(), (*pn)->m_rails.count()));
118	119	}
119	120	}
120	121
r30664	r30665
268	269	return next_time_step;
269	270	}
270	271
271		void netlist_matrix_solver_t::log_stats()
	272	template <int m_N, int _storage_N>
	273	void netlist_matrix_solver_gauss_seidel_t<m_N, _storage_N>::log_stats()
272	274	{
273	275	#if 0
274	276	printf("==============================================\n");
275		printf("Solver %s\n", name().cstr());
276		printf("       ==> %d nets\n", m_nets.count()); //, (*(*groups[i].first())->m_core_terms.first())->name().cstr());
277		printf("       has %s elements\n", is_dynamic() ? "dynamic" : "no dynamic");
278		printf("       has %s elements\n", is_timestep() ? "timestep" : "no timestep");
	277	printf("Solver %s\n", this->name().cstr());
	278	printf("       ==> %d nets\n", this->N()); //, (*(*groups[i].first())->m_core_terms.first())->name().cstr());
	279	printf("       has %s elements\n", this->is_dynamic() ? "dynamic" : "no dynamic");
	280	printf("       has %s elements\n", this->is_timestep() ? "timestep" : "no timestep");
279	281	printf("       %10d invocations (%6d Hz)  %10d gs fails (%6.2f%%) %4.1f average\n",
280		m_calculations,
281		m_calculations * 10 / (int) (netlist().time().as_double() * 10.0),
282		m_gs_fail,
283		100.0 * (double) m_gs_fail / (double) m_calculations,
284		(double) m_gs_total / (double) m_calculations);
	282	this->m_calculations,
	283	this->m_calculations * 10 / (int) (this->netlist().time().as_double() * 10.0),
	284	this->m_gs_fail,
	285	100.0 * (double) this->m_gs_fail / (double) this->m_calculations,
	286	(double) this->m_gs_total / (double) this->m_calculations);
285	287	#endif
286	288	}
287	289
r30664	r30665
294	296	ATTR_COLD int netlist_matrix_solver_direct_t<m_N, _storage_N>::get_net_idx(netlist_net_t *net)
295	297	{
296	298	for (int k = 0; k < N(); k++)
297		if (m_nets[k] == net)
	299	if (m_nets[k].m_net == net)
298	300	return k;
299	301	return -1;
300	302	}
r30664	r30665
302	304	template <int m_N, int _storage_N>
303	305	ATTR_COLD void netlist_matrix_solver_direct_t<m_N, _storage_N>::vsetup(netlist_analog_net_t::list_t &nets, NETLIB_NAME(solver) &owner)
304	306	{
305		netlist_matrix_solver_t::vsetup(nets, owner);
	307	m_dim = nets.count();
	308	netlist_matrix_solver_t::setup(nets, owner, m_nets);
306	309
307		m_term_num = 0;
	310	m_terms.clear();
308	311	m_rail_start = 0;
309	312	for (int k = 0; k < N(); k++)
310	313	{
311		netlist_analog_net_t *net = m_nets[k];
312		const netlist_terminal_t::list_t &terms = net->m_terms;
	314	const netlist_terminal_t::list_t &terms = m_nets[k].m_terms;
313	315	for (int i = 0; i < terms.count(); i++)
314	316	{
315		m_terms[m_term_num].net_this = k;
316	317	int ot = get_net_idx(&terms[i]->m_otherterm->net());
317		m_terms[m_term_num].net_other = ot;
318		m_terms[m_term_num].term = terms[i];
319	318	if (ot>=0)
320	319	{
321		m_term_num++;
	320	m_terms.add(terms_t(terms[i], k, ot));
322	321	SOLVER_VERBOSE_OUT(("Net %d Term %s %f %f\n", k, terms[i]->name().cstr(), terms[i]->m_gt, terms[i]->m_go));
323	322	}
324	323	}
325	324	}
326		m_rail_start = m_term_num;
	325	m_rail_start = m_terms.count();
327	326	for (int k = 0; k < N(); k++)
328	327	{
329		netlist_analog_net_t *net = m_nets[k];
330		const netlist_terminal_t::list_t &terms = net->m_terms;
331		const netlist_terminal_t::list_t &rails = net->m_rails;
	328	const netlist_terminal_t::list_t &terms = m_nets[k].m_terms;
	329	const netlist_terminal_t::list_t &rails = m_nets[k].m_rails;
332	330	for (int i = 0; i < terms.count(); i++)
333	331	{
334		m_terms[m_term_num].net_this = k;
335		int ot = get_net_idx(&terms[i]->m_otherterm->net());
336		m_terms[m_term_num].net_other = ot;
337		m_terms[m_term_num].term = terms[i];
	332	int ot = get_net_idx(&terms[i]->m_otherterm->net());
338	333	if (ot<0)
339	334	{
340		m_term_num++;
341		SOLVER_VERBOSE_OUT(("found term with missing othernet %s\n", terms[i]->name().cstr()));
	335	m_terms.add(terms_t(terms[i], k, ot));
	336	netlist().warning("found term with missing othernet %s\n", terms[i]->name().cstr());
342	337	}
343	338	}
344	339	for (int i = 0; i < rails.count(); i++)
345	340	{
346		m_terms[m_term_num].net_this = k;
347		m_terms[m_term_num].net_other = -1; //get_net_idx(&rails[i]->m_otherterm->net());
348		m_terms[m_term_num].term = rails[i];
349		m_term_num++;
	341	m_terms.add(terms_t(rails[i], k, -1));
350	342	SOLVER_VERBOSE_OUT(("Net %d Rail %s %f %f\n", k, rails[i]->name().cstr(), rails[i]->m_gt, rails[i]->m_go));
351	343	}
352	344	}
353	345	}
354	346
355	347	template <int m_N, int _storage_N>
356		ATTR_HOT void netlist_matrix_solver_direct_t<m_N, _storage_N>::build_LE(
357		double (* RESTRICT A)[_storage_N],
358		double (* RESTRICT RHS))
	348	ATTR_HOT void netlist_matrix_solver_direct_t<m_N, _storage_N>::build_LE()
359	349	{
	350	for (int k=0; k < _storage_N; k++)
	351	for (int i=0; i < _storage_N; i++)
	352	m_A[k][i] = 0.0;
	353
	354	for (int k=0; k < _storage_N; k++)
	355	m_RHS[k] = 0.0;
360	356	#if 0
	357
361	358	for (int i = 0; i < m_term_num; i++)
362	359	{
363	360	terms_t &t = m_terms[i];
364		RHS[t.net_this] += t.term->m_Idr;
365		A[t.net_this][t.net_this] += t.term->m_gt;
366		if (t.net_other >= 0)
	361	m_RHS[t.m_net_this] += t.m_term->m_Idr;
	362	m_A[t.m_net_this][t.m_net_this] += t.m_term->m_gt;
	363	if (t.m_net_other >= 0)
367	364	{
368	365	//m_A[t.net_other][t.net_other] += t.term->m_otherterm->m_gt;
369		A[t.net_this][t.net_other] += -t.term->m_go;
	366	m_A[t.m_net_this][t.m_net_other] += -t.m_term->m_go;
370	367	//m_A[t.net_other][t.net_this] += -t.term->m_otherterm->m_go;
371	368	}
372	369	else
373		RHS[t.net_this] += t.term->m_go * t.term->m_otherterm->net().as_analog().Q_Analog();
	370	m_RHS[t.m_net_this] += t.m_term->m_go * t.m_term->m_otherterm->net().as_analog().Q_Analog();
374	371	}
375	372	#else
376	373	for (int i = 0; i < m_rail_start; i++)
r30664	r30665
378	375	const terms_t &t = m_terms[i];
379	376	//printf("A %d %d %s %f %f\n",t.net_this, t.net_other, t.term->name().cstr(), t.term->m_gt, t.term->m_go);
380	377
381		RHS[t.net_this] += t.term->m_Idr;
382		A[t.net_this][t.net_this] += t.term->m_gt;
383		A[t.net_this][t.net_other] += -t.term->m_go;
	378	m_RHS[t.m_net_this] += t.m_term->m_Idr;
	379	m_A[t.m_net_this][t.m_net_this] += t.m_term->m_gt;
	380	m_A[t.m_net_this][t.m_net_other] += -t.m_term->m_go;
384	381	}
385		for (int i = m_rail_start; i < m_term_num; i++)
	382	for (int i = m_rail_start; i < m_terms.count(); i++)
386	383	{
387	384	const terms_t &t = m_terms[i];
388	385
389		RHS[t.net_this] += t.term->m_Idr;
390		A[t.net_this][t.net_this] += t.term->m_gt;
391		RHS[t.net_this] += t.term->m_go * t.term->m_otherterm->net().as_analog().Q_Analog();
	386	m_RHS[t.m_net_this] += t.m_term->m_Idr;
	387	m_A[t.m_net_this][t.m_net_this] += t.m_term->m_gt;
	388	m_RHS[t.m_net_this] += t.m_term->m_go * t.m_term->m_otherterm->net().as_analog().Q_Analog();
392	389	}
393	390	#endif
394	391	}
395	392
396	393	template <int m_N, int _storage_N>
397	394	ATTR_HOT void netlist_matrix_solver_direct_t<m_N, _storage_N>::gauss_LE(
398		double (* RESTRICT A)[_storage_N],
399		double (* RESTRICT RHS),
400	395	double (* RESTRICT x))
401	396	{
402	397	#if 0
403	398	for (int i = 0; i < N(); i++)
404	399	{
405	400	for (int k = 0; k < N(); k++)
406		printf("%f ", A[i][k]);
407		printf("| %f = %f \n", x[i], RHS[i]);
	401	printf("%f ", m_A[i][k]);
	402	printf("| %f = %f \n", x[i], m_RHS[i]);
408	403	}
409	404	printf("\n");
410	405	#endif
r30664	r30665
419	414	int maxrow = i;
420	415	for (int j = i + 1; j < kN; j++)
421	416	{
422		if (fabs(A[j][i]) > fabs(A[maxrow][i]))
	417	if (fabs(m_A[j][i]) > fabs(m_A[maxrow][i]))
423	418	maxrow = j;
424	419	}
425	420
r30664	r30665
427	422	{
428	423	/* Swap the maxrow and ith row */
429	424	for (int k = i; k < kN; k++) {
430		std::swap(A[i][k], A[maxrow][k]);
	425	std::swap(m_A[i][k], m_A[maxrow][k]);
431	426	}
432		std::swap(RHS[i], RHS[maxrow]);
	427	std::swap(m_RHS[i], m_RHS[maxrow]);
433	428	}
434	429	}
435	430
436	431	/* Singular matrix? */
437		double f = A[i][i];
	432	double f = m_A[i][i];
438	433	//if (fabs(f) < 1e-20) printf("Singular!");
439	434	f = 1.0 / f;
440	435
r30664	r30665
443	438	for (int j = i + 1; j < kN; j++)
444	439	{
445	440	//__builtin_prefetch(&A[j+1][i], 1);
446		const double f1 = A[j][i] * f;
	441	const double f1 = m_A[j][i] * f;
447	442	if (f1 != 0.0)
448	443	{
449	444	for (int k = i; k < kN; k++)
450		A[j][k] -= A[i][k] * f1;
451		RHS[j] -= RHS[i] * f1;
	445	m_A[j][k] -= m_A[i][k] * f1;
	446	m_RHS[j] -= m_RHS[i] * f1;
452	447	}
453	448	}
454	449	}
r30664	r30665
458	453	//__builtin_prefetch(&A[j-1][j], 0);
459	454	double tmp = 0;
460	455	for (int k = j + 1; k < kN; k++)
461		tmp += A[j][k] * x[k];
462		x[j] = (RHS[j] - tmp) / A[j][j];
	456	tmp += m_A[j][k] * x[k];
	457	x[j] = (m_RHS[j] - tmp) / m_A[j][j];
463	458	}
464	459	#if 0
465	460	printf("Solution:\n");
466	461	for (int i = 0; i < N(); i++)
467	462	{
468	463	for (int k = 0; k < N(); k++)
469		printf("%f ", A[i][k]);
470		printf("| %f = %f \n", x[i], RHS[i]);
	464	printf("%f ", m_A[i][k]);
	465	printf("| %f = %f \n", x[i], m_RHS[i]);
471	466	}
472	467	printf("\n");
473	468	#endif
r30664	r30665
476	471
477	472	template <int m_N, int _storage_N>
478	473	ATTR_HOT double netlist_matrix_solver_direct_t<m_N, _storage_N>::delta(
479		const double (* RESTRICT RHS),
480	474	const double (* RESTRICT V))
481	475	{
482	476	double cerr = 0;
483	477	double cerr2 = 0;
484	478	for (int i = 0; i < this->N(); i++)
485	479	{
486		double e = (V[i] - this->m_nets[i]->m_cur_Analog);
487		double e2 = (RHS[i] - this->m_last_RHS[i]);
	480	double e = (V[i] - this->m_nets[i].m_net->m_cur_Analog);
	481	double e2 = (m_RHS[i] - this->m_last_RHS[i]);
488	482	cerr += e * e;
489	483	cerr2 += e2 * e2;
490	484	}
r30664	r30665
493	487
494	488	template <int m_N, int _storage_N>
495	489	ATTR_HOT void netlist_matrix_solver_direct_t<m_N, _storage_N>::store(
496		const double (* RESTRICT RHS),
497		const double (* RESTRICT V))
	490	const double (* RESTRICT V), bool store_RHS)
498	491	{
499	492	for (int i = 0; i < this->N(); i++)
500	493	{
501		this->m_nets[i]->m_cur_Analog = this->m_nets[i]->m_new_Analog = V[i];
	494	this->m_nets[i].m_net->m_cur_Analog = this->m_nets[i].m_net->m_new_Analog = V[i];
502	495	}
503		if (RHS != NULL)
	496	if (store_RHS)
504	497	{
505	498	for (int i = 0; i < this->N(); i++)
506	499	{
507		this->m_last_RHS[i] = RHS[i];
	500	this->m_last_RHS[i] = m_RHS[i];
508	501	}
509	502	}
510	503	}
511	504
512	505	template <int m_N, int _storage_N>
513		ATTR_HOT int netlist_matrix_solver_direct_t<m_N, _storage_N>::solve_non_dynamic(double (* RESTRICT A)[_storage_N], double (* RESTRICT RHS))
	506	ATTR_HOT int netlist_matrix_solver_direct_t<m_N, _storage_N>::solve_non_dynamic()
514	507	{
515	508	double new_v[_storage_N] = { 0.0 };
516	509
517		this->gauss_LE(A, RHS, new_v);
	510	this->gauss_LE(new_v);
518	511
519	512	if (this->is_dynamic())
520	513	{
521		double err = delta(RHS, new_v);
	514	double err = delta(new_v);
522	515
523		store(RHS, new_v);
	516	store(new_v, true);
524	517
525	518	if (err > this->m_params.m_accuracy * this->m_params.m_accuracy)
526	519	{
r30664	r30665
528	521	}
529	522	return 1;
530	523	}
531		store(NULL, new_v);  // ==> No need to store RHS
	524	store(new_v, false);  // ==> No need to store RHS
532	525	return 1;
533	526	}
534	527
535	528	template <int m_N, int _storage_N>
536	529	ATTR_HOT int netlist_matrix_solver_direct_t<m_N, _storage_N>::vsolve_non_dynamic()
537	530	{
538		double A[_storage_N][_storage_N] = { { 0.0 } };
539		double RHS[_storage_N] = { 0.0 };
	531	this->build_LE();
540	532
541		this->build_LE(A, RHS);
542
543		return this->solve_non_dynamic(A, RHS);
	533	return this->solve_non_dynamic();
544	534	}
545	535
546	536
r30664	r30665
551	541	ATTR_HOT int netlist_matrix_solver_direct1_t::vsolve_non_dynamic()
552	542	{
553	543
554		netlist_analog_net_t *net = m_nets[0];
555		double m_A[1][1] = { {0.0} };
556		double m_RHS[1] = { 0.0 };
557		build_LE(m_A, m_RHS);
	544	netlist_analog_net_t *net = m_nets[0].m_net;
	545	this->build_LE();
558	546	//NL_VERBOSE_OUT(("%f %f\n", new_val, m_RHS[0] / m_A[0][0]);
559	547
560	548	double new_val =  m_RHS[0] / m_A[0][0];
r30664	r30665
581	569
582	570	ATTR_HOT int netlist_matrix_solver_direct2_t::vsolve_non_dynamic()
583	571	{
584		double A[2][2] = { { 0.0 } };
585		double RHS[2] = { 0.0 };
586	572
587		build_LE(A, RHS);
	573	build_LE();
588	574
589		const double a = A[0][0];
590		const double b = A[0][1];
591		const double c = A[1][0];
592		const double d = A[1][1];
	575	const double a = m_A[0][0];
	576	const double b = m_A[0][1];
	577	const double c = m_A[1][0];
	578	const double d = m_A[1][1];
593	579
594	580	double new_val[2];
595		new_val[1] = a / (a * d - b * c) * (RHS[1] - c / a * RHS[0]);
596		new_val[0] = (RHS[0] - b * new_val[1]) / a;
	581	new_val[1] = a / (a * d - b * c) * (m_RHS[1] - c / a * m_RHS[0]);
	582	new_val[0] = (m_RHS[0] - b * new_val[1]) / a;
597	583
598	584	if (is_dynamic())
599	585	{
600		double err = delta(RHS, new_val);
601		store(RHS, new_val);
	586	double err = delta(new_val);
	587	store(new_val, true);
602	588	if (err > m_params.m_accuracy * m_params.m_accuracy)
603	589	return 2;
604	590	else
605	591	return 1;
606	592	}
607		store(NULL, new_val);
	593	store(new_val, false);
608	594	return 1;
609	595	}
610	596
r30664	r30665
621	607	*/
622	608
623	609	#if 0
624		double A[_storage_N][_storage_N] = { { 0.0 } };
625		double RHS[_storage_N] = { 0.0 };
626	610	double new_v[_storage_N] = { 0.0 };
627	611	const int iN = this->N();
628	612
r30664	r30665
630	614
631	615	int  resched_cnt = 0;
632	616
633		this->build_LE(A, RHS);
	617	this->build_LE();
634	618
635	619	for (int k = 0; k < iN; k++)
636	620	{
637		new_v[k] = this->m_nets[k]->m_cur_Analog;
	621	new_v[k] = this->m_nets[k].m_net->m_cur_Analog;
638	622	}
639	623	do {
640	624	resched = false;
r30664	r30665
647	631
648	632	// loop auto-vectorized
649	633	for (int i = 0; i < iN; i++)
650		Idrive -= A[pk][i] * new_v[i];
	634	Idrive -= this->m_A[pk][i] * new_v[i];
651	635
652		const double new_val = (RHS[pk] + Idrive + A[pk][pk] * new_v[pk]) / A[pk][pk];
	636	const double new_val = (this->m_RHS[pk] + Idrive + this->m_A[pk][pk] * new_v[pk]) / this->m_A[pk][pk];
653	637
654	638	const double e = (new_val - new_v[k]);
655	639	cerr += e * e;
r30664	r30665
668	652	{
669	653	//this->netlist().warning("Falling back to direct solver .. Consider increasing RESCHED_LOOPS");
670	654	this->m_gs_fail++;
671		int tmp = netlist_matrix_solver_direct_t<m_N, _storage_N>::solve_non_dynamic(A, RHS);
	655	int tmp = netlist_matrix_solver_direct_t<m_N, _storage_N>::solve_non_dynamic();
672	656	this->m_calculations++;
673	657	return tmp;
674	658	}
675	659	else {
676	660	this->m_calculations++;
677	661
678		this->store(NULL, new_v);
	662	this->store(new_v, false);
679	663
680	664	return resched_cnt;
681	665	}
r30664	r30665
697	681
698	682	for (int k = 0; k < iN; k++)
699	683	{
700		this->m_nets[k]->m_new_Analog = this->m_nets[k]->m_cur_Analog;
	684	this->m_nets[k].m_net->m_new_Analog = this->m_nets[k].m_net->m_cur_Analog;
701	685	}
702	686
703	687	for (int k = 0; k < iN; k++)
r30664	r30665
706	690	double gabs_t = 0.0;
707	691	double RHS_t = 0.0;
708	692
709		const netlist_analog_net_t &net = *this->m_nets[k];
710		const netlist_terminal_t::list_t &terms = net.m_terms;
711		const netlist_terminal_t::list_t &rails = net.m_rails;
	693	//const netlist_analog_net_t &net = *this->m_nets[k];
	694	const netlist_terminal_t::list_t &terms = this->m_nets[k].m_terms;
	695	const netlist_terminal_t::list_t &rails = this->m_nets[k].m_rails;
712	696	const int term_count = terms.count();
713	697	const int rail_count = rails.count();
714	698
r30664	r30665
717	701	gtot_t += rails[i]->m_gt;
718	702	gabs_t += fabs(rails[i]->m_go);
719	703	RHS_t += rails[i]->m_Idr;
720		//RHS_t += rails[i]->m_go * rails[i]->m_otherterm->net().as_analog().Q_Analog();
721		RHS_t += rails[i]->m_go * *rails[i]->m_otherterm_ptr;
	704	RHS_t += rails[i]->m_go * rails[i]->m_otherterm->net().as_analog().Q_Analog();
722	705	}
723	706
724	707	for (int i = 0; i < term_count; i++)
r30664	r30665
731	714	gabs_t *= 1.0;
732	715	if (gabs_t > gtot_t)
733	716	{
734		//this->netlist().warning("Falling back to direct solver .. Consider increasing RESCHED_LOOPS");
735	717	w[k] = 1.0 / (gtot_t + gabs_t);
736	718	one_m_w[k] = 1.0 - 1.0 * gtot_t / (gtot_t + gabs_t);
737	719	}
738	720	else
739	721	{
740		w[k] = 1.0/ gtot_t;
741		one_m_w[k] = 1.0 - 1.0;
	722	const double ws = 1.0;
	723	w[k] = ws / gtot_t;
	724	one_m_w[k] = 1.0 - ws;
742	725	}
743	726
744	727	RHS[k] = RHS_t;
r30664	r30665
750	733
751	734	for (int k = 0; k < iN; k++)
752	735	{
753		netlist_analog_net_t &net = *this->m_nets[k];
754		const netlist_terminal_t::list_t &terms = net.m_terms;
	736	netlist_analog_net_t & RESTRICT net = *this->m_nets[k].m_net;
	737	const netlist_terminal_t::list_t &terms = this->m_nets[k].m_terms;
755	738	const int term_count = terms.count();
756	739	double Idrive = 0;
757	740
758		for (int i = 0; i < term_count; i++)
759		Idrive += terms[i]->m_go * *(terms[i]->m_otherterm_ptr);
	741	for (int i = 0; i < term_count; i++)
	742	Idrive += terms[i]->m_go * *(terms[i]->m_new_analog_ptr);
760	743
761		//double new_val = (net->m_cur_Analog * gabs[k] + iIdr) / (gtot[k]);
	744	//double new_val = (net->m_cur_Analog * gabs[k] + iIdr) / (gtot[k]);
762	745	const double new_val = net.m_new_Analog * one_m_w[k] + (Idrive + RHS[k]) * w[k];
763	746
764	747	const double e = (new_val - net.m_new_Analog);
r30664	r30665
786	769	this->m_calculations++;
787	770
788	771	for (int k = 0; k < this->N(); k++)
789		this->m_nets[k]->m_cur_Analog = this->m_nets[k]->m_new_Analog;
	772	this->m_nets[k].m_net->m_cur_Analog = this->m_nets[k].m_net->m_new_Analog;
790	773
791	774	return resched_cnt;
792	775	}

trunk/src/emu/netlist/plists.h

r30664	r30665
61	61	{
62	62	if (m_list != NULL)
63	63	delete[] m_list;
	64	m_list = NULL;
64	65	}
65	66
66	67	ATTR_HOT inline void add(const _ListClass &elem)

Previous

199869 Revisions

Next

---

© 1997-2024 The MAME Team