1 /* Linear_System class implementation: non-inline template functions.
2 Copyright (C) 2001-2010 Roberto Bagnara <bagnara@cs.unipr.it>
3 Copyright (C) 2010-2016 BUGSENG srl (http://bugseng.com)
4
5 This file is part of the Parma Polyhedra Library (PPL).
6
7 The PPL is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11
12 The PPL is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software Foundation,
19 Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA.
20
21 For the most up-to-date information see the Parma Polyhedra Library
22 site: http://bugseng.com/products/ppl/ . */
23
24 #ifndef PPL_Linear_System_templates_hh
25 #define PPL_Linear_System_templates_hh 1
26
27 #include "Bit_Matrix_defs.hh"
28 #include "Linear_Expression_defs.hh"
29 #include "Scalar_Products_defs.hh"
30 #include "Scalar_Products_inlines.hh"
31 #include "swapping_sort_templates.hh"
32 #include <algorithm>
33 #include <iostream>
34 #include <string>
35 #include <deque>
36
37 namespace Parma_Polyhedra_Library {
38
39 template <typename Row>
40 dimension_type
num_lines_or_equalities() const41 Linear_System<Row>::num_lines_or_equalities() const {
42 PPL_ASSERT(num_pending_rows() == 0);
43 const Linear_System& x = *this;
44 dimension_type n = 0;
45 for (dimension_type i = num_rows(); i-- > 0; ) {
46 if (x[i].is_line_or_equality()) {
47 ++n;
48 }
49 }
50 return n;
51 }
52
53 template <typename Row>
54 void
merge_rows_assign(const Linear_System & y)55 Linear_System<Row>::merge_rows_assign(const Linear_System& y) {
56 PPL_ASSERT(space_dimension() >= y.space_dimension());
57 // Both systems have to be sorted and have no pending rows.
58 PPL_ASSERT(check_sorted() && y.check_sorted());
59 PPL_ASSERT(num_pending_rows() == 0 && y.num_pending_rows() == 0);
60
61 Linear_System& x = *this;
62
63 // A temporary vector...
64 Swapping_Vector<Row> tmp;
65 // ... with enough capacity not to require any reallocations.
66 tmp.reserve(compute_capacity(x.rows.size() + y.rows.size(),
67 tmp.max_num_rows()));
68
69 dimension_type xi = 0;
70 const dimension_type x_num_rows = x.num_rows();
71 dimension_type yi = 0;
72 const dimension_type y_num_rows = y.num_rows();
73
74 while (xi < x_num_rows && yi < y_num_rows) {
75 const int comp = compare(x[xi], y[yi]);
76 if (comp <= 0) {
77 // Elements that can be taken from `x' are actually _stolen_ from `x'
78 tmp.resize(tmp.size() + 1);
79 swap(tmp.back(), x.rows[xi++]);
80 tmp.back().set_representation(representation());
81 if (comp == 0) {
82 // A duplicate element.
83 ++yi;
84 }
85 }
86 else {
87 // (comp > 0)
88 tmp.resize(tmp.size() + 1);
89 Row copy(y[yi++], space_dimension(), representation());
90 swap(tmp.back(), copy);
91 }
92 }
93 // Insert what is left.
94 if (xi < x_num_rows) {
95 while (xi < x_num_rows) {
96 tmp.resize(tmp.size() + 1);
97 swap(tmp.back(), x.rows[xi++]);
98 tmp.back().set_representation(representation());
99 }
100 }
101 else {
102 while (yi < y_num_rows) {
103 tmp.resize(tmp.size() + 1);
104 Row copy(y[yi++], space_dimension(), representation());
105 swap(tmp.back(), copy);
106 }
107 }
108
109 // We get the result matrix and let the old one be destroyed.
110 swap(tmp, rows);
111 // There are no pending rows.
112 unset_pending_rows();
113 PPL_ASSERT(check_sorted());
114 PPL_ASSERT(OK());
115 }
116
117 template <typename Row>
118 void
ascii_dump(std::ostream & s) const119 Linear_System<Row>::ascii_dump(std::ostream& s) const {
120 // Prints the topology, the number of rows, the number of columns
121 // and the sorted flag. The specialized methods provided by
122 // Constraint_System and Generator_System take care of properly
123 // printing the contents of the system.
124 s << "topology " << (is_necessarily_closed()
125 ? "NECESSARILY_CLOSED"
126 : "NOT_NECESSARILY_CLOSED")
127 << "\n"
128 << num_rows() << " x " << space_dimension() << " ";
129 Parma_Polyhedra_Library::ascii_dump(s, representation());
130 s << " " << (sorted ? "(sorted)" : "(not_sorted)")
131 << "\n"
132 << "index_first_pending " << first_pending_row()
133 << "\n";
134 for (dimension_type i = 0; i < rows.size(); ++i) {
135 rows[i].ascii_dump(s);
136 }
137 }
138
PPL_OUTPUT_TEMPLATE_DEFINITIONS_ASCII_ONLY(Row,Linear_System<Row>)139 PPL_OUTPUT_TEMPLATE_DEFINITIONS_ASCII_ONLY(Row, Linear_System<Row>)
140
141 template <typename Row>
142 bool
143 Linear_System<Row>::ascii_load(std::istream& s) {
144 std::string str;
145 if (!(s >> str) || str != "topology") {
146 return false;
147 }
148 if (!(s >> str)) {
149 return false;
150 }
151
152 clear();
153
154 Topology t;
155 if (str == "NECESSARILY_CLOSED") {
156 t = NECESSARILY_CLOSED;
157 }
158 else {
159 if (str != "NOT_NECESSARILY_CLOSED") {
160 return false;
161 }
162 t = NOT_NECESSARILY_CLOSED;
163 }
164
165 set_topology(t);
166
167 dimension_type nrows;
168 dimension_type space_dims;
169 if (!(s >> nrows)) {
170 return false;
171 }
172 if (!(s >> str) || str != "x") {
173 return false;
174 }
175 if (!(s >> space_dims)) {
176 return false;
177 }
178
179 space_dimension_ = space_dims;
180
181 if (!Parma_Polyhedra_Library::ascii_load(s, representation_)) {
182 return false;
183 }
184
185 if (!(s >> str) || (str != "(sorted)" && str != "(not_sorted)")) {
186 return false;
187 }
188 const bool sortedness = (str == "(sorted)");
189 dimension_type index;
190 if (!(s >> str) || str != "index_first_pending") {
191 return false;
192 }
193 if (!(s >> index)) {
194 return false;
195 }
196
197 Row row;
198 for (dimension_type i = 0; i < nrows; ++i) {
199 if (!row.ascii_load(s)) {
200 return false;
201 }
202 insert(row, Recycle_Input());
203 }
204 index_first_pending = index;
205 sorted = sortedness;
206
207 // Check invariants.
208 PPL_ASSERT(OK());
209 return true;
210 }
211
212 template <typename Row>
213 void
insert(const Row & r)214 Linear_System<Row>::insert(const Row& r) {
215 Row tmp(r, representation());
216 insert(tmp, Recycle_Input());
217 }
218
219 template <typename Row>
220 void
insert(Row & r,Recycle_Input)221 Linear_System<Row>::insert(Row& r, Recycle_Input) {
222 insert_no_ok(r, Recycle_Input());
223 PPL_ASSERT(OK());
224 }
225
226 template <typename Row>
227 void
insert_no_ok(Row & r,Recycle_Input)228 Linear_System<Row>::insert_no_ok(Row& r, Recycle_Input) {
229 PPL_ASSERT(topology() == r.topology());
230 // This method is only used when the system has no pending rows.
231 PPL_ASSERT(num_pending_rows() == 0);
232
233 const bool was_sorted = is_sorted();
234
235 insert_pending_no_ok(r, Recycle_Input());
236
237 if (was_sorted) {
238 const dimension_type nrows = num_rows();
239 // The added row may have caused the system to be not sorted anymore.
240 if (nrows > 1) {
241 // If the system is not empty and the inserted row is the
242 // greatest one, the system is set to be sorted.
243 // If it is not the greatest one then the system is no longer sorted.
244 sorted = (compare(rows[nrows-2], rows[nrows-1]) <= 0);
245 }
246 else {
247 // A system having only one row is sorted.
248 sorted = true;
249 }
250 }
251
252 unset_pending_rows();
253 }
254
255 template <typename Row>
256 void
insert_pending_no_ok(Row & r,Recycle_Input)257 Linear_System<Row>::insert_pending_no_ok(Row& r, Recycle_Input) {
258 // TODO: A Grid_Generator_System may contain non-normalized lines that
259 // represent parameters, so this check is disabled. Consider re-enabling it
260 // when it's possibile.
261 #if 0
262 // The added row must be strongly normalized and have the same
263 // number of elements as the existing rows of the system.
264 PPL_ASSERT(r.check_strong_normalized());
265 #endif
266
267 PPL_ASSERT(r.topology() == topology());
268
269 r.set_representation(representation());
270
271 if (space_dimension() < r.space_dimension()) {
272 set_space_dimension_no_ok(r.space_dimension());
273 }
274 else {
275 r.set_space_dimension_no_ok(space_dimension());
276 }
277
278 rows.resize(rows.size() + 1);
279 swap(rows.back(), r);
280 }
281
282 template <typename Row>
283 void
insert_pending(const Row & r)284 Linear_System<Row>::insert_pending(const Row& r) {
285 Row tmp(r, representation());
286 insert_pending(tmp, Recycle_Input());
287 }
288
289 template <typename Row>
290 void
insert_pending(Row & r,Recycle_Input)291 Linear_System<Row>::insert_pending(Row& r, Recycle_Input) {
292 insert_pending_no_ok(r, Recycle_Input());
293 PPL_ASSERT(OK());
294 }
295
296 template <typename Row>
297 void
insert_pending(const Linear_System & y)298 Linear_System<Row>::insert_pending(const Linear_System& y) {
299 Linear_System tmp(y, representation(), With_Pending());
300 insert_pending(tmp, Recycle_Input());
301 }
302
303 template <typename Row>
304 void
insert_pending(Linear_System & y,Recycle_Input)305 Linear_System<Row>::insert_pending(Linear_System& y, Recycle_Input) {
306 Linear_System& x = *this;
307 PPL_ASSERT(x.space_dimension() == y.space_dimension());
308
309 // Steal the rows of `y'.
310 // This loop must use an increasing index (instead of a decreasing one) to
311 // preserve the row ordering.
312 for (dimension_type i = 0; i < y.num_rows(); ++i) {
313 x.insert_pending(y.rows[i], Recycle_Input());
314 }
315
316 y.clear();
317
318 PPL_ASSERT(x.OK());
319 }
320
321 template <typename Row>
322 void
insert(const Linear_System & y)323 Linear_System<Row>::insert(const Linear_System& y) {
324 Linear_System tmp(y, representation(), With_Pending());
325 insert(tmp, Recycle_Input());
326 }
327
328 template <typename Row>
329 void
insert(Linear_System & y,Recycle_Input)330 Linear_System<Row>::insert(Linear_System& y, Recycle_Input) {
331 PPL_ASSERT(num_pending_rows() == 0);
332
333 // Adding no rows is a no-op.
334 if (y.has_no_rows()) {
335 return;
336 }
337
338 // Check if sortedness is preserved.
339 if (is_sorted()) {
340 if (!y.is_sorted() || y.num_pending_rows() > 0) {
341 sorted = false;
342 }
343 else {
344 // `y' is sorted and has no pending rows.
345 const dimension_type n_rows = num_rows();
346 if (n_rows > 0) {
347 sorted = (compare(rows[n_rows-1], y[0]) <= 0);
348 }
349 }
350 }
351
352 // Add the rows of `y' as if they were pending.
353 insert_pending(y, Recycle_Input());
354
355 // TODO: May y have pending rows? Should they remain pending?
356
357 // There are no pending_rows.
358 unset_pending_rows();
359
360 PPL_ASSERT(OK());
361 }
362
363 template <typename Row>
364 void
remove_space_dimensions(const Variables_Set & vars)365 Linear_System<Row>::remove_space_dimensions(const Variables_Set& vars) {
366 // Dimension-compatibility assertion.
367 PPL_ASSERT(space_dimension() >= vars.space_dimension());
368
369 // The removal of no dimensions from any system is a no-op. This
370 // case also captures the only legal removal of dimensions from a
371 // 0-dim system.
372 if (vars.empty()) {
373 return;
374 }
375
376 // NOTE: num_rows() is *not* constant, because it may be decreased by
377 // remove_row_no_ok().
378 for (dimension_type i = 0; i < num_rows(); ) {
379 const bool valid = rows[i].remove_space_dimensions(vars);
380 if (!valid) {
381 // Remove the current row.
382 // We can't call remove_row(i) here, because the system is not OK as
383 // some rows already have the new space dimension and others still have
384 // the old one.
385 remove_row_no_ok(i, false);
386 }
387 else {
388 ++i;
389 }
390 }
391
392 space_dimension_ -= vars.size();
393
394 PPL_ASSERT(OK());
395 }
396
397 template <typename Row>
398 void
shift_space_dimensions(Variable v,dimension_type n)399 Linear_System<Row>::shift_space_dimensions(Variable v, dimension_type n) {
400 // NOTE: v.id() may be equal to the space dimension of the system
401 // (when no space dimension need to be shifted).
402 PPL_ASSERT(v.id() <= space_dimension());
403 for (dimension_type i = rows.size(); i-- > 0; ) {
404 rows[i].shift_space_dimensions(v, n);
405 }
406 space_dimension_ += n;
407 PPL_ASSERT(OK());
408 }
409
410 template <typename Row>
411 void
sort_rows()412 Linear_System<Row>::sort_rows() {
413 // We sort the non-pending rows only.
414 sort_rows(0, first_pending_row());
415 sorted = true;
416 PPL_ASSERT(OK());
417 }
418
419 template <typename Row>
420 void
sort_rows(const dimension_type first_row,const dimension_type last_row)421 Linear_System<Row>::sort_rows(const dimension_type first_row,
422 const dimension_type last_row) {
423 PPL_ASSERT(first_row <= last_row && last_row <= num_rows());
424 // We cannot mix pending and non-pending rows.
425 PPL_ASSERT(first_row >= first_pending_row()
426 || last_row <= first_pending_row());
427
428 const bool sorting_pending = (first_row >= first_pending_row());
429 const dimension_type old_num_pending = num_pending_rows();
430
431 const dimension_type num_elems = last_row - first_row;
432 if (num_elems < 2) {
433 return;
434 }
435
436 // Build the function objects implementing indirect sort comparison,
437 // indirect unique comparison and indirect swap operation.
438 using namespace Implementation;
439 typedef Swapping_Vector<Row> Cont;
440 typedef Indirect_Sort_Compare<Cont, Row_Less_Than> Sort_Compare;
441 typedef Indirect_Swapper<Cont> Swapper;
442 const dimension_type num_duplicates
443 = indirect_sort_and_unique(num_elems,
444 Sort_Compare(rows, first_row),
445 Unique_Compare(rows, first_row),
446 Swapper(rows, first_row));
447
448 if (num_duplicates > 0) {
449 typedef typename Cont::iterator Iter;
450 typedef typename std::iterator_traits<Iter>::difference_type diff_t;
451 Iter last = rows.begin() + static_cast<diff_t>(last_row);
452 Iter first = last - + static_cast<diff_t>(num_duplicates);
453 rows.erase(first, last);
454 }
455
456 if (sorting_pending) {
457 PPL_ASSERT(old_num_pending >= num_duplicates);
458 index_first_pending = num_rows() - (old_num_pending - num_duplicates);
459 }
460 else {
461 index_first_pending = num_rows() - old_num_pending;
462 }
463
464 PPL_ASSERT(OK());
465 }
466
467 template <typename Row>
468 void
strong_normalize()469 Linear_System<Row>::strong_normalize() {
470 const dimension_type nrows = rows.size();
471 // We strongly normalize also the pending rows.
472 for (dimension_type i = nrows; i-- > 0; ) {
473 rows[i].strong_normalize();
474 }
475 sorted = (nrows <= 1);
476 PPL_ASSERT(OK());
477 }
478
479 template <typename Row>
480 void
sign_normalize()481 Linear_System<Row>::sign_normalize() {
482 const dimension_type nrows = rows.size();
483 // We sign-normalize also the pending rows.
484 for (dimension_type i = nrows; i-- > 0; ) {
485 rows[i].sign_normalize();
486 }
487 sorted = (nrows <= 1);
488 PPL_ASSERT(OK());
489 }
490
491 /*! \relates Parma_Polyhedra_Library::Linear_System */
492 template <typename Row>
493 bool
operator ==(const Linear_System<Row> & x,const Linear_System<Row> & y)494 operator==(const Linear_System<Row>& x, const Linear_System<Row>& y) {
495 if (x.space_dimension() != y.space_dimension()) {
496 return false;
497 }
498 const dimension_type x_num_rows = x.num_rows();
499 const dimension_type y_num_rows = y.num_rows();
500 if (x_num_rows != y_num_rows) {
501 return false;
502 }
503 if (x.first_pending_row() != y.first_pending_row()) {
504 return false;
505 }
506 // TODO: Check if the following comment is up to date.
507 // Notice that calling operator==(const Swapping_Vector<Row>&,
508 // const Swapping_Vector<Row>&)
509 // would be wrong here, as equality of the type fields would
510 // not be checked.
511 for (dimension_type i = x_num_rows; i-- > 0; ) {
512 if (x[i] != y[i]) {
513 return false;
514 }
515 }
516 return true;
517 }
518
519 template <typename Row>
520 void
sort_and_remove_with_sat(Bit_Matrix & sat)521 Linear_System<Row>::sort_and_remove_with_sat(Bit_Matrix& sat) {
522 // We can only sort the non-pending part of the system.
523 PPL_ASSERT(first_pending_row() == sat.num_rows());
524 if (first_pending_row() <= 1) {
525 set_sorted(true);
526 return;
527 }
528
529 const dimension_type num_elems = sat.num_rows();
530 // Build the function objects implementing indirect sort comparison,
531 // indirect unique comparison and indirect swap operation.
532 typedef Swapping_Vector<Row> Cont;
533 const Implementation::Indirect_Sort_Compare<Cont, Row_Less_Than>
534 sort_cmp(rows);
535 const Unique_Compare unique_cmp(rows);
536 const Implementation::Indirect_Swapper2<Cont, Bit_Matrix> swapper(rows, sat);
537
538 const dimension_type num_duplicates
539 = Implementation::indirect_sort_and_unique(num_elems, sort_cmp,
540 unique_cmp, swapper);
541
542 const dimension_type new_first_pending_row
543 = first_pending_row() - num_duplicates;
544
545 if (num_pending_rows() > 0) {
546 // In this case, we must put the duplicates after the pending rows.
547 const dimension_type n_rows = num_rows() - 1;
548 for (dimension_type i = 0; i < num_duplicates; ++i) {
549 swap(rows[new_first_pending_row + i], rows[n_rows - i]);
550 }
551 }
552
553 // Erasing the duplicated rows...
554 rows.resize(rows.size() - num_duplicates);
555 index_first_pending = new_first_pending_row;
556 // ... and the corresponding rows of the saturation matrix.
557 sat.remove_trailing_rows(num_duplicates);
558
559 // Now the system is sorted.
560 sorted = true;
561
562 PPL_ASSERT(OK());
563 }
564
565 template <typename Row>
566 dimension_type
gauss(const dimension_type n_lines_or_equalities)567 Linear_System<Row>::gauss(const dimension_type n_lines_or_equalities) {
568 // This method is only applied to a linear system having no pending rows and
569 // exactly `n_lines_or_equalities' lines or equalities, all of which occur
570 // before the rays or points or inequalities.
571 PPL_ASSERT(num_pending_rows() == 0);
572 PPL_ASSERT(n_lines_or_equalities == num_lines_or_equalities());
573 #ifndef NDEBUG
574 for (dimension_type i = n_lines_or_equalities; i-- > 0; ) {
575 PPL_ASSERT((*this)[i].is_line_or_equality());
576 }
577 #endif
578
579 dimension_type rank = 0;
580 // Will keep track of the variations on the system of equalities.
581 bool changed = false;
582 // TODO: Don't use the number of columns.
583 const dimension_type num_cols
584 = is_necessarily_closed() ? space_dimension() + 1 : space_dimension() + 2;
585 // TODO: Consider exploiting the row (possible) sparseness of rows in the
586 // following loop, if needed. It would probably make it more cache-efficient
587 // for dense rows, too.
588 for (dimension_type j = num_cols; j-- > 0; ) {
589 for (dimension_type i = rank; i < n_lines_or_equalities; ++i) {
590 // Search for the first row having a non-zero coefficient
591 // (the pivot) in the j-th column.
592 if ((*this)[i].expr.get(j) == 0) {
593 continue;
594 }
595 // Pivot found: if needed, swap rows so that this one becomes
596 // the rank-th row in the linear system.
597 if (i > rank) {
598 swap(rows[i], rows[rank]);
599 // After swapping the system is no longer sorted.
600 changed = true;
601 }
602 // Combine the row containing the pivot with all the lines or
603 // equalities following it, so that all the elements on the j-th
604 // column in these rows become 0.
605 for (dimension_type k = i + 1; k < n_lines_or_equalities; ++k) {
606 if (rows[k].expr.get(Variable(j - 1)) != 0) {
607 rows[k].linear_combine(rows[rank], j);
608 changed = true;
609 }
610 }
611 // Already dealt with the rank-th row.
612 ++rank;
613 // Consider another column index `j'.
614 break;
615 }
616 }
617 if (changed) {
618 sorted = false;
619 }
620
621 PPL_ASSERT(OK());
622 return rank;
623 }
624
625 template <typename Row>
626 void
627 Linear_System<Row>
back_substitute(const dimension_type n_lines_or_equalities)628 ::back_substitute(const dimension_type n_lines_or_equalities) {
629 // This method is only applied to a system having no pending rows and
630 // exactly `n_lines_or_equalities' lines or equalities, all of which occur
631 // before the first ray or point or inequality.
632 PPL_ASSERT(num_pending_rows() == 0);
633 PPL_ASSERT(n_lines_or_equalities <= num_lines_or_equalities());
634 #ifndef NDEBUG
635 for (dimension_type i = n_lines_or_equalities; i-- > 0; ) {
636 PPL_ASSERT((*this)[i].is_line_or_equality());
637 }
638 #endif
639
640 const dimension_type nrows = num_rows();
641 // Trying to keep sortedness.
642 bool still_sorted = is_sorted();
643 // This deque of Booleans will be used to flag those rows that,
644 // before exiting, need to be re-checked for sortedness.
645 std::deque<bool> check_for_sortedness;
646 if (still_sorted) {
647 check_for_sortedness.insert(check_for_sortedness.end(), nrows, false);
648 }
649
650 for (dimension_type k = n_lines_or_equalities; k-- > 0; ) {
651 // For each line or equality, starting from the last one,
652 // looks for the last non-zero element.
653 // `j' will be the index of such a element.
654 Row& row_k = rows[k];
655 const dimension_type j = row_k.expr.last_nonzero();
656 // TODO: Check this.
657 PPL_ASSERT(j != 0);
658
659 // Go through the equalities above `row_k'.
660 for (dimension_type i = k; i-- > 0; ) {
661 Row& row_i = rows[i];
662 if (row_i.expr.get(Variable(j - 1)) != 0) {
663 // Combine linearly `row_i' with `row_k'
664 // so that `row_i[j]' becomes zero.
665 row_i.linear_combine(row_k, j);
666 if (still_sorted) {
667 // Trying to keep sortedness: remember which rows
668 // have to be re-checked for sortedness at the end.
669 if (i > 0) {
670 check_for_sortedness[i-1] = true;
671 }
672 check_for_sortedness[i] = true;
673 }
674 }
675 }
676
677 // Due to strong normalization during previous iterations,
678 // the pivot coefficient `row_k[j]' may now be negative.
679 // Since an inequality (or ray or point) cannot be multiplied
680 // by a negative factor, the coefficient of the pivot must be
681 // forced to be positive.
682 const bool have_to_negate = (row_k.expr.get(Variable(j - 1)) < 0);
683 if (have_to_negate) {
684 neg_assign(row_k.expr);
685 }
686
687 // NOTE: Here row_k will *not* be ok if we have negated it.
688
689 // Note: we do not mark index `k' in `check_for_sortedness',
690 // because we will later negate back the row.
691
692 // Go through all the other rows of the system.
693 for (dimension_type i = n_lines_or_equalities; i < nrows; ++i) {
694 Row& row_i = rows[i];
695 if (row_i.expr.get(Variable(j - 1)) != 0) {
696 // Combine linearly the `row_i' with `row_k'
697 // so that `row_i[j]' becomes zero.
698 row_i.linear_combine(row_k, j);
699 if (still_sorted) {
700 // Trying to keep sortedness: remember which rows
701 // have to be re-checked for sortedness at the end.
702 if (i > n_lines_or_equalities) {
703 check_for_sortedness[i-1] = true;
704 }
705 check_for_sortedness[i] = true;
706 }
707 }
708 }
709 if (have_to_negate) {
710 // Negate `row_k' to restore strong-normalization.
711 neg_assign(row_k.expr);
712 }
713
714 PPL_ASSERT(row_k.OK());
715 }
716
717 // Trying to keep sortedness.
718 for (dimension_type i = 0; still_sorted && i+1 < nrows; ++i) {
719 if (check_for_sortedness[i]) {
720 // Have to check sortedness of `(*this)[i]' with respect to `(*this)[i+1]'.
721 still_sorted = (compare((*this)[i], (*this)[i+1]) <= 0);
722 }
723 }
724
725 // Set the sortedness flag.
726 sorted = still_sorted;
727
728 PPL_ASSERT(OK());
729 }
730
731 template <typename Row>
732 void
simplify()733 Linear_System<Row>::simplify() {
734 // This method is only applied to a system having no pending rows.
735 PPL_ASSERT(num_pending_rows() == 0);
736
737 // Partially sort the linear system so that all lines/equalities come first.
738 const dimension_type old_nrows = num_rows();
739 dimension_type nrows = old_nrows;
740 dimension_type n_lines_or_equalities = 0;
741 for (dimension_type i = 0; i < nrows; ++i) {
742 if ((*this)[i].is_line_or_equality()) {
743 if (n_lines_or_equalities < i) {
744 swap(rows[i], rows[n_lines_or_equalities]);
745 // The system was not sorted.
746 PPL_ASSERT(!sorted);
747 }
748 ++n_lines_or_equalities;
749 }
750 }
751 // Apply Gaussian elimination to the subsystem of lines/equalities.
752 const dimension_type rank = gauss(n_lines_or_equalities);
753 // Eliminate any redundant line/equality that has been detected.
754 if (rank < n_lines_or_equalities) {
755 const dimension_type
756 n_rays_or_points_or_inequalities = nrows - n_lines_or_equalities;
757 const dimension_type
758 num_swaps = std::min(n_lines_or_equalities - rank,
759 n_rays_or_points_or_inequalities);
760 for (dimension_type i = num_swaps; i-- > 0; ) {
761 swap(rows[--nrows], rows[rank + i]);
762 }
763 remove_trailing_rows(old_nrows - nrows);
764 if (n_rays_or_points_or_inequalities > num_swaps) {
765 set_sorted(false);
766 }
767 unset_pending_rows();
768 n_lines_or_equalities = rank;
769 }
770 // Apply back-substitution to the system of rays/points/inequalities.
771 back_substitute(n_lines_or_equalities);
772
773 PPL_ASSERT(OK());
774 }
775
776 template <typename Row>
777 void
778 Linear_System<Row>
add_universe_rows_and_space_dimensions(const dimension_type n)779 ::add_universe_rows_and_space_dimensions(const dimension_type n) {
780 PPL_ASSERT(n > 0);
781 const bool was_sorted = is_sorted();
782 const dimension_type old_n_rows = num_rows();
783 const dimension_type old_space_dim
784 = is_necessarily_closed() ? space_dimension() : space_dimension() + 1;
785 set_space_dimension(space_dimension() + n);
786 rows.resize(rows.size() + n);
787 // The old system is moved to the bottom.
788 for (dimension_type i = old_n_rows; i-- > 0; ) {
789 swap(rows[i], rows[i + n]);
790 }
791 for (dimension_type i = n, c = old_space_dim; i-- > 0; ) {
792 // The top right-hand sub-system (i.e., the system made of new
793 // rows and columns) is set to the specular image of the identity
794 // matrix.
795 if (Variable(c).space_dimension() <= space_dimension()) {
796 // Variable(c) is a user variable.
797 Linear_Expression le(representation());
798 le.set_space_dimension(space_dimension());
799 le += Variable(c);
800 Row r(le, Row::LINE_OR_EQUALITY, row_topology);
801 swap(r, rows[i]);
802 }
803 else {
804 // Variable(c) is the epsilon dimension.
805 PPL_ASSERT(row_topology == NOT_NECESSARILY_CLOSED);
806 Linear_Expression le(Variable(c), representation());
807 Row r(le, Row::LINE_OR_EQUALITY, NECESSARILY_CLOSED);
808 r.mark_as_not_necessarily_closed();
809 swap(r, rows[i]);
810 // Note: `r' is strongly normalized.
811 }
812 ++c;
813 }
814 // If the old system was empty, the last row added is either
815 // a positivity constraint or a point.
816 if (was_sorted) {
817 sorted = (compare(rows[n-1], rows[n]) <= 0);
818 }
819
820 // If the system is not necessarily closed, move the epsilon coefficients to
821 // the last column.
822 if (!is_necessarily_closed()) {
823 // Try to preserve sortedness of `gen_sys'.
824 PPL_ASSERT(old_space_dim != 0);
825 if (!is_sorted()) {
826 for (dimension_type i = n; i-- > 0; ) {
827 rows[i].expr.swap_space_dimensions(Variable(old_space_dim - 1),
828 Variable(old_space_dim - 1 + n));
829 PPL_ASSERT(rows[i].OK());
830 }
831 }
832 else {
833 dimension_type old_eps_index = old_space_dim - 1;
834 // The upper-right corner of `rows' contains the J matrix:
835 // swap coefficients to preserve sortedness.
836 for (dimension_type i = n; i-- > 0; ++old_eps_index) {
837 rows[i].expr.swap_space_dimensions(Variable(old_eps_index),
838 Variable(old_eps_index + 1));
839 PPL_ASSERT(rows[i].OK());
840 }
841
842 sorted = true;
843 }
844 }
845 // NOTE: this already checks for OK().
846 set_index_first_pending_row(index_first_pending + n);
847 }
848
849 template <typename Row>
850 void
sort_pending_and_remove_duplicates()851 Linear_System<Row>::sort_pending_and_remove_duplicates() {
852 PPL_ASSERT(num_pending_rows() > 0);
853 PPL_ASSERT(is_sorted());
854
855 // The non-pending part of the system is already sorted.
856 // Now sorting the pending part..
857 const dimension_type first_pending = first_pending_row();
858 sort_rows(first_pending, num_rows());
859 // Recompute the number of rows, because we may have removed
860 // some rows occurring more than once in the pending part.
861 const dimension_type old_num_rows = num_rows();
862 dimension_type num_rows = old_num_rows;
863
864 dimension_type k1 = 0;
865 dimension_type k2 = first_pending;
866 dimension_type num_duplicates = 0;
867 // In order to erase them, put at the end of the system
868 // those pending rows that also occur in the non-pending part.
869 while (k1 < first_pending && k2 < num_rows) {
870 const int cmp = compare(rows[k1], rows[k2]);
871 if (cmp == 0) {
872 // We found the same row.
873 ++num_duplicates;
874 --num_rows;
875 // By initial sortedness, we can increment index `k1'.
876 ++k1;
877 // Do not increment `k2'; instead, swap there the next pending row.
878 if (k2 < num_rows) {
879 swap(rows[k2], rows[k2 + num_duplicates]);
880 }
881 }
882 else if (cmp < 0) {
883 // By initial sortedness, we can increment `k1'.
884 ++k1;
885 }
886 else {
887 // Here `cmp > 0'.
888 // Increment `k2' and, if we already found any duplicate,
889 // swap the next pending row in position `k2'.
890 ++k2;
891 if (num_duplicates > 0 && k2 < num_rows) {
892 swap(rows[k2], rows[k2 + num_duplicates]);
893 }
894 }
895 }
896 // If needed, swap any duplicates found past the pending rows
897 // that has not been considered yet; then erase the duplicates.
898 if (num_duplicates > 0) {
899 if (k2 < num_rows) {
900 for (++k2; k2 < num_rows; ++k2) {
901 swap(rows[k2], rows[k2 + num_duplicates]);
902 }
903 }
904 rows.resize(num_rows);
905 }
906 sorted = true;
907 PPL_ASSERT(OK());
908 }
909
910 template <typename Row>
911 bool
check_sorted() const912 Linear_System<Row>::check_sorted() const {
913 for (dimension_type i = first_pending_row(); i-- > 1; ) {
914 if (compare(rows[i], rows[i-1]) < 0) {
915 return false;
916 }
917 }
918 return true;
919 }
920
921 template <typename Row>
922 bool
OK() const923 Linear_System<Row>::OK() const {
924 #ifndef NDEBUG
925 using std::endl;
926 using std::cerr;
927 #endif
928
929 for (dimension_type i = rows.size(); i-- > 0; ) {
930 if (rows[i].representation() != representation()) {
931 #ifndef NDEBUG
932 cerr << "Linear_System has a row with the wrong representation!"
933 << endl;
934 #endif
935 return false;
936 }
937 if (rows[i].space_dimension() != space_dimension()) {
938 #ifndef NDEBUG
939 cerr << "Linear_System has a row with the wrong number of space dimensions!"
940 << endl;
941 #endif
942 return false;
943 }
944 }
945
946 for (dimension_type i = rows.size(); i-- > 0; ) {
947 if (rows[i].topology() != topology()) {
948 #ifndef NDEBUG
949 cerr << "Linear_System has a row with the wrong topology!"
950 << endl;
951 #endif
952 return false;
953 }
954
955 }
956 // `index_first_pending' must be less than or equal to `num_rows()'.
957 if (first_pending_row() > num_rows()) {
958 #ifndef NDEBUG
959 cerr << "Linear_System has a negative number of pending rows!"
960 << endl;
961 #endif
962 return false;
963 }
964
965 // Check for topology mismatches.
966 const dimension_type n_rows = num_rows();
967 for (dimension_type i = 0; i < n_rows; ++i) {
968 if (topology() != rows[i].topology()) {
969 #ifndef NDEBUG
970 cerr << "Topology mismatch between the system "
971 << "and one of its rows!"
972 << endl;
973 #endif
974 return false;
975 }
976
977 }
978 if (sorted && !check_sorted()) {
979 #ifndef NDEBUG
980 cerr << "The system declares itself to be sorted but it is not!"
981 << endl;
982 #endif
983 return false;
984 }
985
986 // All checks passed.
987 return true;
988 }
989
990 } // namespace Parma_Polyhedra_Library
991
992 #endif // !defined(PPL_Linear_System_templates_hh)
993