matrix.hpp

Go to the documentation of this file.

#pragma once
 
#include <cmr/config.h>
#include <cmr/export.h>
 
#include <boost/numeric/ublas/matrix.hpp>
 
#include "matrix_transposed.hpp"
#include <iomanip>
 
namespace tu
{
  namespace detail
  {
 
    template <typename Iterator>
    class Range
    {
    public:
      Range(Iterator begin, Iterator end)
        : _begin(begin), _end(end)
      {
 
      }
 
      Range(const Range<Iterator>& other)
        : _begin(other._begin), _end(other._end)
      {
 
      }
 
      Iterator begin()
      {
        return _begin;
      }
 
      Iterator end()
      {
        return _end;
      }
 
    private:
      Iterator _begin;
      Iterator _end;
    };
 
  } /* namespace tu::detail */
 
  template <typename V>
  class Matrix
  {
  public:
    typedef V Value; 
    typedef std::size_t Index; 
 
    struct Nonzero
    {
      Index row;
      Index column;
      const Value& value;
 
      Nonzero(Index r, Index c, const Value& v)
        : row(r), column(c), value(v)
      {
 
      }
    };
  };
 
  template <typename V>
  class DenseMatrix : public Matrix<V>
  {
  public:
    typedef typename Matrix<V>::Value Value;
    typedef typename Matrix<V>::Index Index;
    typedef typename Matrix<V>::Nonzero Nonzero;
 
    template <bool Row>
    class NonzeroIterator
    {
    public:
      NonzeroIterator(const DenseMatrix<V>& matrix, Index major)
        : _matrix(matrix), _major(major), _minor(0)
      {
        if (Row)
        {
          while (_minor < _matrix.numColumns() && _matrix.get(_major, _minor) == 0)
            ++_minor;
        }
        else
        {
          while (_minor < _matrix.numRows() && _matrix.get(_minor, _major) == 0)
            ++_minor;
        }
      }
 
      NonzeroIterator(const DenseMatrix<V>& matrix, Index major, int dummy)
        : _matrix(matrix), _major(major)
      {
        CMR_UNUSED(dummy);
 
        if (Row)
          _minor = _matrix.numColumns();
        else
          _minor = _matrix.numRows();
      }
 
      NonzeroIterator<Row>& operator++()
      {
        ++_minor;
        if (Row)
        {
          while (_minor < _matrix.numColumns() && _matrix.get(_major, _minor) == 0)
            ++_minor;
        }
        else
        {
          while (_minor < _matrix.numRows() && _matrix.get(_minor, _major) == 0)
            ++_minor;
        }
      }
 
      Nonzero operator*() const
      {
        if (Row)
        {
          assert(_major < _matrix.numRows());
          assert(_minor < _matrix.numColumns());
          return Nonzero(_major, _minor, _matrix.get(_major, _minor));
        }
        else
        {
          assert(_minor < _matrix.numRows());
          assert(_major < _matrix.numColumns());
          return Nonzero(_minor, _major, _matrix.get(_minor, _major));
        }
      }
 
      bool operator!=(const NonzeroIterator<Row>& other) const
      {
        assert(&_matrix == &other._matrix);
        assert(_major == other._major);
        return _minor != other._minor;
      }
 
    private:
      const DenseMatrix<V>& _matrix;
      std::size_t _major, _minor;
    };
 
    typedef detail::Range<NonzeroIterator<true>> NonzeroRowRange;
    typedef detail::Range<NonzeroIterator<false>> NonzeroColumnRange;
 
    DenseMatrix()
      : _data(nullptr), _numRows(0), _numColumns(0)
    {
 
    }
 
    DenseMatrix(DenseMatrix<V>&& other)
      : _data(other._data), _numRows(other._numRows), _numColumns(other._numColumns)
    {
      other._data = nullptr;
      other._numRows = 0;
      other._numColumns = 0;
    }
 
    template <typename V2>
    DenseMatrix(const DenseMatrix<V2>& other)
      : _numRows(other._numRows), _numColumns(other._numColumns)
    {
      std::size_t size = other._numRows * other._numColumns;
      _data = new Value[size];
      for (std::size_t i = 0; i < size; ++i)
        _data[i] = other._data[i];
    }
 
    DenseMatrix<V>& operator=(DenseMatrix<V>&& other)
    {
      _data = other._data;
      _numRows = other._numRows;
      _numColumns = other._numColumns;
      other._data = nullptr;
      other._numRows = 0;
      other._numColumns = 0;
      return *this;
    }
 
    template <typename V2>
    DenseMatrix<V>& operator=(const DenseMatrix<V2>& other)
    {
      _numRows = other._numRows;
      _numColumns = other._numColumns;
      std::size_t size = other._numRows * other._numColumns;
      _data = new Value[size];
      for (std::size_t i = 0; i < size; ++i)
        _data[i] = other._data[i];
      return *this;
    }
 
    ~DenseMatrix()
    {
      if (_data != nullptr)
        delete[] _data;
    }
 
    DenseMatrix(Index numRows, Index numColumns)
      : _numRows(numRows), _numColumns(numColumns)
    {
      std::size_t size = numRows * numColumns;
      _data = new Value[size];
      for (std::size_t i = 0; i < size; ++i)
        _data[i] = 0;
    }
 
    Index numRows() const
    {
      return _numRows;
    }
 
    Index numColumns() const
    {
      return _numColumns;
    }
 
    const Value& get(Index row, Index column) const
    {
      assert(row < _numRows);
      assert(column < _numColumns);
      return _data[_numColumns * row + column];
    }
 
    template <typename V2>
    void set(Index row, Index column, const V2& value)
    {
      assert(row < _numRows);
      assert(column < _numColumns);
      _data[_numColumns * row + column] = value;
    }
 
    std::size_t countRowNonzeros(Index row) const
    {
      std::size_t begin = _numColumns * row;
      std::size_t end = _numColumns * (row + 1);
      std::size_t count = 0;
      for (std::size_t i = begin; i < end; ++i)
      {
        if (_data[i] != 0)
          ++count;
      }
      return count;
    }
 
    std::size_t countColumnNonzeros(Index column) const
    {
      std::size_t begin = column;
      std::size_t end = _numRows * _numColumns;
      std::size_t count = 0;
      for (std::size_t i = begin; i < end; i += _numColumns)
      {
        if (_data[i] != 0)
          ++count;
      }
      return count;
    }
 
    NonzeroRowRange iterateRowNonzeros(Index row) const
    {
      return NonzeroRowRange(NonzeroIterator<true>(*this, row),
        NonzeroIterator<true>(*this, row, 0));
    }
 
    NonzeroColumnRange iterateColumnNonzeros(Index column) const
    {
      return NonzeroColumnRange(NonzeroIterator<false>(*this, column),
        NonzeroIterator<false>(*this, column, 0));
    }
 
    void ensureConsistency() const
    {
 
    }
 
  protected:
    Value* _data; 
    Index _numRows; 
    Index _numColumns; 
  };
 
  template<typename V>
  class SparseMatrix : public Matrix<V>
  {
  public:
    typedef typename Matrix<V>::Value Value;
    typedef typename Matrix<V>::Index Index;
    typedef typename Matrix<V>::Nonzero Nonzero;
 
    struct Data
    {
      std::vector<std::pair<Index, Index>> range;
      std::vector<std::pair<Index, Value>> entries;
      bool isSorted;
 
      Data()
        : isSorted(true)
      {
 
      }
 
      Data(Data&& other)
        : range(std::move(other.range)), entries(std::move(other.entries)), isSorted(other.isSorted)
      {
 
      }
 
      Data(const Data& other)
        : range(other.range), entries(other.entries), isSorted(other.isSorted)
      {
 
      }
 
      Data& operator=(Data&& other)
      {
        range = std::move(other.range);
        entries = std::move(other.entries);
        isSorted = other.isSorted;
        return *this;
      }
 
      Data& operator=(const Data& other)
      {
        range = other.range;
        entries = other.entries;
        isSorted = other.isSorted;
        return *this;
      }
 
      ~Data() = default;
 
      void constructFromTransposed(const Data& transposedData, Index numMinor)
      {
        // Iterate over entries and count how many in each major.
        range.resize(numMinor);
        for (auto& r : range)
          r.second = 0;
 
        for (Index i = 0; i < Index(transposedData.range.size()); ++i)
        {
          for (Index j = transposedData.range[i].first; j < transposedData.range[i].second; ++j)
          {
            ++range[transposedData.entries[j].first].second;
          }
        }
 
        // Initialize start of each major. This will be incremented when copying.
        range[0].first = 0;
        for (Index i = 1; i < Index(range.size()); ++i)
        {
          range[i].first = range[i-1].first + range[i-1].second;
          range[i-1].second = range[i].first;
        }
        range.back().second = transposedData.entries.size();
 
        // Copy entries major-wise into slots indicated by first[i], incrementing the latter.
        entries.resize(transposedData.entries.size());
        for (Index transposedMajor = 0; transposedMajor < Index(transposedData.range.size());
          ++transposedMajor)
        {
          for (Index j = transposedData.range[transposedMajor].first;
            j < transposedData.range[transposedMajor].second; ++j)
          {
            auto transposedEntry = transposedData.entries[j];
            entries[range[transposedEntry.first].first] = std::make_pair(transposedMajor,
              transposedEntry.second);
            ++range[transposedEntry.first].first;
          }
        }
 
        // Re-initialize start of each major.
        range[0].first = 0;
        for (Index i = 1; i < Index(range.size()); ++i)
          range[i].first = range[i-1].second;
 
        isSorted = true;
      }
 
      void ensureConsistency(Index numMinor = std::numeric_limits<Index>::max()) const
      {
        for (Index i = 0; i < Index(range.size()); ++i)
        {
          if (range[i].first > entries.size())
            throw std::runtime_error(
                "Inconsistent SparseMatrix::Data: range.first contains index out of range.");
          if (range[i].second > entries.size())
            throw std::runtime_error(
                "Inconsistent SparseMatrix::Data: range.second contains index out of range.");
          if (range[i].first > range[i].second)
            throw std::runtime_error(
                "Inconsistent SparseMatrix::Data: range.first > range.beyond.");
        }
        for (Index i = 0; i < Index(entries.size()); ++i)
        {
          if (entries[i].first >= numMinor)
            throw std::runtime_error("Inconsistent SparseMatrix::Data: minor entry exceeds bound.");
          if (entries[i].second == 0)
            throw std::runtime_error("Inconsistent SparseMatrix::Data: zero entry found.");
        }
        for (Index i = 0; i < Index(range.size()); ++i)
        {
          for (Index j = range[i].first + 1; j < range[i].second; ++j)
          {
            if (entries[j-1].first == entries[j].first)
            {
              std::stringstream ss;
              ss << "Inconsistent SparseMatrix::Data: minor entries of major " << i
                << " contain duplicate indices " << (j-1) << "->" << entries[j-1].first << " and "
                << j << "->" << entries[j].first << ".";
              throw std::runtime_error(ss.str());
            }
            else if (entries[j-1].first > entries[j].first && isSorted)
            {
              std::stringstream ss;
              ss << "Inconsistent SparseMatrix::Data: minor entries of major " << i
                << " are not sorted for indices " << (j-1) << "->" << entries[j-1].first << " and "
                << j << "->" << entries[j].first << ".";
              throw std::runtime_error(ss.str());
            }
          }
        }
      }
 
      void swap(Data& other)
      {
        range.swap(other.range);
        entries.swap(other.entries);
        std::swap(isSorted, other.isSorted);
      }
 
      std::size_t find(Index major, Index minor) const
      {
        if (isSorted)
        {
          // Binary search on interval [lower, upper)
 
          std::size_t lower = range[major].first;
          std::size_t upper = range[major].second;
          std::size_t mid;
          while (lower < upper)
          {
            mid = (lower + upper) / 2;
            if (entries[mid].first < minor)
              lower = mid + 1;
            else if (entries[mid].first > minor)
              upper = mid;
            else
              return mid;
          }
        }
        else
        {
          // Linear scan.
 
          for (std::size_t i = range[major].first; i < range[major].second; ++i)
          {
            if (entries[i].first == minor)
              return i;
          }
        }
 
        return std::numeric_limits<std::size_t>::max();
      }
 
      void sort()
      {
        if (!isSorted)
          return;
 
        for (std::size_t i = 0; i < entries.size(); ++i)
        {
          auto compare = [] (const std::pair<Index, Value>& a,
            const std::pair<Index, Value>& b)
            {
              return a.first < b.first;
            };
 
          std::sort(entries.begin() + range[i].first, entries.begin() + range[i].second, compare);
        }
 
        isSorted = true;
      }
    };
 
    template <bool Row>
    struct NonzeroIterator
    {
    public:
      NonzeroIterator(const SparseMatrix<V>& matrix, Index major)
        : _data(Row ? matrix._rowData : matrix._columnData), _major(major),
        _index(_data.range[major].first)
      {
 
      }
 
      NonzeroIterator(const SparseMatrix<V>& matrix, Index major, int dummy)
        : _data(Row ? matrix._rowData : matrix._columnData), _major(major),
        _index(_data.range[major].second)
      {
        CMR_UNUSED(dummy);
      }
 
      NonzeroIterator(const NonzeroIterator& other)
        : _data(other._data), _major(other._major), _index(other._index)
      {
 
      }
 
      inline const Nonzero operator*() const
      {
        // Statically known, so one branch will be optimized away.
        if (Row)
          return Nonzero(_major, _data.entries[_index].first, _data.entries[_index].second);
        else
          return Nonzero(_data.entries[_index].first, _major, _data.entries[_index].second);
      }
 
      inline void operator++()
      {
        ++_index;
        assert(_index <= _data.range[_major].second);
      }
 
      inline bool operator!=(const NonzeroIterator<Row>& other) const 
      {
        assert(&_data == &other._data);
        return _major != other._major || _index != other._index;
      }
 
      private:
        const Data& _data;
        Index _major;
        Value _index;
    };
 
    typedef detail::Range<NonzeroIterator<true>> NonzeroRowRange;
    typedef detail::Range<NonzeroIterator<false>> NonzeroColumnRange;
 
    SparseMatrix()
      : _zero(0)
    {
 
    }
 
    SparseMatrix(SparseMatrix&& other)
      : _rowData(std::move(other._rowData)), _columnData(std::move(other._columnData)), _zero(0)
    {
 
    }
 
    SparseMatrix& operator=(SparseMatrix&& other)
    {
      _rowData = std::move(other._rowData);
      _columnData = std::move(other._columnData);
      return *this;
    }
 
    SparseMatrix(bool majorRow, Index numMajor, Index numMinor, Index numEntries, Index* first,
      Index* beyond, Index* entryMinors, Value* entryValues, bool mayContainZeros = true,
      bool isSorted = false)
      : _zero(0)
    {
      set(majorRow, numMajor, numMinor, numEntries, first, beyond, entryMinors, entryValues,
        mayContainZeros, isSorted);
    }
 
    SparseMatrix(bool majorRow, const Data& data, Index numMinor, bool majorIsSorted = false)
      : _zero(0)
    {
      set(majorRow, data, numMinor, majorIsSorted);
    }
 
    SparseMatrix(bool majorRow, Data&& data, Index numMinor)
      : _zero(0)
    {
      set(majorRow, data, numMinor);
    }
 
    SparseMatrix(const Data& rowData, const Data& columnData)
      : _zero(0)
    {
      set(rowData, columnData);
    }
 
    SparseMatrix(Data&& rowData, Data&& columnData)
      : _zero(0)
    {
      set(std::move(rowData), std::move(columnData));
    }
 
    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)
    {
      if (mayContainZeros)
      {
        // Use first run over entries to count nonzeros in each row (resp. column).
        _rowData.range.resize(numMajor);
        for (auto& range : _rowData.range)
          range.first = 0;
        Index current = 0;
        for (Index i = 0; i < numMajor; ++i)
        {
          _rowData.range[i].first = current;
          Index b;
          if (beyond != nullptr)
            b = beyond[i];
          else if (i + 1 < numMajor)
            b = first[i+1];
          else
            b = numEntries;
          for (Index j = first[i]; j < b; ++j)
          {
            if (entryValues[j] != 0)
              ++current;
          }
          _rowData.range[i].second = current;
        }
 
        // Use second run to copy the nonzeros.
        _rowData.entries.resize(current);
        current = 0;
        for (Index i = 0; i < numMajor; ++i)
        {
          Index b;
          if (beyond != nullptr)
            b = beyond[i];
          else if (i + 1 < numMajor)
            b = first[i+1];
          else
            b = numEntries;
          for (Index j = first[i]; j < b; ++j)
          {
            if (entryValues[j] != 0)
            {
              _rowData.entries[current] = std::make_pair(entryMinors[j], entryValues[j]);
              ++current;
            }
          }
        }
      }
      else
      {
        // If we can trust the input then we can copy directly.
        _rowData.range.resize(numMajor);
        if (beyond != nullptr)
        {
          for (Index i = 0; i < numMajor; ++i)
            _rowData.range[i] = std::make_pair(first[i], beyond[i]);
        }
        else
        {
          for (Index i = 0; i + 1 < numMajor; ++i)
            _rowData.range[i] = std::make_pair(first[i], first[i+1]);
          _rowData.range.back() = std::make_pair(first[numMajor - 1], Index(numEntries));
        }
        _rowData.entries.resize(numEntries);
        for (Index i = 0; i < numEntries; ++i)
          _rowData.entries[i] = std::make_pair(entryMinors[i], entryValues[i]);
      }
 
      _rowData.isSorted = isSorted;
 
      // Construct column-wise (resp. row-wise) representation.
      _columnData.constructFromTransposed(_rowData, numMinor);
 
      if (!majorRow)
        _rowData.swap(_columnData);
 
  #if !defined(NDEBUG)
      ensureConsistency();
  #endif /* !NDEBUG */
    }
 
    void set(bool majorRow, const Data& data, Index numMinor)
    {
      // We can trust the input and thus we can copy directly.
      _rowData.range = data.range;
      _rowData.entries = data.entries;
      _rowData.isSorted = data.isSorted;
 
      // Construct column-wise (resp. row-wise) representation.
      _columnData.constructFromTransposed(_rowData, numMinor);
 
      if (!majorRow)
        _rowData.swap(_columnData);
 
  #if !defined(NDEBUG)
      ensureConsistency();
  #endif /* !NDEBUG */
    }
 
    void set(bool majorRow, Data&& data, Index numMinor)
    {
      _rowData = std::move(data);
      _columnData.constructFromTransposed(_rowData, numMinor);
      if (!majorRow)
        _rowData.swap(_columnData);
 
#if !defined(NDEBUG)
      ensureConsistency();
#endif /* !NDEBUG */
    }
 
    void set(const Data& rowData, const Data& columnData)
    {
      _rowData.range = rowData.range;
      _rowData.entries = rowData.entries;
      _rowData.isSorted = rowData.isSorted;
      _rowData.ensureConsistency();
 
      _columnData.range = columnData.range;
      _columnData.entries = columnData.entries;
      _columnData.isSorted = columnData.isSorted;
      _columnData.ensureConsistency();
 
#if !defined(NDEBUG)
      ensureConsistency();
#endif /* !NDEBUG */
    }
 
    void set(Data&& rowData, Data&& columnData)
    {
      _rowData = std::move(rowData);
      _columnData = std::move(columnData);
#if !defined(NDEBUG)
      ensureConsistency();
#endif /* !NDEBUG */
    }
 
    ~SparseMatrix() = default;
 
    std::size_t numRows() const
    {
      return _rowData.range.size();
    }
 
    std::size_t numColumns() const
    {
      return _columnData.range.size();
    }
 
    bool hasSortedRows() const
    {
      return _rowData.isSorted;
    }
 
    bool hasSortedColumns() const
    {
      return _columnData.isSorted;
    }
 
    const Value& get(Index row, Index column) const
    {
      assert(row < numRows());
      assert(column < numColumns());
 
      if (_rowData.isSorted)
      {
        std::size_t columnIndex = _rowData.find(row, column);
        if (columnIndex < std::numeric_limits<std::size_t>::max())
          return _rowData.entries[columnIndex].second;
      }
      else
      {
        std::size_t rowIndex = _columnData.find(column, row);
        if (rowIndex < std::numeric_limits<std::size_t>::max())
          return _columnData.entries[rowIndex].second;
      }
 
      return _zero;
    }
 
    std::size_t numNonzeros() const
    {
      return _rowData.entries.size();
    }
 
    std::size_t countRowNonzeros(Index row) const
    {
      return _rowData.range[row].second - _rowData.range[row].first;
    }
 
    std::size_t countColumnNonzeros(Index column) const
    {
      return _columnData.range[column].second - _columnData.range[column].first;
    }
 
    NonzeroRowRange iterateRowNonzeros(Index row) const
    {
      return NonzeroRowRange(NonzeroIterator<true>(*this, row),
        NonzeroIterator<true>(*this, row, 0));
    }
 
    NonzeroColumnRange iterateColumnNonzeros(Index column) const
    {
      return NonzeroColumnRange(NonzeroIterator<false>(*this, column),
        NonzeroIterator<false>(*this, column, 0));
    }
 
    SparseMatrix transposed() const
    {
      return SparseMatrix(_columnData, _rowData);
    }
 
    void ensureConsistency() const
    {
      // First ensure individual consistency of row and column data.
      _rowData.ensureConsistency();
      _columnData.ensureConsistency();
 
      // Copy of a nonzero.
      struct NZ
      {
        Index row;
        Index column;
        Value value;
      };
 
      // Comparison for entries.
      auto compare = [](const NZ& a, const NZ& b)
      {
        if (a.row != b.row)
          return a.row < b.row;
        if (a.column != b.column)
          return a.column < b.column;
        return a.value < b.value;
      };
 
      // Construct sorted vector of entries from row data.
      std::vector<NZ> rowNonzeros;
      for (Index i = 0; i < Index(_rowData.range.size()); ++i)
      {
        for (Index j = _rowData.range[i].first; j < _rowData.range[i].second; ++j)
        {
          NZ nz = { i, _rowData.entries[j].first, _rowData.entries[j].second };
          rowNonzeros.push_back(nz);
        }
      }
      std::sort(rowNonzeros.begin(), rowNonzeros.end(), compare);
 
      // Construct sorted vector of entries from column data.
      std::vector<NZ> columnNonzeros;
      for (Index i = 0; i < Index(_columnData.range.size()); ++i)
      {
        for (Index j = _columnData.range[i].first; j < _columnData.range[i].second; ++j)
        {
          NZ nz = { _columnData.entries[j].first, i, _columnData.entries[j].second };
          columnNonzeros.push_back(nz);
        }
      }
      std::sort(columnNonzeros.begin(), columnNonzeros.end(), compare);
 
      for (std::size_t i = 0; i < rowNonzeros.size(); ++i)
      {
        if (rowNonzeros[i].row != columnNonzeros[i].row
          || rowNonzeros[i].column != columnNonzeros[i].column
          || rowNonzeros[i].value != columnNonzeros[i].value)
        {
          throw std::runtime_error("Inconsistent Matrix: row and column data differ.");
        }
      }
    }
 
    void transpose()
    {
      _rowData.swap(_columnData);
    }
 
    
    void sortRows()
    {
      _rowData.sort();
    }
 
    void sortColumns()
    {
      _columnData.sort();
    }
 
    void sort()
    {
      sortRows();
      sortColumns();
    }
 
  private:
    Data _rowData;
    Data _columnData;
    Value _zero;
    bool _isSorted;
  };
 
  template <typename M>
  class TransposedMatrix : public Matrix<typename M::Value>
  {
  public:
    typedef typename M::Value Value;
    typedef typename M::Index Index;
    typedef typename M::Nonzero Nonzero;
 
    template <bool Row>
    class NonzeroIterator
    {
    public:
      NonzeroIterator(const TransposedMatrix<M>& matrix, Index major)
        : _matrix(matrix._matrix), _iterator(typename M::template NonzeroIterator<!Row>(matrix._matrix, major))
      {
 
      }
 
      NonzeroIterator(const TransposedMatrix<M>& matrix, Index major, int dummy)
        : _matrix(matrix._matrix), _iterator(typename M::template NonzeroIterator<!Row>(matrix._matrix, major, 0))
      {
        CMR_UNUSED(dummy);
      }
 
      NonzeroIterator<Row>& operator++()
      {
        ++_iterator;
      }
 
      Nonzero operator*() const
      {
        Nonzero nz = *_iterator;
        std::swap(nz.row, nz.column);
        return nz;
      }
 
      bool operator!=(const NonzeroIterator<Row>& other) const
      {
        assert(&_matrix == &other._matrix);
        return _iterator != other._iterator;
      }
 
    private:
      const M& _matrix;
      typename M::template NonzeroIterator<!Row> _iterator;
    };
 
    typedef detail::Range<NonzeroIterator<true>> NonzeroRowRange;
    typedef detail::Range<NonzeroIterator<false>> NonzeroColumnRange;
 
    TransposedMatrix(M& matrix)
      : _matrix(matrix)
    {
 
    }
 
    TransposedMatrix(TransposedMatrix&& other)
      : _matrix(other._matrix)
    {
 
    }
 
    std::size_t numRows() const
    {
      return _matrix.numColumns();
    }
 
    std::size_t numColumns() const
    {
      return _matrix.numRows();
    }
 
    bool hasSortedRows() const
    {
      return _matrix.hasSortedColumns();
    }
 
    bool hasSortedColumns() const
    {
      return _matrix.hasSortedRows();
    }
 
    const Value& get(Index row, Index column) const
    {
      assert(row < numRows());
      assert(column < numColumns());
 
      return _matrix.get(column, row);
    }
 
    std::size_t numNonzeros() const
    {
      return _matrix.numNonzeros();
    }
 
    std::size_t countRowNonzeros(Index row) const
    {
      return _matrix.countColumnNonzeros(row);
    }
 
    std::size_t countColumnNonzeros(Index column) const
    {
      return _matrix.countRowNonzeros(column);
    }
 
    NonzeroRowRange iterateRowNonzeros(Index row) const
    {
      return NonzeroRowRange(NonzeroIterator<true>(_matrix, row),
        NonzeroIterator<true>(_matrix, row, 0));
    }
 
    NonzeroColumnRange iterateColumnNonzeros(Index column) const
    {
      return NonzeroColumnRange(NonzeroIterator<false>(_matrix, column),
        NonzeroIterator<false>(_matrix, column, 0));
    }
 
    void ensureConsistency() const
    {
      _matrix.ensureConsistency();
    }
 
  protected:
    M& _matrix;
  };
 
  template <typename M>
  TransposedMatrix<M> transposedMatrix(M& matrix)
  {
    return TransposedMatrix<M>(matrix);
  }
 
  class SubmatrixIndices
  {
  public:
    CMR_EXPORT
    SubmatrixIndices(const SubmatrixIndices& other);
 
    CMR_EXPORT
    SubmatrixIndices(SubmatrixIndices&& other);
 
    CMR_EXPORT
    SubmatrixIndices(const std::vector<std::size_t>& rows, const std::vector<std::size_t>& columns);
 
    CMR_EXPORT
    SubmatrixIndices(std::size_t row, std::size_t column);
 
    inline
    std::size_t numRows() const
    {
      return _rows.size();
    }
 
    inline
    std::size_t numColumns() const
    {
      return _columns.size();
    }
 
    inline
    const std::vector<std::size_t>& rows() const
    {
      return _rows;
    }
 
    inline
    const std::vector<std::size_t>& columns() const
    {
      return _columns;
    }
 
  private:
    std::vector<std::size_t> _rows;
 
    std::vector<std::size_t> _columns;
  };
 
  template <typename M>
  class Submatrix : public Matrix<typename M::Value>
  {
  public:
    typedef typename M::Index Index;
    typedef typename M::Value Value;
    typedef typename M::Nonzero Nonzero;
 
    Submatrix(const M& matrix, const SubmatrixIndices& indices)
      : _matrix(matrix), _indices(indices)
    {
 
    }
 
    Submatrix(Submatrix && other)
      : _matrix(other._matrix), _indices(other._indices)
    {
 
    }
 
    std::size_t numRows() const
    {
      return _indices.numRows();
    }
 
    std::size_t numColumns() const
    {
      return _indices.numColumns();
    }
 
    bool hasSortedRows() const
    {
      return false;
    }
 
    bool hasSortedColumns() const
    {
      return false;
    }
 
    const Value& get(Index row, Index column) const
    {
      assert(row < numRows());
      assert(column < numColumns());
 
      return _matrix.get(_indices.rows()[row], _indices.columns()[column]);
    }
 
    void ensureConsistency() const
    {
      for (auto row : _indices.rows())
      {
        if (row >= _matrix.numRows())
          throw std::runtime_error("Inconsistent Submatrix: row index too large.");
      }
      for (auto column : _indices.columns())
      {
        if (column >= _matrix.numColumns())
          throw std::runtime_error("Inconsistent Submatrix: column index too large.");
      }
      _matrix.ensureConsistency();
    }
 
  private:
    const M& _matrix;
    const SubmatrixIndices& _indices;
  };
 
  
  template <typename Matrix>
  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)
  {
    bool result = false;
    int current_value = 0;
    for (size_t r = row_first; r != row_beyond; ++r)
    {
      for (size_t c = column_first; c != column_beyond; ++c)
      {
        int value = matrix(r, c);
        if (value == 0)
          continue;
 
        value = value >= 0 ? value : -value;
 
        if (!result || value < current_value)
        {
          result = true;
          row = r;
          column = c;
          current_value = value;
        }
      }
    }
    return result;
  }
  
  template <typename Matrix>
  bool matrix_row_zero(const Matrix& matrix, size_t row, size_t column_first, size_t column_beyond)
  {
    for (size_t c = column_first; c != column_beyond; ++c)
      if (matrix(row, c) != 0)
        return false;
    return true;
  }
 
  template <typename Matrix>
  bool matrix_column_zero(const Matrix& matrix, size_t column, size_t row_first, size_t row_beyond)
  {
    return matrix_row_zero(make_transposed_matrix(matrix), column, row_first, row_beyond);
  }
 
  
  template <typename Matrix1, typename Matrix2>
  bool equals(const Matrix1& first, const Matrix2& second)
  {
    if (first.size1() != second.size1())
      return false;
    if (first.size2() != second.size2())
      return false;
    for (size_t r = 0; r < first.size1(); ++r)
    {
      for (size_t c = 0; c < first.size2(); ++c)
      {
        if (first(r, c) != second(r, c))
          return false;
      }
    }
    return true;
  }
 
  
  template <typename MatrixType>
  inline void matrix_permute1(MatrixType& matrix, size_t index1, size_t index2)
  {
    for (size_t index = 0; index < matrix.size2(); ++index)
    {
      std::swap(matrix(index1, index), matrix(index2, index));
    }
  }
 
  template <typename MatrixType>
  inline void matrix_permute2(MatrixType& matrix, size_t index1, size_t index2)
  {
    for (size_t index = 0; index < matrix.size1(); ++index)
    {
      std::swap(matrix(index, index1), matrix(index, index2));
    }
  }
 
  
} /* namespace tu */