Fix crash in stat output when no calculations were run

nl_examples/kidniki.c

@@ -14,9 +14,10 @@ NETLIST_START(dummy)
 	PARAM(Solver.NR_LOOPS, 300)
 	PARAM(Solver.GS_LOOPS, 1)
 	PARAM(Solver.GS_THRESHOLD, 6)
+	PARAM(Solver.ITERATIVE, "SM")
 	//PARAM(Solver.ITERATIVE, "MAT")
 	//PARAM(Solver.ITERATIVE, "GMRES")
-	PARAM(Solver.ITERATIVE, "SOR")
+	//PARAM(Solver.ITERATIVE, "SOR")
 	PARAM(Solver.DYNAMIC_TS, 0)
 	PARAM(Solver.DYNAMIC_LTE, 5e-3)
 	PARAM(Solver.MIN_TIMESTEP, 10e-6)

src/lib/netlist/solver/nld_matrix_solver.h

@@ -46,7 +46,7 @@ class terms_t
 	unsigned m_railstart;
 
 	pvector_t<unsigned> m_nz;   /* all non zero for multiplication */
-	pvector_t<unsigned> m_nzrd; /* non zero right of the diagonal for elimination, includes RHS element */
+	pvector_t<unsigned> m_nzrd; /* non zero right of the diagonal for elimination, may include RHS element */
 	pvector_t<unsigned> m_nzbd; /* non zero below of the diagonal for elimination */
 private:
 	pvector_t<terminal_t *> m_term;

src/lib/netlist/solver/nld_ms_direct.h

@@ -269,8 +269,6 @@ ATTR_COLD void matrix_solver_direct_t<m_N, _storage_N>::vsetup(analog_net_t::lis
 			}
 		}
 
-		for (unsigned j = 0; j < N(); j++)
-		{
 		for (unsigned i = 0; i < t->m_railstart; i++)
 		{
 			if (!t->m_nzrd.contains(other[i]) && other[i] >= (int) (k + 1))
@@ -278,7 +276,7 @@ ATTR_COLD void matrix_solver_direct_t<m_N, _storage_N>::vsetup(analog_net_t::lis
 			if (!t->m_nz.contains(other[i]))
 				t->m_nz.push_back(other[i]);
 		}
-		}
+
 		/* Add RHS element */
 		if (!t->m_nzrd.contains(N()))
 			t->m_nzrd.push_back(N());

src/lib/netlist/solver/nld_ms_sm.h

@@ -13,11 +13,6 @@
 #include "solver/nld_solver.h"
 #include "solver/vector_base.h"
 
-/* Disabling dynamic allocation gives a ~10% boost in performance
- * This flag has been added to support continuous storage for arrays
- * going forward in case we implement cuda solvers in the future.
- */
-
 NETLIB_NAMESPACE_DEVICES_START()
 
 //#define nl_ext_double __float128 // slow, very slow
@@ -273,8 +268,6 @@ ATTR_COLD void matrix_solver_sm_t<m_N, _storage_N>::vsetup(analog_net_t::list_t
 			}
 		}
 
-		for (unsigned j = 0; j < N(); j++)
-		{
 		for (unsigned i = 0; i < t->m_railstart; i++)
 		{
 			if (!t->m_nzrd.contains(other[i]) && other[i] >= (int) (k + 1))
@@ -282,7 +275,7 @@ ATTR_COLD void matrix_solver_sm_t<m_N, _storage_N>::vsetup(analog_net_t::list_t
 			if (!t->m_nz.contains(other[i]))
 				t->m_nz.push_back(other[i]);
 		}
-		}
+
 		/* and sort */
 		psort_list(t->m_nzrd);
 
@@ -476,8 +469,10 @@ void matrix_solver_sm_t<m_N, _storage_N>::LE_compute_x(
 
 	for (int k=0; k<kN; k++)
 	{
+		const nl_double f = RHS(k);
+
 		for (int i=0; i<kN; i++)
-			x[i] += Ainv(i,k) * RHS(k);
+			x[i] += Ainv(i,k) * f;
 	}
 }
 
@@ -514,11 +509,12 @@ void matrix_solver_sm_t<m_N, _storage_N>::store(
 template <unsigned m_N, unsigned _storage_N>
 int matrix_solver_sm_t<m_N, _storage_N>::solve_non_dynamic(ATTR_UNUSED const bool newton_raphson)
 {
+	static const bool incremental = true;
 	static uint cnt = 0;
 
 	nl_double new_V[_storage_N]; // = { 0.0 };
 
-	if (0 || (cnt % 100 == 0))
+	if (0 || ((cnt % 200) == 0))
 	{
 		/* complete calculation */
 		this->LE_invert();
@@ -526,40 +522,57 @@ int matrix_solver_sm_t<m_N, _storage_N>::solve_non_dynamic(ATTR_UNUSED const boo
 	else
 	{
 		const auto iN = N();
+
+		if (not incremental)
+		{
 			for (int row = 0; row < iN; row ++)
 				for (int k = 0; k < iN; k++)
 					Ainv(row,k) = lAinv(row, k);
+		}
+		for (int row = 0; row < iN; row ++)
+		{
+			nl_double v[m_pitch] = {0};
+			unsigned cols[m_pitch];
+			unsigned colcount = 0;
 
-		for (int row = 0; row < N(); row ++)
+			auto &nz = m_terms[row]->m_nz;
+			for (auto & col : nz)
 			{
-			nl_double v[m_pitch];
+				v[col] = A(row,col) - lA(row,col);
-			bool changed = false;
+				if (incremental)
-			for (int k = 0; k < N(); k++)
+					lA(row,col) = A(row,col);
-				v[k] = A(row,k) - lA(row,k);
+				if (v[col] != 0.0)
-			for (int k = 0; k < N(); k++)
+					cols[colcount++] = col;
-				if (v[k] != 0.0)
-				{
-					changed = true;
-					break;
 			}
 
-			if (changed)
+			if (colcount > 0)
 			{
 				nl_double lamba = 0.0;
 				nl_double w[m_pitch] = {0};
-				nl_double z[m_pitch] = {0};
+				nl_double z[m_pitch];
 				/* compute w and lamba */
-				for (int k = 0; k < N(); k++)
+				for (unsigned i = 0; i < iN; i++)
+					z[i] = Ainv(i, row); /* u is row'th column */
+
+				for (unsigned j = 0; j < colcount; j++)
+					lamba += v[cols[j]] * z[cols[j]];
+
+				for (unsigned j=0; j<colcount; j++)
 				{
-					for (int j=0; j<N(); j++)
+					auto col = cols[j];
-						w[k] += Ainv(j,k) * v[j]; /* Transpose(Ainv) * v */
+					auto f = v[col];
-					z[k] = Ainv(k, row); /* u is row'th column */
+					for (unsigned k = 0; k < iN; k++)
-					lamba += v[k] * z[k];
+						w[k] += Ainv(col,k) * f; /* Transpose(Ainv) * v */
 				}
+
 				lamba = -1.0 / (1.0 + lamba);
-				for (int i=0; i<N(); i++)
+				for (int i=0; i<iN; i++)
-					for (int k = 0; k < N(); k++)
+				{
-						Ainv(i,k) += lamba * z[i] * w[k];
+					const nl_double f = lamba * z[i];
+					if (f != 0.0)
+						for (int k = 0; k < iN; k++)
+							Ainv(i,k) += f * w[k];
+				}
 			}
 
 		}

src/lib/netlist/solver/nld_solver.cpp

@@ -309,7 +309,7 @@ ATTR_COLD int matrix_solver_t::get_net_idx(net_t *net)
 
 void matrix_solver_t::log_stats()
 {
-	//if (this->m_stat_calculations != 0 && this->m_params.m_log_stats)
+	if (this->m_stat_calculations != 0 && this->m_stat_vsolver_calls && this->m_params.m_log_stats)
 	{
 		log().verbose("==============================================");
 		log().verbose("Solver {1}", this->name());
@@ -462,6 +462,12 @@ matrix_solver_t * NETLIB_NAME(solver)::create_solver(int size, const bool use_sp
 				typedef matrix_solver_sm_t<m_N,_storage_N> solver_mat;
 				return palloc(solver_mat(&m_params, size));
 			}
+			else if (pstring("SM").equals(m_iterative_solver))
+			{
+				/* Sherman-Morrison Formula */
+				typedef matrix_solver_sm_t<m_N,_storage_N> solver_mat;
+				return palloc(solver_mat(&m_params, size));
+			}
 			else if (pstring("SOR").equals(m_iterative_solver))
 			{
 				typedef matrix_solver_SOR_t<m_N,_storage_N> solver_GS;