CMR  1.3.0
matrix.hpp
Go to the documentation of this file.
1 #pragma once
2 
3 #include <cmr/config.h>
4 #include <cmr/export.h>
5 
6 #include <boost/numeric/ublas/matrix.hpp>
7 
8 #include "matrix_transposed.hpp"
9 #include <iomanip>
10 
11 namespace tu
12 {
13  namespace detail
14  {
15 
16  template <typename Iterator>
17  class Range
18  {
19  public:
20  Range(Iterator begin, Iterator end)
21  : _begin(begin), _end(end)
22  {
23 
24  }
25 
26  Range(const Range<Iterator>& other)
27  : _begin(other._begin), _end(other._end)
28  {
29 
30  }
31 
32  Iterator begin()
33  {
34  return _begin;
35  }
36 
37  Iterator end()
38  {
39  return _end;
40  }
41 
42  private:
43  Iterator _begin;
44  Iterator _end;
45  };
46 
47  } /* namespace tu::detail */
48 
53  template <typename V>
54  class Matrix
55  {
56  public:
57  typedef V Value;
58  typedef std::size_t Index;
59 
60  struct Nonzero
61  {
64  const Value& value;
65 
66  Nonzero(Index r, Index c, const Value& v)
67  : row(r), column(c), value(v)
68  {
69 
70  }
71  };
72  };
73 
78  template <typename V>
79  class DenseMatrix : public Matrix<V>
80  {
81  public:
82  typedef typename Matrix<V>::Value Value;
83  typedef typename Matrix<V>::Index Index;
84  typedef typename Matrix<V>::Nonzero Nonzero;
85 
86  template <bool Row>
88  {
89  public:
90  NonzeroIterator(const DenseMatrix<V>& matrix, Index major)
91  : _matrix(matrix), _major(major), _minor(0)
92  {
93  if (Row)
94  {
95  while (_minor < _matrix.numColumns() && _matrix.get(_major, _minor) == 0)
96  ++_minor;
97  }
98  else
99  {
100  while (_minor < _matrix.numRows() && _matrix.get(_minor, _major) == 0)
101  ++_minor;
102  }
103  }
104 
105  NonzeroIterator(const DenseMatrix<V>& matrix, Index major, int dummy)
106  : _matrix(matrix), _major(major)
107  {
108  CMR_UNUSED(dummy);
109 
110  if (Row)
111  _minor = _matrix.numColumns();
112  else
113  _minor = _matrix.numRows();
114  }
115 
117  {
118  ++_minor;
119  if (Row)
120  {
121  while (_minor < _matrix.numColumns() && _matrix.get(_major, _minor) == 0)
122  ++_minor;
123  }
124  else
125  {
126  while (_minor < _matrix.numRows() && _matrix.get(_minor, _major) == 0)
127  ++_minor;
128  }
129  }
130 
132  {
133  if (Row)
134  {
135  assert(_major < _matrix.numRows());
136  assert(_minor < _matrix.numColumns());
137  return Nonzero(_major, _minor, _matrix.get(_major, _minor));
138  }
139  else
140  {
141  assert(_minor < _matrix.numRows());
142  assert(_major < _matrix.numColumns());
143  return Nonzero(_minor, _major, _matrix.get(_minor, _major));
144  }
145  }
146 
147  bool operator!=(const NonzeroIterator<Row>& other) const
148  {
149  assert(&_matrix == &other._matrix);
150  assert(_major == other._major);
151  return _minor != other._minor;
152  }
153 
154  private:
155  const DenseMatrix<V>& _matrix;
156  std::size_t _major, _minor;
157  };
158 
161 
167  : _data(nullptr), _numRows(0), _numColumns(0)
168  {
169 
170  }
171 
177  : _data(other._data), _numRows(other._numRows), _numColumns(other._numColumns)
178  {
179  other._data = nullptr;
180  other._numRows = 0;
181  other._numColumns = 0;
182  }
183 
188  template <typename V2>
190  : _numRows(other._numRows), _numColumns(other._numColumns)
191  {
192  std::size_t size = other._numRows * other._numColumns;
193  _data = new Value[size];
194  for (std::size_t i = 0; i < size; ++i)
195  _data[i] = other._data[i];
196  }
197 
203  {
204  _data = other._data;
205  _numRows = other._numRows;
206  _numColumns = other._numColumns;
207  other._data = nullptr;
208  other._numRows = 0;
209  other._numColumns = 0;
210  return *this;
211  }
212 
217  template <typename V2>
219  {
220  _numRows = other._numRows;
221  _numColumns = other._numColumns;
222  std::size_t size = other._numRows * other._numColumns;
223  _data = new Value[size];
224  for (std::size_t i = 0; i < size; ++i)
225  _data[i] = other._data[i];
226  return *this;
227  }
228 
234  {
235  if (_data != nullptr)
236  delete[] _data;
237  }
238 
245  {
246  std::size_t size = numRows * numColumns;
247  _data = new Value[size];
248  for (std::size_t i = 0; i < size; ++i)
249  _data[i] = 0;
250  }
251 
256  Index numRows() const
257  {
258  return _numRows;
259  }
260 
266  {
267  return _numColumns;
268  }
269 
274  const Value& get(Index row, Index column) const
275  {
276  assert(row < _numRows);
277  assert(column < _numColumns);
278  return _data[_numColumns * row + column];
279  }
280 
285  template <typename V2>
286  void set(Index row, Index column, const V2& value)
287  {
288  assert(row < _numRows);
289  assert(column < _numColumns);
290  _data[_numColumns * row + column] = value;
291  }
292 
297  std::size_t countRowNonzeros(Index row) const
298  {
299  std::size_t begin = _numColumns * row;
300  std::size_t end = _numColumns * (row + 1);
301  std::size_t count = 0;
302  for (std::size_t i = begin; i < end; ++i)
303  {
304  if (_data[i] != 0)
305  ++count;
306  }
307  return count;
308  }
309 
314  std::size_t countColumnNonzeros(Index column) const
315  {
316  std::size_t begin = column;
317  std::size_t end = _numRows * _numColumns;
318  std::size_t count = 0;
319  for (std::size_t i = begin; i < end; i += _numColumns)
320  {
321  if (_data[i] != 0)
322  ++count;
323  }
324  return count;
325  }
326 
332  {
333  return NonzeroRowRange(NonzeroIterator<true>(*this, row),
334  NonzeroIterator<true>(*this, row, 0));
335  }
336 
342  {
343  return NonzeroColumnRange(NonzeroIterator<false>(*this, column),
344  NonzeroIterator<false>(*this, column, 0));
345  }
346 
351  void ensureConsistency() const
352  {
353 
354  }
355 
356  protected:
360  };
361 
366  template<typename V>
367  class SparseMatrix : public Matrix<V>
368  {
369  public:
370  typedef typename Matrix<V>::Value Value;
371  typedef typename Matrix<V>::Index Index;
372  typedef typename Matrix<V>::Nonzero Nonzero;
373 
381  struct Data
382  {
384  std::vector<std::pair<Index, Index>> range;
386  std::vector<std::pair<Index, Value>> entries;
388  bool isSorted;
389 
395  : isSorted(true)
396  {
397 
398  }
399 
404  Data(Data&& other)
405  : range(std::move(other.range)), entries(std::move(other.entries)), isSorted(other.isSorted)
406  {
407 
408  }
409 
414  Data(const Data& other)
415  : range(other.range), entries(other.entries), isSorted(other.isSorted)
416  {
417 
418  }
419 
424  Data& operator=(Data&& other)
425  {
426  range = std::move(other.range);
427  entries = std::move(other.entries);
428  isSorted = other.isSorted;
429  return *this;
430  }
431 
436  Data& operator=(const Data& other)
437  {
438  range = other.range;
439  entries = other.entries;
440  isSorted = other.isSorted;
441  return *this;
442  }
443 
448  ~Data() = default;
449 
455  void constructFromTransposed(const Data& transposedData, Index numMinor)
456  {
457  // Iterate over entries and count how many in each major.
458  range.resize(numMinor);
459  for (auto& r : range)
460  r.second = 0;
461 
462  for (Index i = 0; i < Index(transposedData.range.size()); ++i)
463  {
464  for (Index j = transposedData.range[i].first; j < transposedData.range[i].second; ++j)
465  {
466  ++range[transposedData.entries[j].first].second;
467  }
468  }
469 
470  // Initialize start of each major. This will be incremented when copying.
471  range[0].first = 0;
472  for (Index i = 1; i < Index(range.size()); ++i)
473  {
474  range[i].first = range[i-1].first + range[i-1].second;
475  range[i-1].second = range[i].first;
476  }
477  range.back().second = transposedData.entries.size();
478 
479  // Copy entries major-wise into slots indicated by first[i], incrementing the latter.
480  entries.resize(transposedData.entries.size());
481  for (Index transposedMajor = 0; transposedMajor < Index(transposedData.range.size());
482  ++transposedMajor)
483  {
484  for (Index j = transposedData.range[transposedMajor].first;
485  j < transposedData.range[transposedMajor].second; ++j)
486  {
487  auto transposedEntry = transposedData.entries[j];
488  entries[range[transposedEntry.first].first] = std::make_pair(transposedMajor,
489  transposedEntry.second);
490  ++range[transposedEntry.first].first;
491  }
492  }
493 
494  // Re-initialize start of each major.
495  range[0].first = 0;
496  for (Index i = 1; i < Index(range.size()); ++i)
497  range[i].first = range[i-1].second;
498 
499  isSorted = true;
500  }
501 
506  void ensureConsistency(Index numMinor = std::numeric_limits<Index>::max()) const
507  {
508  for (Index i = 0; i < Index(range.size()); ++i)
509  {
510  if (range[i].first > entries.size())
511  throw std::runtime_error(
512  "Inconsistent SparseMatrix::Data: range.first contains index out of range.");
513  if (range[i].second > entries.size())
514  throw std::runtime_error(
515  "Inconsistent SparseMatrix::Data: range.second contains index out of range.");
516  if (range[i].first > range[i].second)
517  throw std::runtime_error(
518  "Inconsistent SparseMatrix::Data: range.first > range.beyond.");
519  }
520  for (Index i = 0; i < Index(entries.size()); ++i)
521  {
522  if (entries[i].first >= numMinor)
523  throw std::runtime_error("Inconsistent SparseMatrix::Data: minor entry exceeds bound.");
524  if (entries[i].second == 0)
525  throw std::runtime_error("Inconsistent SparseMatrix::Data: zero entry found.");
526  }
527  for (Index i = 0; i < Index(range.size()); ++i)
528  {
529  for (Index j = range[i].first + 1; j < range[i].second; ++j)
530  {
531  if (entries[j-1].first == entries[j].first)
532  {
533  std::stringstream ss;
534  ss << "Inconsistent SparseMatrix::Data: minor entries of major " << i
535  << " contain duplicate indices " << (j-1) << "->" << entries[j-1].first << " and "
536  << j << "->" << entries[j].first << ".";
537  throw std::runtime_error(ss.str());
538  }
539  else if (entries[j-1].first > entries[j].first && isSorted)
540  {
541  std::stringstream ss;
542  ss << "Inconsistent SparseMatrix::Data: minor entries of major " << i
543  << " are not sorted for indices " << (j-1) << "->" << entries[j-1].first << " and "
544  << j << "->" << entries[j].first << ".";
545  throw std::runtime_error(ss.str());
546  }
547  }
548  }
549  }
550 
555  void swap(Data& other)
556  {
557  range.swap(other.range);
558  entries.swap(other.entries);
559  std::swap(isSorted, other.isSorted);
560  }
561 
569  std::size_t find(Index major, Index minor) const
570  {
571  if (isSorted)
572  {
573  // Binary search on interval [lower, upper)
574 
575  std::size_t lower = range[major].first;
576  std::size_t upper = range[major].second;
577  std::size_t mid;
578  while (lower < upper)
579  {
580  mid = (lower + upper) / 2;
581  if (entries[mid].first < minor)
582  lower = mid + 1;
583  else if (entries[mid].first > minor)
584  upper = mid;
585  else
586  return mid;
587  }
588  }
589  else
590  {
591  // Linear scan.
592 
593  for (std::size_t i = range[major].first; i < range[major].second; ++i)
594  {
595  if (entries[i].first == minor)
596  return i;
597  }
598  }
599 
600  return std::numeric_limits<std::size_t>::max();
601  }
602 
607  void sort()
608  {
609  if (!isSorted)
610  return;
611 
612  for (std::size_t i = 0; i < entries.size(); ++i)
613  {
614  auto compare = [] (const std::pair<Index, Value>& a,
615  const std::pair<Index, Value>& b)
616  {
617  return a.first < b.first;
618  };
619 
620  std::sort(entries.begin() + range[i].first, entries.begin() + range[i].second, compare);
621  }
622 
623  isSorted = true;
624  }
625  };
626 
631  template <bool Row>
633  {
634  public:
635  NonzeroIterator(const SparseMatrix<V>& matrix, Index major)
636  : _data(Row ? matrix._rowData : matrix._columnData), _major(major),
637  _index(_data.range[major].first)
638  {
639 
640  }
641 
642  NonzeroIterator(const SparseMatrix<V>& matrix, Index major, int dummy)
643  : _data(Row ? matrix._rowData : matrix._columnData), _major(major),
644  _index(_data.range[major].second)
645  {
646  CMR_UNUSED(dummy);
647  }
648 
654  : _data(other._data), _major(other._major), _index(other._index)
655  {
656 
657  }
658 
663  inline const Nonzero operator*() const
664  {
665  // Statically known, so one branch will be optimized away.
666  if (Row)
667  return Nonzero(_major, _data.entries[_index].first, _data.entries[_index].second);
668  else
669  return Nonzero(_data.entries[_index].first, _major, _data.entries[_index].second);
670  }
671 
676  inline void operator++()
677  {
678  ++_index;
679  assert(_index <= _data.range[_major].second);
680  }
681 
686  inline bool operator!=(const NonzeroIterator<Row>& other) const
687  {
688  assert(&_data == &other._data);
689  return _major != other._major || _index != other._index;
690  }
691 
692  private:
694  const Data& _data;
696  Index _major;
698  Value _index;
699  };
700 
703 
709  : _zero(0)
710  {
711 
712  }
713 
719  : _rowData(std::move(other._rowData)), _columnData(std::move(other._columnData)), _zero(0)
720  {
721 
722  }
723 
729  {
730  _rowData = std::move(other._rowData);
731  _columnData = std::move(other._columnData);
732  return *this;
733  }
734 
754  SparseMatrix(bool majorRow, Index numMajor, Index numMinor, Index numEntries, Index* first,
755  Index* beyond, Index* entryMinors, Value* entryValues, bool mayContainZeros = true,
756  bool isSorted = false)
757  : _zero(0)
758  {
759  set(majorRow, numMajor, numMinor, numEntries, first, beyond, entryMinors, entryValues,
760  mayContainZeros, isSorted);
761  }
762 
773  SparseMatrix(bool majorRow, const Data& data, Index numMinor, bool majorIsSorted = false)
774  : _zero(0)
775  {
776  set(majorRow, data, numMinor, majorIsSorted);
777  }
778 
788  SparseMatrix(bool majorRow, Data&& data, Index numMinor)
789  : _zero(0)
790  {
791  set(majorRow, data, numMinor);
792  }
793 
803  SparseMatrix(const Data& rowData, const Data& columnData)
804  : _zero(0)
805  {
806  set(rowData, columnData);
807  }
808 
818  SparseMatrix(Data&& rowData, Data&& columnData)
819  : _zero(0)
820  {
821  set(std::move(rowData), std::move(columnData));
822  }
823 
843  void set(bool majorRow, Index numMajor, Index numMinor, Index numEntries, const Index* first,
844  const Index* beyond, const Index* entryMinors, const Value* entryValues,
845  bool mayContainZeros = true, bool isSorted = false)
846  {
847  if (mayContainZeros)
848  {
849  // Use first run over entries to count nonzeros in each row (resp. column).
850  _rowData.range.resize(numMajor);
851  for (auto& range : _rowData.range)
852  range.first = 0;
853  Index current = 0;
854  for (Index i = 0; i < numMajor; ++i)
855  {
856  _rowData.range[i].first = current;
857  Index b;
858  if (beyond != nullptr)
859  b = beyond[i];
860  else if (i + 1 < numMajor)
861  b = first[i+1];
862  else
863  b = numEntries;
864  for (Index j = first[i]; j < b; ++j)
865  {
866  if (entryValues[j] != 0)
867  ++current;
868  }
869  _rowData.range[i].second = current;
870  }
871 
872  // Use second run to copy the nonzeros.
873  _rowData.entries.resize(current);
874  current = 0;
875  for (Index i = 0; i < numMajor; ++i)
876  {
877  Index b;
878  if (beyond != nullptr)
879  b = beyond[i];
880  else if (i + 1 < numMajor)
881  b = first[i+1];
882  else
883  b = numEntries;
884  for (Index j = first[i]; j < b; ++j)
885  {
886  if (entryValues[j] != 0)
887  {
888  _rowData.entries[current] = std::make_pair(entryMinors[j], entryValues[j]);
889  ++current;
890  }
891  }
892  }
893  }
894  else
895  {
896  // If we can trust the input then we can copy directly.
897  _rowData.range.resize(numMajor);
898  if (beyond != nullptr)
899  {
900  for (Index i = 0; i < numMajor; ++i)
901  _rowData.range[i] = std::make_pair(first[i], beyond[i]);
902  }
903  else
904  {
905  for (Index i = 0; i + 1 < numMajor; ++i)
906  _rowData.range[i] = std::make_pair(first[i], first[i+1]);
907  _rowData.range.back() = std::make_pair(first[numMajor - 1], Index(numEntries));
908  }
909  _rowData.entries.resize(numEntries);
910  for (Index i = 0; i < numEntries; ++i)
911  _rowData.entries[i] = std::make_pair(entryMinors[i], entryValues[i]);
912  }
913 
914  _rowData.isSorted = isSorted;
915 
916  // Construct column-wise (resp. row-wise) representation.
917  _columnData.constructFromTransposed(_rowData, numMinor);
918 
919  if (!majorRow)
920  _rowData.swap(_columnData);
921 
922  #if !defined(NDEBUG)
924  #endif /* !NDEBUG */
925  }
926 
937  void set(bool majorRow, const Data& data, Index numMinor)
938  {
939  // We can trust the input and thus we can copy directly.
940  _rowData.range = data.range;
941  _rowData.entries = data.entries;
942  _rowData.isSorted = data.isSorted;
943 
944  // Construct column-wise (resp. row-wise) representation.
945  _columnData.constructFromTransposed(_rowData, numMinor);
946 
947  if (!majorRow)
948  _rowData.swap(_columnData);
949 
950  #if !defined(NDEBUG)
952  #endif /* !NDEBUG */
953  }
954 
965  void set(bool majorRow, Data&& data, Index numMinor)
966  {
967  _rowData = std::move(data);
968  _columnData.constructFromTransposed(_rowData, numMinor);
969  if (!majorRow)
970  _rowData.swap(_columnData);
971 
972 #if !defined(NDEBUG)
974 #endif /* !NDEBUG */
975  }
976 
986  void set(const Data& rowData, const Data& columnData)
987  {
988  _rowData.range = rowData.range;
989  _rowData.entries = rowData.entries;
990  _rowData.isSorted = rowData.isSorted;
991  _rowData.ensureConsistency();
992 
993  _columnData.range = columnData.range;
994  _columnData.entries = columnData.entries;
995  _columnData.isSorted = columnData.isSorted;
996  _columnData.ensureConsistency();
997 
998 #if !defined(NDEBUG)
1000 #endif /* !NDEBUG */
1001  }
1002 
1012  void set(Data&& rowData, Data&& columnData)
1013  {
1014  _rowData = std::move(rowData);
1015  _columnData = std::move(columnData);
1016 #if !defined(NDEBUG)
1018 #endif /* !NDEBUG */
1019  }
1020 
1025  ~SparseMatrix() = default;
1026 
1031  std::size_t numRows() const
1032  {
1033  return _rowData.range.size();
1034  }
1035 
1040  std::size_t numColumns() const
1041  {
1042  return _columnData.range.size();
1043  }
1044 
1049  bool hasSortedRows() const
1050  {
1051  return _rowData.isSorted;
1052  }
1053 
1058  bool hasSortedColumns() const
1059  {
1060  return _columnData.isSorted;
1061  }
1062 
1070  const Value& get(Index row, Index column) const
1071  {
1072  assert(row < numRows());
1073  assert(column < numColumns());
1074 
1075  if (_rowData.isSorted)
1076  {
1077  std::size_t columnIndex = _rowData.find(row, column);
1078  if (columnIndex < std::numeric_limits<std::size_t>::max())
1079  return _rowData.entries[columnIndex].second;
1080  }
1081  else
1082  {
1083  std::size_t rowIndex = _columnData.find(column, row);
1084  if (rowIndex < std::numeric_limits<std::size_t>::max())
1085  return _columnData.entries[rowIndex].second;
1086  }
1087 
1088  return _zero;
1089  }
1090 
1095  std::size_t numNonzeros() const
1096  {
1097  return _rowData.entries.size();
1098  }
1099 
1104  std::size_t countRowNonzeros(Index row) const
1105  {
1106  return _rowData.range[row].second - _rowData.range[row].first;
1107  }
1108 
1113  std::size_t countColumnNonzeros(Index column) const
1114  {
1115  return _columnData.range[column].second - _columnData.range[column].first;
1116  }
1117 
1123  {
1124  return NonzeroRowRange(NonzeroIterator<true>(*this, row),
1125  NonzeroIterator<true>(*this, row, 0));
1126  }
1127 
1133  {
1134  return NonzeroColumnRange(NonzeroIterator<false>(*this, column),
1135  NonzeroIterator<false>(*this, column, 0));
1136  }
1137 
1143  {
1144  return SparseMatrix(_columnData, _rowData);
1145  }
1146 
1154  void ensureConsistency() const
1155  {
1156  // First ensure individual consistency of row and column data.
1157  _rowData.ensureConsistency();
1158  _columnData.ensureConsistency();
1159 
1160  // Copy of a nonzero.
1161  struct NZ
1162  {
1163  Index row;
1164  Index column;
1165  Value value;
1166  };
1167 
1168  // Comparison for entries.
1169  auto compare = [](const NZ& a, const NZ& b)
1170  {
1171  if (a.row != b.row)
1172  return a.row < b.row;
1173  if (a.column != b.column)
1174  return a.column < b.column;
1175  return a.value < b.value;
1176  };
1177 
1178  // Construct sorted vector of entries from row data.
1179  std::vector<NZ> rowNonzeros;
1180  for (Index i = 0; i < Index(_rowData.range.size()); ++i)
1181  {
1182  for (Index j = _rowData.range[i].first; j < _rowData.range[i].second; ++j)
1183  {
1184  NZ nz = { i, _rowData.entries[j].first, _rowData.entries[j].second };
1185  rowNonzeros.push_back(nz);
1186  }
1187  }
1188  std::sort(rowNonzeros.begin(), rowNonzeros.end(), compare);
1189 
1190  // Construct sorted vector of entries from column data.
1191  std::vector<NZ> columnNonzeros;
1192  for (Index i = 0; i < Index(_columnData.range.size()); ++i)
1193  {
1194  for (Index j = _columnData.range[i].first; j < _columnData.range[i].second; ++j)
1195  {
1196  NZ nz = { _columnData.entries[j].first, i, _columnData.entries[j].second };
1197  columnNonzeros.push_back(nz);
1198  }
1199  }
1200  std::sort(columnNonzeros.begin(), columnNonzeros.end(), compare);
1201 
1202  for (std::size_t i = 0; i < rowNonzeros.size(); ++i)
1203  {
1204  if (rowNonzeros[i].row != columnNonzeros[i].row
1205  || rowNonzeros[i].column != columnNonzeros[i].column
1206  || rowNonzeros[i].value != columnNonzeros[i].value)
1207  {
1208  throw std::runtime_error("Inconsistent Matrix: row and column data differ.");
1209  }
1210  }
1211  }
1212 
1219  void transpose()
1220  {
1221  _rowData.swap(_columnData);
1222  }
1223 
1224 
1229  void sortRows()
1230  {
1231  _rowData.sort();
1232  }
1233 
1239  {
1240  _columnData.sort();
1241  }
1242 
1247  void sort()
1248  {
1249  sortRows();
1250  sortColumns();
1251  }
1252 
1253  private:
1255  Data _rowData;
1257  Data _columnData;
1259  Value _zero;
1261  bool _isSorted;
1262  };
1263 
1271  template <typename M>
1272  class TransposedMatrix : public Matrix<typename M::Value>
1273  {
1274  public:
1275  typedef typename M::Value Value;
1276  typedef typename M::Index Index;
1277  typedef typename M::Nonzero Nonzero;
1278 
1279  template <bool Row>
1281  {
1282  public:
1284  : _matrix(matrix._matrix), _iterator(typename M::template NonzeroIterator<!Row>(matrix._matrix, major))
1285  {
1286 
1287  }
1288 
1289  NonzeroIterator(const TransposedMatrix<M>& matrix, Index major, int dummy)
1290  : _matrix(matrix._matrix), _iterator(typename M::template NonzeroIterator<!Row>(matrix._matrix, major, 0))
1291  {
1292  CMR_UNUSED(dummy);
1293  }
1294 
1296  {
1297  ++_iterator;
1298  }
1299 
1301  {
1302  Nonzero nz = *_iterator;
1303  std::swap(nz.row, nz.column);
1304  return nz;
1305  }
1306 
1307  bool operator!=(const NonzeroIterator<Row>& other) const
1308  {
1309  assert(&_matrix == &other._matrix);
1310  return _iterator != other._iterator;
1311  }
1312 
1313  private:
1314  const M& _matrix;
1315  typename M::template NonzeroIterator<!Row> _iterator;
1316  };
1317 
1320 
1325  TransposedMatrix(M& matrix)
1326  : _matrix(matrix)
1327  {
1328 
1329  }
1330 
1336  : _matrix(other._matrix)
1337  {
1338 
1339  }
1340 
1345  std::size_t numRows() const
1346  {
1347  return _matrix.numColumns();
1348  }
1349 
1354  std::size_t numColumns() const
1355  {
1356  return _matrix.numRows();
1357  }
1358 
1363  bool hasSortedRows() const
1364  {
1365  return _matrix.hasSortedColumns();
1366  }
1367 
1372  bool hasSortedColumns() const
1373  {
1374  return _matrix.hasSortedRows();
1375  }
1376 
1384  const Value& get(Index row, Index column) const
1385  {
1386  assert(row < numRows());
1387  assert(column < numColumns());
1388 
1389  return _matrix.get(column, row);
1390  }
1391 
1396  std::size_t numNonzeros() const
1397  {
1398  return _matrix.numNonzeros();
1399  }
1400 
1405  std::size_t countRowNonzeros(Index row) const
1406  {
1407  return _matrix.countColumnNonzeros(row);
1408  }
1409 
1414  std::size_t countColumnNonzeros(Index column) const
1415  {
1416  return _matrix.countRowNonzeros(column);
1417  }
1418 
1424  {
1426  NonzeroIterator<true>(_matrix, row, 0));
1427  }
1428 
1434  {
1436  NonzeroIterator<false>(_matrix, column, 0));
1437  }
1438 
1443  void ensureConsistency() const
1444  {
1445  _matrix.ensureConsistency();
1446  }
1447 
1448  protected:
1450  };
1451 
1456  template <typename M>
1458  {
1459  return TransposedMatrix<M>(matrix);
1460  }
1461 
1467  {
1468  public:
1473  CMR_EXPORT
1474  SubmatrixIndices(const SubmatrixIndices& other);
1475 
1480  CMR_EXPORT
1482 
1487  CMR_EXPORT
1488  SubmatrixIndices(const std::vector<std::size_t>& rows, const std::vector<std::size_t>& columns);
1489 
1496  CMR_EXPORT
1497  SubmatrixIndices(std::size_t row, std::size_t column);
1498 
1503  inline
1504  std::size_t numRows() const
1505  {
1506  return _rows.size();
1507  }
1508 
1513  inline
1514  std::size_t numColumns() const
1515  {
1516  return _columns.size();
1517  }
1518 
1523  inline
1524  const std::vector<std::size_t>& rows() const
1525  {
1526  return _rows;
1527  }
1528 
1533  inline
1534  const std::vector<std::size_t>& columns() const
1535  {
1536  return _columns;
1537  }
1538 
1539  private:
1541  std::vector<std::size_t> _rows;
1542 
1544  std::vector<std::size_t> _columns;
1545  };
1546 
1547  template <typename M>
1548  class Submatrix : public Matrix<typename M::Value>
1549  {
1550  public:
1551  typedef typename M::Index Index;
1552  typedef typename M::Value Value;
1553  typedef typename M::Nonzero Nonzero;
1554 
1559  Submatrix(const M& matrix, const SubmatrixIndices& indices)
1560  : _matrix(matrix), _indices(indices)
1561  {
1562 
1563  }
1564 
1570  : _matrix(other._matrix), _indices(other._indices)
1571  {
1572 
1573  }
1574 
1579  std::size_t numRows() const
1580  {
1581  return _indices.numRows();
1582  }
1583 
1588  std::size_t numColumns() const
1589  {
1590  return _indices.numColumns();
1591  }
1592 
1597  bool hasSortedRows() const
1598  {
1599  return false;
1600  }
1601 
1606  bool hasSortedColumns() const
1607  {
1608  return false;
1609  }
1610 
1618  const Value& get(Index row, Index column) const
1619  {
1620  assert(row < numRows());
1621  assert(column < numColumns());
1622 
1623  return _matrix.get(_indices.rows()[row], _indices.columns()[column]);
1624  }
1625 
1630  void ensureConsistency() const
1631  {
1632  for (auto row : _indices.rows())
1633  {
1634  if (row >= _matrix.numRows())
1635  throw std::runtime_error("Inconsistent Submatrix: row index too large.");
1636  }
1637  for (auto column : _indices.columns())
1638  {
1639  if (column >= _matrix.numColumns())
1640  throw std::runtime_error("Inconsistent Submatrix: column index too large.");
1641  }
1642  _matrix.ensureConsistency();
1643  }
1644 
1645  private:
1646  const M& _matrix;
1647  const SubmatrixIndices& _indices;
1648  };
1649 
1650 
1651  template <typename Matrix>
1652  bool find_smallest_nonzero_matrix_entry(const Matrix& matrix, size_t row_first, size_t row_beyond, size_t column_first, size_t column_beyond,
1653  size_t& row, size_t& column)
1654  {
1655  bool result = false;
1656  int current_value = 0;
1657  for (size_t r = row_first; r != row_beyond; ++r)
1658  {
1659  for (size_t c = column_first; c != column_beyond; ++c)
1660  {
1661  int value = matrix(r, c);
1662  if (value == 0)
1663  continue;
1664 
1665  value = value >= 0 ? value : -value;
1666 
1667  if (!result || value < current_value)
1668  {
1669  result = true;
1670  row = r;
1671  column = c;
1672  current_value = value;
1673  }
1674  }
1675  }
1676  return result;
1677  }
1678 
1679  template <typename Matrix>
1680  bool matrix_row_zero(const Matrix& matrix, size_t row, size_t column_first, size_t column_beyond)
1681  {
1682  for (size_t c = column_first; c != column_beyond; ++c)
1683  if (matrix(row, c) != 0)
1684  return false;
1685  return true;
1686  }
1687 
1688  template <typename Matrix>
1689  bool matrix_column_zero(const Matrix& matrix, size_t column, size_t row_first, size_t row_beyond)
1690  {
1691  return matrix_row_zero(make_transposed_matrix(matrix), column, row_first, row_beyond);
1692  }
1693 
1694 
1695  template <typename Matrix1, typename Matrix2>
1696  bool equals(const Matrix1& first, const Matrix2& second)
1697  {
1698  if (first.size1() != second.size1())
1699  return false;
1700  if (first.size2() != second.size2())
1701  return false;
1702  for (size_t r = 0; r < first.size1(); ++r)
1703  {
1704  for (size_t c = 0; c < first.size2(); ++c)
1705  {
1706  if (first(r, c) != second(r, c))
1707  return false;
1708  }
1709  }
1710  return true;
1711  }
1712 
1713 
1722  template <typename MatrixType>
1723  inline void matrix_permute1(MatrixType& matrix, size_t index1, size_t index2)
1724  {
1725  for (size_t index = 0; index < matrix.size2(); ++index)
1726  {
1727  std::swap(matrix(index1, index), matrix(index2, index));
1728  }
1729  }
1730 
1739  template <typename MatrixType>
1740  inline void matrix_permute2(MatrixType& matrix, size_t index1, size_t index2)
1741  {
1742  for (size_t index = 0; index < matrix.size1(); ++index)
1743  {
1744  std::swap(matrix(index, index1), matrix(index, index2));
1745  }
1746  }
1747 
1748 
1749 } /* namespace tu */
Definition: matrix.hpp:88
NonzeroIterator(const DenseMatrix< V > &matrix, Index major, int dummy)
Definition: matrix.hpp:105
bool operator!=(const NonzeroIterator< Row > &other) const
Definition: matrix.hpp:147
NonzeroIterator(const DenseMatrix< V > &matrix, Index major)
Definition: matrix.hpp:90
NonzeroIterator< Row > & operator++()
Definition: matrix.hpp:116
Nonzero operator*() const
Definition: matrix.hpp:131
Dense matrix with entries of type V.
Definition: matrix.hpp:80
DenseMatrix(Index numRows, Index numColumns)
Constructs zero matrix of given size.
Definition: matrix.hpp:243
std::size_t countRowNonzeros(Index row) const
Returns the number of nonzeros in row.
Definition: matrix.hpp:297
DenseMatrix< V > & operator=(const DenseMatrix< V2 > &other)
Assignment operator.
Definition: matrix.hpp:218
DenseMatrix()
Default constructor for 0x0 matrix.
Definition: matrix.hpp:166
detail::Range< NonzeroIterator< false > > NonzeroColumnRange
Definition: matrix.hpp:160
Matrix< V >::Nonzero Nonzero
Definition: matrix.hpp:84
Index _numColumns
Number of rows.
Definition: matrix.hpp:359
DenseMatrix(const DenseMatrix< V2 > &other)
Copy constructor.
Definition: matrix.hpp:189
Index _numRows
Matrix entries with row-major ordering.
Definition: matrix.hpp:358
Index numColumns() const
Returns the number of columns.
Definition: matrix.hpp:265
Index numRows() const
Returns the number of rows.
Definition: matrix.hpp:256
detail::Range< NonzeroIterator< true > > NonzeroRowRange
Definition: matrix.hpp:159
void ensureConsistency() const
Checks for consistency, raising a std::runtime_error if inconsistent.
Definition: matrix.hpp:351
std::size_t countColumnNonzeros(Index column) const
Returns the number of nonzeros in column.
Definition: matrix.hpp:314
Matrix< V >::Value Value
Definition: matrix.hpp:82
Value * _data
Definition: matrix.hpp:357
~DenseMatrix()
Destructor.
Definition: matrix.hpp:233
NonzeroRowRange iterateRowNonzeros(Index row) const
Returns a range for iterating over the nonzeros of row.
Definition: matrix.hpp:331
DenseMatrix(DenseMatrix< V > &&other)
Move constructor that steals the data from other.
Definition: matrix.hpp:176
const Value & get(Index row, Index column) const
Returns entry at row, column.
Definition: matrix.hpp:274
Matrix< V >::Index Index
Definition: matrix.hpp:83
NonzeroColumnRange iterateColumnNonzeros(Index column) const
Returns a range for iterating over the nonzeros of column.
Definition: matrix.hpp:341
void set(Index row, Index column, const V2 &value)
Sets entry at row, column to copy of value.
Definition: matrix.hpp:286
DenseMatrix< V > & operator=(DenseMatrix< V > &&other)
Move assignment operator that steals the data from other.
Definition: matrix.hpp:202
Base class for matrices whose entries have type V.
Definition: matrix.hpp:55
V Value
Definition: matrix.hpp:57
std::size_t Index
Type of entries.
Definition: matrix.hpp:58
Sparse matrix with entries of type V.
Definition: matrix.hpp:368
SparseMatrix & operator=(SparseMatrix &&other)
Move assignment operator.
Definition: matrix.hpp:728
detail::Range< NonzeroIterator< true > > NonzeroRowRange
Definition: matrix.hpp:701
Matrix< V >::Index Index
Definition: matrix.hpp:371
std::size_t numRows() const
Returns the number of rows.
Definition: matrix.hpp:1031
std::size_t numColumns() const
Returns the number of columns.
Definition: matrix.hpp:1040
void ensureConsistency() const
Checks for consistency, raising a std::runtime_error if inconsistent.
Definition: matrix.hpp:1154
bool hasSortedRows() const
Indicates if the row data is sorted by column.
Definition: matrix.hpp:1049
void sortRows()
Sort row data (by column).
Definition: matrix.hpp:1229
~SparseMatrix()=default
Destructor.
SparseMatrix(Data &&rowData, Data &&columnData)
Constructs a matrix.
Definition: matrix.hpp:818
bool hasSortedColumns() const
Indicates if the column data is sorted by row.
Definition: matrix.hpp:1058
SparseMatrix(bool majorRow, Data &&data, Index numMinor)
Constructs a matrix.
Definition: matrix.hpp:788
std::size_t numNonzeros() const
Returns the number of nonzeros.
Definition: matrix.hpp:1095
void set(bool majorRow, Data &&data, Index numMinor)
Sets the contents of the matrix.
Definition: matrix.hpp:965
Matrix< V >::Value Value
Definition: matrix.hpp:370
NonzeroRowRange iterateRowNonzeros(Index row) const
Returns a range for iterating over the nonzeros of row.
Definition: matrix.hpp:1122
std::size_t countColumnNonzeros(Index column) const
Returns the number of nonzeros in column.
Definition: matrix.hpp:1113
const Value & get(Index row, Index column) const
Returns entry at row, column.
Definition: matrix.hpp:1070
void sort()
Sort row and column data (by column and row, respectively).
Definition: matrix.hpp:1247
void set(const Data &rowData, const Data &columnData)
Sets the contents of the matrix.
Definition: matrix.hpp:986
Matrix< V >::Nonzero Nonzero
Definition: matrix.hpp:372
SparseMatrix(SparseMatrix &&other)
Move constructor.
Definition: matrix.hpp:718
SparseMatrix transposed() const
Returns a copy of the matrix' transpose.
Definition: matrix.hpp:1142
detail::Range< NonzeroIterator< false > > NonzeroColumnRange
Definition: matrix.hpp:702
SparseMatrix(const Data &rowData, const Data &columnData)
Constructs a matrix.
Definition: matrix.hpp:803
NonzeroColumnRange iterateColumnNonzeros(Index column) const
Returns a range for iterating over the nonzeros of column.
Definition: matrix.hpp:1132
void set(bool majorRow, const Data &data, Index numMinor)
Sets the contents of the matrix.
Definition: matrix.hpp:937
SparseMatrix(bool majorRow, const Data &data, Index numMinor, bool majorIsSorted=false)
Constructs a matrix.
Definition: matrix.hpp:773
void set(Data &&rowData, Data &&columnData)
Sets the contents of the matrix.
Definition: matrix.hpp:1012
void set(bool majorRow, Index numMajor, Index numMinor, Index numEntries, const Index *first, const Index *beyond, const Index *entryMinors, const Value *entryValues, bool mayContainZeros=true, bool isSorted=false)
Sets the contents of the matrix.
Definition: matrix.hpp:843
void sortColumns()
Sort column data (by row).
Definition: matrix.hpp:1238
SparseMatrix()
Constructs a 0x0 matrix.
Definition: matrix.hpp:708
void transpose()
Transposes the matrix.
Definition: matrix.hpp:1219
std::size_t countRowNonzeros(Index row) const
Returns the number of nonzeros in row.
Definition: matrix.hpp:1104
SparseMatrix(bool majorRow, Index numMajor, Index numMinor, Index numEntries, Index *first, Index *beyond, Index *entryMinors, Value *entryValues, bool mayContainZeros=true, bool isSorted=false)
Constructs a matrix.
Definition: matrix.hpp:754
Indexes a submatrix.
Definition: matrix.hpp:1467
CMR_EXPORT SubmatrixIndices(std::size_t row, std::size_t column)
Constructor for a 1x1 submatrix.
const std::vector< std::size_t > & columns() const
Returns the vector of column indices.
Definition: matrix.hpp:1534
CMR_EXPORT SubmatrixIndices(const SubmatrixIndices &other)
Copy constructor.
Definition: matrix.cpp:7
const std::vector< std::size_t > & rows() const
Returns the vector of row indices.
Definition: matrix.hpp:1524
std::size_t numRows() const
Returns the number of row indices.
Definition: matrix.hpp:1504
std::size_t numColumns() const
Returns the number of column indices.
Definition: matrix.hpp:1514
Definition: matrix.hpp:1549
M::Value Value
Definition: matrix.hpp:1552
bool hasSortedRows() const
Indicates if the row data is sorted by column.
Definition: matrix.hpp:1597
M::Index Index
Definition: matrix.hpp:1551
Submatrix(Submatrix &&other)
Move constructor.
Definition: matrix.hpp:1569
bool hasSortedColumns() const
Indicates if the column data is sorted by row.
Definition: matrix.hpp:1606
M::Nonzero Nonzero
Definition: matrix.hpp:1553
std::size_t numRows() const
Returns the number of rows.
Definition: matrix.hpp:1579
std::size_t numColumns() const
Returns the number of columns.
Definition: matrix.hpp:1588
Submatrix(const M &matrix, const SubmatrixIndices &indices)
Constructs submatrix of matrix indexed by indices.
Definition: matrix.hpp:1559
void ensureConsistency() const
Checks for consistency, raising a std::runtime_error if inconsistent.
Definition: matrix.hpp:1630
const Value & get(Index row, Index column) const
Returns entry at row, column.
Definition: matrix.hpp:1618
Definition: matrix.hpp:1281
Nonzero operator*() const
Definition: matrix.hpp:1300
NonzeroIterator< Row > & operator++()
Definition: matrix.hpp:1295
bool operator!=(const NonzeroIterator< Row > &other) const
Definition: matrix.hpp:1307
NonzeroIterator(const TransposedMatrix< M > &matrix, Index major)
Definition: matrix.hpp:1283
NonzeroIterator(const TransposedMatrix< M > &matrix, Index major, int dummy)
Definition: matrix.hpp:1289
Matrix proxy for the transpose of a matrix.
Definition: matrix.hpp:1273
std::size_t numColumns() const
Returns the number of columns.
Definition: matrix.hpp:1354
bool hasSortedColumns() const
Indicates if the column data is sorted by row.
Definition: matrix.hpp:1372
M::Nonzero Nonzero
Definition: matrix.hpp:1277
std::size_t countRowNonzeros(Index row) const
Returns the number of nonzeros in row.
Definition: matrix.hpp:1405
bool hasSortedRows() const
Indicates if the row data is sorted by column.
Definition: matrix.hpp:1363
M::Value Value
Definition: matrix.hpp:1275
detail::Range< NonzeroIterator< false > > NonzeroColumnRange
Definition: matrix.hpp:1319
std::size_t countColumnNonzeros(Index column) const
Returns the number of nonzeros in column.
Definition: matrix.hpp:1414
void ensureConsistency() const
Checks for consistency, raising a std::runtime_error if inconsistent.
Definition: matrix.hpp:1443
M & _matrix
Definition: matrix.hpp:1449
TransposedMatrix(M &matrix)
Constructs from a matrix of type M.
Definition: matrix.hpp:1325
std::size_t numNonzeros() const
Returns the number of nonzeros.
Definition: matrix.hpp:1396
M::Index Index
Definition: matrix.hpp:1276
std::size_t numRows() const
Returns the number of rows.
Definition: matrix.hpp:1345
NonzeroRowRange iterateRowNonzeros(Index row) const
Returns a range for iterating over the nonzeros of row.
Definition: matrix.hpp:1423
const Value & get(Index row, Index column) const
Returns entry at row, column.
Definition: matrix.hpp:1384
NonzeroColumnRange iterateColumnNonzeros(Index column) const
Returns a range for iterating over the nonzeros of column.
Definition: matrix.hpp:1433
detail::Range< NonzeroIterator< true > > NonzeroRowRange
Definition: matrix.hpp:1318
TransposedMatrix(TransposedMatrix &&other)
Move constructor.
Definition: matrix.hpp:1335
Definition: matrix.hpp:18
Iterator begin()
Definition: matrix.hpp:32
Range(Iterator begin, Iterator end)
Definition: matrix.hpp:20
Range(const Range< Iterator > &other)
Definition: matrix.hpp:26
Iterator end()
Definition: matrix.hpp:37
#define CMR_UNUSED(x)
Definition: env.h:24
Definition: algorithm.hpp:14
void sort(permutation &permutation, size_t first, size_t beyond, Less &less)
Definition: permutations.hpp:583
bool equals(const Matrix1 &first, const Matrix2 &second)
Definition: matrix.hpp:1696
void matrix_permute1(MatrixType &matrix, size_t index1, size_t index2)
Definition: matrix.hpp:1723
bool find_smallest_nonzero_matrix_entry(const Matrix &matrix, size_t row_first, size_t row_beyond, size_t column_first, size_t column_beyond, size_t &row, size_t &column)
Definition: matrix.hpp:1652
TransposedMatrix< M > transposedMatrix(M &matrix)
Returns a TransposedMatrix for matrix.
Definition: matrix.hpp:1457
matrix_transposed< MatrixType > make_transposed_matrix(MatrixType &matrix)
Definition: matrix_transposed.hpp:148
bool matrix_column_zero(const Matrix &matrix, size_t column, size_t row_first, size_t row_beyond)
Definition: matrix.hpp:1689
bool matrix_row_zero(const Matrix &matrix, size_t row, size_t column_first, size_t column_beyond)
Definition: matrix.hpp:1680
void matrix_permute2(MatrixType &matrix, size_t index1, size_t index2)
Definition: matrix.hpp:1740
Row/column index type.
Definition: matrix.hpp:61
Index row
Definition: matrix.hpp:62
Index column
Definition: matrix.hpp:63
const Value & value
Definition: matrix.hpp:64
Nonzero(Index r, Index c, const Value &v)
Definition: matrix.hpp:66
Matrix data, row-wise or column-wise.
Definition: matrix.hpp:382
~Data()=default
Destructor.
std::vector< std::pair< Index, Value > > entries
Array that maps the nonzero indices to pairs of minor index and value.
Definition: matrix.hpp:386
void swap(Data &other)
Swaps data with other.
Definition: matrix.hpp:555
Data & operator=(const Data &other)
Assignment operator.
Definition: matrix.hpp:436
std::vector< std::pair< Index, Index > > range
Array that maps a major to its nonzero range specified by first and beyond entries.
Definition: matrix.hpp:384
std::size_t find(Index major, Index minor) const
Finds entry major,minor.
Definition: matrix.hpp:569
Data & operator=(Data &&other)
Move operator.
Definition: matrix.hpp:424
void sort()
Sort data by minor.
Definition: matrix.hpp:607
bool isSorted
Whether the data is sorted by minor.
Definition: matrix.hpp:388
Data(const Data &other)
Copy constructor.
Definition: matrix.hpp:414
Data()
Default constructor.
Definition: matrix.hpp:394
Data(Data &&other)
Move constructor.
Definition: matrix.hpp:404
void constructFromTransposed(const Data &transposedData, Index numMinor)
Computes the transpose version of data. numMinor is the new number of minor indices.
Definition: matrix.hpp:455
void ensureConsistency(Index numMinor=std::numeric_limits< Index >::max()) const
Checks for consistency, raising a std::runtime_error if inconsistent.
Definition: matrix.hpp:506
Iterator for row- or column-wise iteration over the entries.
Definition: matrix.hpp:633
const Nonzero operator*() const
Returns the current entry.
Definition: matrix.hpp:663
NonzeroIterator(const SparseMatrix< V > &matrix, Index major, int dummy)
Definition: matrix.hpp:642
bool operator!=(const NonzeroIterator< Row > &other) const
Compares two iterators for being not equal.
Definition: matrix.hpp:686
void operator++()
Advances to the next entry in the same major.
Definition: matrix.hpp:676
NonzeroIterator(const NonzeroIterator &other)
Copy constructor.
Definition: matrix.hpp:653
NonzeroIterator(const SparseMatrix< V > &matrix, Index major)
Definition: matrix.hpp:635