Google OR-Tools v9.12
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
sparse.h
Go to the documentation of this file.
1// Copyright 2010-2025 Google LLC
2// Licensed under the Apache License, Version 2.0 (the "License");
3// you may not use this file except in compliance with the License.
4// You may obtain a copy of the License at
5//
6// http://www.apache.org/licenses/LICENSE-2.0
7//
8// Unless required by applicable law or agreed to in writing, software
9// distributed under the License is distributed on an "AS IS" BASIS,
10// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11// See the License for the specific language governing permissions and
12// limitations under the License.
13
14//
15// The following are very good references for terminology, data structures,
16// and algorithms:
17//
18// I.S. Duff, A.M. Erisman and J.K. Reid, "Direct Methods for Sparse Matrices",
19// Clarendon, Oxford, UK, 1987, ISBN 0-19-853421-3,
20// http://www.amazon.com/dp/0198534213.
21//
22//
23// T.A. Davis, "Direct methods for Sparse Linear Systems", SIAM, Philadelphia,
24// 2006, ISBN-13: 978-0-898716-13, http://www.amazon.com/dp/0898716136.
25//
26//
27// Both books also contain a wealth of references.
28
29#ifndef OR_TOOLS_LP_DATA_SPARSE_H_
30#define OR_TOOLS_LP_DATA_SPARSE_H_
31
32#include <algorithm>
33#include <cstdint>
34#include <string>
35#include <vector>
36
37#include "absl/log/check.h"
38#include "absl/types/span.h"
39#include "ortools/base/types.h"
45#include "ortools/util/bitset.h"
47
48namespace operations_research {
49namespace glop {
50
52
53// --------------------------------------------------------
54// SparseMatrix
55// --------------------------------------------------------
56// SparseMatrix is a class for sparse matrices suitable for computation.
57// Data is represented using the so-called compressed-column storage scheme.
58// Entries (row, col, value) are stored by column using a SparseColumn.
59//
60// Citing [Duff et al, 1987], a matrix is sparse if many of its coefficients are
61// zero and if there is an advantage in exploiting its zeros.
62// For practical reasons, not all zeros are exploited (for example those that
63// result from calculations.) The term entry refers to those coefficients that
64// are handled explicitly. All non-zeros are entries while some zero
65// coefficients may also be entries.
66//
67// Note that no special ordering of entries is assumed.
68class SparseMatrix {
69 public:
71
72 // Useful for testing. This makes it possible to write:
73 // SparseMatrix matrix {
74 // {1, 2, 3},
75 // {4, 5, 6},
76 // {7, 8, 9}};
77#if (!defined(_MSC_VER) || _MSC_VER >= 1800)
79 std::initializer_list<std::initializer_list<Fractional>> init_list);
80
81 // This type is neither copyable nor movable.
82 SparseMatrix(const SparseMatrix&) = delete;
83 SparseMatrix& operator=(const SparseMatrix&) = delete;
85#endif
86 // Clears internal data structure, i.e. erases all the columns and set
87 // the number of rows to zero.
88 void Clear();
89
90 // Returns true if the matrix is empty.
91 // That is if num_rows() OR num_cols() are zero.
92 bool IsEmpty() const;
93
94 // Cleans the columns, i.e. removes zero-values entries, removes duplicates
95 // entries and sorts remaining entries in increasing row order.
96 // Call with care: Runs in O(num_cols * column_cleanup), with each column
97 // cleanup running in O(num_entries * log(num_entries)).
98 void CleanUp();
99
100 // Call CheckNoDuplicates() on all columns, useful for doing a DCHECK.
101 bool CheckNoDuplicates() const;
102
103 // Call IsCleanedUp() on all columns, useful for doing a DCHECK.
104 bool IsCleanedUp() const;
105
106 // Change the number of row of this matrix.
107 void SetNumRows(RowIndex num_rows);
108
109 // Appends an empty column and returns its index.
110 ColIndex AppendEmptyColumn();
111
112 // Appends a unit vector defined by the single entry (row, value).
113 // Note that the row should be smaller than the number of rows of the matrix.
114 void AppendUnitVector(RowIndex row, Fractional value);
115
116 // Swaps the content of this SparseMatrix with the one passed as argument.
117 // Works in O(1).
118 void Swap(SparseMatrix* matrix);
119
120 // Populates the matrix with num_cols columns of zeros. As the number of rows
121 // is specified by num_rows, the matrix is not necessarily square.
122 // Previous columns/values are deleted.
123 void PopulateFromZero(RowIndex num_rows, ColIndex num_cols);
124
125 // Populates the matrix from the Identity matrix of size num_cols.
126 // Previous columns/values are deleted.
127 void PopulateFromIdentity(ColIndex num_cols);
128
129 // Populates the matrix from the transposed of the given matrix.
130 // Note that this preserve the property of lower/upper triangular matrix
131 // to have the diagonal coefficients first/last in each columns. It actually
132 // sorts the entries in each columns by their indices.
133 template <typename Matrix>
134 void PopulateFromTranspose(const Matrix& input);
135
136 // Populates a SparseMatrix from another one (copy), note that this run in
137 // O(number of entries in the matrix).
138 void PopulateFromSparseMatrix(const SparseMatrix& matrix);
139
140 // Populates a SparseMatrix from the image of a matrix A through the given
141 // row_perm and inverse_col_perm. See permutation.h for more details.
142 template <typename Matrix>
143 void PopulateFromPermutedMatrix(const Matrix& a,
144 const RowPermutation& row_perm,
145 const ColumnPermutation& inverse_col_perm);
146
147 // Populates a SparseMatrix from the result of alpha * A + beta * B,
148 // where alpha and beta are Fractionals, A and B are sparse matrices.
150 Fractional beta, const SparseMatrix& b);
151
152 // Multiplies SparseMatrix a by SparseMatrix b.
153 void PopulateFromProduct(const SparseMatrix& a, const SparseMatrix& b);
154
155 // Removes the marked columns from the matrix and adjust its size.
156 // This runs in O(num_cols).
157 void DeleteColumns(const DenseBooleanRow& columns_to_delete);
158
159 // Applies the given row permutation and deletes the rows for which
160 // permutation[row] is kInvalidRow. Sets the new number of rows to num_rows.
161 // This runs in O(num_entries).
162 void DeleteRows(RowIndex num_rows, const RowPermutation& permutation);
163
164 // Appends all rows from the given matrix to the calling object after the last
165 // row of the calling object. Both matrices must have the same number of
166 // columns. The method returns true if the rows were added successfully and
167 // false if it can't add the rows because the number of columns of the
168 // matrices are different.
169 bool AppendRowsFromSparseMatrix(const SparseMatrix& matrix);
170
171 // Applies the row permutation.
172 void ApplyRowPermutation(const RowPermutation& row_perm);
173
174 // Returns the coefficient at position row in column col.
175 // Call with care: runs in O(num_entries_in_col) as entries may not be sorted.
176 Fractional LookUpValue(RowIndex row, ColIndex col) const;
177
178 // Returns true if the matrix equals a (with a maximum error smaller than
179 // given the tolerance).
180 bool Equals(const SparseMatrix& a, Fractional tolerance) const;
181
182 // Returns, in min_magnitude and max_magnitude, the minimum and maximum
183 // magnitudes of the non-zero coefficients of the calling object.
184 void ComputeMinAndMaxMagnitudes(Fractional* min_magnitude,
185 Fractional* max_magnitude) const;
186
187 // Return the matrix dimension.
188 RowIndex num_rows() const { return num_rows_; }
189 ColIndex num_cols() const { return ColIndex(columns_.size()); }
191 // Access the underlying sparse columns.
192 const SparseColumn& column(ColIndex col) const { return columns_[col]; }
193 SparseColumn* mutable_column(ColIndex col) { return &(columns_[col]); }
195 // Returns the total numbers of entries in the matrix.
196 // Runs in O(num_cols).
197 EntryIndex num_entries() const;
198
199 // Computes the 1-norm of the matrix.
200 // The 1-norm |A| is defined as max_j sum_i |a_ij| or
201 // max_col sum_row |a(row,col)|.
203
204 // Computes the oo-norm (infinity-norm) of the matrix.
205 // The oo-norm |A| is defined as max_i sum_j |a_ij| or
206 // max_row sum_col |a(row,col)|.
208
209 // Returns a dense representation of the matrix.
210 std::string Dump() const;
211
212 private:
213 // Resets the internal data structure and create an empty rectangular
214 // matrix of size num_rows x num_cols.
215 void Reset(ColIndex num_cols, RowIndex num_rows);
216
217 // Vector of sparse columns.
219
220 // Number of rows. This is needed as sparse columns don't have a maximum
221 // number of rows.
222 RowIndex num_rows_;
223};
224
225// A matrix constructed from a list of already existing SparseColumn. This class
226// does not take ownership of the underlying columns, and thus they must outlive
227// this class (and keep the same address in memory).
228class MatrixView {
229 public:
230 MatrixView() = default;
231 explicit MatrixView(const SparseMatrix& matrix) {
234
235 // Takes all the columns of the given matrix.
236 void PopulateFromMatrix(const SparseMatrix& matrix) {
237 const ColIndex num_cols = matrix.num_cols();
238 columns_.resize(num_cols, nullptr);
239 for (ColIndex col(0); col < num_cols; ++col) {
240 columns_[col] = &matrix.column(col);
241 }
242 num_rows_ = matrix.num_rows();
243 }
244
245 // Takes all the columns of the first matrix followed by the columns of the
246 // second matrix.
247 void PopulateFromMatrixPair(const SparseMatrix& matrix_a,
248 const SparseMatrix& matrix_b) {
249 const ColIndex num_cols = matrix_a.num_cols() + matrix_b.num_cols();
250 columns_.resize(num_cols, nullptr);
251 for (ColIndex col(0); col < matrix_a.num_cols(); ++col) {
252 columns_[col] = &matrix_a.column(col);
253 }
254 for (ColIndex col(0); col < matrix_b.num_cols(); ++col) {
255 columns_[matrix_a.num_cols() + col] = &matrix_b.column(col);
256 }
257 num_rows_ = std::max(matrix_a.num_rows(), matrix_b.num_rows());
258 }
259
260 // Takes only the columns of the given matrix that belongs to the given basis.
261 void PopulateFromBasis(const MatrixView& matrix,
262 const RowToColMapping& basis) {
263 columns_.resize(RowToColIndex(basis.size()), nullptr);
264 for (RowIndex row(0); row < basis.size(); ++row) {
265 columns_[RowToColIndex(row)] = &matrix.column(basis[row]);
266 }
267 num_rows_ = matrix.num_rows();
268 }
269
270 // Same behavior as the SparseMatrix functions above.
271 bool IsEmpty() const { return columns_.empty(); }
272 RowIndex num_rows() const { return num_rows_; }
273 ColIndex num_cols() const { return columns_.size(); }
274 const SparseColumn& column(ColIndex col) const { return *columns_[col]; }
275 EntryIndex num_entries() const;
278
279 private:
280 RowIndex num_rows_;
282};
283
285 const SparseMatrix& input);
287 const SparseMatrix& a, const RowPermutation& row_perm,
288 const ColumnPermutation& inverse_col_perm);
289extern template void
291 const CompactSparseMatrixView& a, const RowPermutation& row_perm,
292 const ColumnPermutation& inverse_col_perm);
293
294// Another matrix representation which is more efficient than a SparseMatrix but
295// doesn't allow matrix modification. It is faster to construct, uses less
296// memory and provides a better cache locality when iterating over the non-zeros
297// of the matrix columns.
299 public:
300 // When iteration performance matter, getting a ConstView allows the compiler
301 // to do better aliasing analysis and not re-read vectors address all the
302 // time.
303 class ConstView {
304 public:
305 explicit ConstView(const CompactSparseMatrix* matrix)
306 : coefficients_(matrix->coefficients_.data()),
307 rows_(matrix->rows_.data()),
308 starts_(matrix->starts_.data()) {}
309
310 // Functions to iterate on the entries of a given column:
311 // const auto view = compact_matrix.view();
312 // for (const EntryIndex i : view.Column(col)) {
313 // const RowIndex row = view.EntryRow(i);
314 // const Fractional coefficient = view.EntryCoefficient(i);
315 // }
316 ::util::IntegerRange<EntryIndex> Column(ColIndex col) const {
317 return ::util::IntegerRange<EntryIndex>(starts_[col.value()],
318 starts_[col.value() + 1]);
319 }
320 Fractional EntryCoefficient(EntryIndex i) const {
321 return coefficients_[i.value()];
323 RowIndex EntryRow(EntryIndex i) const { return rows_[i.value()]; }
324
325 EntryIndex ColumnNumEntries(ColIndex col) const {
326 return starts_[col.value() + 1] - starts_[col.value()];
328
329 // Returns the scalar product of the given row vector with the column of
330 // index col of this matrix.
331 Fractional ColumnScalarProduct(ColIndex col,
332 DenseRow::ConstView vector) const;
333
334 private:
335 const Fractional* const coefficients_;
336 const RowIndex* const rows_;
337 const EntryIndex* const starts_;
338 };
339
340 CompactSparseMatrix() = default;
341 ConstView view() const { return ConstView(this); }
343 // Convenient constructors for tests.
344 // TODO(user): If this is needed in production code, it can be done faster.
345 explicit CompactSparseMatrix(const SparseMatrix& matrix) {
348
349 // This type is neither copyable nor movable.
351 CompactSparseMatrix& operator=(const CompactSparseMatrix&) = delete;
353 // Creates a CompactSparseMatrix from the given MatrixView. The matrices are
354 // the same, only the representation differ. Note that the entry order in
355 // each column is preserved.
357
358 // Creates a CompactSparseMatrix by copying the input and adding an identity
359 // matrix to the left of it.
361
362 // Creates a CompactSparseMatrix from the transpose of the given
363 // CompactSparseMatrix. Note that the entries in each columns will be ordered
364 // by row indices.
366
367 // Clears the matrix and sets its number of rows. If none of the Populate()
368 // function has been called, Reset() must be called before calling any of the
369 // Add*() functions below.
370 void Reset(RowIndex num_rows);
371
372 // Api to add columns one at the time.
373 void AddEntryToCurrentColumn(RowIndex row, Fractional coeff);
374 void CloseCurrentColumn();
375
376 // Adds a dense column to the CompactSparseMatrix (only the non-zero will be
377 // actually stored). This work in O(input.size()) and returns the index of the
378 // added column.
379 ColIndex AddDenseColumn(const DenseColumn& dense_column);
380
381 // Same as AddDenseColumn(), but only adds the non-zero from the given start.
382 ColIndex AddDenseColumnPrefix(DenseColumn::ConstView dense_column,
383 RowIndex start);
384
385 // Same as AddDenseColumn(), but uses the given non_zeros pattern of input.
386 // If non_zeros is empty, this actually calls AddDenseColumn().
387 ColIndex AddDenseColumnWithNonZeros(const DenseColumn& dense_column,
388 absl::Span<const RowIndex> non_zeros);
389
390 // Adds a dense column for which we know the non-zero positions and clears it.
391 // Note that this function supports duplicate indices in non_zeros. The
392 // complexity is in O(non_zeros.size()). Only the indices present in non_zeros
393 // will be cleared. Returns the index of the added column.
395 std::vector<RowIndex>* non_zeros);
396
397 // Returns the number of entries (i.e. degree) of the given column.
398 EntryIndex ColumnNumEntries(ColIndex col) const {
399 return starts_[col + 1] - starts_[col];
401
402 // Returns the matrix dimensions. See same functions in SparseMatrix.
403 EntryIndex num_entries() const {
404 DCHECK_EQ(coefficients_.size(), rows_.size());
405 return coefficients_.size();
406 }
407 RowIndex num_rows() const { return num_rows_; }
408 ColIndex num_cols() const { return num_cols_; }
410 // Returns whether or not this matrix contains any non-zero entries.
411 bool IsEmpty() const {
412 DCHECK_EQ(coefficients_.size(), rows_.size());
413 return coefficients_.empty();
414 }
415
416 // Alternative iteration API compatible with the one from SparseMatrix.
417 // The ConstView alternative should be faster.
418 ColumnView column(ColIndex col) const {
419 DCHECK_LT(col, num_cols_);
421 // Note that the start may be equal to row.size() if the last columns
422 // are empty, it is why we don't use &row[start].
423 const EntryIndex start = starts_[col];
424 return ColumnView(starts_[col + 1] - start, rows_.data() + start.value(),
425 coefficients_.data() + start.value());
426 }
427
428 // Returns true if the given column is empty. Note that for triangular matrix
429 // this does not include the diagonal coefficient (see below).
430 bool ColumnIsEmpty(ColIndex col) const {
431 return starts_[col + 1] == starts_[col];
433
434 // Returns the scalar product of the given row vector with the column of index
435 // col of this matrix.
436 Fractional ColumnScalarProduct(ColIndex col, const DenseRow& vector) const {
437 return view().ColumnScalarProduct(col, vector.const_view());
439
440 // Adds a multiple of the given column of this matrix to the given
441 // dense_column. If multiplier is 0.0, this function does nothing. This
442 // function is declared in the .h for efficiency.
443 void ColumnAddMultipleToDenseColumn(ColIndex col, Fractional multiplier,
444 DenseColumn::View dense_column) const {
445 if (multiplier == 0.0) return;
446 const auto entry_rows = rows_.view();
447 const auto entry_coeffs = coefficients_.view();
448 for (const EntryIndex i : Column(col)) {
449 dense_column[entry_rows[i]] += multiplier * entry_coeffs[i];
450 }
451 }
452 void ColumnAddMultipleToDenseColumn(ColIndex col, Fractional multiplier,
453 DenseColumn* dense_column) const {
454 return ColumnAddMultipleToDenseColumn(col, multiplier,
455 dense_column->view());
456 }
457
458 // Same as ColumnAddMultipleToDenseColumn() but also adds the new non-zeros to
459 // the non_zeros vector. A non-zero is "new" if is_non_zero[row] was false,
460 // and we update dense_column[row]. This function also updates is_non_zero.
461 void ColumnAddMultipleToSparseScatteredColumn(ColIndex col,
462 Fractional multiplier,
465 if (multiplier == 0.0) return;
466 const auto entry_rows = rows_.view();
467 const auto entry_coeffs = coefficients_.view();
468 for (const EntryIndex i : Column(col)) {
469 column->Add(entry_rows[i], multiplier * entry_coeffs[i]);
470 }
471 }
472
473 // Copies the given column of this matrix into the given dense_column.
474 // This function is declared in the .h for efficiency.
475 void ColumnCopyToDenseColumn(ColIndex col, DenseColumn* dense_column) const {
476 RETURN_IF_NULL(dense_column);
477 dense_column->AssignToZero(num_rows_);
478 ColumnCopyToClearedDenseColumn(col, dense_column);
479 }
480
481 // Same as ColumnCopyToDenseColumn() but assumes the column to be initially
482 // all zero.
483 void ColumnCopyToClearedDenseColumn(ColIndex col,
484 DenseColumn* dense_column) const {
485 RETURN_IF_NULL(dense_column);
486 dense_column->resize(num_rows_, 0.0);
487 const auto entry_rows = rows_.view();
488 const auto entry_coeffs = coefficients_.view();
489 for (const EntryIndex i : Column(col)) {
490 (*dense_column)[entry_rows[i]] = entry_coeffs[i];
491 }
492 }
493
494 // Same as ColumnCopyToClearedDenseColumn() but also fills non_zeros.
495 void ColumnCopyToClearedDenseColumnWithNonZeros(
496 ColIndex col, DenseColumn* dense_column,
497 RowIndexVector* non_zeros) const {
498 RETURN_IF_NULL(dense_column);
499 dense_column->resize(num_rows_, 0.0);
500 non_zeros->clear();
501 const auto entry_rows = rows_.view();
502 const auto entry_coeffs = coefficients_.view();
503 for (const EntryIndex i : Column(col)) {
504 const RowIndex row = entry_rows[i];
505 (*dense_column)[row] = entry_coeffs[i];
506 non_zeros->push_back(row);
507 }
508 }
509
510 void Swap(CompactSparseMatrix* other);
511
512 protected:
513 // Functions to iterate on the entries of a given column.
514 ::util::IntegerRange<EntryIndex> Column(ColIndex col) const {
515 return ::util::IntegerRange<EntryIndex>(starts_[col], starts_[col + 1]);
517
518 // The matrix dimensions, properly updated by full and incremental builders.
519 RowIndex num_rows_;
520 ColIndex num_cols_;
522 // Holds the columns non-zero coefficients and row positions.
523 // The entries for the column of index col are stored in the entries
524 // [starts_[col], starts_[col + 1]).
531 ColIndex col, DenseRow::ConstView vector) const {
532 // We expand ourselves since we don't really care about the floating
533 // point order of operation and this seems faster.
534 int i = starts_[col.value()].value();
535 const int end = starts_[col.value() + 1].value();
536 const int shifted_end = end - 3;
537 Fractional result1 = 0.0;
538 Fractional result2 = 0.0;
539 Fractional result3 = 0.0;
540 Fractional result4 = 0.0;
541 for (; i < shifted_end; i += 4) {
542 result1 += coefficients_[i] * vector[RowToColIndex(rows_[i])];
543 result2 += coefficients_[i + 1] * vector[RowToColIndex(rows_[i + 1])];
544 result3 += coefficients_[i + 2] * vector[RowToColIndex(rows_[i + 2])];
545 result4 += coefficients_[i + 3] * vector[RowToColIndex(rows_[i + 3])];
546 }
547 Fractional result = result1 + result2 + result3 + result4;
548 if (i < end) {
549 result += coefficients_[i] * vector[RowToColIndex(rows_[i])];
550 if (i + 1 < end) {
551 result += coefficients_[i + 1] * vector[RowToColIndex(rows_[i + 1])];
552 if (i + 2 < end) {
553 result += coefficients_[i + 2] * vector[RowToColIndex(rows_[i + 2])];
554 }
555 }
556 }
557 return result;
558}
559
560// A matrix view of the basis columns of a CompactSparseMatrix, with basis
561// specified as a RowToColMapping. This class does not take ownership of the
562// underlying matrix or basis, and thus they must outlive this class (and keep
563// the same address in memory).
564class CompactSparseMatrixView {
565 public:
567 const RowToColMapping* basis)
568 : compact_matrix_(*compact_matrix),
569 columns_(basis->data(), basis->size().value()) {}
570 CompactSparseMatrixView(const CompactSparseMatrix* compact_matrix,
571 const std::vector<ColIndex>* columns)
572 : compact_matrix_(*compact_matrix), columns_(*columns) {}
573
574 // Same behavior as the SparseMatrix functions above.
575 bool IsEmpty() const { return compact_matrix_.IsEmpty(); }
576 RowIndex num_rows() const { return compact_matrix_.num_rows(); }
577 ColIndex num_cols() const { return ColIndex(columns_.size()); }
578 ColumnView column(ColIndex col) const {
579 return compact_matrix_.column(columns_[col.value()]);
581 EntryIndex num_entries() const;
582 Fractional ComputeOneNorm() const;
583 Fractional ComputeInfinityNorm() const;
584
585 private:
586 // We require that the underlying CompactSparseMatrix and RowToColMapping
587 // continue to own the (potentially large) data accessed via this view.
588 const CompactSparseMatrix& compact_matrix_;
589 const absl::Span<const ColIndex> columns_;
590};
591
592// Specialization of a CompactSparseMatrix used for triangular matrices.
593// To be able to solve triangular systems as efficiently as possible, the
594// diagonal entries are stored in a separate vector and not in the underlying
595// CompactSparseMatrix.
596//
597// Advanced usage: this class also support matrices that can be permuted into a
598// triangular matrix and some functions work directly on such matrices.
599class TriangularMatrix : private CompactSparseMatrix {
600 public:
601 TriangularMatrix() : all_diagonal_coefficients_are_one_(true) {}
602
603 // This type is neither copyable nor movable.
604 TriangularMatrix(const TriangularMatrix&) = delete;
607 // Only a subset of the functions from CompactSparseMatrix are exposed (note
608 // the private inheritance). They are extended to deal with diagonal
609 // coefficients properly.
611 void Swap(TriangularMatrix* other);
612 bool IsEmpty() const { return diagonal_coefficients_.empty(); }
613 RowIndex num_rows() const { return num_rows_; }
614 ColIndex num_cols() const { return num_cols_; }
615 EntryIndex num_entries() const {
616 return EntryIndex(num_cols_.value()) + coefficients_.size();
618
619 // On top of the CompactSparseMatrix functionality, TriangularMatrix::Reset()
620 // also pre-allocates space of size col_size for a number of internal vectors.
621 // This helps reduce costly push_back operations for large problems.
622 //
623 // WARNING: Reset() must be called with a sufficiently large col_capacity
624 // prior to any Add* calls (e.g., AddTriangularColumn).
625 void Reset(RowIndex num_rows, ColIndex col_capacity);
626
627 // Constructs a triangular matrix from the given SparseMatrix. The input is
628 // assumed to be lower or upper triangular without any permutations. This is
629 // checked in debug mode.
630 void PopulateFromTriangularSparseMatrix(const SparseMatrix& input);
631
632 // Functions to create a triangular matrix incrementally, column by column.
633 // A client needs to call Reset(num_rows) first, and then each column must be
634 // added by calling one of the 3 functions below.
635 //
636 // Note that the row indices of the columns are allowed to be permuted: the
637 // diagonal entry of the column #col not being necessarily on the row #col.
638 // This is why these functions require the 'diagonal_row' parameter. The
639 // permutation can be fixed at the end by a call to
640 // ApplyRowPermutationToNonDiagonalEntries() or accounted directly in the case
641 // of PermutedLowerSparseSolve().
642 void AddTriangularColumn(const ColumnView& column, RowIndex diagonal_row);
643 void AddTriangularColumnWithGivenDiagonalEntry(const SparseColumn& column,
644 RowIndex diagonal_row,
645 Fractional diagonal_value);
646 void AddDiagonalOnlyColumn(Fractional diagonal_value);
647
648 // Adds the given sparse column divided by diagonal_coefficient.
649 // The diagonal_row is assumed to be present and its value should be the
650 // same as the one given in diagonal_coefficient. Note that this function
651 // tests for zero coefficients in the input column and removes them.
652 void AddAndNormalizeTriangularColumn(const SparseColumn& column,
653 RowIndex diagonal_row,
654 Fractional diagonal_coefficient);
655
656 // Applies the given row permutation to all entries except the diagonal ones.
657 void ApplyRowPermutationToNonDiagonalEntries(const RowPermutation& row_perm);
658
659 // Copy a triangular column with its diagonal entry to the given SparseColumn.
660 void CopyColumnToSparseColumn(ColIndex col, SparseColumn* output) const;
661
662 // Copy a triangular matrix to the given SparseMatrix.
663 void CopyToSparseMatrix(SparseMatrix* output) const;
664
665 // Returns the index of the first column which is not an identity column (i.e.
666 // a column j with only one entry of value 1 at the j-th row). This is always
667 // zero if the matrix is not triangular.
668 ColIndex GetFirstNonIdentityColumn() const {
669 return first_non_identity_column_;
671
672 // Returns the diagonal coefficient of the given column.
673 Fractional GetDiagonalCoefficient(ColIndex col) const {
674 return diagonal_coefficients_[col];
676
677 // Returns true iff the column contains no non-diagonal entries.
678 bool ColumnIsDiagonalOnly(ColIndex col) const {
681
682 // --------------------------------------------------------------------------
683 // Triangular solve functions.
684 //
685 // All the functions containing the word Lower (resp. Upper) require the
686 // matrix to be lower (resp. upper_) triangular without any permutation.
687 // --------------------------------------------------------------------------
688
689 // Solve the system L.x = rhs for a lower triangular matrix.
690 // The result overwrite rhs.
691 void LowerSolve(DenseColumn* rhs) const;
692
693 // Solves the system U.x = rhs for an upper triangular matrix.
694 void UpperSolve(DenseColumn* rhs) const;
695
696 // Solves the system Transpose(U).x = rhs where U is upper triangular.
697 // This can be used to do a left-solve for a row vector (i.e. y.Y = rhs).
698 void TransposeUpperSolve(DenseColumn* rhs) const;
699
700 // This assumes that the rhs is all zero before the given position.
701 void LowerSolveStartingAt(ColIndex start, DenseColumn* rhs) const;
702
703 // Solves the system Transpose(L).x = rhs, where L is lower triangular.
704 // This can be used to do a left-solve for a row vector (i.e., y.Y = rhs).
705 void TransposeLowerSolve(DenseColumn* rhs) const;
706
707 // Hyper-sparse version of the triangular solve functions. The passed
708 // non_zero_rows should contain the positions of the symbolic non-zeros of the
709 // result in the order in which they need to be accessed (or in the reverse
710 // order for the Reverse*() versions).
711 //
712 // The non-zero vector is mutable so that the symbolic non-zeros that are
713 // actually zero because of numerical cancellations can be removed.
714 //
715 // The non-zeros can be computed by one of these two methods:
716 // - ComputeRowsToConsiderWithDfs() which will give them in the reverse order
717 // of the one they need to be accessed in. This is only a topological order,
718 // and it will not necessarily be "sorted".
719 // - ComputeRowsToConsiderInSortedOrder() which will always give them in
720 // increasing order.
721 //
722 // Note that if the non-zeros are given in a sorted order, then the
723 // hyper-sparse functions will return EXACTLY the same results as the non
724 // hyper-sparse version above.
725 //
726 // For a given solve, here is the required order:
727 // - For a lower solve, increasing non-zeros order.
728 // - For an upper solve, decreasing non-zeros order.
729 // - for a transpose lower solve, decreasing non-zeros order.
730 // - for a transpose upper solve, increasing non_zeros order.
731 //
732 // For a general discussion of hyper-sparsity in LP, see:
733 // J.A.J. Hall, K.I.M. McKinnon, "Exploiting hyper-sparsity in the revised
734 // simplex method", December 1999, MS 99-014.
735 // http://www.maths.ed.ac.uk/hall/MS-99/MS9914.pdf
736 void HyperSparseSolve(DenseColumn* rhs, RowIndexVector* non_zero_rows) const;
737 void HyperSparseSolveWithReversedNonZeros(
738 DenseColumn* rhs, RowIndexVector* non_zero_rows) const;
739 void TransposeHyperSparseSolve(DenseColumn* rhs,
740 RowIndexVector* non_zero_rows) const;
741 void TransposeHyperSparseSolveWithReversedNonZeros(
742 DenseColumn* rhs, RowIndexVector* non_zero_rows) const;
743
744 // Given the positions of the non-zeros of a vector, computes the non-zero
745 // positions of the vector after a solve by this triangular matrix. The order
746 // of the returned non-zero positions will be in the REVERSE elimination
747 // order. If the function detects that there are too many non-zeros, then it
748 // aborts early and non_zero_rows is cleared.
749 void ComputeRowsToConsiderWithDfs(RowIndexVector* non_zero_rows) const;
750
751 // Same as TriangularComputeRowsToConsider() but always returns the non-zeros
752 // sorted by rows. It is up to the client to call the direct or reverse
753 // hyper-sparse solve function depending if the matrix is upper or lower
754 // triangular.
755 void ComputeRowsToConsiderInSortedOrder(RowIndexVector* non_zero_rows) const;
756
757 // This is currently only used for testing. It achieves the same result as
758 // PermutedLowerSparseSolve() below, but the latter exploits the sparsity of
759 // rhs and is thus faster for our use case.
760 //
761 // Note that partial_inverse_row_perm only permutes the first k rows, where k
762 // is the same as partial_inverse_row_perm.size(). It is the inverse
763 // permutation of row_perm which only permutes k rows into is [0, k), the
764 // other row images beeing kInvalidRow. The other arguments are the same as
765 // for PermutedLowerSparseSolve() and described there.
766 //
767 // IMPORTANT: lower will contain all the "symbolic" non-zero entries.
768 // A "symbolic" zero entry is one that will be zero whatever the coefficients
769 // of the rhs entries. That is it only depends on the position of its
770 // entries, not on their values. Thus, some of its coefficients may be zero.
771 // This fact is exploited by the LU factorization code. The zero coefficients
772 // of upper will be cleaned, however.
773 void PermutedLowerSolve(const SparseColumn& rhs,
774 const RowPermutation& row_perm,
775 const RowMapping& partial_inverse_row_perm,
776 SparseColumn* lower, SparseColumn* upper) const;
777
778 // This solves a lower triangular system with only ones on the diagonal where
779 // the matrix and the input rhs are permuted by the inverse of row_perm. Note
780 // that the output will also be permuted by the inverse of row_perm. The
781 // function also supports partial permutation. That is if row_perm[i] < 0 then
782 // column row_perm[i] is assumed to be an identity column.
783 //
784 // The output is given as follow:
785 // - lower is cleared, and receives the rows for which row_perm[row] < 0
786 // meaning not yet examined as a pivot (see markowitz.cc).
787 // - upper is NOT cleared, and the other rows (row_perm[row] >= 0) are
788 // appended to it.
789 // - Note that lower and upper can point to the same SparseColumn.
790 //
791 // Note: This function is non-const because ComputeRowsToConsider() also
792 // prunes the underlying dependency graph of the lower matrix while doing a
793 // solve. See marked_ and pruned_ends_ below.
794 void PermutedLowerSparseSolve(const ColumnView& rhs,
795 const RowPermutation& row_perm,
796 SparseColumn* lower, SparseColumn* upper);
797
798 // This is used to compute the deterministic time of a matrix factorization.
799 int64_t NumFpOperationsInLastPermutedLowerSparseSolve() const {
800 return num_fp_operations_;
801 }
802
803 // To be used in DEBUG mode by the client code. This check that the matrix is
804 // lower- (resp. upper-) triangular without any permutation and that there is
805 // no zero on the diagonal. We can't do that on each Solve() that require so,
806 // otherwise it will be too slow in debug.
807 bool IsLowerTriangular() const;
808 bool IsUpperTriangular() const;
809
810 // Visible for testing. This is used by PermutedLowerSparseSolve() to compute
811 // the non-zero indices of the result. The output is as follow:
812 // - lower_column_rows will contains the rows for which row_perm[row] < 0.
813 // - upper_column_rows will contains the other rows in the reverse topological
814 // order in which they should be considered in PermutedLowerSparseSolve().
815 //
816 // This function is non-const because it prunes the underlying dependency
817 // graph of the lower matrix while doing a solve. See marked_ and pruned_ends_
818 // below.
819 //
820 // Pruning the graph at the same time is slower but not by too much (< 2x) and
821 // seems worth doing. Note that when the lower matrix is dense, most of the
822 // graph will likely be pruned. As a result, the symbolic phase will be
823 // negligible compared to the numerical phase so we don't really need a dense
824 // version of PermutedLowerSparseSolve().
826 const RowPermutation& row_perm,
827 RowIndexVector* lower_column_rows,
828 RowIndexVector* upper_column_rows);
829
830 // The upper bound is computed using one of the algorithm presented in
831 // "A Survey of Condition Number Estimation for Triangular Matrices"
832 // https:epubs.siam.org/doi/pdf/10.1137/1029112/
835
836 private:
837 // Internal versions of some Solve() functions to avoid code duplication.
838 template <bool diagonal_of_ones>
839 void LowerSolveStartingAtInternal(ColIndex start,
840 DenseColumn::View rhs) const;
841 template <bool diagonal_of_ones>
842 void UpperSolveInternal(DenseColumn::View rhs) const;
843 template <bool diagonal_of_ones>
844 void TransposeLowerSolveInternal(DenseColumn::View rhs) const;
845 template <bool diagonal_of_ones>
846 void TransposeUpperSolveInternal(DenseColumn::View rhs) const;
847 template <bool diagonal_of_ones>
848 void HyperSparseSolveInternal(DenseColumn::View rhs,
849 RowIndexVector* non_zero_rows) const;
850 template <bool diagonal_of_ones>
851 void HyperSparseSolveWithReversedNonZerosInternal(
852 DenseColumn::View rhs, RowIndexVector* non_zero_rows) const;
853 template <bool diagonal_of_ones>
854 void TransposeHyperSparseSolveInternal(DenseColumn::View rhs,
855 RowIndexVector* non_zero_rows) const;
856 template <bool diagonal_of_ones>
857 void TransposeHyperSparseSolveWithReversedNonZerosInternal(
858 DenseColumn::View rhs, RowIndexVector* non_zero_rows) const;
859
860 // Internal function used by the Add*() functions to finish adding
861 // a new column to a triangular matrix.
862 void CloseCurrentColumn(Fractional diagonal_value);
863
864 // Extra data for "triangular" matrices. The diagonal coefficients are
865 // stored in a separate vector instead of beeing stored in each column.
866 StrictITIVector<ColIndex, Fractional> diagonal_coefficients_;
867
868 // Index of the first column which is not a diagonal only column with a
869 // coefficient of 1. This is used to optimize the solves.
870 ColIndex first_non_identity_column_;
871
872 // This common case allows for more efficient Solve() functions.
873 // TODO(user): Do not even construct diagonal_coefficients_ in this case?
874 bool all_diagonal_coefficients_are_one_;
875
876 // For the hyper-sparse version. These are used to implement a DFS, see
877 // TriangularComputeRowsToConsider() for more details.
878 mutable Bitset64<RowIndex> stored_;
879 mutable std::vector<RowIndex> nodes_to_explore_;
880
881 // For PermutedLowerSparseSolve().
882 int64_t num_fp_operations_;
883 mutable std::vector<RowIndex> lower_column_rows_;
884 mutable std::vector<RowIndex> upper_column_rows_;
885 mutable DenseColumn initially_all_zero_scratchpad_;
886
887 // This boolean vector is used to detect entries that can be pruned during
888 // the DFS used for the symbolic phase of ComputeRowsToConsider().
889 //
890 // Problem: We have a DAG where each node has outgoing arcs towards other
891 // nodes (this adjacency list is NOT sorted by any order). We want to compute
892 // the reachability of a set of nodes S and its topological order. While doing
893 // this, we also want to prune the adjacency lists to exploit the simple fact
894 // that if a -> (b, c) and b -> (c) then c can be removed from the adjacency
895 // list of a since it will be implied through b. Note that this doesn't change
896 // the reachability of any set nor a valid topological ordering of such a set.
897 //
898 // The concept is known as the transitive reduction of a DAG, see
899 // http://en.wikipedia.org/wiki/Transitive_reduction.
900 //
901 // Heuristic algorithm: While doing the DFS to compute Reach(S) and its
902 // topological order, each time we process a node, we mark all its adjacent
903 // node while going down in the DFS, and then we unmark all of them when we go
904 // back up. During the un-marking, if a node is already un-marked, it means
905 // that it was implied by some other path starting at the current node and we
906 // can prune it and remove it from the adjacency list of the current node.
907 //
908 // Note(user): I couldn't find any reference for this algorithm, even though
909 // I suspect I am not the first one to need something similar.
910 mutable DenseBooleanColumn marked_;
911
912 // This is used to represent a pruned sub-matrix of the current matrix that
913 // corresponds to the pruned DAG as described in the comment above for
914 // marked_. This vector is used to encode the sub-matrix as follow:
915 // - Both the rows and the coefficients of the pruned matrix are still stored
916 // in rows_ and coefficients_.
917 // - The data of column 'col' is still stored starting at starts_[col].
918 // - But, its end is given by pruned_ends_[col] instead of starts_[col + 1].
919 //
920 // The idea of using a smaller graph for the symbolic phase is well known in
921 // sparse linear algebra. See:
922 // - John R. Gilbert and Joseph W. H. Liu, "Elimination structures for
923 // unsymmetric sparse LU factors", Tech. Report CS-90-11. Departement of
924 // Computer Science, York University, North York. Ontario, Canada, 1990.
925 // - Stanley C. Eisenstat and Joseph W. H. Liu, "Exploiting structural
926 // symmetry in a sparse partial pivoting code". SIAM J. Sci. Comput. Vol
927 // 14, No 1, pp. 253-257, January 1993.
928 //
929 // Note that we use an original algorithm and prune the graph while performing
930 // the symbolic phase. Hence the pruning will only benefit the next symbolic
931 // phase. This is different from Eisenstat-Liu's symmetric pruning. It is
932 // still a heuristic and will not necessarily find the minimal graph that
933 // has the same result for the symbolic phase though.
934 //
935 // TODO(user): Use this during the "normal" hyper-sparse solves so that
936 // we can benefit from the pruned lower matrix there?
938};
939
940} // namespace glop
941} // namespace operations_research
942
943#endif // OR_TOOLS_LP_DATA_SPARSE_H_
CompactSparseMatrixView(const CompactSparseMatrix *compact_matrix, const RowToColMapping *basis)
Definition sparse.h:568
EntryIndex ColumnNumEntries(ColIndex col) const
Definition sparse.h:327
ConstView(const CompactSparseMatrix *matrix)
Definition sparse.h:307
Fractional ColumnScalarProduct(ColIndex col, DenseRow::ConstView vector) const
Definition sparse.h:532
void PopulateFromMatrixView(const MatrixView &input)
Definition sparse.cc:447
CompactSparseMatrix & operator=(const CompactSparseMatrix &)=delete
ColIndex AddDenseColumnPrefix(DenseColumn::ConstView dense_column, RowIndex start)
Same as AddDenseColumn(), but only adds the non-zero from the given start.
Definition sparse.cc:597
StrictITIVector< EntryIndex, Fractional > coefficients_
Definition sparse.h:527
void ColumnAddMultipleToDenseColumn(ColIndex col, Fractional multiplier, DenseColumn *dense_column) const
Definition sparse.h:454
ColIndex AddDenseColumn(const DenseColumn &dense_column)
Definition sparse.cc:593
void PopulateFromSparseMatrixAndAddSlacks(const SparseMatrix &input)
Definition sparse.cc:466
StrictITIVector< ColIndex, EntryIndex > starts_
Definition sparse.h:529
EntryIndex num_entries() const
Returns the matrix dimensions. See same functions in SparseMatrix.
Definition sparse.h:405
ColIndex AddAndClearColumnWithNonZeros(DenseColumn *column, std::vector< RowIndex > *non_zeros)
Definition sparse.cc:626
bool IsEmpty() const
Returns whether or not this matrix contains any non-zero entries.
Definition sparse.h:413
ColIndex AddDenseColumnWithNonZeros(const DenseColumn &dense_column, absl::Span< const RowIndex > non_zeros)
Definition sparse.cc:611
void AddEntryToCurrentColumn(RowIndex row, Fractional coeff)
Api to add columns one at the time.
Definition sparse.cc:582
bool ColumnIsEmpty(ColIndex col) const
Definition sparse.h:432
EntryIndex ColumnNumEntries(ColIndex col) const
Returns the number of entries (i.e. degree) of the given column.
Definition sparse.h:400
StrictITIVector< EntryIndex, RowIndex > rows_
Definition sparse.h:528
void Swap(CompactSparseMatrix *other)
Definition sparse.cc:642
void ColumnCopyToClearedDenseColumn(ColIndex col, DenseColumn *dense_column) const
Definition sparse.h:485
Fractional ColumnScalarProduct(ColIndex col, const DenseRow &vector) const
Definition sparse.h:438
void PopulateFromTranspose(const CompactSparseMatrix &input)
Definition sparse.cc:495
RowIndex num_rows_
The matrix dimensions, properly updated by full and incremental builders.
Definition sparse.h:521
::util::IntegerRange< EntryIndex > Column(ColIndex col) const
Functions to iterate on the entries of a given column.
Definition sparse.h:516
ColumnView column(ColIndex col) const
Definition sparse.h:420
void PopulateFromMatrixPair(const SparseMatrix &matrix_a, const SparseMatrix &matrix_b)
Definition sparse.h:249
Fractional ComputeInfinityNorm() const
Definition sparse.cc:433
void PopulateFromMatrix(const SparseMatrix &matrix)
Takes all the columns of the given matrix.
Definition sparse.h:238
bool IsEmpty() const
Same behavior as the SparseMatrix functions above.
Definition sparse.h:273
Fractional ComputeOneNorm() const
Definition sparse.cc:430
void PopulateFromBasis(const MatrixView &matrix, const RowToColMapping &basis)
Takes only the columns of the given matrix that belongs to the given basis.
Definition sparse.h:263
const SparseColumn & column(ColIndex col) const
Definition sparse.h:276
void ComputeMinAndMaxMagnitudes(Fractional *min_magnitude, Fractional *max_magnitude) const
Definition sparse.cc:379
void PopulateFromProduct(const SparseMatrix &a, const SparseMatrix &b)
Multiplies SparseMatrix a by SparseMatrix b.
Definition sparse.cc:260
std::string Dump() const
Returns a dense representation of the matrix.
Definition sparse.cc:409
void PopulateFromPermutedMatrix(const Matrix &a, const RowPermutation &row_perm, const ColumnPermutation &inverse_col_perm)
Definition sparse.cc:222
const SparseColumn & column(ColIndex col) const
Access the underlying sparse columns.
Definition sparse.h:194
SparseMatrix & operator=(const SparseMatrix &)=delete
void PopulateFromTranspose(const Matrix &input)
Instantiate needed templates.
Definition sparse.cc:191
bool CheckNoDuplicates() const
Call CheckNoDuplicates() on all columns, useful for doing a DCHECK.
Definition sparse.cc:136
void PopulateFromLinearCombination(Fractional alpha, const SparseMatrix &a, Fractional beta, const SparseMatrix &b)
Definition sparse.cc:235
SparseColumn * mutable_column(ColIndex col)
Definition sparse.h:195
void AppendUnitVector(RowIndex row, Fractional value)
Definition sparse.cc:161
void Swap(SparseMatrix *matrix)
Definition sparse.cc:168
bool IsCleanedUp() const
Call IsCleanedUp() on all columns, useful for doing a DCHECK.
Definition sparse.cc:145
void ApplyRowPermutation(const RowPermutation &row_perm)
Applies the row permutation.
Definition sparse.cc:326
void DeleteRows(RowIndex num_rows, const RowPermutation &permutation)
Definition sparse.cc:299
void PopulateFromIdentity(ColIndex num_cols)
Definition sparse.cc:182
RowIndex num_rows() const
Return the matrix dimension.
Definition sparse.h:190
Fractional ComputeInfinityNorm() const
Definition sparse.cc:405
void PopulateFromZero(RowIndex num_rows, ColIndex num_cols)
Definition sparse.cc:174
bool AppendRowsFromSparseMatrix(const SparseMatrix &matrix)
Definition sparse.cc:312
void DeleteColumns(const DenseBooleanRow &columns_to_delete)
Definition sparse.cc:286
void PopulateFromSparseMatrix(const SparseMatrix &matrix)
Definition sparse.cc:216
bool Equals(const SparseMatrix &a, Fractional tolerance) const
Definition sparse.cc:337
ColIndex AppendEmptyColumn()
Appends an empty column and returns its index.
Definition sparse.cc:155
void SetNumRows(RowIndex num_rows)
Change the number of row of this matrix.
Definition sparse.cc:153
Fractional LookUpValue(RowIndex row, ColIndex col) const
Definition sparse.cc:333
StrictITISpan< ColIndex, const Fractional > ConstView
Definition lp_types.h:291
Fractional ComputeInverseInfinityNormUpperBound() const
Definition sparse.cc:1516
void PermutedComputeRowsToConsider(const ColumnView &rhs, const RowPermutation &row_perm, RowIndexVector *lower_column_rows, RowIndexVector *upper_column_rows)
Definition sparse.cc:1273
StrictITIVector< RowIndex, ColIndex > RowToColMapping
Definition lp_types.h:394
std::vector< RowIndex > RowIndexVector
Definition lp_types.h:361
StrictITIVector< RowIndex, bool > DenseBooleanColumn
Column of booleans.
Definition lp_types.h:383
Permutation< ColIndex > ColumnPermutation
Permutation< RowIndex > RowPermutation
ColIndex RowToColIndex(RowIndex row)
Get the ColIndex corresponding to the column # row.
Definition lp_types.h:54
StrictITIVector< RowIndex, RowIndex > RowMapping
Column of row indices. Used to represent mappings between rows.
Definition lp_types.h:389
StrictITIVector< RowIndex, Fractional > DenseColumn
Column-vector types. Column-vector types are indexed by a row index.
Definition lp_types.h:380
StrictITIVector< ColIndex, Fractional > DenseRow
Row-vector types. Row-vector types are indexed by a column index.
Definition lp_types.h:351
StrictITIVector< ColIndex, bool > DenseBooleanRow
Row of booleans.
Definition lp_types.h:354
In SWIG mode, we don't want anything besides these top-level includes.
static int input(yyscan_t yyscanner)
#define RETURN_IF_NULL(x)