00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #ifndef PPL_DB_Matrix_inlines_hh
00024 #define PPL_DB_Matrix_inlines_hh 1
00025
00026 #include "globals.defs.hh"
00027 #include "Checked_Number.defs.hh"
00028 #include <cassert>
00029 #include <iostream>
00030
00031 namespace Parma_Polyhedra_Library {
00032
00033 template <typename T>
00034 inline void
00035 DB_Matrix<T>::swap(DB_Matrix& y) {
00036 std::swap(rows, y.rows);
00037 std::swap(row_size, y.row_size);
00038 std::swap(row_capacity, y.row_capacity);
00039 }
00040
00041 template <typename T>
00042 inline dimension_type
00043 DB_Matrix<T>::max_num_rows() {
00044 return std::vector<DB_Row<T> >().max_size();
00045 }
00046
00047 template <typename T>
00048 inline dimension_type
00049 DB_Matrix<T>::max_num_columns() {
00050 return DB_Row<T>::max_size();
00051 }
00052
00053 template <typename T>
00054 inline
00055 DB_Matrix<T>::const_iterator::const_iterator()
00056 : i(Iter()) {
00057 }
00058
00059 template <typename T>
00060 inline
00061 DB_Matrix<T>::const_iterator::const_iterator(const Iter& b)
00062 : i(b) {
00063 }
00064
00065 template <typename T>
00066 inline
00067 DB_Matrix<T>::const_iterator::const_iterator(const const_iterator& y)
00068 : i(y.i) {
00069 }
00070
00071 template <typename T>
00072 inline typename DB_Matrix<T>::const_iterator&
00073 DB_Matrix<T>::const_iterator::operator=(const const_iterator& y) {
00074 i = y.i;
00075 return *this;
00076 }
00077
00078 template <typename T>
00079 inline typename DB_Matrix<T>::const_iterator::reference
00080 DB_Matrix<T>::const_iterator::operator*() const {
00081 return *i;
00082 }
00083
00084 template <typename T>
00085 inline typename DB_Matrix<T>::const_iterator::pointer
00086 DB_Matrix<T>::const_iterator::operator->() const {
00087 return &*i;
00088 }
00089
00090 template <typename T>
00091 inline typename DB_Matrix<T>::const_iterator&
00092 DB_Matrix<T>::const_iterator::operator++() {
00093 ++i;
00094 return *this;
00095 }
00096
00097 template <typename T>
00098 inline typename DB_Matrix<T>::const_iterator
00099 DB_Matrix<T>::const_iterator::operator++(int) {
00100 return const_iterator(i++);
00101 }
00102
00103 template <typename T>
00104 inline bool
00105 DB_Matrix<T>::const_iterator::operator==(const const_iterator& y) const {
00106 return i == y.i;
00107 }
00108
00109 template <typename T>
00110 inline bool
00111 DB_Matrix<T>::const_iterator::operator!=(const const_iterator& y) const {
00112 return !operator==(y);
00113 }
00114
00115 template <typename T>
00116 inline typename DB_Matrix<T>::const_iterator
00117 DB_Matrix<T>::begin() const {
00118 return const_iterator(rows.begin());
00119 }
00120
00121 template <typename T>
00122 inline typename DB_Matrix<T>::const_iterator
00123 DB_Matrix<T>::end() const {
00124 return const_iterator(rows.end());
00125 }
00126
00127 template <typename T>
00128 inline
00129 DB_Matrix<T>::DB_Matrix()
00130 : rows(),
00131 row_size(0),
00132 row_capacity(0) {
00133 }
00134
00135 template <typename T>
00136 inline
00137 DB_Matrix<T>::~DB_Matrix() {
00138 }
00139
00140 template <typename T>
00141 inline DB_Row<T>&
00142 DB_Matrix<T>::operator[](const dimension_type k) {
00143 assert(k < rows.size());
00144 return rows[k];
00145 }
00146
00147 template <typename T>
00148 inline const DB_Row<T>&
00149 DB_Matrix<T>::operator[](const dimension_type k) const {
00150 assert(k < rows.size());
00151 return rows[k];
00152 }
00153
00154 template <typename T>
00155 inline dimension_type
00156 DB_Matrix<T>::num_rows() const {
00157 return rows.size();
00158 }
00159
00160 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00161
00162 #endif // PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00163 template <typename T>
00164 inline bool
00165 operator!=(const DB_Matrix<T>& x, const DB_Matrix<T>& y) {
00166 return !(x == y);
00167 }
00168
00169 template <typename T>
00170 inline
00171 DB_Matrix<T>::DB_Matrix(const dimension_type n_rows)
00172 : rows(n_rows),
00173 row_size(n_rows),
00174 row_capacity(compute_capacity(n_rows, max_num_columns())) {
00175
00176 for (dimension_type i = 0; i < n_rows; ++i)
00177 rows[i].construct(n_rows, row_capacity);
00178 assert(OK());
00179 }
00180
00181 template <typename T>
00182 inline
00183 DB_Matrix<T>::DB_Matrix(const DB_Matrix& y)
00184 : rows(y.rows),
00185 row_size(y.row_size),
00186 row_capacity(compute_capacity(y.row_size, max_num_columns())) {
00187 }
00188
00189 template <typename T>
00190 template <typename U>
00191 inline
00192 DB_Matrix<T>::DB_Matrix(const DB_Matrix<U>& y)
00193 : rows(y.rows.size()),
00194 row_size(y.row_size),
00195 row_capacity(compute_capacity(y.row_size, max_num_columns())) {
00196
00197 for (dimension_type i = 0, n_rows = rows.size(); i < n_rows; ++i)
00198 rows[i].construct_upward_approximation(y[i], row_capacity);
00199 assert(OK());
00200 }
00201
00202 template <typename T>
00203 inline DB_Matrix<T>&
00204 DB_Matrix<T>::operator=(const DB_Matrix& y) {
00205
00206
00207
00208
00209 if (this != &y) {
00210
00211 rows = y.rows;
00212 row_size = y.row_size;
00213
00214
00215 row_capacity = compute_capacity(y.row_size, max_num_columns());
00216 }
00217 return *this;
00218 }
00219
00220 template <typename T>
00221 void
00222 DB_Matrix<T>::grow(const dimension_type new_n_rows) {
00223 const dimension_type old_n_rows = rows.size();
00224 assert(new_n_rows >= old_n_rows);
00225
00226 if (new_n_rows > old_n_rows) {
00227 if (new_n_rows <= row_capacity) {
00228
00229 if (rows.capacity() < new_n_rows) {
00230
00231 std::vector<DB_Row<T> > new_rows;
00232 new_rows.reserve(compute_capacity(new_n_rows, max_num_rows()));
00233 new_rows.insert(new_rows.end(), new_n_rows, DB_Row<T>());
00234
00235 dimension_type i = new_n_rows;
00236 while (i-- > old_n_rows)
00237 new_rows[i].construct(new_n_rows, row_capacity);
00238
00239 ++i;
00240 while (i-- > 0)
00241 new_rows[i].swap(rows[i]);
00242
00243 std::swap(rows, new_rows);
00244 }
00245 else {
00246
00247 rows.insert(rows.end(), new_n_rows - old_n_rows, DB_Row<T>());
00248 for (dimension_type i = new_n_rows; i-- > old_n_rows; )
00249 rows[i].construct(new_n_rows, row_capacity);
00250 }
00251 }
00252 else {
00253
00254 DB_Matrix new_matrix;
00255 new_matrix.rows.reserve(compute_capacity(new_n_rows, max_num_rows()));
00256 new_matrix.rows.insert(new_matrix.rows.end(), new_n_rows, DB_Row<T>());
00257
00258 new_matrix.row_size = new_n_rows;
00259 new_matrix.row_capacity = compute_capacity(new_n_rows,
00260 max_num_columns());
00261 dimension_type i = new_n_rows;
00262 while (i-- > old_n_rows)
00263 new_matrix.rows[i].construct(new_matrix.row_size,
00264 new_matrix.row_capacity);
00265
00266 ++i;
00267 while (i-- > 0) {
00268 DB_Row<T> new_row(rows[i],
00269 new_matrix.row_size,
00270 new_matrix.row_capacity);
00271 std::swap(new_matrix.rows[i], new_row);
00272 }
00273
00274 swap(new_matrix);
00275 return;
00276 }
00277 }
00278
00279 if (new_n_rows > row_size) {
00280
00281 if (new_n_rows <= row_capacity)
00282
00283 for (dimension_type i = old_n_rows; i-- > 0; )
00284 rows[i].expand_within_capacity(new_n_rows);
00285 else {
00286
00287
00288 const dimension_type new_row_capacity
00289 = compute_capacity(new_n_rows, max_num_columns());
00290 for (dimension_type i = old_n_rows; i-- > 0; ) {
00291 DB_Row<T> new_row(rows[i], new_n_rows, new_row_capacity);
00292 std::swap(rows[i], new_row);
00293 }
00294 row_capacity = new_row_capacity;
00295 }
00296
00297 row_size = new_n_rows;
00298 }
00299 }
00300
00301 template <typename T>
00302 void
00303 DB_Matrix<T>::resize_no_copy(const dimension_type new_n_rows) {
00304 dimension_type old_n_rows = rows.size();
00305
00306 if (new_n_rows > old_n_rows) {
00307
00308 if (new_n_rows <= row_capacity) {
00309
00310 if (rows.capacity() < new_n_rows) {
00311
00312 std::vector<DB_Row<T> > new_rows;
00313 new_rows.reserve(compute_capacity(new_n_rows, max_num_rows()));
00314 new_rows.insert(new_rows.end(), new_n_rows, DB_Row<T>());
00315
00316
00317 dimension_type i = new_n_rows;
00318 while (i-- > old_n_rows)
00319 new_rows[i].construct(new_n_rows, row_capacity);
00320
00321 ++i;
00322 while (i-- > 0)
00323 new_rows[i].swap(rows[i]);
00324
00325 std::swap(rows, new_rows);
00326 }
00327 else {
00328
00329 rows.insert(rows.end(), new_n_rows - old_n_rows, DB_Row<T>());
00330
00331
00332 for (dimension_type i = new_n_rows; i-- > old_n_rows; )
00333 rows[i].construct(new_n_rows, row_capacity);
00334 }
00335 }
00336 else {
00337
00338 DB_Matrix new_matrix(new_n_rows);
00339 swap(new_matrix);
00340 return;
00341 }
00342 }
00343 else if (new_n_rows < old_n_rows) {
00344
00345 rows.erase(rows.begin() + new_n_rows, rows.end());
00346
00347 for (dimension_type i = new_n_rows; i-- > 0; )
00348 rows[i].shrink(new_n_rows);
00349 old_n_rows = new_n_rows;
00350 }
00351
00352 if (new_n_rows > row_size) {
00353
00354 if (new_n_rows <= row_capacity)
00355
00356 for (dimension_type i = old_n_rows; i-- > 0; )
00357 rows[i].expand_within_capacity(new_n_rows);
00358 else {
00359
00360
00361 const dimension_type new_row_capacity
00362 = compute_capacity(new_n_rows, max_num_columns());
00363 for (dimension_type i = old_n_rows; i-- > 0; ) {
00364 DB_Row<T> new_row(new_n_rows, new_row_capacity);
00365 std::swap(rows[i], new_row);
00366 }
00367 row_capacity = new_row_capacity;
00368 }
00369 }
00370
00371 row_size = new_n_rows;
00372 }
00373
00374 template <typename T>
00375 void
00376 DB_Matrix<T>::ascii_dump(std::ostream& s) const {
00377 const DB_Matrix<T>& x = *this;
00378 const char separator = ' ';
00379 const dimension_type nrows = x.num_rows();
00380 s << nrows << separator << "\n";
00381 for (dimension_type i = 0; i < nrows; ++i) {
00382 for (dimension_type j = 0; j < nrows; ++j) {
00383 using namespace IO_Operators;
00384 s << x[i][j] << separator;
00385 }
00386 s << "\n";
00387 }
00388 }
00389
00390 PPL_OUTPUT_TEMPLATE_DEFINITIONS(T, DB_Matrix<T>);
00391
00392 template <typename T>
00393 bool
00394 DB_Matrix<T>::ascii_load(std::istream& s) {
00395 dimension_type nrows;
00396 if (!(s >> nrows))
00397 return false;
00398 resize_no_copy(nrows);
00399 DB_Matrix& x = *this;
00400 for (dimension_type i = 0; i < nrows; ++i)
00401 for (dimension_type j = 0; j < nrows; ++j) {
00402 Result r = input(x[i][j], s, ROUND_UP);
00403
00404 if (!s || r == V_CVT_STR_UNK)
00405 return false;
00406 }
00407
00408 assert(OK());
00409 return true;
00410 }
00411
00412
00413 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00414
00415 #endif // PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00416 template <typename T>
00417 inline bool
00418 operator==(const DB_Matrix<T>& x, const DB_Matrix<T>& y) {
00419 const dimension_type x_num_rows = x.num_rows();
00420 if (x_num_rows != y.num_rows())
00421 return false;
00422 for (dimension_type i = x_num_rows; i-- > 0; )
00423 if (x[i] != y[i])
00424 return false;
00425 return true;
00426 }
00427
00428 template <typename To, typename From>
00429 struct maybe_assign_struct {
00430 static inline Result
00431 function(const To*& top, To& tmp, const From& from, Rounding_Dir dir) {
00432
00433
00434 top = &tmp;
00435 return assign_r(tmp, from, dir);
00436 }
00437 };
00438
00439 template <typename Type>
00440 struct maybe_assign_struct<Type, Type> {
00441 static inline Result
00442 function(const Type*& top, Type&, const Type& from, Rounding_Dir) {
00443
00444 top = &from;
00445 return V_EQ;
00446 }
00447 };
00448
00449 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00450
00456 #endif // PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00457 template <typename To, typename From>
00458 inline Result
00459 maybe_assign(const To*& top, To& tmp, const From& from, Rounding_Dir dir) {
00460 return maybe_assign_struct<To, From>::function(top, tmp, from, dir);
00461 }
00462
00463 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00464
00465 #endif // PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00466 template <typename Specialization, typename Temp, typename To, typename T>
00467 inline bool
00468 l_m_distance_assign(Checked_Number<To, Extended_Number_Policy>& r,
00469 const DB_Matrix<T>& x,
00470 const DB_Matrix<T>& y,
00471 const Rounding_Dir dir,
00472 Temp& tmp0,
00473 Temp& tmp1,
00474 Temp& tmp2) {
00475 const dimension_type x_num_rows = x.num_rows();
00476 if (x_num_rows != y.num_rows())
00477 return false;
00478 assign_r(tmp0, 0, ROUND_NOT_NEEDED);
00479 for (dimension_type i = x_num_rows; i-- > 0; ) {
00480 const DB_Row<T>& x_i = x[i];
00481 const DB_Row<T>& y_i = y[i];
00482 for (dimension_type j = x_num_rows; j-- > 0; ) {
00483 const T& x_i_j = x_i[j];
00484 const T& y_i_j = y_i[j];
00485 if (is_plus_infinity(x_i_j)) {
00486 if (is_plus_infinity(y_i_j))
00487 continue;
00488 else {
00489 pinf:
00490 r = PLUS_INFINITY;
00491 return true;
00492 }
00493 }
00494 else if (is_plus_infinity(y_i_j))
00495 goto pinf;
00496
00497 const Temp* tmp1p;
00498 const Temp* tmp2p;
00499 if (x_i_j > y_i_j) {
00500 maybe_assign(tmp1p, tmp1, x_i_j, dir);
00501 maybe_assign(tmp2p, tmp2, y_i_j, inverse(dir));
00502 }
00503 else {
00504 maybe_assign(tmp1p, tmp1, y_i_j, dir);
00505 maybe_assign(tmp2p, tmp2, x_i_j, inverse(dir));
00506 }
00507 sub_assign_r(tmp1, *tmp1p, *tmp2p, dir);
00508 assert(tmp1 >= 0);
00509 Specialization::combine(tmp0, tmp1, dir);
00510 }
00511 }
00512 Specialization::finalize(tmp0, dir);
00513 assign_r(r, tmp0, dir);
00514 return true;
00515 }
00516
00517 template <typename Temp>
00518 struct Rectilinear_Distance_Specialization {
00519 static inline void
00520 combine(Temp& running, const Temp& current, Rounding_Dir dir) {
00521 add_assign_r(running, running, current, dir);
00522 }
00523
00524 static inline void
00525 finalize(Temp&, Rounding_Dir) {
00526 }
00527 };
00528
00529 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00530
00531 #endif // PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00532 template <typename Temp, typename To, typename T>
00533 inline bool
00534 rectilinear_distance_assign(Checked_Number<To, Extended_Number_Policy>& r,
00535 const DB_Matrix<T>& x,
00536 const DB_Matrix<T>& y,
00537 const Rounding_Dir dir,
00538 Temp& tmp0,
00539 Temp& tmp1,
00540 Temp& tmp2) {
00541 return
00542 l_m_distance_assign<Rectilinear_Distance_Specialization<Temp> >(r, x, y,
00543 dir,
00544 tmp0,
00545 tmp1,
00546 tmp2);
00547 }
00548
00549
00550 template <typename Temp>
00551 struct Euclidean_Distance_Specialization {
00552 static inline void
00553 combine(Temp& running, Temp& current, Rounding_Dir dir) {
00554 mul_assign_r(current, current, current, dir);
00555 add_assign_r(running, running, current, dir);
00556 }
00557
00558 static inline void
00559 finalize(Temp& running, Rounding_Dir dir) {
00560 sqrt_assign_r(running, running, dir);
00561 }
00562 };
00563
00564 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00565
00566 #endif // PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00567 template <typename Temp, typename To, typename T>
00568 inline bool
00569 euclidean_distance_assign(Checked_Number<To, Extended_Number_Policy>& r,
00570 const DB_Matrix<T>& x,
00571 const DB_Matrix<T>& y,
00572 const Rounding_Dir dir,
00573 Temp& tmp0,
00574 Temp& tmp1,
00575 Temp& tmp2) {
00576 return
00577 l_m_distance_assign<Euclidean_Distance_Specialization<Temp> >(r, x, y,
00578 dir,
00579 tmp0,
00580 tmp1,
00581 tmp2);
00582 }
00583
00584
00585 template <typename Temp>
00586 struct L_Infinity_Distance_Specialization {
00587 static inline void
00588 combine(Temp& running, const Temp& current, Rounding_Dir) {
00589 if (current > running)
00590 running = current;
00591 }
00592
00593 static inline void
00594 finalize(Temp&, Rounding_Dir) {
00595 }
00596 };
00597
00598 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00599
00600 #endif // PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00601 template <typename Temp, typename To, typename T>
00602 inline bool
00603 l_infinity_distance_assign(Checked_Number<To, Extended_Number_Policy>& r,
00604 const DB_Matrix<T>& x,
00605 const DB_Matrix<T>& y,
00606 const Rounding_Dir dir,
00607 Temp& tmp0,
00608 Temp& tmp1,
00609 Temp& tmp2) {
00610 return
00611 l_m_distance_assign<L_Infinity_Distance_Specialization<Temp> >(r, x, y,
00612 dir,
00613 tmp0,
00614 tmp1,
00615 tmp2);
00616 }
00617
00618 template <typename T>
00619 bool
00620 DB_Matrix<T>::OK() const {
00621 #ifndef NDEBUG
00622 using std::endl;
00623 using std::cerr;
00624 #endif
00625
00626
00627 if (num_rows() != row_size) {
00628 #ifndef NDEBUG
00629 cerr << "DB_Matrix has fewer columns than rows:\n"
00630 << "row_size is " << row_size
00631 << ", num_rows() is " << num_rows() << "!"
00632 << endl;
00633 #endif
00634 return false;
00635 }
00636
00637 const DB_Matrix& x = *this;
00638 const dimension_type n_rows = x.num_rows();
00639 for (dimension_type i = 0; i < n_rows; ++i) {
00640 if (!x[i].OK(row_size, row_capacity))
00641 return false;
00642 }
00643
00644
00645 return true;
00646 }
00647
00648 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00649
00650 #endif // PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00651 template <typename T>
00652 std::ostream&
00653 IO_Operators::operator<<(std::ostream& s, const DB_Matrix<T>& c) {
00654 const dimension_type n = c.num_rows();
00655 for (dimension_type i = 0; i < n; ++i) {
00656 for (dimension_type j = 0; j < n; ++j)
00657 s << c[i][j] << " ";
00658 s << "\n";
00659 }
00660 return s;
00661 }
00662
00663 }
00664
00665 namespace std {
00666
00667 #ifdef PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00668
00669 #endif // PPL_DOXYGEN_INCLUDE_IMPLEMENTATION_DETAILS
00670 template <typename T>
00671 inline void
00672 swap(Parma_Polyhedra_Library::DB_Matrix<T>& x,
00673 Parma_Polyhedra_Library::DB_Matrix<T>& y) {
00674 x.swap(y);
00675 }
00676
00677 }
00678
00679 #endif // !defined(PPL_DB_Matrix_inlines_hh)