25 : is_identity_factorization_(true),
29 inverse_row_perm_() {}
33 lower_.
Reset(RowIndex(0), ColIndex(0));
34 upper_.
Reset(RowIndex(0), ColIndex(0));
35 transpose_upper_.
Reset(RowIndex(0), ColIndex(0));
36 transpose_lower_.
Reset(RowIndex(0), ColIndex(0));
37 is_identity_factorization_ =
true;
40 inverse_row_perm_.
clear();
41 inverse_col_perm_.
clear();
53 markowitz_.
ComputeLU(matrix, &row_perm_, &col_perm_, &lower_, &upper_));
56 ComputeTransposeUpper();
57 ComputeTransposeLower();
59 is_identity_factorization_ =
false;
62 stats_.basis_num_entries.Add(matrix.
num_entries().value());
70 if (is_identity_factorization_)
return;
80 if (is_identity_factorization_)
return;
95 const std::vector<RowIndex>& non_zeros,
DenseColumn* column) {
97 if (non_zeros.empty()) {
101 for (
const RowIndex
row : non_zeros) {
103 (*column)[
row] = 0.0;
114 non_zero_rows_.clear();
119 const RowIndex permuted_row = row_perm_[e.row()];
120 dense_zero_scratchpad_[permuted_row] = e.coefficient();
121 non_zero_rows_.push_back(permuted_row);
125 if (non_zero_rows_.empty()) {
131 if (non_zero_rows_.empty()) {
137 return ComputeSquaredNormAndResetToZero(non_zero_rows_,
138 &dense_zero_scratchpad_);
142 if (is_identity_factorization_)
return 1.0;
144 const RowIndex permuted_row =
147 non_zero_rows_.clear();
150 dense_zero_scratchpad_[permuted_row] = 1.0;
151 non_zero_rows_.push_back(permuted_row);
154 if (non_zero_rows_.empty()) {
156 &dense_zero_scratchpad_);
158 transpose_upper_.
HyperSparseSolve(&dense_zero_scratchpad_, &non_zero_rows_);
161 if (non_zero_rows_.empty()) {
162 transpose_lower_.
UpperSolve(&dense_zero_scratchpad_);
165 &dense_zero_scratchpad_, &non_zero_rows_);
167 return ComputeSquaredNormAndResetToZero(non_zero_rows_,
168 &dense_zero_scratchpad_);
176 const RowIndex num_rows = perm.
size();
177 for (RowIndex
row(0);
row < num_rows; ++
row) {
178 if (
a[
row] !=
b[perm[
row]])
return false;
187 if (!is_identity_factorization_) {
198 template <
typename Column>
199 void LuFactorization::RightSolveLInternal(
const Column&
b,
208 for (
const auto e :
b) {
209 const RowIndex permuted_row = row_perm_[e.row()];
210 (*x)[permuted_row] = e.coefficient();
220 first_column_to_consider =
std::min(first_column_to_consider,
col);
237 if (is_identity_factorization_) {
239 (*x)[e.row()] = e.coefficient();
245 RightSolveLInternal(
b, x);
249 if (is_identity_factorization_)
return;
273 if (is_identity_factorization_) {
278 if (
b.non_zeros.empty()) {
283 RightSolveLInternal(
b, x);
289 if (is_identity_factorization_)
return;
305 if (is_identity_factorization_)
return;
323 if (is_identity_factorization_) {
339 if (result_before_permutation ==
nullptr) {
351 for (
const RowIndex
row : *nz) {
364 for (RowIndex
row(0);
row < inverse_row_perm_.
size(); ++
row) {
367 const RowIndex permuted_row = inverse_row_perm_[
row];
368 (*x)[permuted_row] =
value;
372 nz->swap(result_before_permutation->
non_zeros);
373 nz->reserve(result_before_permutation->
non_zeros.size());
374 for (
const RowIndex
row : result_before_permutation->
non_zeros) {
376 const RowIndex permuted_row = inverse_row_perm_[
row];
377 (*x)[permuted_row] =
value;
378 nz->push_back(permuted_row);
394 if (is_identity_factorization_) {
399 const ColIndex permuted_col = col_perm_.
empty() ?
col : col_perm_[
col];
400 (*y)[permuted_col] = 1.0;
423 if (is_identity_factorization_) {
424 column_of_upper_.
Clear();
426 return column_of_upper_;
430 return column_of_upper_;
435 const int initial_num_entries = matrix.
num_entries().value();
436 const int lu_num_entries =
438 if (is_identity_factorization_ || initial_num_entries == 0)
return 1.0;
439 return static_cast<double>(lu_num_entries) /
440 static_cast<double>(initial_num_entries);
444 return is_identity_factorization_
450 if (is_identity_factorization_)
return 1.0;
461 if (is_identity_factorization_)
return 1.0;
462 const RowIndex num_rows = lower_.
num_rows();
463 const ColIndex num_cols = lower_.
num_cols();
465 for (ColIndex
col(0);
col < num_cols; ++
col) {
472 for (RowIndex
row(0);
row < num_rows; ++
row) {
473 column_norm += std::abs(right_hand_side[
row]);
482 if (is_identity_factorization_)
return 1.0;
483 const RowIndex num_rows = lower_.
num_rows();
484 const ColIndex num_cols = lower_.
num_cols();
486 for (ColIndex
col(0);
col < num_cols; ++
col) {
492 for (RowIndex
row(0);
row < num_rows; ++
row) {
493 row_sum[
row] += std::abs(right_hand_side[
row]);
498 for (RowIndex
row(0);
row < num_rows; ++
row) {
506 if (is_identity_factorization_)
return 1.0;
512 if (is_identity_factorization_)
return 1.0;
524 double density = 0.0;
526 if (row_perm[e.row()] !=
kNonPivotal && e.coefficient() != 0.0) {
530 const RowIndex num_rows = row_perm.
size();
531 return density / num_rows.value();
535 void LuFactorization::ComputeTransposeUpper() {
540 void LuFactorization::ComputeTransposeLower()
const {
545 bool LuFactorization::CheckFactorization(
const CompactSparseMatrixView& matrix,
547 if (is_identity_factorization_)
return true;
551 paq.PopulateFromPermutedMatrix(matrix, row_perm_, inverse_col_perm_);
552 if (!row_perm_.
Check()) {
555 if (!inverse_col_perm_.
Check()) {
559 SparseMatrix should_be_zero;
560 should_be_zero.PopulateFromLinearCombination(
Fractional(1.0), paq,
563 for (ColIndex
col(0);
col < should_be_zero.num_cols(); ++
col) {
565 const Fractional magnitude = std::abs(e.coefficient());
566 if (magnitude > tolerance) {
567 VLOG(2) << magnitude <<
" != 0, at column " <<
col;
#define DCHECK_NE(val1, val2)
#define DCHECK(condition)
#define DCHECK_EQ(val1, val2)
#define VLOG(verboselevel)
void swap(StrongVector &x)
Fractional ComputeInfinityNorm() const
ColIndex num_cols() const
Fractional ComputeOneNorm() const
RowIndex num_rows() const
EntryIndex num_entries() const
Fractional ComputeInverseOneNorm() const
void LeftSolveUWithNonZeros(ScatteredRow *y) const
const SparseColumn & GetColumnOfU(ColIndex col) const
Fractional ComputeInverseInfinityNorm() const
void RightSolveLForColumnView(const ColumnView &b, ScatteredColumn *x) const
void RightSolveLWithPermutedInput(const DenseColumn &a, ScatteredColumn *x) const
double GetFillInPercentage(const CompactSparseMatrixView &matrix) const
Fractional RightSolveSquaredNorm(const ColumnView &a) const
void RightSolveUWithNonZeros(ScatteredColumn *x) const
void LeftSolve(DenseRow *y) const
bool LeftSolveLWithNonZeros(ScatteredRow *y, ScatteredColumn *result_before_permutation) const
void RightSolve(DenseColumn *x) const
ColIndex LeftSolveUForUnitRow(ColIndex col, ScatteredRow *y) const
Fractional DualEdgeSquaredNorm(RowIndex row) const
void RightSolveLForScatteredColumn(const ScatteredColumn &b, ScatteredColumn *x) const
void RightSolveLWithNonZeros(ScatteredColumn *x) const
Fractional ComputeDeterminant() const
Fractional ComputeInfinityNormConditionNumber(const CompactSparseMatrixView &matrix) const
void ComputeLowerTimesUpper(SparseMatrix *product) const
ABSL_MUST_USE_RESULT Status ComputeFactorization(const CompactSparseMatrixView &compact_matrix)
EntryIndex NumberOfEntries() const
Fractional ComputeInverseInfinityNormUpperBound() const
Fractional ComputeOneNormConditionNumber(const CompactSparseMatrixView &matrix) const
ABSL_MUST_USE_RESULT Status ComputeLU(const CompactSparseMatrixView &basis_matrix, RowPermutation *row_perm, ColumnPermutation *col_perm, TriangularMatrix *lower, TriangularMatrix *upper)
int ComputeSignature() const
void PopulateFromInverse(const Permutation &inverse)
typename Iterator::Entry Entry
void SetCoefficient(Index index, Fractional value)
void resize(IntType size)
void TransposeHyperSparseSolve(DenseColumn *rhs, RowIndexVector *non_zero_rows) const
void UpperSolve(DenseColumn *rhs) const
void HyperSparseSolve(DenseColumn *rhs, RowIndexVector *non_zero_rows) const
ColIndex GetFirstNonIdentityColumn() const
ColIndex num_cols() const
Fractional GetDiagonalCoefficient(ColIndex col) const
void LowerSolve(DenseColumn *rhs) const
void TransposeHyperSparseSolveWithReversedNonZeros(DenseColumn *rhs, RowIndexVector *non_zero_rows) const
void LowerSolveStartingAt(ColIndex start, DenseColumn *rhs) const
void PopulateFromTranspose(const TriangularMatrix &input)
void CopyColumnToSparseColumn(ColIndex col, SparseColumn *output) const
RowIndex num_rows() const
void ComputeRowsToConsiderInSortedOrder(RowIndexVector *non_zero_rows, Fractional sparsity_ratio, Fractional num_ops_ratio) const
void TransposeLowerSolve(DenseColumn *rhs) const
void TransposeUpperSolve(DenseColumn *rhs) const
Fractional ComputeInverseInfinityNormUpperBound() const
void Reset(RowIndex num_rows, ColIndex col_capacity)
bool ColumnIsDiagonalOnly(ColIndex col) const
void HyperSparseSolveWithReversedNonZeros(DenseColumn *rhs, RowIndexVector *non_zero_rows) const
EntryIndex num_entries() const
void PermuteWithScratchpad(const Permutation< PermutationIndexType > &permutation, StrictITIVector< IndexType, Fractional > *zero_scratchpad, StrictITIVector< IndexType, Fractional > *input_output)
Fractional Square(Fractional f)
Fractional SquaredNorm(const SparseColumn &v)
void ApplyInversePermutation(const Permutation< IndexType > &perm, const ITIVectorType &b, ITIVectorType *result)
bool IsAllZero(const Container &input)
void PermuteWithKnownNonZeros(const Permutation< IndexType > &permutation, StrictITIVector< IndexType, Fractional > *zero_scratchpad, StrictITIVector< IndexType, Fractional > *output, std::vector< IndexType > *non_zeros)
ColIndex RowToColIndex(RowIndex row)
void ClearAndResizeVectorWithNonZeros(IndexType size, ScatteredRowOrCol *v)
RowIndex ColToRowIndex(ColIndex col)
const RowIndex kNonPivotal(-1)
std::vector< RowIndex > RowIndexVector
void ApplyPermutation(const Permutation< IndexType > &perm, const ITIVectorType &b, ITIVectorType *result)
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
#define IF_STATS_ENABLED(instructions)
#define SCOPED_TIME_STAT(stats)
#define GLOP_RETURN_IF_ERROR(function_call)
#define GLOP_RETURN_AND_LOG_ERROR(error_code, message)
bool non_zeros_are_sorted
std::vector< Index > non_zeros
StrictITIVector< Index, Fractional > values