Move fillmatrix to to matrix solver base class. (nw)

src/lib/netlist/plib/mat_cr.h

@@ -67,13 +67,20 @@ namespace plib
 		, m_size(n)
 		{
 			for (index_type i=0; i<n+1; i++)
-				A[i] = 0;
+				row_idx[i] = 0;
 		}
 
 		~matrix_compressed_rows_t() = default;
 
 		constexpr index_type size() const { return static_cast<index_type>((N>0) ? N : m_size); }
 
+		void clear()
+		{
+			nz_num = 0;
+			for (index_type i=0; i < size() + 1; i++)
+				row_idx[i] = 0;
+		}
+
 		void set_scalar(const T scalar)
 		{
 			for (index_type i=0, e=nz_num; i<e; i++)
@@ -86,7 +93,7 @@ namespace plib
 			while (ri < row_idx[r+1] && col_idx[ri] < c)
 			  ri++;
 			// we have the position now;
-			if (nz_num > 0 && col_idx[ri] == c)
+			if (ri < row_idx[r+1] && col_idx[ri] == c)
 				A[ri] = val;
 			else
 			{
@@ -97,7 +104,7 @@ namespace plib
 				}
 				A[ri] = val;
 				col_idx[ri] = c;
-				for (C i = row_idx[r]; i < size()+1;i++)
+				for (C i = r + 1; i < size() + 1; i++)
 					row_idx[i]++;
 				nz_num++;
 				if (c==r)

src/lib/netlist/plib/parray.h

@@ -60,6 +60,7 @@ namespace plib {
 		using size_type = typename base_type::size_type;
 		using reference = typename base_type::reference;
 		using const_reference = typename base_type::const_reference;
+		using value_type = typename base_type::value_type;
 
 		template <int X = SIZE >
 		parray(size_type size, typename std::enable_if<X==0, int>::type = 0)

src/lib/netlist/solver/nld_matrix_solver.h

@@ -53,42 +53,6 @@ namespace devices
 
 		void set_pointers();
 
-		/* FIXME: this works a bit better for larger matrices */
-		template <typename AP, typename FT>
-		void fill_matrix/*_larger*/(AP &tcr, FT &RHS)
-		{
-
-			const std::size_t term_count = this->count();
-			const std::size_t railstart = this->m_railstart;
-
-			for (std::size_t i = 0; i < railstart; i++)
-				*tcr[i]       -= m_go[i];
-
-	#if 1
-			FT gtot_t = 0.0;
-			FT RHS_t = 0.0;
-
-			for (std::size_t i = 0; i < term_count; i++)
-			{
-				gtot_t        += m_gt[i];
-				RHS_t         += m_Idr[i];
-			}
-			// FIXME: Code above is faster than vec_sum - Check this
-	#else
-			auto gtot_t = plib::vec_sum<FT>(term_count, m_gt);
-			auto RHS_t = plib::vec_sum<FT>(term_count, m_Idr);
-	#endif
-
-			for (std::size_t i = railstart; i < term_count; i++)
-			{
-				RHS_t += (/*m_Idr[i]*/ + m_go[i] * *m_connected_net_V[i]);
-			}
-
-			RHS = RHS_t;
-			// update diagonal element ...
-			*tcr[railstart] += gtot_t; //mat.A[mat.diag[k]] += gtot_t;
-		}
-
 		std::size_t m_railstart;
 
 		std::vector<unsigned> m_nz;   /* all non zero for multiplication */
@@ -207,6 +171,46 @@ namespace devices
 		template <typename T>
 		void build_LE_RHS(T &child);
 
+		template <typename AP, typename FT>
+		void fill_matrix(std::size_t N, AP &tcr, FT &RHS)
+		{
+			for (std::size_t k = 0; k < N; k++)
+			{
+				auto *net = m_terms[k].get();
+				auto **tcr_r = tcr[k].data();
+
+				const std::size_t term_count = net->count();
+				const std::size_t railstart = net->m_railstart;
+
+				for (std::size_t i = 0; i < railstart; i++)
+					*tcr_r[i]       -= net->m_go[i];
+
+				typename FT::value_type gtot_t = 0.0;
+				typename FT::value_type RHS_t = 0.0;
+
+				for (std::size_t i = 0; i < term_count; i++)
+				{
+					gtot_t        += net->m_gt[i];
+					RHS_t         += net->m_Idr[i];
+				}
+				// FIXME: Code above is faster than vec_sum - Check this
+		#if 0
+				auto gtot_t = plib::vec_sum<FT>(term_count, m_gt);
+				auto RHS_t = plib::vec_sum<FT>(term_count, m_Idr);
+		#endif
+
+				for (std::size_t i = railstart; i < term_count; i++)
+				{
+					RHS_t += (/*m_Idr[i]*/ + net->m_go[i] * *net->m_connected_net_V[i]);
+				}
+
+				RHS[k] = RHS_t;
+				// update diagonal element ...
+				*tcr_r[railstart] += gtot_t; //mat.A[mat.diag[k]] += gtot_t;
+			}
+
+		}
+
 		std::vector<plib::unique_ptr<terms_for_net_t>> m_terms;
 		std::vector<analog_net_t *> m_nets;
 		std::vector<poolptr<proxied_analog_output_t>> m_inps;

src/lib/netlist/solver/nld_ms_gcr.h

@@ -296,8 +296,8 @@ unsigned matrix_solver_GCR_t<FT, SIZE>::vsolve_non_dynamic(const bool newton_rap
 	mat.set_scalar(0.0);
 
 	/* populate matrix */
-	for (std::size_t k = 0; k < iN; k++)
+
-		this->m_terms[k]->fill_matrix(m_term_cr[k], RHS[k]);
+	this->fill_matrix(iN, m_term_cr, RHS);
 
 	/* now solve it */
 

src/lib/netlist/solver/nld_ms_gmres.h

@@ -111,9 +111,10 @@ namespace devices
 		m_ops.m_mat.set_scalar(0.0);
 
 		/* populate matrix and V for first estimate */
+		this->fill_matrix(iN, m_term_cr, RHS);
+
 		for (std::size_t k = 0; k < iN; k++)
 		{
-			this->m_terms[k]->fill_matrix(m_term_cr[k], RHS[k]);
 			this->m_new_V[k] = this->m_nets[k]->Q_Analog();
 		}