27 const Fractional EtaMatrix::kSparseThreshold = 0.5;
35 eta_coeff_ = direction.
values;
40 kSparseThreshold * eta_coeff_.
size().value()) {
54 if (!sparse_eta_coeff_.
IsEmpty()) {
55 LeftSolveWithSparseEta(y);
57 LeftSolveWithDenseEta(y);
68 if (!sparse_eta_coeff_.
IsEmpty()) {
69 RightSolveWithSparseEta(d);
71 RightSolveWithDenseEta(d);
80 bool is_eta_col_in_pos =
false;
81 const int size = pos->size();
82 for (
int i = 0; i < size; ++i) {
83 const ColIndex
col = (*pos)[i];
85 if (
col == eta_col_) {
86 is_eta_col_in_pos =
true;
89 y_value -= (*y)[
col] * eta_coeff_[
row];
92 (*y)[eta_col_] = y_value / eta_col_coefficient_;
95 if (!is_eta_col_in_pos) pos->push_back(eta_col_);
98 void EtaMatrix::LeftSolveWithDenseEta(
DenseRow* y)
const {
100 const RowIndex num_rows(eta_coeff_.
size());
101 for (RowIndex
row(0);
row < num_rows; ++
row) {
104 (*y)[eta_col_] = y_value / eta_col_coefficient_;
107 void EtaMatrix::LeftSolveWithSparseEta(
DenseRow* y)
const {
112 (*y)[eta_col_] = y_value / eta_col_coefficient_;
115 void EtaMatrix::RightSolveWithDenseEta(
DenseColumn* d)
const {
117 const Fractional coeff = (*d)[eta_row] / eta_col_coefficient_;
118 const RowIndex num_rows(eta_coeff_.
size());
119 for (RowIndex
row(0);
row < num_rows; ++
row) {
120 (*d)[
row] -= eta_coeff_[
row] * coeff;
122 (*d)[eta_row] = coeff;
125 void EtaMatrix::RightSolveWithSparseEta(
DenseColumn* d)
const {
127 const Fractional coeff = (*d)[eta_row] / eta_col_coefficient_;
129 (*d)[e.row()] -= e.coefficient() * coeff;
131 (*d)[eta_row] = coeff;
144 RowIndex leaving_variable_row,
146 const ColIndex leaving_variable_col =
RowToColIndex(leaving_variable_row);
148 new EtaMatrix(leaving_variable_col, direction);
149 eta_matrix_.push_back(eta_factorization);
154 for (
int i = eta_matrix_.size() - 1; i >= 0; --i) {
155 eta_matrix_[i]->LeftSolve(y);
161 for (
int i = eta_matrix_.size() - 1; i >= 0; --i) {
162 eta_matrix_[i]->SparseLeftSolve(y, pos);
168 const size_t num_eta_matrices = eta_matrix_.size();
169 for (
int i = 0; i < num_eta_matrices; ++i) {
170 eta_matrix_[i]->RightSolve(d);
180 compact_matrix_(*compact_matrix),
182 tau_is_computed_(false),
185 eta_factorization_(),
187 deterministic_time_(0.0) {
196 tau_computation_can_be_optimized_ =
false;
197 eta_factorization_.
Clear();
198 lu_factorization_.
Clear();
199 rank_one_factorization_.
Clear();
223 stats_.refactorization_interval.Add(num_updates_);
228 const double kLuComplexityFactor = 10;
229 deterministic_time_ +=
243 Status BasisFactorization::MiddleProductFormUpdate(
244 ColIndex entering_col, RowIndex leaving_variable_row) {
245 const ColIndex right_index = right_pool_mapping_[entering_col];
246 const ColIndex left_index =
249 LOG(
INFO) <<
"One update vector is missing!!!";
257 for (
const EntryIndex i : right_storage_.
Column(right_index)) {
260 scratchpad_non_zeros_.push_back(
row);
263 const SparseColumn& column_of_u =
266 scratchpad_[e.row()] -= e.coefficient();
267 scratchpad_non_zeros_.push_back(e.row());
274 &scratchpad_, &scratchpad_non_zeros_);
275 RankOneUpdateElementaryMatrix elementary_update_matrix(
276 &storage_, u_index, left_index, scalar_product);
277 if (elementary_update_matrix.IsSingular()) {
280 rank_one_factorization_.
Update(elementary_update_matrix);
285 RowIndex leaving_variable_row,
287 if (num_updates_ < max_num_updates_) {
294 if (use_middle_product_form_update_) {
296 MiddleProductFormUpdate(entering_col, leaving_variable_row));
298 eta_factorization_.
Update(entering_col, leaving_variable_row, direction);
300 tau_computation_can_be_optimized_ =
false;
309 BumpDeterministicTimeForSolve(compact_matrix_.
num_rows().value());
310 if (use_middle_product_form_update_) {
325 BumpDeterministicTimeForSolve(d->
non_zeros.size());
326 if (use_middle_product_form_update_) {
341 BumpDeterministicTimeForSolve(compact_matrix_.
num_rows().value());
342 if (use_middle_product_form_update_) {
343 if (tau_computation_can_be_optimized_) {
346 tau_computation_can_be_optimized_ =
false;
360 tau_is_computed_ =
true;
368 BumpDeterministicTimeForSolve(1);
371 if (!use_middle_product_form_update_) {
404 if (tau_is_computed_) {
405 tau_computation_can_be_optimized_ =
408 tau_computation_can_be_optimized_ =
false;
411 tau_is_computed_ =
false;
420 BumpDeterministicTimeForSolve(1);
432 BumpDeterministicTimeForSolve(
436 if (!use_middle_product_form_update_) {
447 if (
col >= right_pool_mapping_.
size()) {
459 right_pool_mapping_[
col] =
470 BumpDeterministicTimeForSolve(
a.num_entries().value());
477 BumpDeterministicTimeForSolve(1);
481 bool BasisFactorization::IsIdentityBasis()
const {
482 const RowIndex num_rows = compact_matrix_.
num_rows();
483 for (RowIndex
row(0);
row < num_rows; ++
row) {
484 const ColIndex
col = basis_[
row];
488 if (entry_row !=
row || coeff != 1.0)
return false;
494 if (IsIdentityBasis())
return 1.0;
500 if (IsIdentityBasis())
return 1.0;
509 if (IsIdentityBasis())
return 1.0;
510 const RowIndex num_rows = compact_matrix_.
num_rows();
513 for (ColIndex
col(0);
col < num_cols; ++
col) {
521 for (RowIndex
row(0);
row < num_rows; ++
row) {
522 column_norm += std::abs(right_hand_side[
row]);
531 if (IsIdentityBasis())
return 1.0;
532 const RowIndex num_rows = compact_matrix_.
num_rows();
535 for (ColIndex
col(0);
col < num_cols; ++
col) {
542 for (RowIndex
row(0);
row < num_rows; ++
row) {
543 row_sum[
row] += std::abs(right_hand_side[
row]);
548 for (RowIndex
row(0);
row < num_rows; ++
row) {
555 if (IsIdentityBasis())
return 1.0;
560 if (IsIdentityBasis())
return 1.0;
566 if (IsIdentityBasis())
return 1.0;
567 BumpDeterministicTimeForSolve(compact_matrix_.
num_rows().value());
573 return deterministic_time_;
576 void BasisFactorization::BumpDeterministicTimeForSolve(
int num_entries)
const {
578 if (compact_matrix_.
num_rows().value() == 0)
return;
579 const double density =
581 static_cast<double>(compact_matrix_.
num_rows().value());
582 deterministic_time_ +=
#define DCHECK_NE(val1, val2)
#define DCHECK(condition)
#define DCHECK_EQ(val1, val2)
Fractional ComputeInverseOneNorm() const
BasisFactorization(const CompactSparseMatrix *compact_matrix, const RowToColMapping *basis)
Fractional ComputeInfinityNormConditionNumberUpperBound() const
ABSL_MUST_USE_RESULT Status Refactorize()
Fractional ComputeInverseInfinityNorm() const
Fractional ComputeInfinityNorm() const
const DenseColumn & RightSolveForTau(const ScatteredColumn &a) const
void LeftSolveForUnitRow(ColIndex j, ScatteredRow *y) const
Fractional RightSolveSquaredNorm(const ColumnView &a) const
Fractional ComputeOneNorm() const
Fractional ComputeOneNormConditionNumber() const
ABSL_MUST_USE_RESULT Status Initialize()
bool IsRefactorized() const
void TemporaryLeftSolveForUnitRow(ColIndex j, ScatteredRow *y) const
ABSL_MUST_USE_RESULT Status Update(ColIndex entering_col, RowIndex leaving_variable_row, const ScatteredColumn &direction)
virtual ~BasisFactorization()
Fractional DualEdgeSquaredNorm(RowIndex row) const
void LeftSolve(ScatteredRow *y) const
void RightSolveForProblemColumn(ColIndex col, ScatteredColumn *d) const
Fractional ComputeInfinityNormConditionNumber() const
void SetParameters(const GlopParameters ¶meters)
void RightSolve(ScatteredColumn *d) const
double DeterministicTime() const
ABSL_MUST_USE_RESULT Status ForceRefactorization()
RowIndex GetFirstRow() const
Fractional GetFirstCoefficient() const
EntryIndex num_entries() const
ColIndex AddDenseColumn(const DenseColumn &dense_column)
ColIndex num_cols() const
ColIndex AddDenseColumnWithNonZeros(const DenseColumn &dense_column, const std::vector< RowIndex > &non_zeros)
::util::IntegerRange< EntryIndex > Column(ColIndex col) const
ColIndex AddAndClearColumnWithNonZeros(DenseColumn *column, std::vector< RowIndex > *non_zeros)
void ColumnCopyToClearedDenseColumnWithNonZeros(ColIndex col, DenseColumn *dense_column, RowIndexVector *non_zeros) const
ColIndex AddDenseColumnPrefix(const DenseColumn &dense_column, RowIndex start)
RowIndex num_rows() const
Fractional EntryCoefficient(EntryIndex i) const
void ColumnCopyToClearedDenseColumn(ColIndex col, DenseColumn *dense_column) const
Fractional ColumnScalarProduct(ColIndex col, const DenseRow &vector) const
void Reset(RowIndex num_rows)
RowIndex EntryRow(EntryIndex i) const
ColumnView column(ColIndex col) const
Fractional ComputeInfinityNorm() const
Fractional ComputeOneNorm() const
void RightSolve(DenseColumn *d) const
void LeftSolve(DenseRow *y) const
virtual ~EtaFactorization()
void SparseLeftSolve(DenseRow *y, ColIndexVector *pos) const
void Update(ColIndex entering_col, RowIndex leaving_variable_row, const ScatteredColumn &direction)
void RightSolve(DenseColumn *d) const
void LeftSolve(DenseRow *y) const
EtaMatrix(ColIndex eta_col, const ScatteredColumn &direction)
void SparseLeftSolve(DenseRow *y, ColIndexVector *pos) const
void LeftSolveUWithNonZeros(ScatteredRow *y) const
const SparseColumn & GetColumnOfU(ColIndex col) const
void RightSolveLForColumnView(const ColumnView &b, ScatteredColumn *x) const
void RightSolveLWithPermutedInput(const DenseColumn &a, ScatteredColumn *x) 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
ABSL_MUST_USE_RESULT Status ComputeFactorization(const CompactSparseMatrixView &compact_matrix)
EntryIndex NumberOfEntries() const
Fractional ComputeInverseInfinityNormUpperBound() const
void RightSolveWithNonZeros(ScatteredColumn *d) const
void LeftSolveWithNonZeros(ScatteredRow *y) const
void Update(const RankOneUpdateElementaryMatrix &update_matrix)
EntryIndex num_entries() const
bool CheckNoDuplicates() const
typename Iterator::Entry Entry
void SetCoefficient(Index index, Fractional value)
void AssignToZero(IntType size)
void resize(IntType size)
void assign(IntType size, const T &v)
void STLDeleteElements(T *container)
static double DeterministicTimeForFpOperations(int64 n)
std::vector< ColIndex > ColIndexVector
bool IsAllZero(const Container &input)
StrictITIVector< ColIndex, Fractional > DenseRow
ColIndex RowToColIndex(RowIndex row)
void ClearAndResizeVectorWithNonZeros(IndexType size, ScatteredRowOrCol *v)
const DenseRow & Transpose(const DenseColumn &col)
RowIndex ColToRowIndex(ColIndex col)
std::vector< RowIndex > RowIndexVector
StrictITIVector< RowIndex, Fractional > DenseColumn
const ColIndex kInvalidCol(-1)
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
#define RETURN_IF_NULL(x)
#define SCOPED_TIME_STAT(stats)
#define GLOP_RETURN_IF_ERROR(function_call)
#define GLOP_RETURN_AND_LOG_ERROR(error_code, message)
void SortNonZerosIfNeeded()
std::vector< Index > non_zeros
StrictITIVector< Index, Fractional > values