fixes

libsolutil/LP.cpp

@@ -51,6 +51,7 @@ using namespace solidity::util;
 
 using rational = boost::rational<bigint>;
 
+#define DEBUG
 
 namespace
 {
@@ -279,7 +280,7 @@ void LPSolver::combineSubProblems(size_t _combineInto, size_t _combineFrom)
 
 	for (size_t row = 0; row < combineFrom.factors.rows(); row++)
 		for (auto&& entry: combineFrom.factors.iterateRow(row))
-			combineInto.factors.insert(entry.row + rowShift, entry.col + colShift, move(entry.value));
+			combineInto.factors.entry(entry.row + rowShift, entry.col + colShift).value = move(entry.value);
 #else
 	size_t rowShift = combineInto.factors.size();
 	size_t newRowLength = combineInto.variables.size() + combineFrom.variables.size();
@@ -362,7 +363,9 @@ void LPSolver::addConstraintToSubProblem(
 	// Introduce the slack variable.
 	size_t slackIndex = addNewVariableToSubProblem(_subProblem);
 	// Name is only needed for printing
-	//problem.variables[slackIndex].name = "_s" + to_string(m_slackVariableCounter++);
+#ifdef DEBUG
+	problem.variables[slackIndex].name = "_s" + to_string(m_slackVariableCounter++);
+#endif
 #ifdef SPARSE
 	problem.basicVariables[slackIndex] = problem.factors.rows();
 #else
@@ -398,8 +401,7 @@ void LPSolver::addConstraintToSubProblem(
 					row,
 					entry
 				);
-				// TODO could remove
+				problem.factors.remove(problem.factors.entry(row, innerIndex));
-				problem.factors.entry(row, innerIndex).value = 0;
 			}
 		}
 	for (auto const& [outerIndex, entry]: _constraint.data.enumerateTail())
@@ -407,7 +409,12 @@ void LPSolver::addConstraintToSubProblem(
 		{
 			size_t innerIndex = problem.varMapping.at(outerIndex);
 			if (!problem.basicVariables.count(innerIndex))
-				problem.factors.entry(row, innerIndex).value += entry;
+			{
+				SparseMatrix::Entry& e = problem.factors.entry(row, innerIndex);
+				e.value += entry;
+				if (!e.value)
+					problem.factors.remove(e);
+			}
 			valueForSlack += problem.variables[innerIndex].value * entry;
 		}
 
@@ -506,6 +513,7 @@ LPResult LPSolver::SubProblem::check()
 			{
 #ifdef DEBUG
 				cerr << "Replacing by " << variables[*replacementVar].name << endl;
+				cerr << "Setting basic var to to " << basicVar.bounds.lower->m_main << endl;
 #endif
 
 				pivotAndUpdate(*bvi, *basicVar.bounds.lower, *replacementVar);
@@ -575,6 +583,7 @@ string LPSolver::SubProblem::toString() const
 		for (auto&& entry: const_cast<SparseMatrix&>(factors).iterateRow(rowIndex))
 		{
 			rational const& f = entry.value;
+			solAssert(!!f);
 			size_t i = entry.col;
 #else
 		for (auto&& [i, f]: row.enumerate())
@@ -755,9 +764,14 @@ void LPSolver::SubProblem::pivot(size_t _old, size_t _new)
 	if (pivot != -1)
 		factors.multiplyRowByFactor(pivotRow, rational{-1} / pivot);
 
-	for (auto&& entry: factors.iterateColumn(_new))
+	for (auto it = factors.iterateColumn(_new).begin(); it != factors.iterateColumn(_new).end(); )
+	{
+		SparseMatrix::Entry& entry = *it;
+		// Increment becasue "addMultipleOfRow" might invalidate the iterator
+		++it;
 		if (entry.row != pivotRow)
 			factors.addMultipleOfRow(pivotRow, entry.row, entry.value);
+	}
 #else
 	LinearExpression& pivotRowData = factors[pivotRow];
 
@@ -794,7 +808,7 @@ void LPSolver::SubProblem::pivotAndUpdate(
 )
 {
 #ifdef SPARSE
-	RationalWithDelta theta = (_newValue - variables[_oldBasicVar].value) / factors.entry(_oldBasicVar, _newBasicVar).value;
+	RationalWithDelta theta = (_newValue - variables[_oldBasicVar].value) / factors.entry(basicVariables[_oldBasicVar], _newBasicVar).value;
 #else
 	RationalWithDelta theta = (_newValue - variables[_oldBasicVar].value) / factors[basicVariables[_oldBasicVar]][_newBasicVar];
 #endif

libsolutil/LP.h

@@ -19,6 +19,7 @@
 
 // use sparse matrices
 #define SPARSE 1
+#define DEBUG
 
 #include <libsolutil/Numeric.h>
 #include <libsolutil/LinearExpression.h>

libsolutil/LinearExpression.cpp

@@ -71,13 +71,10 @@ void SparseMatrix::multiplyRowByFactor(size_t _row, rational const& _factor)
 	}
 }
 
-void SparseMatrix::addMultipleOfRow(size_t _sourceRow, size_t _targetRow, rational const& _factor)
+void SparseMatrix::addMultipleOfRow(size_t _sourceRow, size_t _targetRow, rational const _factor)
 {
-	if (_targetRow >= m_row_start.size())
+	ensureSize(_targetRow, 0);
-	{
+	solAssert(_sourceRow != _targetRow);
-		m_row_start.resize(_targetRow + 1);
-		m_row_end.resize(_targetRow + 1);
-	}
 	Entry* source = m_row_start[_sourceRow];
 	Entry* target = m_row_start[_targetRow];
 
@@ -94,8 +91,13 @@ void SparseMatrix::addMultipleOfRow(size_t _sourceRow, size_t _targetRow, ration
 				remove(*target);
 				target = next;
 			}
+			else
+				target = target->next_in_row;
 		}
-		target = prependInRow(target, _targetRow, source->col, _factor * source->value)->next_in_row;
+		else if (rational newValue = _factor * source->value)
+			target = prependInRow(target, _targetRow, source->col, newValue)->next_in_row;
+		else if (target)
+			target = target->next_in_row;
 
 		source = source->next_in_row;
 	}
@@ -111,18 +113,35 @@ SparseMatrix::Entry& SparseMatrix::entry(size_t _row, size_t _column)
 		return *prependInRow(successor, _row, _column, {});
 }
 
-void SparseMatrix::insert(size_t _row, size_t _column, rational _value)
+void SparseMatrix::remove(SparseMatrix::Entry& _e)
 {
-	ensureSize(_row, _column);
+	// TODO this does not deallocate the entry.
-	prependInRow(entryOrSuccessorInRow(_row, _column), _row, _column, move(_value));
+	// At some point we should perform cleanup
+	// and swap entries in the vector
+	if (_e.prev_in_row)
+		_e.prev_in_row->next_in_row = _e.next_in_row;
+	else
+		m_row_start[_e.row] = _e.next_in_row;
+	if (_e.next_in_row)
+		_e.next_in_row->prev_in_row = _e.prev_in_row;
+	else
+		m_row_end[_e.row] = _e.prev_in_row;
+	if (_e.prev_in_col)
+		_e.prev_in_col->next_in_col = _e.next_in_col;
+	else
+		m_col_start[_e.col] = _e.next_in_col;
+	if (_e.next_in_col)
+		_e.next_in_col->prev_in_col = _e.prev_in_col;
+	else
+		m_col_end[_e.col] = _e.prev_in_col;
 }
 
 void SparseMatrix::appendRow(LinearExpression const& _entries)
 {
-	m_row_start.push_back(nullptr);
+	size_t row_nr = rows();
-	m_row_end.push_back(nullptr);
+	ensureSize(row_nr, 0);
-	size_t row_nr = m_row_start.size() - 1;
+	for (auto&& [i, v]: _entries.enumerate())
-	for (auto&& [i, v]: _entries.enumerate()) {
+	{
 		if (!v)
 			continue;
 		prependInRow(nullptr, row_nr, i, move(v));
@@ -156,26 +175,6 @@ SparseMatrix::Entry* SparseMatrix::entryOrSuccessorInRow(size_t _row, size_t _co
 	return successor;
 }
 
-void SparseMatrix::remove(SparseMatrix::Entry& _e)
-{
-	if (_e.prev_in_row)
-		_e.prev_in_row->next_in_row = _e.next_in_row;
-	else
-		m_row_start[_e.row] = _e.next_in_row;
-	if (_e.next_in_row)
-		_e.next_in_row->prev_in_row = _e.prev_in_row;
-	else
-		m_row_end[_e.row] = _e.prev_in_row;
-	if (_e.prev_in_col)
-		_e.prev_in_col->next_in_col = _e.next_in_col;
-	else
-		m_col_start[_e.col] = _e.next_in_col;
-	if (_e.next_in_col)
-		_e.next_in_col->prev_in_col = _e.prev_in_col;
-	else
-		m_col_end[_e.col] = _e.prev_in_col;
-}
-
 SparseMatrix::Entry* SparseMatrix::prependInRow(Entry* _successor, size_t _row, size_t _column, rational _value)
 {
 	m_elements.emplace_back(make_unique<Entry>(Entry{
@@ -211,11 +210,7 @@ void SparseMatrix::adjustColumnProperties(Entry& _entry)
 {
 	size_t column = _entry.col;
 
-	if (column >= m_col_start.size())
+	solAssert(column < m_col_start.size());
-	{
-		m_col_start.resize(column + 1);
-		m_col_end.resize(column + 1);
-	}
 	Entry* c = nullptr;
 	if (m_col_end[column] && m_col_end[column]->row > _entry.row)
 	{

libsolutil/LinearExpression.h

@@ -247,11 +247,9 @@ public:
 	/// @returns Entry for all non-zero v in the row _row
 	IteratorCombiner iterateRow(size_t _row);
 	void multiplyRowByFactor(size_t _row, rational const& _factor);
-	void addMultipleOfRow(size_t _sourceRow, size_t _targetRow, rational const& _factor);
+	void addMultipleOfRow(size_t _sourceRow, size_t _targetRow, rational _factor);
 	Entry& entry(size_t _row, size_t _column);
-	/// Inserts the value at the row/rolumn.
+	void remove(Entry& _entry);
-	/// Assumes the entry does not exist yet.
-	void insert(size_t _row, size_t _column, rational _value);
 
 	void appendRow(LinearExpression const& _entries);
 
@@ -260,7 +258,6 @@ private:
 	/// @returns the entry at the row/column if it exists or its successor in the row.
 	Entry* entryOrSuccessorInRow(size_t _row, size_t _column);
 
-	void remove(Entry& _entry);
 	/// Prepends a new entry before the given element or at end of row if nullptr.
 	Entry* prependInRow(Entry* _successor, size_t _row, size_t _column, rational _value);
 	void adjustColumnProperties(Entry& _entry);