Bug Summary

File:mlir/lib/Analysis/Presburger/Simplex.cpp
Warning:line 1014, column 17
2nd function call argument is an uninitialized value

Annotated Source Code

Press '?' to see keyboard shortcuts

clang -cc1 -cc1 -triple x86_64-pc-linux-gnu -analyze -disable-free -clear-ast-before-backend -disable-llvm-verifier -discard-value-names -main-file-name Simplex.cpp -analyzer-store=region -analyzer-opt-analyze-nested-blocks -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=cplusplus -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -analyzer-config-compatibility-mode=true -mrelocation-model pic -pic-level 2 -mframe-pointer=none -fmath-errno -ffp-contract=on -fno-rounding-math -mconstructor-aliases -funwind-tables=2 -target-cpu x86-64 -tune-cpu generic -debugger-tuning=gdb -ffunction-sections -fdata-sections -fcoverage-compilation-dir=/build/llvm-toolchain-snapshot-14~++20220125101009+ceec4383681c/build-llvm/tools/clang/stage2-bins -resource-dir /usr/lib/llvm-14/lib/clang/14.0.0 -D MLIR_CUDA_CONVERSIONS_ENABLED=1 -D MLIR_ROCM_CONVERSIONS_ENABLED=1 -D _DEBUG -D _GNU_SOURCE -D __STDC_CONSTANT_MACROS -D __STDC_FORMAT_MACROS -D __STDC_LIMIT_MACROS -I tools/mlir/lib/Analysis/Presburger -I /build/llvm-toolchain-snapshot-14~++20220125101009+ceec4383681c/mlir/lib/Analysis/Presburger -I include -I /build/llvm-toolchain-snapshot-14~++20220125101009+ceec4383681c/llvm/include -I /build/llvm-toolchain-snapshot-14~++20220125101009+ceec4383681c/mlir/include -I tools/mlir/include -D _FORTIFY_SOURCE=2 -D NDEBUG -U NDEBUG -internal-isystem /usr/lib/gcc/x86_64-linux-gnu/10/../../../../include/c++/10 -internal-isystem /usr/lib/gcc/x86_64-linux-gnu/10/../../../../include/x86_64-linux-gnu/c++/10 -internal-isystem /usr/lib/gcc/x86_64-linux-gnu/10/../../../../include/c++/10/backward -internal-isystem /usr/lib/llvm-14/lib/clang/14.0.0/include -internal-isystem /usr/local/include -internal-isystem /usr/lib/gcc/x86_64-linux-gnu/10/../../../../x86_64-linux-gnu/include -internal-externc-isystem /usr/include/x86_64-linux-gnu -internal-externc-isystem /include -internal-externc-isystem /usr/include -fmacro-prefix-map=/build/llvm-toolchain-snapshot-14~++20220125101009+ceec4383681c/build-llvm/tools/clang/stage2-bins=build-llvm/tools/clang/stage2-bins -fmacro-prefix-map=/build/llvm-toolchain-snapshot-14~++20220125101009+ceec4383681c/= -fcoverage-prefix-map=/build/llvm-toolchain-snapshot-14~++20220125101009+ceec4383681c/build-llvm/tools/clang/stage2-bins=build-llvm/tools/clang/stage2-bins -fcoverage-prefix-map=/build/llvm-toolchain-snapshot-14~++20220125101009+ceec4383681c/= -O3 -Wno-unused-command-line-argument -Wno-unused-parameter -Wwrite-strings -Wno-missing-field-initializers -Wno-long-long -Wno-maybe-uninitialized -Wno-class-memaccess -Wno-redundant-move -Wno-pessimizing-move -Wno-noexcept-type -Wno-comment -std=c++14 -fdeprecated-macro -fdebug-compilation-dir=/build/llvm-toolchain-snapshot-14~++20220125101009+ceec4383681c/build-llvm/tools/clang/stage2-bins -fdebug-prefix-map=/build/llvm-toolchain-snapshot-14~++20220125101009+ceec4383681c/build-llvm/tools/clang/stage2-bins=build-llvm/tools/clang/stage2-bins -fdebug-prefix-map=/build/llvm-toolchain-snapshot-14~++20220125101009+ceec4383681c/= -ferror-limit 19 -fvisibility-inlines-hidden -stack-protector 2 -fgnuc-version=4.2.1 -fcolor-diagnostics -vectorize-loops -vectorize-slp -analyzer-output=html -analyzer-config stable-report-filename=true -faddrsig -D__GCC_HAVE_DWARF2_CFI_ASM=1 -o /tmp/scan-build-2022-01-25-232935-20746-1 -x c++ /build/llvm-toolchain-snapshot-14~++20220125101009+ceec4383681c/mlir/lib/Analysis/Presburger/Simplex.cpp
1//===- Simplex.cpp - MLIR Simplex Class -----------------------------------===//
2//
3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4// See https://llvm.org/LICENSE.txt for license information.
5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6//
7//===----------------------------------------------------------------------===//
8
9#include "mlir/Analysis/Presburger/Simplex.h"
10#include "mlir/Analysis/Presburger/Matrix.h"
11#include "mlir/Support/MathExtras.h"
12#include "llvm/ADT/Optional.h"
13
14namespace mlir {
15using Direction = Simplex::Direction;
16
17const int nullIndex = std::numeric_limits<int>::max();
18
19/// Construct a Simplex object with `nVar` variables.
20SimplexBase::SimplexBase(unsigned nVar)
21 : nRow(0), nCol(2), nRedundant(0), tableau(0, 2 + nVar), empty(false) {
22 colUnknown.push_back(nullIndex);
23 colUnknown.push_back(nullIndex);
24 for (unsigned i = 0; i < nVar; ++i) {
25 var.emplace_back(Orientation::Column, /*restricted=*/false, /*pos=*/nCol);
26 colUnknown.push_back(i);
27 nCol++;
28 }
29}
30
31SimplexBase::SimplexBase(const IntegerPolyhedron &constraints)
32 : SimplexBase(constraints.getNumIds()) {
33 for (unsigned i = 0, numIneqs = constraints.getNumInequalities();
34 i < numIneqs; ++i)
35 addInequality(constraints.getInequality(i));
36 for (unsigned i = 0, numEqs = constraints.getNumEqualities(); i < numEqs; ++i)
37 addEquality(constraints.getEquality(i));
38}
39
40const Simplex::Unknown &SimplexBase::unknownFromIndex(int index) const {
41 assert(index != nullIndex && "nullIndex passed to unknownFromIndex")(static_cast <bool> (index != nullIndex && "nullIndex passed to unknownFromIndex"
) ? void (0) : __assert_fail ("index != nullIndex && \"nullIndex passed to unknownFromIndex\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 41, __extension__
__PRETTY_FUNCTION__))
;
42 return index >= 0 ? var[index] : con[~index];
43}
44
45const Simplex::Unknown &SimplexBase::unknownFromColumn(unsigned col) const {
46 assert(col < nCol && "Invalid column")(static_cast <bool> (col < nCol && "Invalid column"
) ? void (0) : __assert_fail ("col < nCol && \"Invalid column\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 46, __extension__
__PRETTY_FUNCTION__))
;
47 return unknownFromIndex(colUnknown[col]);
48}
49
50const Simplex::Unknown &SimplexBase::unknownFromRow(unsigned row) const {
51 assert(row < nRow && "Invalid row")(static_cast <bool> (row < nRow && "Invalid row"
) ? void (0) : __assert_fail ("row < nRow && \"Invalid row\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 51, __extension__
__PRETTY_FUNCTION__))
;
52 return unknownFromIndex(rowUnknown[row]);
53}
54
55Simplex::Unknown &SimplexBase::unknownFromIndex(int index) {
56 assert(index != nullIndex && "nullIndex passed to unknownFromIndex")(static_cast <bool> (index != nullIndex && "nullIndex passed to unknownFromIndex"
) ? void (0) : __assert_fail ("index != nullIndex && \"nullIndex passed to unknownFromIndex\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 56, __extension__
__PRETTY_FUNCTION__))
;
57 return index >= 0 ? var[index] : con[~index];
58}
59
60Simplex::Unknown &SimplexBase::unknownFromColumn(unsigned col) {
61 assert(col < nCol && "Invalid column")(static_cast <bool> (col < nCol && "Invalid column"
) ? void (0) : __assert_fail ("col < nCol && \"Invalid column\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 61, __extension__
__PRETTY_FUNCTION__))
;
62 return unknownFromIndex(colUnknown[col]);
63}
64
65Simplex::Unknown &SimplexBase::unknownFromRow(unsigned row) {
66 assert(row < nRow && "Invalid row")(static_cast <bool> (row < nRow && "Invalid row"
) ? void (0) : __assert_fail ("row < nRow && \"Invalid row\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 66, __extension__
__PRETTY_FUNCTION__))
;
67 return unknownFromIndex(rowUnknown[row]);
68}
69
70/// Add a new row to the tableau corresponding to the given constant term and
71/// list of coefficients. The coefficients are specified as a vector of
72/// (variable index, coefficient) pairs.
73unsigned SimplexBase::addRow(ArrayRef<int64_t> coeffs) {
74 assert(coeffs.size() == 1 + var.size() &&(static_cast <bool> (coeffs.size() == 1 + var.size() &&
"Incorrect number of coefficients!") ? void (0) : __assert_fail
("coeffs.size() == 1 + var.size() && \"Incorrect number of coefficients!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 75, __extension__
__PRETTY_FUNCTION__))
75 "Incorrect number of coefficients!")(static_cast <bool> (coeffs.size() == 1 + var.size() &&
"Incorrect number of coefficients!") ? void (0) : __assert_fail
("coeffs.size() == 1 + var.size() && \"Incorrect number of coefficients!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 75, __extension__
__PRETTY_FUNCTION__))
;
76
77 ++nRow;
78 // If the tableau is not big enough to accomodate the extra row, we extend it.
79 if (nRow >= tableau.getNumRows())
80 tableau.resizeVertically(nRow);
81 rowUnknown.push_back(~con.size());
82 con.emplace_back(Orientation::Row, false, nRow - 1);
83
84 tableau(nRow - 1, 0) = 1;
85 tableau(nRow - 1, 1) = coeffs.back();
86 for (unsigned col = 2; col < nCol; ++col)
87 tableau(nRow - 1, col) = 0;
88
89 // Process each given variable coefficient.
90 for (unsigned i = 0; i < var.size(); ++i) {
91 unsigned pos = var[i].pos;
92 if (coeffs[i] == 0)
93 continue;
94
95 if (var[i].orientation == Orientation::Column) {
96 // If a variable is in column position at column col, then we just add the
97 // coefficient for that variable (scaled by the common row denominator) to
98 // the corresponding entry in the new row.
99 tableau(nRow - 1, pos) += coeffs[i] * tableau(nRow - 1, 0);
100 continue;
101 }
102
103 // If the variable is in row position, we need to add that row to the new
104 // row, scaled by the coefficient for the variable, accounting for the two
105 // rows potentially having different denominators. The new denominator is
106 // the lcm of the two.
107 int64_t lcm = mlir::lcm(tableau(nRow - 1, 0), tableau(pos, 0));
108 int64_t nRowCoeff = lcm / tableau(nRow - 1, 0);
109 int64_t idxRowCoeff = coeffs[i] * (lcm / tableau(pos, 0));
110 tableau(nRow - 1, 0) = lcm;
111 for (unsigned col = 1; col < nCol; ++col)
112 tableau(nRow - 1, col) =
113 nRowCoeff * tableau(nRow - 1, col) + idxRowCoeff * tableau(pos, col);
114 }
115
116 normalizeRow(nRow - 1);
117 // Push to undo log along with the index of the new constraint.
118 undoLog.push_back(UndoLogEntry::RemoveLastConstraint);
119 return con.size() - 1;
120}
121
122/// Normalize the row by removing factors that are common between the
123/// denominator and all the numerator coefficients.
124void SimplexBase::normalizeRow(unsigned row) {
125 int64_t gcd = 0;
126 for (unsigned col = 0; col < nCol; ++col) {
127 gcd = llvm::greatestCommonDivisor(gcd, std::abs(tableau(row, col)));
128 // If the gcd becomes 1 then the row is already normalized.
129 if (gcd == 1)
130 return;
131 }
132
133 // Note that the gcd can never become zero since the first element of the row,
134 // the denominator, is non-zero.
135 assert(gcd != 0)(static_cast <bool> (gcd != 0) ? void (0) : __assert_fail
("gcd != 0", "mlir/lib/Analysis/Presburger/Simplex.cpp", 135
, __extension__ __PRETTY_FUNCTION__))
;
136 for (unsigned col = 0; col < nCol; ++col)
137 tableau(row, col) /= gcd;
138}
139
140namespace {
141bool signMatchesDirection(int64_t elem, Direction direction) {
142 assert(elem != 0 && "elem should not be 0")(static_cast <bool> (elem != 0 && "elem should not be 0"
) ? void (0) : __assert_fail ("elem != 0 && \"elem should not be 0\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 142, __extension__
__PRETTY_FUNCTION__))
;
143 return direction == Direction::Up ? elem > 0 : elem < 0;
144}
145
146Direction flippedDirection(Direction direction) {
147 return direction == Direction::Up ? Direction::Down : Simplex::Direction::Up;
148}
149} // namespace
150
151/// Find a pivot to change the sample value of the row in the specified
152/// direction. The returned pivot row will involve `row` if and only if the
153/// unknown is unbounded in the specified direction.
154///
155/// To increase (resp. decrease) the value of a row, we need to find a live
156/// column with a non-zero coefficient. If the coefficient is positive, we need
157/// to increase (decrease) the value of the column, and if the coefficient is
158/// negative, we need to decrease (increase) the value of the column. Also,
159/// we cannot decrease the sample value of restricted columns.
160///
161/// If multiple columns are valid, we break ties by considering a lexicographic
162/// ordering where we prefer unknowns with lower index.
163Optional<SimplexBase::Pivot> SimplexBase::findPivot(int row,
164 Direction direction) const {
165 Optional<unsigned> col;
166 for (unsigned j = 2; j < nCol; ++j) {
167 int64_t elem = tableau(row, j);
168 if (elem == 0)
169 continue;
170
171 if (unknownFromColumn(j).restricted &&
172 !signMatchesDirection(elem, direction))
173 continue;
174 if (!col || colUnknown[j] < colUnknown[*col])
175 col = j;
176 }
177
178 if (!col)
179 return {};
180
181 Direction newDirection =
182 tableau(row, *col) < 0 ? flippedDirection(direction) : direction;
183 Optional<unsigned> maybePivotRow = findPivotRow(row, newDirection, *col);
184 return Pivot{maybePivotRow.getValueOr(row), *col};
185}
186
187/// Swap the associated unknowns for the row and the column.
188///
189/// First we swap the index associated with the row and column. Then we update
190/// the unknowns to reflect their new position and orientation.
191void SimplexBase::swapRowWithCol(unsigned row, unsigned col) {
192 std::swap(rowUnknown[row], colUnknown[col]);
193 Unknown &uCol = unknownFromColumn(col);
194 Unknown &uRow = unknownFromRow(row);
195 uCol.orientation = Orientation::Column;
196 uRow.orientation = Orientation::Row;
197 uCol.pos = col;
198 uRow.pos = row;
199}
200
201void SimplexBase::pivot(Pivot pair) { pivot(pair.row, pair.column); }
202
203/// Pivot pivotRow and pivotCol.
204///
205/// Let R be the pivot row unknown and let C be the pivot col unknown.
206/// Since initially R = a*C + sum b_i * X_i
207/// (where the sum is over the other column's unknowns, x_i)
208/// C = (R - (sum b_i * X_i))/a
209///
210/// Let u be some other row unknown.
211/// u = c*C + sum d_i * X_i
212/// So u = c*(R - sum b_i * X_i)/a + sum d_i * X_i
213///
214/// This results in the following transform:
215/// pivot col other col pivot col other col
216/// pivot row a b -> pivot row 1/a -b/a
217/// other row c d other row c/a d - bc/a
218///
219/// Taking into account the common denominators p and q:
220///
221/// pivot col other col pivot col other col
222/// pivot row a/p b/p -> pivot row p/a -b/a
223/// other row c/q d/q other row cp/aq (da - bc)/aq
224///
225/// The pivot row transform is accomplished be swapping a with the pivot row's
226/// common denominator and negating the pivot row except for the pivot column
227/// element.
228void SimplexBase::pivot(unsigned pivotRow, unsigned pivotCol) {
229 assert(pivotCol >= 2 && "Refusing to pivot invalid column")(static_cast <bool> (pivotCol >= 2 && "Refusing to pivot invalid column"
) ? void (0) : __assert_fail ("pivotCol >= 2 && \"Refusing to pivot invalid column\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 229, __extension__
__PRETTY_FUNCTION__))
;
230
231 swapRowWithCol(pivotRow, pivotCol);
232 std::swap(tableau(pivotRow, 0), tableau(pivotRow, pivotCol));
233 // We need to negate the whole pivot row except for the pivot column.
234 if (tableau(pivotRow, 0) < 0) {
235 // If the denominator is negative, we negate the row by simply negating the
236 // denominator.
237 tableau(pivotRow, 0) = -tableau(pivotRow, 0);
238 tableau(pivotRow, pivotCol) = -tableau(pivotRow, pivotCol);
239 } else {
240 for (unsigned col = 1; col < nCol; ++col) {
241 if (col == pivotCol)
242 continue;
243 tableau(pivotRow, col) = -tableau(pivotRow, col);
244 }
245 }
246 normalizeRow(pivotRow);
247
248 for (unsigned row = 0; row < nRow; ++row) {
249 if (row == pivotRow)
250 continue;
251 if (tableau(row, pivotCol) == 0) // Nothing to do.
252 continue;
253 tableau(row, 0) *= tableau(pivotRow, 0);
254 for (unsigned j = 1; j < nCol; ++j) {
255 if (j == pivotCol)
256 continue;
257 // Add rather than subtract because the pivot row has been negated.
258 tableau(row, j) = tableau(row, j) * tableau(pivotRow, 0) +
259 tableau(row, pivotCol) * tableau(pivotRow, j);
260 }
261 tableau(row, pivotCol) *= tableau(pivotRow, pivotCol);
262 normalizeRow(row);
263 }
264}
265
266/// Perform pivots until the unknown has a non-negative sample value or until
267/// no more upward pivots can be performed. Return success if we were able to
268/// bring the row to a non-negative sample value, and failure otherwise.
269LogicalResult SimplexBase::restoreRow(Unknown &u) {
270 assert(u.orientation == Orientation::Row &&(static_cast <bool> (u.orientation == Orientation::Row &&
"unknown should be in row position") ? void (0) : __assert_fail
("u.orientation == Orientation::Row && \"unknown should be in row position\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 271, __extension__
__PRETTY_FUNCTION__))
271 "unknown should be in row position")(static_cast <bool> (u.orientation == Orientation::Row &&
"unknown should be in row position") ? void (0) : __assert_fail
("u.orientation == Orientation::Row && \"unknown should be in row position\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 271, __extension__
__PRETTY_FUNCTION__))
;
272
273 while (tableau(u.pos, 1) < 0) {
274 Optional<Pivot> maybePivot = findPivot(u.pos, Direction::Up);
275 if (!maybePivot)
276 break;
277
278 pivot(*maybePivot);
279 if (u.orientation == Orientation::Column)
280 return success(); // the unknown is unbounded above.
281 }
282 return success(tableau(u.pos, 1) >= 0);
283}
284
285/// Find a row that can be used to pivot the column in the specified direction.
286/// This returns an empty optional if and only if the column is unbounded in the
287/// specified direction (ignoring skipRow, if skipRow is set).
288///
289/// If skipRow is set, this row is not considered, and (if it is restricted) its
290/// restriction may be violated by the returned pivot. Usually, skipRow is set
291/// because we don't want to move it to column position unless it is unbounded,
292/// and we are either trying to increase the value of skipRow or explicitly
293/// trying to make skipRow negative, so we are not concerned about this.
294///
295/// If the direction is up (resp. down) and a restricted row has a negative
296/// (positive) coefficient for the column, then this row imposes a bound on how
297/// much the sample value of the column can change. Such a row with constant
298/// term c and coefficient f for the column imposes a bound of c/|f| on the
299/// change in sample value (in the specified direction). (note that c is
300/// non-negative here since the row is restricted and the tableau is consistent)
301///
302/// We iterate through the rows and pick the row which imposes the most
303/// stringent bound, since pivoting with a row changes the row's sample value to
304/// 0 and hence saturates the bound it imposes. We break ties between rows that
305/// impose the same bound by considering a lexicographic ordering where we
306/// prefer unknowns with lower index value.
307Optional<unsigned> SimplexBase::findPivotRow(Optional<unsigned> skipRow,
308 Direction direction,
309 unsigned col) const {
310 Optional<unsigned> retRow;
311 // Initialize these to zero in order to silence a warning about retElem and
312 // retConst being used uninitialized in the initialization of `diff` below. In
313 // reality, these are always initialized when that line is reached since these
314 // are set whenever retRow is set.
315 int64_t retElem = 0, retConst = 0;
316 for (unsigned row = nRedundant; row < nRow; ++row) {
317 if (skipRow && row == *skipRow)
318 continue;
319 int64_t elem = tableau(row, col);
320 if (elem == 0)
321 continue;
322 if (!unknownFromRow(row).restricted)
323 continue;
324 if (signMatchesDirection(elem, direction))
325 continue;
326 int64_t constTerm = tableau(row, 1);
327
328 if (!retRow) {
329 retRow = row;
330 retElem = elem;
331 retConst = constTerm;
332 continue;
333 }
334
335 int64_t diff = retConst * elem - constTerm * retElem;
336 if ((diff == 0 && rowUnknown[row] < rowUnknown[*retRow]) ||
337 (diff != 0 && !signMatchesDirection(diff, direction))) {
338 retRow = row;
339 retElem = elem;
340 retConst = constTerm;
341 }
342 }
343 return retRow;
344}
345
346bool SimplexBase::isEmpty() const { return empty; }
347
348void SimplexBase::swapRows(unsigned i, unsigned j) {
349 if (i == j)
350 return;
351 tableau.swapRows(i, j);
352 std::swap(rowUnknown[i], rowUnknown[j]);
353 unknownFromRow(i).pos = i;
354 unknownFromRow(j).pos = j;
355}
356
357void SimplexBase::swapColumns(unsigned i, unsigned j) {
358 assert(i < nCol && j < nCol && "Invalid columns provided!")(static_cast <bool> (i < nCol && j < nCol
&& "Invalid columns provided!") ? void (0) : __assert_fail
("i < nCol && j < nCol && \"Invalid columns provided!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 358, __extension__
__PRETTY_FUNCTION__))
;
359 if (i == j)
360 return;
361 tableau.swapColumns(i, j);
362 std::swap(colUnknown[i], colUnknown[j]);
363 unknownFromColumn(i).pos = i;
364 unknownFromColumn(j).pos = j;
365}
366
367/// Mark this tableau empty and push an entry to the undo stack.
368void SimplexBase::markEmpty() {
369 // If the set is already empty, then we shouldn't add another UnmarkEmpty log
370 // entry, since in that case the Simplex will be erroneously marked as
371 // non-empty when rolling back past this point.
372 if (empty)
373 return;
374 undoLog.push_back(UndoLogEntry::UnmarkEmpty);
375 empty = true;
376}
377
378/// Add an inequality to the tableau. If coeffs is c_0, c_1, ... c_n, where n
379/// is the current number of variables, then the corresponding inequality is
380/// c_n + c_0*x_0 + c_1*x_1 + ... + c_{n-1}*x_{n-1} >= 0.
381///
382/// We add the inequality and mark it as restricted. We then try to make its
383/// sample value non-negative. If this is not possible, the tableau has become
384/// empty and we mark it as such.
385void SimplexBase::addInequality(ArrayRef<int64_t> coeffs) {
386 unsigned conIndex = addRow(coeffs);
387 Unknown &u = con[conIndex];
388 u.restricted = true;
389 LogicalResult result = restoreRow(u);
390 if (failed(result))
391 markEmpty();
392}
393
394/// Add an equality to the tableau. If coeffs is c_0, c_1, ... c_n, where n
395/// is the current number of variables, then the corresponding equality is
396/// c_n + c_0*x_0 + c_1*x_1 + ... + c_{n-1}*x_{n-1} == 0.
397///
398/// We simply add two opposing inequalities, which force the expression to
399/// be zero.
400void SimplexBase::addEquality(ArrayRef<int64_t> coeffs) {
401 addInequality(coeffs);
402 SmallVector<int64_t, 8> negatedCoeffs;
403 for (int64_t coeff : coeffs)
404 negatedCoeffs.emplace_back(-coeff);
405 addInequality(negatedCoeffs);
406}
407
408unsigned SimplexBase::getNumVariables() const { return var.size(); }
409unsigned SimplexBase::getNumConstraints() const { return con.size(); }
410
411/// Return a snapshot of the current state. This is just the current size of the
412/// undo log.
413unsigned SimplexBase::getSnapshot() const { return undoLog.size(); }
414
415void SimplexBase::undo(UndoLogEntry entry) {
416 if (entry == UndoLogEntry::RemoveLastConstraint) {
417 Unknown &constraint = con.back();
418 if (constraint.orientation == Orientation::Column) {
419 unsigned column = constraint.pos;
420 Optional<unsigned> row;
421
422 // Try to find any pivot row for this column that preserves tableau
423 // consistency (except possibly the column itself, which is going to be
424 // deallocated anyway).
425 //
426 // If no pivot row is found in either direction, then the unknown is
427 // unbounded in both directions and we are free to
428 // perform any pivot at all. To do this, we just need to find any row with
429 // a non-zero coefficient for the column.
430 if (Optional<unsigned> maybeRow =
431 findPivotRow({}, Direction::Up, column)) {
432 row = *maybeRow;
433 } else if (Optional<unsigned> maybeRow =
434 findPivotRow({}, Direction::Down, column)) {
435 row = *maybeRow;
436 } else {
437 // The loop doesn't find a pivot row only if the column has zero
438 // coefficients for every row. But the unknown is a constraint,
439 // so it was added initially as a row. Such a row could never have been
440 // pivoted to a column. So a pivot row will always be found.
441 for (unsigned i = nRedundant; i < nRow; ++i) {
442 if (tableau(i, column) != 0) {
443 row = i;
444 break;
445 }
446 }
447 }
448 assert(row.hasValue() && "No pivot row found!")(static_cast <bool> (row.hasValue() && "No pivot row found!"
) ? void (0) : __assert_fail ("row.hasValue() && \"No pivot row found!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 448, __extension__
__PRETTY_FUNCTION__))
;
449 pivot(*row, column);
450 }
451
452 // Move this unknown to the last row and remove the last row from the
453 // tableau.
454 swapRows(constraint.pos, nRow - 1);
455 // It is not strictly necessary to shrink the tableau, but for now we
456 // maintain the invariant that the tableau has exactly nRow rows.
457 tableau.resizeVertically(nRow - 1);
458 nRow--;
459 rowUnknown.pop_back();
460 con.pop_back();
461 } else if (entry == UndoLogEntry::RemoveLastVariable) {
462 // Whenever we are rolling back the addition of a variable, it is guaranteed
463 // that the variable will be in column position.
464 //
465 // We can see this as follows: any constraint that depends on this variable
466 // was added after this variable was added, so the addition of such
467 // constraints should already have been rolled back by the time we get to
468 // rolling back the addition of the variable. Therefore, no constraint
469 // currently has a component along the variable, so the variable itself must
470 // be part of the basis.
471 assert(var.back().orientation == Orientation::Column &&(static_cast <bool> (var.back().orientation == Orientation
::Column && "Variable to be removed must be in column orientation!"
) ? void (0) : __assert_fail ("var.back().orientation == Orientation::Column && \"Variable to be removed must be in column orientation!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 472, __extension__
__PRETTY_FUNCTION__))
472 "Variable to be removed must be in column orientation!")(static_cast <bool> (var.back().orientation == Orientation
::Column && "Variable to be removed must be in column orientation!"
) ? void (0) : __assert_fail ("var.back().orientation == Orientation::Column && \"Variable to be removed must be in column orientation!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 472, __extension__
__PRETTY_FUNCTION__))
;
473
474 // Move this variable to the last column and remove the column from the
475 // tableau.
476 swapColumns(var.back().pos, nCol - 1);
477 tableau.resizeHorizontally(nCol - 1);
478 var.pop_back();
479 colUnknown.pop_back();
480 nCol--;
481 } else if (entry == UndoLogEntry::UnmarkEmpty) {
482 empty = false;
483 } else if (entry == UndoLogEntry::UnmarkLastRedundant) {
484 nRedundant--;
485 }
486}
487
488/// Rollback to the specified snapshot.
489///
490/// We undo all the log entries until the log size when the snapshot was taken
491/// is reached.
492void SimplexBase::rollback(unsigned snapshot) {
493 while (undoLog.size() > snapshot) {
494 undo(undoLog.back());
495 undoLog.pop_back();
496 }
497}
498
499void SimplexBase::appendVariable(unsigned count) {
500 if (count == 0)
501 return;
502 var.reserve(var.size() + count);
503 colUnknown.reserve(colUnknown.size() + count);
504 for (unsigned i = 0; i < count; ++i) {
505 nCol++;
506 var.emplace_back(Orientation::Column, /*restricted=*/false,
507 /*pos=*/nCol - 1);
508 colUnknown.push_back(var.size() - 1);
509 }
510 tableau.resizeHorizontally(nCol);
511 undoLog.insert(undoLog.end(), count, UndoLogEntry::RemoveLastVariable);
512}
513
514/// Add all the constraints from the given IntegerPolyhedron.
515void SimplexBase::intersectIntegerPolyhedron(const IntegerPolyhedron &poly) {
516 assert(poly.getNumIds() == getNumVariables() &&(static_cast <bool> (poly.getNumIds() == getNumVariables
() && "IntegerPolyhedron must have same dimensionality as simplex"
) ? void (0) : __assert_fail ("poly.getNumIds() == getNumVariables() && \"IntegerPolyhedron must have same dimensionality as simplex\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 517, __extension__
__PRETTY_FUNCTION__))
517 "IntegerPolyhedron must have same dimensionality as simplex")(static_cast <bool> (poly.getNumIds() == getNumVariables
() && "IntegerPolyhedron must have same dimensionality as simplex"
) ? void (0) : __assert_fail ("poly.getNumIds() == getNumVariables() && \"IntegerPolyhedron must have same dimensionality as simplex\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 517, __extension__
__PRETTY_FUNCTION__))
;
518 for (unsigned i = 0, e = poly.getNumInequalities(); i < e; ++i)
519 addInequality(poly.getInequality(i));
520 for (unsigned i = 0, e = poly.getNumEqualities(); i < e; ++i)
521 addEquality(poly.getEquality(i));
522}
523
524Optional<Fraction> Simplex::computeRowOptimum(Direction direction,
525 unsigned row) {
526 // Keep trying to find a pivot for the row in the specified direction.
527 while (Optional<Pivot> maybePivot = findPivot(row, direction)) {
528 // If findPivot returns a pivot involving the row itself, then the optimum
529 // is unbounded, so we return None.
530 if (maybePivot->row == row)
531 return {};
532 pivot(*maybePivot);
533 }
534
535 // The row has reached its optimal sample value, which we return.
536 // The sample value is the entry in the constant column divided by the common
537 // denominator for this row.
538 return Fraction(tableau(row, 1), tableau(row, 0));
539}
540
541/// Compute the optimum of the specified expression in the specified direction,
542/// or None if it is unbounded.
543Optional<Fraction> Simplex::computeOptimum(Direction direction,
544 ArrayRef<int64_t> coeffs) {
545 assert(!empty && "Simplex should not be empty")(static_cast <bool> (!empty && "Simplex should not be empty"
) ? void (0) : __assert_fail ("!empty && \"Simplex should not be empty\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 545, __extension__
__PRETTY_FUNCTION__))
;
546
547 unsigned snapshot = getSnapshot();
548 unsigned conIndex = addRow(coeffs);
549 unsigned row = con[conIndex].pos;
550 Optional<Fraction> optimum = computeRowOptimum(direction, row);
551 rollback(snapshot);
552 return optimum;
553}
554
555Optional<Fraction> Simplex::computeOptimum(Direction direction, Unknown &u) {
556 assert(!empty && "Simplex should not be empty!")(static_cast <bool> (!empty && "Simplex should not be empty!"
) ? void (0) : __assert_fail ("!empty && \"Simplex should not be empty!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 556, __extension__
__PRETTY_FUNCTION__))
;
557 if (u.orientation == Orientation::Column) {
558 unsigned column = u.pos;
559 Optional<unsigned> pivotRow = findPivotRow({}, direction, column);
560 // If no pivot is returned, the constraint is unbounded in the specified
561 // direction.
562 if (!pivotRow)
563 return {};
564 pivot(*pivotRow, column);
565 }
566
567 unsigned row = u.pos;
568 Optional<Fraction> optimum = computeRowOptimum(direction, row);
569 if (u.restricted && direction == Direction::Down &&
570 (!optimum || *optimum < Fraction(0, 1))) {
571 if (failed(restoreRow(u)))
572 llvm_unreachable("Could not restore row!")::llvm::llvm_unreachable_internal("Could not restore row!", "mlir/lib/Analysis/Presburger/Simplex.cpp"
, 572)
;
573 }
574 return optimum;
575}
576
577bool Simplex::isBoundedAlongConstraint(unsigned constraintIndex) {
578 assert(!empty && "It is not meaningful to ask whether a direction is bounded "(static_cast <bool> (!empty && "It is not meaningful to ask whether a direction is bounded "
"in an empty set.") ? void (0) : __assert_fail ("!empty && \"It is not meaningful to ask whether a direction is bounded \" \"in an empty set.\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 579, __extension__
__PRETTY_FUNCTION__))
579 "in an empty set.")(static_cast <bool> (!empty && "It is not meaningful to ask whether a direction is bounded "
"in an empty set.") ? void (0) : __assert_fail ("!empty && \"It is not meaningful to ask whether a direction is bounded \" \"in an empty set.\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 579, __extension__
__PRETTY_FUNCTION__))
;
580 // The constraint's perpendicular is already bounded below, since it is a
581 // constraint. If it is also bounded above, we can return true.
582 return computeOptimum(Direction::Up, con[constraintIndex]).hasValue();
583}
584
585/// Redundant constraints are those that are in row orientation and lie in
586/// rows 0 to nRedundant - 1.
587bool Simplex::isMarkedRedundant(unsigned constraintIndex) const {
588 const Unknown &u = con[constraintIndex];
589 return u.orientation == Orientation::Row && u.pos < nRedundant;
590}
591
592/// Mark the specified row redundant.
593///
594/// This is done by moving the unknown to the end of the block of redundant
595/// rows (namely, to row nRedundant) and incrementing nRedundant to
596/// accomodate the new redundant row.
597void Simplex::markRowRedundant(Unknown &u) {
598 assert(u.orientation == Orientation::Row &&(static_cast <bool> (u.orientation == Orientation::Row &&
"Unknown should be in row position!") ? void (0) : __assert_fail
("u.orientation == Orientation::Row && \"Unknown should be in row position!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 599, __extension__
__PRETTY_FUNCTION__))
599 "Unknown should be in row position!")(static_cast <bool> (u.orientation == Orientation::Row &&
"Unknown should be in row position!") ? void (0) : __assert_fail
("u.orientation == Orientation::Row && \"Unknown should be in row position!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 599, __extension__
__PRETTY_FUNCTION__))
;
600 assert(u.pos >= nRedundant && "Unknown is already marked redundant!")(static_cast <bool> (u.pos >= nRedundant && "Unknown is already marked redundant!"
) ? void (0) : __assert_fail ("u.pos >= nRedundant && \"Unknown is already marked redundant!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 600, __extension__
__PRETTY_FUNCTION__))
;
601 swapRows(u.pos, nRedundant);
602 ++nRedundant;
603 undoLog.emplace_back(UndoLogEntry::UnmarkLastRedundant);
604}
605
606/// Find a subset of constraints that is redundant and mark them redundant.
607void Simplex::detectRedundant() {
608 // It is not meaningful to talk about redundancy for empty sets.
609 if (empty)
610 return;
611
612 // Iterate through the constraints and check for each one if it can attain
613 // negative sample values. If it can, it's not redundant. Otherwise, it is.
614 // We mark redundant constraints redundant.
615 //
616 // Constraints that get marked redundant in one iteration are not respected
617 // when checking constraints in later iterations. This prevents, for example,
618 // two identical constraints both being marked redundant since each is
619 // redundant given the other one. In this example, only the first of the
620 // constraints that is processed will get marked redundant, as it should be.
621 for (Unknown &u : con) {
622 if (u.orientation == Orientation::Column) {
623 unsigned column = u.pos;
624 Optional<unsigned> pivotRow = findPivotRow({}, Direction::Down, column);
625 // If no downward pivot is returned, the constraint is unbounded below
626 // and hence not redundant.
627 if (!pivotRow)
628 continue;
629 pivot(*pivotRow, column);
630 }
631
632 unsigned row = u.pos;
633 Optional<Fraction> minimum = computeRowOptimum(Direction::Down, row);
634 if (!minimum || *minimum < Fraction(0, 1)) {
635 // Constraint is unbounded below or can attain negative sample values and
636 // hence is not redundant.
637 if (failed(restoreRow(u)))
638 llvm_unreachable("Could not restore non-redundant row!")::llvm::llvm_unreachable_internal("Could not restore non-redundant row!"
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 638)
;
639 continue;
640 }
641
642 markRowRedundant(u);
643 }
644}
645
646bool Simplex::isUnbounded() {
647 if (empty)
648 return false;
649
650 SmallVector<int64_t, 8> dir(var.size() + 1);
651 for (unsigned i = 0; i < var.size(); ++i) {
652 dir[i] = 1;
653
654 Optional<Fraction> maybeMax = computeOptimum(Direction::Up, dir);
655 if (!maybeMax)
656 return true;
657
658 Optional<Fraction> maybeMin = computeOptimum(Direction::Down, dir);
659 if (!maybeMin)
660 return true;
661
662 dir[i] = 0;
663 }
664 return false;
665}
666
667/// Make a tableau to represent a pair of points in the original tableau.
668///
669/// The product constraints and variables are stored as: first A's, then B's.
670///
671/// The product tableau has row layout:
672/// A's redundant rows, B's redundant rows, A's other rows, B's other rows.
673///
674/// It has column layout:
675/// denominator, constant, A's columns, B's columns.
676Simplex Simplex::makeProduct(const Simplex &a, const Simplex &b) {
677 unsigned numVar = a.getNumVariables() + b.getNumVariables();
678 unsigned numCon = a.getNumConstraints() + b.getNumConstraints();
679 Simplex result(numVar);
680
681 result.tableau.resizeVertically(numCon);
682 result.empty = a.empty || b.empty;
683
684 auto concat = [](ArrayRef<Unknown> v, ArrayRef<Unknown> w) {
685 SmallVector<Unknown, 8> result;
686 result.reserve(v.size() + w.size());
687 result.insert(result.end(), v.begin(), v.end());
688 result.insert(result.end(), w.begin(), w.end());
689 return result;
690 };
691 result.con = concat(a.con, b.con);
692 result.var = concat(a.var, b.var);
693
694 auto indexFromBIndex = [&](int index) {
695 return index >= 0 ? a.getNumVariables() + index
696 : ~(a.getNumConstraints() + ~index);
697 };
698
699 result.colUnknown.assign(2, nullIndex);
700 for (unsigned i = 2; i < a.nCol; ++i) {
701 result.colUnknown.push_back(a.colUnknown[i]);
702 result.unknownFromIndex(result.colUnknown.back()).pos =
703 result.colUnknown.size() - 1;
704 }
705 for (unsigned i = 2; i < b.nCol; ++i) {
706 result.colUnknown.push_back(indexFromBIndex(b.colUnknown[i]));
707 result.unknownFromIndex(result.colUnknown.back()).pos =
708 result.colUnknown.size() - 1;
709 }
710
711 auto appendRowFromA = [&](unsigned row) {
712 for (unsigned col = 0; col < a.nCol; ++col)
713 result.tableau(result.nRow, col) = a.tableau(row, col);
714 result.rowUnknown.push_back(a.rowUnknown[row]);
715 result.unknownFromIndex(result.rowUnknown.back()).pos =
716 result.rowUnknown.size() - 1;
717 result.nRow++;
718 };
719
720 // Also fixes the corresponding entry in rowUnknown and var/con (as the case
721 // may be).
722 auto appendRowFromB = [&](unsigned row) {
723 result.tableau(result.nRow, 0) = b.tableau(row, 0);
724 result.tableau(result.nRow, 1) = b.tableau(row, 1);
725
726 unsigned offset = a.nCol - 2;
727 for (unsigned col = 2; col < b.nCol; ++col)
728 result.tableau(result.nRow, offset + col) = b.tableau(row, col);
729 result.rowUnknown.push_back(indexFromBIndex(b.rowUnknown[row]));
730 result.unknownFromIndex(result.rowUnknown.back()).pos =
731 result.rowUnknown.size() - 1;
732 result.nRow++;
733 };
734
735 result.nRedundant = a.nRedundant + b.nRedundant;
736 for (unsigned row = 0; row < a.nRedundant; ++row)
737 appendRowFromA(row);
738 for (unsigned row = 0; row < b.nRedundant; ++row)
739 appendRowFromB(row);
740 for (unsigned row = a.nRedundant; row < a.nRow; ++row)
741 appendRowFromA(row);
742 for (unsigned row = b.nRedundant; row < b.nRow; ++row)
743 appendRowFromB(row);
744
745 return result;
746}
747
748Optional<SmallVector<Fraction, 8>> SimplexBase::getRationalSample() const {
749 if (empty)
750 return {};
751
752 SmallVector<Fraction, 8> sample;
753 sample.reserve(var.size());
754 // Push the sample value for each variable into the vector.
755 for (const Unknown &u : var) {
756 if (u.orientation == Orientation::Column) {
757 // If the variable is in column position, its sample value is zero.
758 sample.emplace_back(0, 1);
759 } else {
760 // If the variable is in row position, its sample value is the entry in
761 // the constant column divided by the entry in the common denominator
762 // column.
763 sample.emplace_back(tableau(u.pos, 1), tableau(u.pos, 0));
764 }
765 }
766 return sample;
767}
768
769Optional<SmallVector<int64_t, 8>>
770SimplexBase::getSamplePointIfIntegral() const {
771 // If the tableau is empty, no sample point exists.
772 if (empty)
773 return {};
774
775 // The value will always exist since the Simplex is non-empty.
776 SmallVector<Fraction, 8> rationalSample = *getRationalSample();
777 SmallVector<int64_t, 8> integerSample;
778 integerSample.reserve(var.size());
779 for (const Fraction &coord : rationalSample) {
780 // If the sample is non-integral, return None.
781 if (coord.num % coord.den != 0)
782 return {};
783 integerSample.push_back(coord.num / coord.den);
784 }
785 return integerSample;
786}
787
788/// Given a simplex for a polytope, construct a new simplex whose variables are
789/// identified with a pair of points (x, y) in the original polytope. Supports
790/// some operations needed for generalized basis reduction. In what follows,
791/// dotProduct(x, y) = x_1 * y_1 + x_2 * y_2 + ... x_n * y_n where n is the
792/// dimension of the original polytope.
793///
794/// This supports adding equality constraints dotProduct(dir, x - y) == 0. It
795/// also supports rolling back this addition, by maintaining a snapshot stack
796/// that contains a snapshot of the Simplex's state for each equality, just
797/// before that equality was added.
798class GBRSimplex {
799 using Orientation = Simplex::Orientation;
800
801public:
802 GBRSimplex(const Simplex &originalSimplex)
803 : simplex(Simplex::makeProduct(originalSimplex, originalSimplex)),
804 simplexConstraintOffset(simplex.getNumConstraints()) {}
805
806 /// Add an equality dotProduct(dir, x - y) == 0.
807 /// First pushes a snapshot for the current simplex state to the stack so
808 /// that this can be rolled back later.
809 void addEqualityForDirection(ArrayRef<int64_t> dir) {
810 assert((static_cast <bool> (std::any_of(dir.begin(), dir.end()
, [](int64_t x) { return x != 0; }) && "Direction passed is the zero vector!"
) ? void (0) : __assert_fail ("std::any_of(dir.begin(), dir.end(), [](int64_t x) { return x != 0; }) && \"Direction passed is the zero vector!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 812, __extension__
__PRETTY_FUNCTION__))
811 std::any_of(dir.begin(), dir.end(), [](int64_t x) { return x != 0; }) &&(static_cast <bool> (std::any_of(dir.begin(), dir.end()
, [](int64_t x) { return x != 0; }) && "Direction passed is the zero vector!"
) ? void (0) : __assert_fail ("std::any_of(dir.begin(), dir.end(), [](int64_t x) { return x != 0; }) && \"Direction passed is the zero vector!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 812, __extension__
__PRETTY_FUNCTION__))
812 "Direction passed is the zero vector!")(static_cast <bool> (std::any_of(dir.begin(), dir.end()
, [](int64_t x) { return x != 0; }) && "Direction passed is the zero vector!"
) ? void (0) : __assert_fail ("std::any_of(dir.begin(), dir.end(), [](int64_t x) { return x != 0; }) && \"Direction passed is the zero vector!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 812, __extension__
__PRETTY_FUNCTION__))
;
813 snapshotStack.push_back(simplex.getSnapshot());
814 simplex.addEquality(getCoeffsForDirection(dir));
815 }
816 /// Compute max(dotProduct(dir, x - y)).
817 Fraction computeWidth(ArrayRef<int64_t> dir) {
818 Optional<Fraction> maybeWidth =
819 simplex.computeOptimum(Direction::Up, getCoeffsForDirection(dir));
820 assert(maybeWidth.hasValue() && "Width should not be unbounded!")(static_cast <bool> (maybeWidth.hasValue() && "Width should not be unbounded!"
) ? void (0) : __assert_fail ("maybeWidth.hasValue() && \"Width should not be unbounded!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 820, __extension__
__PRETTY_FUNCTION__))
;
821 return *maybeWidth;
822 }
823
824 /// Compute max(dotProduct(dir, x - y)) and save the dual variables for only
825 /// the direction equalities to `dual`.
826 Fraction computeWidthAndDuals(ArrayRef<int64_t> dir,
827 SmallVectorImpl<int64_t> &dual,
828 int64_t &dualDenom) {
829 // We can't just call into computeWidth or computeOptimum since we need to
830 // access the state of the tableau after computing the optimum, and these
831 // functions rollback the insertion of the objective function into the
832 // tableau before returning. We instead add a row for the objective function
833 // ourselves, call into computeOptimum, compute the duals from the tableau
834 // state, and finally rollback the addition of the row before returning.
835 unsigned snap = simplex.getSnapshot();
836 unsigned conIndex = simplex.addRow(getCoeffsForDirection(dir));
837 unsigned row = simplex.con[conIndex].pos;
838 Optional<Fraction> maybeWidth =
839 simplex.computeRowOptimum(Simplex::Direction::Up, row);
840 assert(maybeWidth.hasValue() && "Width should not be unbounded!")(static_cast <bool> (maybeWidth.hasValue() && "Width should not be unbounded!"
) ? void (0) : __assert_fail ("maybeWidth.hasValue() && \"Width should not be unbounded!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 840, __extension__
__PRETTY_FUNCTION__))
;
841 dualDenom = simplex.tableau(row, 0);
842 dual.clear();
843
844 // The increment is i += 2 because equalities are added as two inequalities,
845 // one positive and one negative. Each iteration processes one equality.
846 for (unsigned i = simplexConstraintOffset; i < conIndex; i += 2) {
847 // The dual variable for an inequality in column orientation is the
848 // negative of its coefficient at the objective row. If the inequality is
849 // in row orientation, the corresponding dual variable is zero.
850 //
851 // We want the dual for the original equality, which corresponds to two
852 // inequalities: a positive inequality, which has the same coefficients as
853 // the equality, and a negative equality, which has negated coefficients.
854 //
855 // Note that at most one of these inequalities can be in column
856 // orientation because the column unknowns should form a basis and hence
857 // must be linearly independent. If the positive inequality is in column
858 // position, its dual is the dual corresponding to the equality. If the
859 // negative inequality is in column position, the negation of its dual is
860 // the dual corresponding to the equality. If neither is in column
861 // position, then that means that this equality is redundant, and its dual
862 // is zero.
863 //
864 // Note that it is NOT valid to perform pivots during the computation of
865 // the duals. This entire dual computation must be performed on the same
866 // tableau configuration.
867 assert(!(simplex.con[i].orientation == Orientation::Column &&(static_cast <bool> (!(simplex.con[i].orientation == Orientation
::Column && simplex.con[i + 1].orientation == Orientation
::Column) && "Both inequalities for the equality cannot be in column "
"orientation!") ? void (0) : __assert_fail ("!(simplex.con[i].orientation == Orientation::Column && simplex.con[i + 1].orientation == Orientation::Column) && \"Both inequalities for the equality cannot be in column \" \"orientation!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 870, __extension__
__PRETTY_FUNCTION__))
868 simplex.con[i + 1].orientation == Orientation::Column) &&(static_cast <bool> (!(simplex.con[i].orientation == Orientation
::Column && simplex.con[i + 1].orientation == Orientation
::Column) && "Both inequalities for the equality cannot be in column "
"orientation!") ? void (0) : __assert_fail ("!(simplex.con[i].orientation == Orientation::Column && simplex.con[i + 1].orientation == Orientation::Column) && \"Both inequalities for the equality cannot be in column \" \"orientation!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 870, __extension__
__PRETTY_FUNCTION__))
869 "Both inequalities for the equality cannot be in column "(static_cast <bool> (!(simplex.con[i].orientation == Orientation
::Column && simplex.con[i + 1].orientation == Orientation
::Column) && "Both inequalities for the equality cannot be in column "
"orientation!") ? void (0) : __assert_fail ("!(simplex.con[i].orientation == Orientation::Column && simplex.con[i + 1].orientation == Orientation::Column) && \"Both inequalities for the equality cannot be in column \" \"orientation!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 870, __extension__
__PRETTY_FUNCTION__))
870 "orientation!")(static_cast <bool> (!(simplex.con[i].orientation == Orientation
::Column && simplex.con[i + 1].orientation == Orientation
::Column) && "Both inequalities for the equality cannot be in column "
"orientation!") ? void (0) : __assert_fail ("!(simplex.con[i].orientation == Orientation::Column && simplex.con[i + 1].orientation == Orientation::Column) && \"Both inequalities for the equality cannot be in column \" \"orientation!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 870, __extension__
__PRETTY_FUNCTION__))
;
871 if (simplex.con[i].orientation == Orientation::Column)
872 dual.push_back(-simplex.tableau(row, simplex.con[i].pos));
873 else if (simplex.con[i + 1].orientation == Orientation::Column)
874 dual.push_back(simplex.tableau(row, simplex.con[i + 1].pos));
875 else
876 dual.push_back(0);
877 }
878 simplex.rollback(snap);
879 return *maybeWidth;
880 }
881
882 /// Remove the last equality that was added through addEqualityForDirection.
883 ///
884 /// We do this by rolling back to the snapshot at the top of the stack, which
885 /// should be a snapshot taken just before the last equality was added.
886 void removeLastEquality() {
887 assert(!snapshotStack.empty() && "Snapshot stack is empty!")(static_cast <bool> (!snapshotStack.empty() && "Snapshot stack is empty!"
) ? void (0) : __assert_fail ("!snapshotStack.empty() && \"Snapshot stack is empty!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 887, __extension__
__PRETTY_FUNCTION__))
;
888 simplex.rollback(snapshotStack.back());
889 snapshotStack.pop_back();
890 }
891
892private:
893 /// Returns coefficients of the expression 'dot_product(dir, x - y)',
894 /// i.e., dir_1 * x_1 + dir_2 * x_2 + ... + dir_n * x_n
895 /// - dir_1 * y_1 - dir_2 * y_2 - ... - dir_n * y_n,
896 /// where n is the dimension of the original polytope.
897 SmallVector<int64_t, 8> getCoeffsForDirection(ArrayRef<int64_t> dir) {
898 assert(2 * dir.size() == simplex.getNumVariables() &&(static_cast <bool> (2 * dir.size() == simplex.getNumVariables
() && "Direction vector has wrong dimensionality") ? void
(0) : __assert_fail ("2 * dir.size() == simplex.getNumVariables() && \"Direction vector has wrong dimensionality\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 899, __extension__
__PRETTY_FUNCTION__))
899 "Direction vector has wrong dimensionality")(static_cast <bool> (2 * dir.size() == simplex.getNumVariables
() && "Direction vector has wrong dimensionality") ? void
(0) : __assert_fail ("2 * dir.size() == simplex.getNumVariables() && \"Direction vector has wrong dimensionality\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 899, __extension__
__PRETTY_FUNCTION__))
;
900 SmallVector<int64_t, 8> coeffs(dir.begin(), dir.end());
901 coeffs.reserve(2 * dir.size());
902 for (int64_t coeff : dir)
903 coeffs.push_back(-coeff);
904 coeffs.push_back(0); // constant term
905 return coeffs;
906 }
907
908 Simplex simplex;
909 /// The first index of the equality constraints, the index immediately after
910 /// the last constraint in the initial product simplex.
911 unsigned simplexConstraintOffset;
912 /// A stack of snapshots, used for rolling back.
913 SmallVector<unsigned, 8> snapshotStack;
914};
915
916// Return a + scale*b;
917static SmallVector<int64_t, 8> scaleAndAdd(ArrayRef<int64_t> a, int64_t scale,
918 ArrayRef<int64_t> b) {
919 assert(a.size() == b.size())(static_cast <bool> (a.size() == b.size()) ? void (0) :
__assert_fail ("a.size() == b.size()", "mlir/lib/Analysis/Presburger/Simplex.cpp"
, 919, __extension__ __PRETTY_FUNCTION__))
;
920 SmallVector<int64_t, 8> res;
921 res.reserve(a.size());
922 for (unsigned i = 0, e = a.size(); i < e; ++i)
923 res.push_back(a[i] + scale * b[i]);
924 return res;
925}
926
927/// Reduce the basis to try and find a direction in which the polytope is
928/// "thin". This only works for bounded polytopes.
929///
930/// This is an implementation of the algorithm described in the paper
931/// "An Implementation of Generalized Basis Reduction for Integer Programming"
932/// by W. Cook, T. Rutherford, H. E. Scarf, D. Shallcross.
933///
934/// Let b_{level}, b_{level + 1}, ... b_n be the current basis.
935/// Let width_i(v) = max <v, x - y> where x and y are points in the original
936/// polytope such that <b_j, x - y> = 0 is satisfied for all level <= j < i.
937///
938/// In every iteration, we first replace b_{i+1} with b_{i+1} + u*b_i, where u
939/// is the integer such that width_i(b_{i+1} + u*b_i) is minimized. Let dual_i
940/// be the dual variable associated with the constraint <b_i, x - y> = 0 when
941/// computing width_{i+1}(b_{i+1}). It can be shown that dual_i is the
942/// minimizing value of u, if it were allowed to be fractional. Due to
943/// convexity, the minimizing integer value is either floor(dual_i) or
944/// ceil(dual_i), so we just need to check which of these gives a lower
945/// width_{i+1} value. If dual_i turned out to be an integer, then u = dual_i.
946///
947/// Now if width_i(b_{i+1}) < 0.75 * width_i(b_i), we swap b_i and (the new)
948/// b_{i + 1} and decrement i (unless i = level, in which case we stay at the
949/// same i). Otherwise, we increment i.
950///
951/// We keep f values and duals cached and invalidate them when necessary.
952/// Whenever possible, we use them instead of recomputing them. We implement the
953/// algorithm as follows.
954///
955/// In an iteration at i we need to compute:
956/// a) width_i(b_{i + 1})
957/// b) width_i(b_i)
958/// c) the integer u that minimizes width_i(b_{i + 1} + u*b_i)
959///
960/// If width_i(b_i) is not already cached, we compute it.
961///
962/// If the duals are not already cached, we compute width_{i+1}(b_{i+1}) and
963/// store the duals from this computation.
964///
965/// We call updateBasisWithUAndGetFCandidate, which finds the minimizing value
966/// of u as explained before, caches the duals from this computation, sets
967/// b_{i+1} to b_{i+1} + u*b_i, and returns the new value of width_i(b_{i+1}).
968///
969/// Now if width_i(b_{i+1}) < 0.75 * width_i(b_i), we swap b_i and b_{i+1} and
970/// decrement i, resulting in the basis
971/// ... b_{i - 1}, b_{i + 1} + u*b_i, b_i, b_{i+2}, ...
972/// with corresponding f values
973/// ... width_{i-1}(b_{i-1}), width_i(b_{i+1} + u*b_i), width_{i+1}(b_i), ...
974/// The values up to i - 1 remain unchanged. We have just gotten the middle
975/// value from updateBasisWithUAndGetFCandidate, so we can update that in the
976/// cache. The value at width_{i+1}(b_i) is unknown, so we evict this value from
977/// the cache. The iteration after decrementing needs exactly the duals from the
978/// computation of width_i(b_{i + 1} + u*b_i), so we keep these in the cache.
979///
980/// When incrementing i, no cached f values get invalidated. However, the cached
981/// duals do get invalidated as the duals for the higher levels are different.
982void Simplex::reduceBasis(Matrix &basis, unsigned level) {
983 const Fraction epsilon(3, 4);
984
985 if (level == basis.getNumRows() - 1)
12
Assuming the condition is false
13
Taking false branch
986 return;
987
988 GBRSimplex gbrSimplex(*this);
989 SmallVector<Fraction, 8> width;
990 SmallVector<int64_t, 8> dual;
991 int64_t dualDenom;
14
'dualDenom' declared without an initial value
992
993 // Finds the value of u that minimizes width_i(b_{i+1} + u*b_i), caches the
994 // duals from this computation, sets b_{i+1} to b_{i+1} + u*b_i, and returns
995 // the new value of width_i(b_{i+1}).
996 //
997 // If dual_i is not an integer, the minimizing value must be either
998 // floor(dual_i) or ceil(dual_i). We compute the expression for both and
999 // choose the minimizing value.
1000 //
1001 // If dual_i is an integer, we don't need to perform these computations. We
1002 // know that in this case,
1003 // a) u = dual_i.
1004 // b) one can show that dual_j for j < i are the same duals we would have
1005 // gotten from computing width_i(b_{i + 1} + u*b_i), so the correct duals
1006 // are the ones already in the cache.
1007 // c) width_i(b_{i+1} + u*b_i) = min_{alpha} width_i(b_{i+1} + alpha * b_i),
1008 // which
1009 // one can show is equal to width_{i+1}(b_{i+1}). The latter value must
1010 // be in the cache, so we get it from there and return it.
1011 auto updateBasisWithUAndGetFCandidate = [&](unsigned i) -> Fraction {
1012 assert(i < level + dual.size() && "dual_i is not known!")(static_cast <bool> (i < level + dual.size() &&
"dual_i is not known!") ? void (0) : __assert_fail ("i < level + dual.size() && \"dual_i is not known!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1012, __extension__
__PRETTY_FUNCTION__))
;
22
'?' condition is true
1013
1014 int64_t u = floorDiv(dual[i - level], dualDenom);
23
2nd function call argument is an uninitialized value
1015 basis.addToRow(i, i + 1, u);
1016 if (dual[i - level] % dualDenom != 0) {
1017 SmallVector<int64_t, 8> candidateDual[2];
1018 int64_t candidateDualDenom[2];
1019 Fraction widthI[2];
1020
1021 // Initially u is floor(dual) and basis reflects this.
1022 widthI[0] = gbrSimplex.computeWidthAndDuals(
1023 basis.getRow(i + 1), candidateDual[0], candidateDualDenom[0]);
1024
1025 // Now try ceil(dual), i.e. floor(dual) + 1.
1026 ++u;
1027 basis.addToRow(i, i + 1, 1);
1028 widthI[1] = gbrSimplex.computeWidthAndDuals(
1029 basis.getRow(i + 1), candidateDual[1], candidateDualDenom[1]);
1030
1031 unsigned j = widthI[0] < widthI[1] ? 0 : 1;
1032 if (j == 0)
1033 // Subtract 1 to go from u = ceil(dual) back to floor(dual).
1034 basis.addToRow(i, i + 1, -1);
1035
1036 // width_i(b{i+1} + u*b_i) should be minimized at our value of u.
1037 // We assert that this holds by checking that the values of width_i at
1038 // u - 1 and u + 1 are greater than or equal to the value at u. If the
1039 // width is lesser at either of the adjacent values, then our computed
1040 // value of u is clearly not the minimizer. Otherwise by convexity the
1041 // computed value of u is really the minimizer.
1042
1043 // Check the value at u - 1.
1044 assert(gbrSimplex.computeWidth(scaleAndAdd((static_cast <bool> (gbrSimplex.computeWidth(scaleAndAdd
( basis.getRow(i + 1), -1, basis.getRow(i))) >= widthI[j] &&
"Computed u value does not minimize the width!") ? void (0) :
__assert_fail ("gbrSimplex.computeWidth(scaleAndAdd( basis.getRow(i + 1), -1, basis.getRow(i))) >= widthI[j] && \"Computed u value does not minimize the width!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1046, __extension__
__PRETTY_FUNCTION__))
1045 basis.getRow(i + 1), -1, basis.getRow(i))) >= widthI[j] &&(static_cast <bool> (gbrSimplex.computeWidth(scaleAndAdd
( basis.getRow(i + 1), -1, basis.getRow(i))) >= widthI[j] &&
"Computed u value does not minimize the width!") ? void (0) :
__assert_fail ("gbrSimplex.computeWidth(scaleAndAdd( basis.getRow(i + 1), -1, basis.getRow(i))) >= widthI[j] && \"Computed u value does not minimize the width!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1046, __extension__
__PRETTY_FUNCTION__))
1046 "Computed u value does not minimize the width!")(static_cast <bool> (gbrSimplex.computeWidth(scaleAndAdd
( basis.getRow(i + 1), -1, basis.getRow(i))) >= widthI[j] &&
"Computed u value does not minimize the width!") ? void (0) :
__assert_fail ("gbrSimplex.computeWidth(scaleAndAdd( basis.getRow(i + 1), -1, basis.getRow(i))) >= widthI[j] && \"Computed u value does not minimize the width!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1046, __extension__
__PRETTY_FUNCTION__))
;
1047 // Check the value at u + 1.
1048 assert(gbrSimplex.computeWidth(scaleAndAdd((static_cast <bool> (gbrSimplex.computeWidth(scaleAndAdd
( basis.getRow(i + 1), +1, basis.getRow(i))) >= widthI[j] &&
"Computed u value does not minimize the width!") ? void (0) :
__assert_fail ("gbrSimplex.computeWidth(scaleAndAdd( basis.getRow(i + 1), +1, basis.getRow(i))) >= widthI[j] && \"Computed u value does not minimize the width!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1050, __extension__
__PRETTY_FUNCTION__))
1049 basis.getRow(i + 1), +1, basis.getRow(i))) >= widthI[j] &&(static_cast <bool> (gbrSimplex.computeWidth(scaleAndAdd
( basis.getRow(i + 1), +1, basis.getRow(i))) >= widthI[j] &&
"Computed u value does not minimize the width!") ? void (0) :
__assert_fail ("gbrSimplex.computeWidth(scaleAndAdd( basis.getRow(i + 1), +1, basis.getRow(i))) >= widthI[j] && \"Computed u value does not minimize the width!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1050, __extension__
__PRETTY_FUNCTION__))
1050 "Computed u value does not minimize the width!")(static_cast <bool> (gbrSimplex.computeWidth(scaleAndAdd
( basis.getRow(i + 1), +1, basis.getRow(i))) >= widthI[j] &&
"Computed u value does not minimize the width!") ? void (0) :
__assert_fail ("gbrSimplex.computeWidth(scaleAndAdd( basis.getRow(i + 1), +1, basis.getRow(i))) >= widthI[j] && \"Computed u value does not minimize the width!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1050, __extension__
__PRETTY_FUNCTION__))
;
1051
1052 dual = std::move(candidateDual[j]);
1053 dualDenom = candidateDualDenom[j];
1054 return widthI[j];
1055 }
1056
1057 assert(i + 1 - level < width.size() && "width_{i+1} wasn't saved")(static_cast <bool> (i + 1 - level < width.size() &&
"width_{i+1} wasn't saved") ? void (0) : __assert_fail ("i + 1 - level < width.size() && \"width_{i+1} wasn't saved\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1057, __extension__
__PRETTY_FUNCTION__))
;
1058 // f_i(b_{i+1} + dual*b_i) == width_{i+1}(b_{i+1}) when `dual` minimizes the
1059 // LHS. (note: the basis has already been updated, so b_{i+1} + dual*b_i in
1060 // the above expression is equal to basis.getRow(i+1) below.)
1061 assert(gbrSimplex.computeWidth(basis.getRow(i + 1)) ==(static_cast <bool> (gbrSimplex.computeWidth(basis.getRow
(i + 1)) == width[i + 1 - level]) ? void (0) : __assert_fail (
"gbrSimplex.computeWidth(basis.getRow(i + 1)) == width[i + 1 - level]"
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1062, __extension__
__PRETTY_FUNCTION__))
1062 width[i + 1 - level])(static_cast <bool> (gbrSimplex.computeWidth(basis.getRow
(i + 1)) == width[i + 1 - level]) ? void (0) : __assert_fail (
"gbrSimplex.computeWidth(basis.getRow(i + 1)) == width[i + 1 - level]"
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1062, __extension__
__PRETTY_FUNCTION__))
;
1063 return width[i + 1 - level];
1064 };
1065
1066 // In the ith iteration of the loop, gbrSimplex has constraints for directions
1067 // from `level` to i - 1.
1068 unsigned i = level;
1069 while (i < basis.getNumRows() - 1) {
15
Assuming the condition is true
16
Loop condition is true. Entering loop body
1070 if (i >= level + width.size()) {
17
Assuming the condition is false
18
Taking false branch
1071 // We don't even know the value of f_i(b_i), so let's find that first.
1072 // We have to do this first since later we assume that width already
1073 // contains values up to and including i.
1074
1075 assert((i == 0 || i - 1 < level + width.size()) &&(static_cast <bool> ((i == 0 || i - 1 < level + width
.size()) && "We are at level i but we don't know the value of width_{i-1}"
) ? void (0) : __assert_fail ("(i == 0 || i - 1 < level + width.size()) && \"We are at level i but we don't know the value of width_{i-1}\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1076, __extension__
__PRETTY_FUNCTION__))
1076 "We are at level i but we don't know the value of width_{i-1}")(static_cast <bool> ((i == 0 || i - 1 < level + width
.size()) && "We are at level i but we don't know the value of width_{i-1}"
) ? void (0) : __assert_fail ("(i == 0 || i - 1 < level + width.size()) && \"We are at level i but we don't know the value of width_{i-1}\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1076, __extension__
__PRETTY_FUNCTION__))
;
1077
1078 // We don't actually use these duals at all, but it doesn't matter
1079 // because this case should only occur when i is level, and there are no
1080 // duals in that case anyway.
1081 assert(i == level && "This case should only occur when i == level")(static_cast <bool> (i == level && "This case should only occur when i == level"
) ? void (0) : __assert_fail ("i == level && \"This case should only occur when i == level\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1081, __extension__
__PRETTY_FUNCTION__))
;
1082 width.push_back(
1083 gbrSimplex.computeWidthAndDuals(basis.getRow(i), dual, dualDenom));
1084 }
1085
1086 if (i >= level + dual.size()) {
19
Assuming the condition is false
20
Taking false branch
1087 assert(i + 1 >= level + width.size() &&(static_cast <bool> (i + 1 >= level + width.size() &&
"We don't know dual_i but we know width_{i+1}") ? void (0) :
__assert_fail ("i + 1 >= level + width.size() && \"We don't know dual_i but we know width_{i+1}\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1088, __extension__
__PRETTY_FUNCTION__))
1088 "We don't know dual_i but we know width_{i+1}")(static_cast <bool> (i + 1 >= level + width.size() &&
"We don't know dual_i but we know width_{i+1}") ? void (0) :
__assert_fail ("i + 1 >= level + width.size() && \"We don't know dual_i but we know width_{i+1}\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1088, __extension__
__PRETTY_FUNCTION__))
;
1089 // We don't know dual for our level, so let's find it.
1090 gbrSimplex.addEqualityForDirection(basis.getRow(i));
1091 width.push_back(gbrSimplex.computeWidthAndDuals(basis.getRow(i + 1), dual,
1092 dualDenom));
1093 gbrSimplex.removeLastEquality();
1094 }
1095
1096 // This variable stores width_i(b_{i+1} + u*b_i).
1097 Fraction widthICandidate = updateBasisWithUAndGetFCandidate(i);
21
Calling 'operator()'
1098 if (widthICandidate < epsilon * width[i - level]) {
1099 basis.swapRows(i, i + 1);
1100 width[i - level] = widthICandidate;
1101 // The values of width_{i+1}(b_{i+1}) and higher may change after the
1102 // swap, so we remove the cached values here.
1103 width.resize(i - level + 1);
1104 if (i == level) {
1105 dual.clear();
1106 continue;
1107 }
1108
1109 gbrSimplex.removeLastEquality();
1110 i--;
1111 continue;
1112 }
1113
1114 // Invalidate duals since the higher level needs to recompute its own duals.
1115 dual.clear();
1116 gbrSimplex.addEqualityForDirection(basis.getRow(i));
1117 i++;
1118 }
1119}
1120
1121/// Search for an integer sample point using a branch and bound algorithm.
1122///
1123/// Each row in the basis matrix is a vector, and the set of basis vectors
1124/// should span the space. Initially this is the identity matrix,
1125/// i.e., the basis vectors are just the variables.
1126///
1127/// In every level, a value is assigned to the level-th basis vector, as
1128/// follows. Compute the minimum and maximum rational values of this direction.
1129/// If only one integer point lies in this range, constrain the variable to
1130/// have this value and recurse to the next variable.
1131///
1132/// If the range has multiple values, perform generalized basis reduction via
1133/// reduceBasis and then compute the bounds again. Now we try constraining
1134/// this direction in the first value in this range and "recurse" to the next
1135/// level. If we fail to find a sample, we try assigning the direction the next
1136/// value in this range, and so on.
1137///
1138/// If no integer sample is found from any of the assignments, or if the range
1139/// contains no integer value, then of course the polytope is empty for the
1140/// current assignment of the values in previous levels, so we return to
1141/// the previous level.
1142///
1143/// If we reach the last level where all the variables have been assigned values
1144/// already, then we simply return the current sample point if it is integral,
1145/// and go back to the previous level otherwise.
1146///
1147/// To avoid potentially arbitrarily large recursion depths leading to stack
1148/// overflows, this algorithm is implemented iteratively.
1149Optional<SmallVector<int64_t, 8>> Simplex::findIntegerSample() {
1150 if (empty)
1
Assuming field 'empty' is false
2
Taking false branch
1151 return {};
1152
1153 unsigned nDims = var.size();
1154 Matrix basis = Matrix::identity(nDims);
1155
1156 unsigned level = 0;
1157 // The snapshot just before constraining a direction to a value at each level.
1158 SmallVector<unsigned, 8> snapshotStack;
1159 // The maximum value in the range of the direction for each level.
1160 SmallVector<int64_t, 8> upperBoundStack;
1161 // The next value to try constraining the basis vector to at each level.
1162 SmallVector<int64_t, 8> nextValueStack;
1163
1164 snapshotStack.reserve(basis.getNumRows());
1165 upperBoundStack.reserve(basis.getNumRows());
1166 nextValueStack.reserve(basis.getNumRows());
1167 while (level != -1u) {
3
Loop condition is true. Entering loop body
1168 if (level == basis.getNumRows()) {
4
Assuming the condition is false
5
Taking false branch
1169 // We've assigned values to all variables. Return if we have a sample,
1170 // or go back up to the previous level otherwise.
1171 if (auto maybeSample = getSamplePointIfIntegral())
1172 return maybeSample;
1173 level--;
1174 continue;
1175 }
1176
1177 if (level >= upperBoundStack.size()) {
6
Assuming the condition is true
7
Taking true branch
1178 // We haven't populated the stack values for this level yet, so we have
1179 // just come down a level ("recursed"). Find the lower and upper bounds.
1180 // If there is more than one integer point in the range, perform
1181 // generalized basis reduction.
1182 SmallVector<int64_t, 8> basisCoeffs =
1183 llvm::to_vector<8>(basis.getRow(level));
1184 basisCoeffs.push_back(0);
1185
1186 int64_t minRoundedUp, maxRoundedDown;
1187 std::tie(minRoundedUp, maxRoundedDown) =
1188 computeIntegerBounds(basisCoeffs);
1189
1190 // Heuristic: if the sample point is integral at this point, just return
1191 // it.
1192 if (auto maybeSample = getSamplePointIfIntegral())
8
Taking false branch
1193 return *maybeSample;
1194
1195 if (minRoundedUp < maxRoundedDown) {
9
Assuming 'minRoundedUp' is < 'maxRoundedDown'
10
Taking true branch
1196 reduceBasis(basis, level);
11
Calling 'Simplex::reduceBasis'
1197 basisCoeffs = llvm::to_vector<8>(basis.getRow(level));
1198 basisCoeffs.push_back(0);
1199 std::tie(minRoundedUp, maxRoundedDown) =
1200 computeIntegerBounds(basisCoeffs);
1201 }
1202
1203 snapshotStack.push_back(getSnapshot());
1204 // The smallest value in the range is the next value to try.
1205 nextValueStack.push_back(minRoundedUp);
1206 upperBoundStack.push_back(maxRoundedDown);
1207 }
1208
1209 assert((snapshotStack.size() - 1 == level &&(static_cast <bool> ((snapshotStack.size() - 1 == level
&& nextValueStack.size() - 1 == level && upperBoundStack
.size() - 1 == level) && "Mismatched variable stack sizes!"
) ? void (0) : __assert_fail ("(snapshotStack.size() - 1 == level && nextValueStack.size() - 1 == level && upperBoundStack.size() - 1 == level) && \"Mismatched variable stack sizes!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1212, __extension__
__PRETTY_FUNCTION__))
1210 nextValueStack.size() - 1 == level &&(static_cast <bool> ((snapshotStack.size() - 1 == level
&& nextValueStack.size() - 1 == level && upperBoundStack
.size() - 1 == level) && "Mismatched variable stack sizes!"
) ? void (0) : __assert_fail ("(snapshotStack.size() - 1 == level && nextValueStack.size() - 1 == level && upperBoundStack.size() - 1 == level) && \"Mismatched variable stack sizes!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1212, __extension__
__PRETTY_FUNCTION__))
1211 upperBoundStack.size() - 1 == level) &&(static_cast <bool> ((snapshotStack.size() - 1 == level
&& nextValueStack.size() - 1 == level && upperBoundStack
.size() - 1 == level) && "Mismatched variable stack sizes!"
) ? void (0) : __assert_fail ("(snapshotStack.size() - 1 == level && nextValueStack.size() - 1 == level && upperBoundStack.size() - 1 == level) && \"Mismatched variable stack sizes!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1212, __extension__
__PRETTY_FUNCTION__))
1212 "Mismatched variable stack sizes!")(static_cast <bool> ((snapshotStack.size() - 1 == level
&& nextValueStack.size() - 1 == level && upperBoundStack
.size() - 1 == level) && "Mismatched variable stack sizes!"
) ? void (0) : __assert_fail ("(snapshotStack.size() - 1 == level && nextValueStack.size() - 1 == level && upperBoundStack.size() - 1 == level) && \"Mismatched variable stack sizes!\""
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1212, __extension__
__PRETTY_FUNCTION__))
;
1213
1214 // Whether we "recursed" or "returned" from a lower level, we rollback
1215 // to the snapshot of the starting state at this level. (in the "recursed"
1216 // case this has no effect)
1217 rollback(snapshotStack.back());
1218 int64_t nextValue = nextValueStack.back();
1219 nextValueStack.back()++;
1220 if (nextValue > upperBoundStack.back()) {
1221 // We have exhausted the range and found no solution. Pop the stack and
1222 // return up a level.
1223 snapshotStack.pop_back();
1224 nextValueStack.pop_back();
1225 upperBoundStack.pop_back();
1226 level--;
1227 continue;
1228 }
1229
1230 // Try the next value in the range and "recurse" into the next level.
1231 SmallVector<int64_t, 8> basisCoeffs(basis.getRow(level).begin(),
1232 basis.getRow(level).end());
1233 basisCoeffs.push_back(-nextValue);
1234 addEquality(basisCoeffs);
1235 level++;
1236 }
1237
1238 return {};
1239}
1240
1241/// Compute the minimum and maximum integer values the expression can take. We
1242/// compute each separately.
1243std::pair<int64_t, int64_t>
1244Simplex::computeIntegerBounds(ArrayRef<int64_t> coeffs) {
1245 int64_t minRoundedUp;
1246 if (Optional<Fraction> maybeMin =
1247 computeOptimum(Simplex::Direction::Down, coeffs))
1248 minRoundedUp = ceil(*maybeMin);
1249 else
1250 llvm_unreachable("Tableau should not be unbounded")::llvm::llvm_unreachable_internal("Tableau should not be unbounded"
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1250)
;
1251
1252 int64_t maxRoundedDown;
1253 if (Optional<Fraction> maybeMax =
1254 computeOptimum(Simplex::Direction::Up, coeffs))
1255 maxRoundedDown = floor(*maybeMax);
1256 else
1257 llvm_unreachable("Tableau should not be unbounded")::llvm::llvm_unreachable_internal("Tableau should not be unbounded"
, "mlir/lib/Analysis/Presburger/Simplex.cpp", 1257)
;
1258
1259 return {minRoundedUp, maxRoundedDown};
1260}
1261
1262void SimplexBase::print(raw_ostream &os) const {
1263 os << "rows = " << nRow << ", columns = " << nCol << "\n";
1264 if (empty)
1265 os << "Simplex marked empty!\n";
1266 os << "var: ";
1267 for (unsigned i = 0; i < var.size(); ++i) {
1268 if (i > 0)
1269 os << ", ";
1270 var[i].print(os);
1271 }
1272 os << "\ncon: ";
1273 for (unsigned i = 0; i < con.size(); ++i) {
1274 if (i > 0)
1275 os << ", ";
1276 con[i].print(os);
1277 }
1278 os << '\n';
1279 for (unsigned row = 0; row < nRow; ++row) {
1280 if (row > 0)
1281 os << ", ";
1282 os << "r" << row << ": " << rowUnknown[row];
1283 }
1284 os << '\n';
1285 os << "c0: denom, c1: const";
1286 for (unsigned col = 2; col < nCol; ++col)
1287 os << ", c" << col << ": " << colUnknown[col];
1288 os << '\n';
1289 for (unsigned row = 0; row < nRow; ++row) {
1290 for (unsigned col = 0; col < nCol; ++col)
1291 os << tableau(row, col) << '\t';
1292 os << '\n';
1293 }
1294 os << '\n';
1295}
1296
1297void SimplexBase::dump() const { print(llvm::errs()); }
1298
1299bool Simplex::isRationalSubsetOf(const IntegerPolyhedron &poly) {
1300 if (isEmpty())
1301 return true;
1302
1303 for (unsigned i = 0, e = poly.getNumInequalities(); i < e; ++i)
1304 if (!isRedundantInequality(poly.getInequality(i)))
1305 return false;
1306
1307 for (unsigned i = 0, e = poly.getNumEqualities(); i < e; ++i)
1308 if (!isRedundantEquality(poly.getEquality(i)))
1309 return false;
1310
1311 return true;
1312}
1313
1314/// Computes the minimum value `coeffs` can take. If the value is greater than
1315/// or equal to zero, the polytope entirely lies in the half-space defined by
1316/// `coeffs >= 0`.
1317bool Simplex::isRedundantInequality(ArrayRef<int64_t> coeffs) {
1318 Optional<Fraction> minimum = computeOptimum(Direction::Down, coeffs);
1319 return minimum && *minimum >= Fraction(0, 1);
1320}
1321
1322/// Check whether the equality given by `coeffs == 0` is redundant given
1323/// the existing constraints. This is redundant when `coeffs` is already
1324/// always zero under the existing constraints. `coeffs` is always zero
1325/// when the minimum and maximum value that `coeffs` can take are both zero.
1326bool Simplex::isRedundantEquality(ArrayRef<int64_t> coeffs) {
1327 Optional<Fraction> minimum = computeOptimum(Direction::Down, coeffs);
1328 Optional<Fraction> maximum = computeOptimum(Direction::Up, coeffs);
1329 return minimum && maximum && *maximum == Fraction(0, 1) &&
1330 *minimum == Fraction(0, 1);
1331}
1332
1333} // namespace mlir