Google OR-Tools v9.11
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
operations_research::glop::TriangularMatrix Class Reference

#include <sparse.h>

Inheritance diagram for operations_research::glop::TriangularMatrix:
operations_research::glop::CompactSparseMatrix

Public Member Functions

 TriangularMatrix ()
 
 TriangularMatrix (const TriangularMatrix &)=delete
 This type is neither copyable nor movable.
 
TriangularMatrixoperator= (const TriangularMatrix &)=delete
 
void PopulateFromTranspose (const TriangularMatrix &input)
 
void Swap (TriangularMatrix *other)
 
bool IsEmpty () const
 
RowIndex num_rows () const
 
ColIndex num_cols () const
 
EntryIndex num_entries () const
 
void Reset (RowIndex num_rows, ColIndex col_capacity)
 
void PopulateFromTriangularSparseMatrix (const SparseMatrix &input)
 
void AddTriangularColumn (const ColumnView &column, RowIndex diagonal_row)
 
void AddTriangularColumnWithGivenDiagonalEntry (const SparseColumn &column, RowIndex diagonal_row, Fractional diagonal_value)
 
void AddDiagonalOnlyColumn (Fractional diagonal_value)
 
void AddAndNormalizeTriangularColumn (const SparseColumn &column, RowIndex diagonal_row, Fractional diagonal_coefficient)
 
void ApplyRowPermutationToNonDiagonalEntries (const RowPermutation &row_perm)
 Applies the given row permutation to all entries except the diagonal ones.
 
void CopyColumnToSparseColumn (ColIndex col, SparseColumn *output) const
 Copy a triangular column with its diagonal entry to the given SparseColumn.
 
void CopyToSparseMatrix (SparseMatrix *output) const
 Copy a triangular matrix to the given SparseMatrix.
 
ColIndex GetFirstNonIdentityColumn () const
 
Fractional GetDiagonalCoefficient (ColIndex col) const
 Returns the diagonal coefficient of the given column.
 
bool ColumnIsDiagonalOnly (ColIndex col) const
 Returns true iff the column contains no non-diagonal entries.
 
void LowerSolve (DenseColumn *rhs) const
 
void UpperSolve (DenseColumn *rhs) const
 Solves the system U.x = rhs for an upper triangular matrix.
 
void TransposeUpperSolve (DenseColumn *rhs) const
 
void LowerSolveStartingAt (ColIndex start, DenseColumn *rhs) const
 This assumes that the rhs is all zero before the given position.
 
void TransposeLowerSolve (DenseColumn *rhs) const
 
void HyperSparseSolve (DenseColumn *rhs, RowIndexVector *non_zero_rows) const
 
void HyperSparseSolveWithReversedNonZeros (DenseColumn *rhs, RowIndexVector *non_zero_rows) const
 
void TransposeHyperSparseSolve (DenseColumn *rhs, RowIndexVector *non_zero_rows) const
 
void TransposeHyperSparseSolveWithReversedNonZeros (DenseColumn *rhs, RowIndexVector *non_zero_rows) const
 
void ComputeRowsToConsiderWithDfs (RowIndexVector *non_zero_rows) const
 
void ComputeRowsToConsiderInSortedOrder (RowIndexVector *non_zero_rows) const
 
void PermutedLowerSolve (const SparseColumn &rhs, const RowPermutation &row_perm, const RowMapping &partial_inverse_row_perm, SparseColumn *lower, SparseColumn *upper) const
 
void PermutedLowerSparseSolve (const ColumnView &rhs, const RowPermutation &row_perm, SparseColumn *lower, SparseColumn *upper)
 
int64_t NumFpOperationsInLastPermutedLowerSparseSolve () const
 This is used to compute the deterministic time of a matrix factorization.
 
bool IsLowerTriangular () const
 
bool IsUpperTriangular () const
 
void PermutedComputeRowsToConsider (const ColumnView &rhs, const RowPermutation &row_perm, RowIndexVector *lower_column_rows, RowIndexVector *upper_column_rows)
 
Fractional ComputeInverseInfinityNormUpperBound () const
 
Fractional ComputeInverseInfinityNorm () const
 

Detailed Description

Specialization of a CompactSparseMatrix used for triangular matrices. To be able to solve triangular systems as efficiently as possible, the diagonal entries are stored in a separate vector and not in the underlying CompactSparseMatrix.

Advanced usage: this class also support matrices that can be permuted into a triangular matrix and some functions work directly on such matrices.

Definition at line 597 of file sparse.h.

Constructor & Destructor Documentation

◆ TriangularMatrix() [1/2]

operations_research::glop::TriangularMatrix::TriangularMatrix ( )
inline

Definition at line 599 of file sparse.h.

◆ TriangularMatrix() [2/2]

operations_research::glop::TriangularMatrix::TriangularMatrix ( const TriangularMatrix & )
delete

This type is neither copyable nor movable.

Member Function Documentation

◆ AddAndNormalizeTriangularColumn()

void operations_research::glop::TriangularMatrix::AddAndNormalizeTriangularColumn ( const SparseColumn & column,
RowIndex diagonal_row,
Fractional diagonal_coefficient )

Adds the given sparse column divided by diagonal_coefficient. The diagonal_row is assumed to be present and its value should be the same as the one given in diagonal_coefficient. Note that this function tests for zero coefficients in the input column and removes them.

Todo
(user): use division by a constant using multiplication.

Definition at line 697 of file sparse.cc.

◆ AddDiagonalOnlyColumn()

void operations_research::glop::TriangularMatrix::AddDiagonalOnlyColumn ( Fractional diagonal_value)

Definition at line 678 of file sparse.cc.

◆ AddTriangularColumn()

void operations_research::glop::TriangularMatrix::AddTriangularColumn ( const ColumnView & column,
RowIndex diagonal_row )

Functions to create a triangular matrix incrementally, column by column. A client needs to call Reset(num_rows) first, and then each column must be added by calling one of the 3 functions below.

Note
the row indices of the columns are allowed to be permuted: the diagonal entry of the column col not being necessarily on the row col. This is why these functions require the 'diagonal_row' parameter. The permutation can be fixed at the end by a call to ApplyRowPermutationToNonDiagonalEntries() or accounted directly in the case of PermutedLowerSparseSolve().

Definition at line 682 of file sparse.cc.

◆ AddTriangularColumnWithGivenDiagonalEntry()

void operations_research::glop::TriangularMatrix::AddTriangularColumnWithGivenDiagonalEntry ( const SparseColumn & column,
RowIndex diagonal_row,
Fractional diagonal_value )

Definition at line 714 of file sparse.cc.

◆ ApplyRowPermutationToNonDiagonalEntries()

void operations_research::glop::TriangularMatrix::ApplyRowPermutationToNonDiagonalEntries ( const RowPermutation & row_perm)

Applies the given row permutation to all entries except the diagonal ones.

Definition at line 754 of file sparse.cc.

◆ ColumnIsDiagonalOnly()

bool operations_research::glop::TriangularMatrix::ColumnIsDiagonalOnly ( ColIndex col) const
inline

Returns true iff the column contains no non-diagonal entries.

Definition at line 676 of file sparse.h.

◆ ComputeInverseInfinityNorm()

Fractional operations_research::glop::TriangularMatrix::ComputeInverseInfinityNorm ( ) const

Get the col-th column of the matrix inverse.

Compute sum_j inverse_ij.

Compute max_i sum_j inverse_ij.

Definition at line 1521 of file sparse.cc.

◆ ComputeInverseInfinityNormUpperBound()

Fractional operations_research::glop::TriangularMatrix::ComputeInverseInfinityNormUpperBound ( ) const

The upper bound is computed using one of the algorithm presented in "A Survey of Condition Number Estimation for Triangular Matrices" https:epubs.siam.org/doi/pdf/10.1137/1029112/

A known upper bound for the infinity norm of T^{-1} is the infinity norm of y where T'*y = x with:

  • x the all 1s vector.
  • Each entry in T' is the absolute value of the same entry in T.

Identity matrix

Definition at line 1493 of file sparse.cc.

◆ ComputeRowsToConsiderInSortedOrder()

void operations_research::glop::TriangularMatrix::ComputeRowsToConsiderInSortedOrder ( RowIndexVector * non_zero_rows) const

Same as TriangularComputeRowsToConsider() but always returns the non-zeros sorted by rows. It is up to the client to call the direct or reverse hyper-sparse solve function depending if the matrix is upper or lower triangular.

Todo
(user): Investigate the best thresholds.

Definition at line 1448 of file sparse.cc.

◆ ComputeRowsToConsiderWithDfs()

void operations_research::glop::TriangularMatrix::ComputeRowsToConsiderWithDfs ( RowIndexVector * non_zero_rows) const

Given the positions of the non-zeros of a vector, computes the non-zero positions of the vector after a solve by this triangular matrix. The order of the returned non-zero positions will be in the REVERSE elimination order. If the function detects that there are too many non-zeros, then it aborts early and non_zero_rows is cleared.

We don't start the DFS if the initial number of non-zeros is under the sparsity_threshold. During the DFS, we abort it if the number of floating points operations get larger than the num_ops_threshold.

In both cases, we make sure to clear non_zero_rows so that the solving part will use the non-hypersparse version of the code.

Todo
(user): Investigate the best thresholds.

Initialize using the non-zero positions of the input.

Topological sort based on Depth-First-Search. Same remarks as the version implemented in PermutedComputeRowsToConsider().

If the depth-first search from the current node is finished, we store the node. This will store the node in reverse topological order.

If the node is already stored, skip.

Go one level forward in the depth-first search, and store the 'adjacent' node on nodes_to_explore_ for further processing.

We reverse the sign of nodes_to_explore_.back() to detect when the DFS will be back on this node.

Abort if the number of operations is not negligible compared to the number of rows. Note that this test also prevents the code from cycling in case the matrix is actually not triangular.

Clear stored_.

If we aborted, clear the result.

Definition at line 1370 of file sparse.cc.

◆ CopyColumnToSparseColumn()

void operations_research::glop::TriangularMatrix::CopyColumnToSparseColumn ( ColIndex col,
SparseColumn * output ) const

Copy a triangular column with its diagonal entry to the given SparseColumn.

Definition at line 762 of file sparse.cc.

◆ CopyToSparseMatrix()

void operations_research::glop::TriangularMatrix::CopyToSparseMatrix ( SparseMatrix * output) const

Copy a triangular matrix to the given SparseMatrix.

Definition at line 774 of file sparse.cc.

◆ GetDiagonalCoefficient()

Fractional operations_research::glop::TriangularMatrix::GetDiagonalCoefficient ( ColIndex col) const
inline

Returns the diagonal coefficient of the given column.

Definition at line 671 of file sparse.h.

◆ GetFirstNonIdentityColumn()

ColIndex operations_research::glop::TriangularMatrix::GetFirstNonIdentityColumn ( ) const
inline

Returns the index of the first column which is not an identity column (i.e. a column j with only one entry of value 1 at the j-th row). This is always zero if the matrix is not triangular.

Definition at line 666 of file sparse.h.

◆ HyperSparseSolve()

void operations_research::glop::TriangularMatrix::HyperSparseSolve ( DenseColumn * rhs,
RowIndexVector * non_zero_rows ) const

Hyper-sparse version of the triangular solve functions. The passed non_zero_rows should contain the positions of the symbolic non-zeros of the result in the order in which they need to be accessed (or in the reverse order for the Reverse*() versions).

The non-zero vector is mutable so that the symbolic non-zeros that are actually zero because of numerical cancellations can be removed.

The non-zeros can be computed by one of these two methods:

Note
if the non-zeros are given in a sorted order, then the hyper-sparse functions will return EXACTLY the same results as the non hyper-sparse version above.

For a given solve, here is the required order:

  • For a lower solve, increasing non-zeros order.
  • For an upper solve, decreasing non-zeros order.
  • for a transpose lower solve, decreasing non-zeros order.
  • for a transpose upper solve, increasing non_zeros order.

For a general discussion of hyper-sparsity in LP, see: J.A.J. Hall, K.I.M. McKinnon, "Exploiting hyper-sparsity in the revised simplex method", December 1999, MS 99-014. http://www.maths.ed.ac.uk/hall/MS-99/MS9914.pdf

Definition at line 960 of file sparse.cc.

◆ HyperSparseSolveWithReversedNonZeros()

void operations_research::glop::TriangularMatrix::HyperSparseSolveWithReversedNonZeros ( DenseColumn * rhs,
RowIndexVector * non_zero_rows ) const

Definition at line 992 of file sparse.cc.

◆ IsEmpty()

bool operations_research::glop::TriangularMatrix::IsEmpty ( ) const
inline

Definition at line 610 of file sparse.h.

◆ IsLowerTriangular()

bool operations_research::glop::TriangularMatrix::IsLowerTriangular ( ) const

To be used in DEBUG mode by the client code. This check that the matrix is lower- (resp. upper-) triangular without any permutation and that there is no zero on the diagonal. We can't do that on each Solve() that require so, otherwise it will be too slow in debug.

Definition at line 734 of file sparse.cc.

◆ IsUpperTriangular()

bool operations_research::glop::TriangularMatrix::IsUpperTriangular ( ) const

Definition at line 744 of file sparse.cc.

◆ LowerSolve()

void operations_research::glop::TriangularMatrix::LowerSolve ( DenseColumn * rhs) const

Triangular solve functions.

All the functions containing the word Lower (resp. Upper) require the matrix to be lower (resp. upper_) triangular without any permutation. Solve the system L.x = rhs for a lower triangular matrix. The result overwrite rhs.

Definition at line 781 of file sparse.cc.

◆ LowerSolveStartingAt()

void operations_research::glop::TriangularMatrix::LowerSolveStartingAt ( ColIndex start,
DenseColumn * rhs ) const

This assumes that the rhs is all zero before the given position.

Definition at line 785 of file sparse.cc.

◆ num_cols()

ColIndex operations_research::glop::TriangularMatrix::num_cols ( ) const
inline

Definition at line 612 of file sparse.h.

◆ num_entries()

EntryIndex operations_research::glop::TriangularMatrix::num_entries ( ) const
inline

Definition at line 613 of file sparse.h.

◆ num_rows()

RowIndex operations_research::glop::TriangularMatrix::num_rows ( ) const
inline

Definition at line 611 of file sparse.h.

◆ NumFpOperationsInLastPermutedLowerSparseSolve()

int64_t operations_research::glop::TriangularMatrix::NumFpOperationsInLastPermutedLowerSparseSolve ( ) const
inline

This is used to compute the deterministic time of a matrix factorization.

Definition at line 799 of file sparse.h.

◆ operator=()

TriangularMatrix & operations_research::glop::TriangularMatrix::operator= ( const TriangularMatrix & )
delete

◆ PermutedComputeRowsToConsider()

void operations_research::glop::TriangularMatrix::PermutedComputeRowsToConsider ( const ColumnView & rhs,
const RowPermutation & row_perm,
RowIndexVector * lower_column_rows,
RowIndexVector * upper_column_rows )

Visible for testing. This is used by PermutedLowerSparseSolve() to compute the non-zero indices of the result. The output is as follow:

  • lower_column_rows will contains the rows for which row_perm[row] < 0.
  • upper_column_rows will contains the other rows in the reverse topological order in which they should be considered in PermutedLowerSparseSolve().

This function is non-const because it prunes the underlying dependency graph of the lower matrix while doing a solve. See marked_ and pruned_ends_ below.

Pruning the graph at the same time is slower but not by too much (< 2x) and seems worth doing. Note that when the lower matrix is dense, most of the graph will likely be pruned. As a result, the symbolic phase will be negligible compared to the numerical phase so we don't really need a dense version of PermutedLowerSparseSolve().

The goal is to find which rows of the working column we will need to look at in PermutedLowerSparseSolve() when solving P^{-1}.L.P.x = rhs, 'P' being a row permutation, 'L' a lower triangular matrix and 'this' being 'P^{-1}.L'.

Note
the columns of L that are identity columns (this is the case for the ones corresponding to a kNonPivotal in P) can be skipped since they will leave the working column unchanged.

Let G denote the graph G = (V,E) of the column-to-row adjacency of A:

  • 'V' is the set of nodes, one node i corresponds to a both a row and a column (the matrix is square).
  • 'E' is the set of arcs. There is an arc from node i to node j iff the coefficient of i-th column, j-th row of A = P^{-1}.L.P is non zero.

Let S denote the set of nodes i such that rhs_i != 0. Let R denote the set of all accessible nodes from S in G. x_k is possibly non-zero iff k is in R, i.e. if k is not in R then x_k = 0 for sure, and there is no need to look a the row k during the solve.

So, to solve P^{-1}.L.P.x = rhs, only rows corresponding to P.R have to be considered (ignoring the one that map to identity column of L). A topological sort of P.R is used to decide in which order one should iterate on them. This will be given by upper_column_rows_ and it will be populated in reverse order.

Topological sort based on Depth-First-Search. A few notes:

  • By construction, if the matrix can be permuted into a lower triangular form, there is no cycle. This code does nothing to test for cycles, but there is a DCHECK() to detect them during debugging.
  • This version uses sentinels (kInvalidRow) on nodes_to_explore_ to know when a node has been explored (i.e. when the recursive dfs goes back in the call stack). This is faster than an alternate implementation that uses another Boolean array to detect when we go back in the depth-first search.

If the depth-first search from the current node is finished (i.e. there is a sentinel on the stack), we store the node (which is just before on the stack). This will store the nodes in reverse topological order.

Unmark and prune the nodes that are already unmarked. See the header comment on marked_ for the algorithm description.

Complexity note: The only difference with the "normal" DFS doing no pruning is this extra loop here and the marked_[entry_row] = true in the loop later in this function. On an already pruned graph, this is probably between 1 and 2 times slower than the "normal" DFS.

Note
we could keep the pruned row in a separate vector and not touch the triangular matrix. But the current solution seems better cache-wise and memory-wise.

If the node is already stored, skip.

Expand only if we are not on a kNonPivotal row. Otherwise we can store the node right away.

Go one level forward in the depth-first search, and store the 'adjacent' node on nodes_to_explore_ for further processing.

The graph contains cycles? this is not supposed to happen.

Clear stored_.

Definition at line 1256 of file sparse.cc.

◆ PermutedLowerSolve()

void operations_research::glop::TriangularMatrix::PermutedLowerSolve ( const SparseColumn & rhs,
const RowPermutation & row_perm,
const RowMapping & partial_inverse_row_perm,
SparseColumn * lower,
SparseColumn * upper ) const

This is currently only used for testing. It achieves the same result as PermutedLowerSparseSolve() below, but the latter exploits the sparsity of rhs and is thus faster for our use case.

Note
partial_inverse_row_perm only permutes the first k rows, where k is the same as partial_inverse_row_perm.size(). It is the inverse permutation of row_perm which only permutes k rows into is [0, k), the other row images beeing kInvalidRow. The other arguments are the same as for PermutedLowerSparseSolve() and described there.

IMPORTANT: lower will contain all the "symbolic" non-zero entries. A "symbolic" zero entry is one that will be zero whatever the coefficients of the rhs entries. That is it only depends on the position of its entries, not on their values. Thus, some of its coefficients may be zero. This fact is exploited by the LU factorization code. The zero coefficients of upper will be cleaned, however.

Definition at line 1133 of file sparse.cc.

◆ PermutedLowerSparseSolve()

void operations_research::glop::TriangularMatrix::PermutedLowerSparseSolve ( const ColumnView & rhs,
const RowPermutation & row_perm,
SparseColumn * lower,
SparseColumn * upper )

This solves a lower triangular system with only ones on the diagonal where the matrix and the input rhs are permuted by the inverse of row_perm. Note that the output will also be permuted by the inverse of row_perm. The function also supports partial permutation. That is if row_perm[i] < 0 then column row_perm[i] is assumed to be an identity column.

The output is given as follow:

  • lower is cleared, and receives the rows for which row_perm[row] < 0 meaning not yet examined as a pivot (see markowitz.cc).
  • upper is NOT cleared, and the other rows (row_perm[row] >= 0) are appended to it.
  • Note that lower and upper can point to the same SparseColumn.
Note
This function is non-const because ComputeRowsToConsider() also prunes the underlying dependency graph of the lower matrix while doing a solve. See marked_ and pruned_ends_ below.

Compute the set of rows that will be non zero in the result (lower_column, upper_column).

Copy rhs into initially_all_zero_scratchpad_.

We clear lower_column first in case upper_column and lower_column point to the same underlying SparseColumn.

rows_to_consider_ contains the row to process in reverse order. Note in particular that each "permuted_row" will never be touched again and so its value is final. We copy the result in (lower_column, upper_column) and clear initially_all_zero_scratchpad_ at the same time.

Note
permuted_row will not appear in the loop below so we already know the value of the solution at this position.
Todo
(user): The size of lower is exact, so we could be slighly faster here.

Definition at line 1176 of file sparse.cc.

◆ PopulateFromTranspose()

void operations_research::glop::TriangularMatrix::PopulateFromTranspose ( const TriangularMatrix & input)

Only a subset of the functions from CompactSparseMatrix are exposed (note the private inheritance). They are extended to deal with diagonal coefficients properly.

This takes care of the triangular special case.

The elimination structure of the transpose is not the same.

Compute first_non_identity_column_. Note that this is not necessarily the same as input.first_non_identity_column_ for an upper triangular matrix.

Definition at line 533 of file sparse.cc.

◆ PopulateFromTriangularSparseMatrix()

void operations_research::glop::TriangularMatrix::PopulateFromTriangularSparseMatrix ( const SparseMatrix & input)

Constructs a triangular matrix from the given SparseMatrix. The input is assumed to be lower or upper triangular without any permutations. This is checked in debug mode.

Definition at line 725 of file sparse.cc.

◆ Reset()

void operations_research::glop::TriangularMatrix::Reset ( RowIndex num_rows,
ColIndex col_capacity )

On top of the CompactSparseMatrix functionality, TriangularMatrix::Reset() also pre-allocates space of size col_size for a number of internal vectors. This helps reduce costly push_back operations for large problems.

Warning
Reset() must be called with a sufficiently large col_capacity prior to any Add* calls (e.g., AddTriangularColumn).

Non-zero entries in the first column always have an offset of 0.

Definition at line 566 of file sparse.cc.

◆ Swap()

void operations_research::glop::TriangularMatrix::Swap ( TriangularMatrix * other)

Definition at line 635 of file sparse.cc.

◆ TransposeHyperSparseSolve()

void operations_research::glop::TriangularMatrix::TransposeHyperSparseSolve ( DenseColumn * rhs,
RowIndexVector * non_zero_rows ) const

Definition at line 1027 of file sparse.cc.

◆ TransposeHyperSparseSolveWithReversedNonZeros()

void operations_research::glop::TriangularMatrix::TransposeHyperSparseSolveWithReversedNonZeros ( DenseColumn * rhs,
RowIndexVector * non_zero_rows ) const

Definition at line 1079 of file sparse.cc.

◆ TransposeLowerSolve()

void operations_research::glop::TriangularMatrix::TransposeLowerSolve ( DenseColumn * rhs) const

Solves the system Transpose(L).x = rhs, where L is lower triangular. This can be used to do a left-solve for a row vector (i.e., y.Y = rhs).

Definition at line 902 of file sparse.cc.

◆ TransposeUpperSolve()

void operations_research::glop::TriangularMatrix::TransposeUpperSolve ( DenseColumn * rhs) const

Solves the system Transpose(U).x = rhs where U is upper triangular. This can be used to do a left-solve for a row vector (i.e. y.Y = rhs).

Definition at line 851 of file sparse.cc.

◆ UpperSolve()

void operations_research::glop::TriangularMatrix::UpperSolve ( DenseColumn * rhs) const

Solves the system U.x = rhs for an upper triangular matrix.

Definition at line 817 of file sparse.cc.


The documentation for this class was generated from the following files: