impl

libsolutil/LP.cpp

@@ -270,15 +270,16 @@ void LPSolver::combineSubProblems(size_t _combineInto, size_t _combineFrom)
 	//cerr << "Combining\n" << m_subProblems.at(_combineFrom)->toString();
 	//cerr << "\ninto\n" << m_subProblems.at(_combineInto)->toString();
 	SubProblem& combineInto = unseal(_combineInto);
-	SubProblem const& combineFrom = *m_subProblems[_combineFrom];
+	SubProblem& combineFrom = *m_subProblems[_combineFrom];
 
 	size_t varShift = combineInto.variables.size();
 #ifdef SPARSE
 	size_t rowShift = combineInto.factors.rows();
+	size_t colShift = combineInto.factors.columns();
 
 	for (size_t row = 0; row < combineFrom.factors.rows(); row++)
-		for (auto&& [col, v]: combineFrom.factors.iterateRow(row))
+		for (auto&& entry: combineFrom.factors.iterateRow(row))
-			combineinto.factors.insert(row + rowShift, col + colShift, move(v));
+			combineInto.factors.insert(entry.row + rowShift, entry.col + colShift, move(entry.value));
 #else
 	size_t rowShift = combineInto.factors.size();
 	size_t newRowLength = combineInto.variables.size() + combineFrom.variables.size();
@@ -362,7 +363,11 @@ void LPSolver::addConstraintToSubProblem(
 	size_t slackIndex = addNewVariableToSubProblem(_subProblem);
 	// Name is only needed for printing
 	//problem.variables[slackIndex].name = "_s" + to_string(m_slackVariableCounter++);
+#ifdef SPARSE
+	problem.basicVariables[slackIndex] = problem.factors.rows();
+#else
 	problem.basicVariables[slackIndex] = problem.factors.size();
+#endif
 	if (_constraint.kind == Constraint::EQUAL)
 	{
 		problem.variables[slackIndex].bounds.lower = _constraint.data[0];
@@ -375,6 +380,41 @@ void LPSolver::addConstraintToSubProblem(
 	// TODO it is a basic var, so we don't add it, unless we use this for basic vars.
 	//problem.variablesPotentiallyOutOfBounds.insert(slackIndex);
 
+#ifdef SPARSE
+	// Compress the constraint, i.e. turn outer variable indices into
+	// inner variable indices.
+	RationalWithDelta valueForSlack;
+	size_t row = problem.factors.rows();
+	// First, handle the basic variables.
+	LinearExpression basicVarNullifier;
+	for (auto const& [outerIndex, entry]: _constraint.data.enumerateTail())
+		if (entry)
+		{
+			size_t innerIndex = problem.varMapping.at(outerIndex);
+			if (problem.basicVariables.count(innerIndex))
+			{
+				problem.factors.addMultipleOfRow(
+					problem.basicVariables[innerIndex],
+					row,
+					entry
+				);
+				// TODO could remove
+				problem.factors.entry(row, innerIndex).value = 0;
+			}
+		}
+	for (auto const& [outerIndex, entry]: _constraint.data.enumerateTail())
+		if (entry)
+		{
+			size_t innerIndex = problem.varMapping.at(outerIndex);
+			if (!problem.basicVariables.count(innerIndex))
+				problem.factors.entry(row, innerIndex).value += entry;
+			valueForSlack += problem.variables[innerIndex].value * entry;
+		}
+
+	problem.factors.entry(row, slackIndex).value = -1;
+
+	problem.basicVariables[slackIndex] = row;
+#else
 	// Compress the constraint, i.e. turn outer variable indices into
 	// inner variable indices.
 	RationalWithDelta valueForSlack;
@@ -399,7 +439,9 @@ void LPSolver::addConstraintToSubProblem(
 
 	compressedConstraint[slackIndex] = -1;
 	problem.factors.emplace_back(move(compressedConstraint));
+
 	problem.basicVariables[slackIndex] = problem.factors.size() - 1;
+#endif
 	problem.variables[slackIndex].value = valueForSlack;
 }
 
@@ -414,8 +456,10 @@ size_t LPSolver::addNewVariableToSubProblem(size_t _subProblem)
 {
 	SubProblem& problem = unseal(_subProblem);
 	size_t index = problem.variables.size();
+#ifndef SPARSE
 	for (LinearExpression& c: problem.factors)
 		c.resize(index + 1);
+#endif
 	problem.variables.emplace_back();
 	return index;
 }
@@ -519,12 +563,23 @@ string LPSolver::SubProblem::toString() const
 			resultString += "       ";
 		resultString += "   := " + v.value.toString() + "\n";
 	}
+#ifdef SPARSE
+	for (size_t rowIndex = 0; rowIndex < factors.rows(); rowIndex++)
+#else
 	for (auto&& [rowIndex, row]: factors | ranges::views::enumerate)
+#endif
 	{
 		string basicVarPrefix;
 		string rowString;
+#ifdef SPARSE
+		for (auto&& entry: const_cast<SparseMatrix&>(factors).iterateRow(rowIndex))
+		{
+			rational const& f = entry.value;
+			size_t i = entry.col;
+#else
 		for (auto&& [i, f]: row.enumerate())
 		{
+#endif
 			if (basicVariables.count(i) && basicVariables.at(i) == rowIndex)
 			{
 				solAssert(f == -1);
@@ -590,13 +645,24 @@ void LPSolver::SubProblem::update(size_t _varIndex, RationalWithDelta const& _va
 {
 	RationalWithDelta delta = _value - variables[_varIndex].value;
 	variables[_varIndex].value = _value;
+#ifdef SPARSE
+	// TODO can we store that?
+	map<size_t, size_t> basicVarForRow = invertMap(basicVariables);
+	for (auto&& entry: factors.iterateColumn(_varIndex))
+		if (entry.value && basicVarForRow.count(entry.row))
+		{
+			size_t j = basicVarForRow[entry.row];
+			variables[j].value += delta * entry.value;
+			//variablesPotentiallyOutOfBounds.insert(j);
+		}
+#else
 	for (auto&& [j, row]: basicVariables)
 		if (factors[row][_varIndex])
 		{
 			variables[j].value += delta * factors[row][_varIndex];
 			//variablesPotentiallyOutOfBounds.insert(j);
 		}
-
+#endif
 }
 
 optional<size_t> LPSolver::SubProblem::firstConflictingBasicVariable() const
@@ -619,9 +685,16 @@ optional<size_t> LPSolver::SubProblem::firstReplacementVar(
 	bool _increasing
 ) const
 {
+#ifdef SPARSE
+	for (auto&& entry: const_cast<SparseMatrix&>(factors).iterateRow(basicVariables.at(_basicVarToReplace)))
+	{
+		size_t i = entry.col;
+		rational const& factor = entry.value;
+#else
 	LinearExpression const& basicVarEquation = factors[basicVariables.at(_basicVarToReplace)];
 	for (auto const& [i, factor]: basicVarEquation.enumerate())
 	{
+#endif
 		if (i == _basicVarToReplace || !factor)
 			continue;
 		bool positive = factor > 0;
@@ -647,9 +720,16 @@ set<size_t> LPSolver::SubProblem::reasonsForUnsat(
 	else if (!_increasing && variables[_basicVarToReplace].bounds.upperReason)
 		r.insert(*variables[_basicVarToReplace].bounds.upperReason);
 
+#ifdef SPARSE
+	for (auto&& entry: const_cast<SparseMatrix&>(factors).iterateRow(basicVariables.at(_basicVarToReplace)))
+	{
+		size_t i = entry.col;
+		rational const& factor = entry.value;
+#else
 	LinearExpression const& basicVarEquation = factors[basicVariables.at(_basicVarToReplace)];
 	for (auto const& [i, factor]: basicVarEquation.enumerate())
 	{
+#endif
 		if (i == _basicVarToReplace || !factor)
 			continue;
 		bool positive = factor > 0;
@@ -669,6 +749,16 @@ void LPSolver::SubProblem::pivot(size_t _old, size_t _new)
 	// Transform pivotRow such that the coefficient for _new is -1
 	// Then use that to set all other coefficients for _new to zero.
 	size_t pivotRow = basicVariables[_old];
+#ifdef SPARSE
+	rational pivot = factors.entry(pivotRow, _new).value;
+	solAssert(pivot != 0, "");
+	if (pivot != -1)
+		factors.multiplyRowByFactor(pivotRow, rational{-1} / pivot);
+
+	for (auto&& entry: factors.iterateColumn(_new))
+		if (entry.row != pivotRow)
+			factors.addMultipleOfRow(pivotRow, entry.row, entry.value);
+#else
 	LinearExpression& pivotRowData = factors[pivotRow];
 
 	rational pivot = pivotRowData[_new];
@@ -691,6 +781,7 @@ void LPSolver::SubProblem::pivot(size_t _old, size_t _new)
 	for (size_t i = 0; i < factors.size(); ++i)
 		if (i != pivotRow)
 			subtractMultipleOfPivotRow(factors[i]);
+#endif
 
 	basicVariables.erase(_old);
 	basicVariables[_new] = pivotRow;
@@ -702,14 +793,30 @@ void LPSolver::SubProblem::pivotAndUpdate(
 	size_t _newBasicVar
 )
 {
+#ifdef SPARSE
+	RationalWithDelta theta = (_newValue - variables[_oldBasicVar].value) / factors.entry(_oldBasicVar, _newBasicVar).value;
+#else
 	RationalWithDelta theta = (_newValue - variables[_oldBasicVar].value) / factors[basicVariables[_oldBasicVar]][_newBasicVar];
+#endif
 
 	variables[_oldBasicVar].value = _newValue;
 	variables[_newBasicVar].value += theta;
 
+#ifdef SPARSE
+	// TODO can we store that?
+	map<size_t, size_t> basicVarForRow = invertMap(basicVariables);
+	for (auto&& entry: factors.iterateColumn(_newBasicVar))
+		if (basicVarForRow.count(entry.row))
+		{
+			size_t i = basicVarForRow[entry.row];
+			if (i != _oldBasicVar)
+				variables[i].value += theta * entry.value;
+		}
+#else
 	for (auto const& [i, row]: basicVariables)
 		if (i != _oldBasicVar && factors[row][_newBasicVar])
 			variables[i].value += theta * factors[row][_newBasicVar];
+#endif
 
 	pivot(_oldBasicVar, _newBasicVar);
 }

libsolutil/LinearExpression.cpp

@@ -23,10 +23,16 @@ using namespace std;
 
 SparseMatrix::SparseMatrixIterator SparseMatrix::IteratorCombiner::begin()
 {
-	return SparseMatrixIterator(
+	if (
-		m_isRow ? m_matrix.m_row_start[m_RowOrColumn] : m_matrix.m_col_start[m_RowOrColumn],
+		(m_isRow && m_rowOrColumn >= m_matrix.m_row_start.size()) ||
-		m_isRow
+		(!m_isRow && m_rowOrColumn >= m_matrix.m_col_start.size())
-	);
+	)
+		return SparseMatrixIterator(nullptr, m_isRow);
+	else
+		return SparseMatrixIterator(
+			m_isRow ? m_matrix.m_row_start[m_rowOrColumn] : m_matrix.m_col_start[m_rowOrColumn],
+			m_isRow
+		);
 }
 
 SparseMatrix::SparseMatrixIterator SparseMatrix::IteratorCombiner::end()
@@ -34,7 +40,7 @@ SparseMatrix::SparseMatrixIterator SparseMatrix::IteratorCombiner::end()
 	return SparseMatrixIterator(nullptr, m_isRow);
 }
 
-SparseMatrix::IteratorCombiner SparseMatrix::enumerateColumn(size_t _column)
+SparseMatrix::IteratorCombiner SparseMatrix::iterateColumn(size_t _column)
 {
 	return IteratorCombiner{
 		_column,
@@ -43,7 +49,7 @@ SparseMatrix::IteratorCombiner SparseMatrix::enumerateColumn(size_t _column)
 	};
 }
 
-SparseMatrix::IteratorCombiner SparseMatrix::enumerateRow(size_t _row)
+SparseMatrix::IteratorCombiner SparseMatrix::iterateRow(size_t _row)
 {
 	return IteratorCombiner{
 		_row,
@@ -67,6 +73,11 @@ void SparseMatrix::multiplyRowByFactor(size_t _row, rational const& _factor)
 
 void SparseMatrix::addMultipleOfRow(size_t _sourceRow, size_t _targetRow, rational const& _factor)
 {
+	if (_targetRow >= m_row_start.size())
+	{
+		m_row_start.resize(_targetRow + 1);
+		m_row_end.resize(_targetRow + 1);
+	}
 	Entry* source = m_row_start[_sourceRow];
 	Entry* target = m_row_start[_targetRow];
 
@@ -90,6 +101,30 @@ void SparseMatrix::addMultipleOfRow(size_t _sourceRow, size_t _targetRow, ration
 	}
 }
 
+SparseMatrix::Entry& SparseMatrix::entry(size_t _row, size_t _column)
+{
+	Entry* successor = entryOrSuccessorInRow(_row, _column);
+	if (successor && successor->col == _column)
+		return *successor;
+	else
+		return *prependInRow(successor, _row, _column, {});
+}
+
+void SparseMatrix::insert(size_t _row, size_t _column, rational _value)
+{
+	if (_column >= m_col_start.size())
+	{
+		m_col_start.resize(_column + 1);
+		m_col_end.resize(_column + 1);
+	}
+	if (_row >= m_row_start.size())
+	{
+		m_row_start.resize(_row + 1);
+		m_row_end.resize(_row + 1);
+	}
+
+	prependInRow(entryOrSuccessorInRow(_row, _column), _row, _column, move(_value));
+}
 
 void SparseMatrix::appendRow(LinearExpression const& _entries)
 {
@@ -103,6 +138,19 @@ void SparseMatrix::appendRow(LinearExpression const& _entries)
 	}
 }
 
+SparseMatrix::Entry* SparseMatrix::entryOrSuccessorInRow(size_t _row, size_t _column)
+{
+	Entry* successor = nullptr;
+	if (m_row_end[_row] && m_row_end[_row]->col >= _column)
+	{
+		successor = m_row_start[_row];
+		// TODO could choose to search from end
+		while (successor && successor->col < _column)
+			successor = successor->next_in_row;
+	}
+	return successor;
+}
+
 void SparseMatrix::remove(SparseMatrix::Entry& _e)
 {
 	if (_e.prev_in_row)

libsolutil/LinearExpression.h

@@ -213,7 +213,7 @@ public:
 		pointer operator->() { return m_ptr; }
 		SparseMatrixIterator& operator++()
 		{
-			m_ptr = m_isRow ? m_ptr->next_in_row : m_pt->next_in_col;
+			m_ptr = m_isRow ? m_ptr->next_in_row : m_ptr->next_in_col;
 			return *this;
 		}
 		SparseMatrixIterator operator++(int) { SparseMatrixIterator tmp = *this; ++(*this); return tmp; }
@@ -232,7 +232,7 @@ public:
 	};
 	struct IteratorCombiner
 	{
-		size_t m_RowOrColumn;
+		size_t m_rowOrColumn;
 		bool m_isRow;
 		SparseMatrix& m_matrix;
 		SparseMatrixIterator begin();
@@ -242,25 +242,30 @@ public:
 	size_t rows() const { return m_row_start.size(); }
 	size_t columns() const { return m_col_start.size(); }
 
-	/// @returns (i, v) for all non-zero v in the column _column
+	/// @returns Entry for all non-zero v in the column _column
-	void enumerateColumn(size_t _column);
+	IteratorCombiner iterateColumn(size_t _column);
-	/// @returns (i, v) for all non-zero v in the row _row
+	/// @returns Entry for all non-zero v in the row _row
-	void enumerateRow(size_t _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);
-	rational entry(size_t _row, size_t _column) const;
+	Entry& entry(size_t _row, size_t _column);
+	/// Inserts the value at the row/rolumn.
+	/// Assumes the entry does not exist yet.
 	void insert(size_t _row, size_t _column, rational _value);
 
 	void appendRow(LinearExpression const& _entries);
 
 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);
 
-	std::vector<std::unique_ptr<Entry>> m_elements;
+	// TODO unique_ptr?
+	std::vector<std::shared_ptr<Entry>> m_elements;
 	std::vector<Entry*> m_row_start;
 	std::vector<Entry*> m_col_start;
 	std::vector<Entry*> m_row_end;