Matrix Classes

Prerequisits

#include <Matrix.h>
#include <SparseMatrix.h>
#include <ListMatrix.h>

using namespace polymake;

Introduction

The Matrix class family implement matrices in the usual algebraic notion. It contains three persistent matrix types with different data storage organization:

template <typename ElementType> class Matrix;

holds the elements in a contiguous array.

template <typename ElementType> class SparseMatrix;

is a two-dimensional associative array with row and column indices as keys; elements equal to the default value (ElementType(), which is 0 for most numerical types) are not stored, but implicitly encoded by the gaps in the key set. Each row and column is organized as a balanced binary search (AVL) tree.

template <typename VectorType> class ListMatrix;

is a list of row vectors. Based on std::list. Can be parameterized with any persistent vector type.

The following brief analysis might help you when choosing the most efficient representation for a concrete case. r and c are number of rows and columns respectively, and n is the average number of non-zero elements in a row/column.

SparseMatrix deploys an adaptive reallocation strategy similar to std::vector, reserving additional stock memory by every reallocation. This way, the amortized reallocation costs are proportional to the logarithm of the number of appended rows/columns.

Performance comparison
Operation	`Matrix`	`SparseMatrix`	`ListMatrix`
rowwise iteration	O(r) with approx. equal overhead
columnwise iteration	O(c)	O(c)	O(c), but with rather high overhead
random row/column access	O(1)	O(1)	—
random element access	O(1)	O(log(n))	—
append a row	O(r c) + reallocation costs	O(r) + reallocation costs + n O(log(n))	O(1)
append a column	O(r c) + reallocation costs	O(c) + reallocation costs + n O(log(n)) See also the special class `RestrictedSparseMatrix`.	c*(appending costs for a `VectorType`)

Matrix and SparseMatrix keep the data attached to a smart pointer with reference counting.

Generic type

template <typename Top, typename ElementType=typename Top::element_type> class GenericMatrix;

The generic class for all matrix classes. Defines the following types and constants:

`element_type`	`ElementType`, the persistent element type
`top_type`	`Top`, the concrete matrix object behind the generic reference
`generic_type`	`GenericMatrix<Top>` itself
`persistent_type`	`Matrix<ElementType>` for dense matrices, `SparseMatrix<ElementType>` for sparse matrices.
`row_type` `col_type`	Alias classes representing single matrix rows and columns respectively. They always belong to the `GenericVector` family.
`bool is_sparse`	`true` if the element sequence in an arbitrary row or column can contain gaps expressing implicit zero elements.

Methods and friend functions applicable to GenericMatrix directly are scattered over the entire page. They are discernible on the prefix GenericMatrix:: or parameters of type GenericMatrix&. To use other methods via a generic reference (or pointer), you must first move to the concrete object using the top() method.

Construction

Any persistent class: Matrix, SparseMatrix, or ListMatrix.

Any persistent dense class: Matrix, or ListMatrix built with dense row vectors.

Matrix();

Create a matrix with 0 rows and 0 columns.

Matrix (int m, int n);

Create a matrix with m rows and n columns, initialize all elements to 0 (in sparse cases implicitly.)

Matrix (int m, int n, const ElementType& value);

Create a matrix with m rows and n columns, initialize all elements with value.

template <typename Iterator>

Matrix (int m, int n, Iterator src);

Create a matrix with m rows and n columns, initialize the elements from a data sequence. src should iterate either over m*n scalar values, corresponding to the elements in the row order (the column index changes first,) or over m vectors of dimension n, corresponding to the matrix rows. In the sparse case, zero input elements are filtered out.

template <typename OtherMatrix>

Matrix (const GenericMatrix<OtherMatrix, ElementType>&);

Create a matrix as a copy of another matrix of the same element type.

template <typename OtherMatrix, typename OtherElementType>

explicit Matrix (const GenericMatrix<OtherMatrix, OtherElementType>&);

Copying with element conversion is also possible, but requires an explicit constructor call.

GenericMatrix convert_to<OtherElementType> (const GenericMatrix&);

Create a matrix with another element type, using explicit type conversion.

Modification

Any persistent matrix: Matrix, SparseMatrix, ListMatrix, or an alias object referring to a mutable (non-const) martix.

Any sparse matrix: SparseMatrix, ListMatrix<SparseVector>, or an alias object referring to a mutable sparse matrix.

This assertion is active only in the debugging compilation mode.

template <typename OtherMatrix>

GenericMatrix::operator= (const GenericMatrix<OtherMatrix, ElementType>&);

Persistent matrix objects have after the assignment the same dimensions as the right hand side operand. Alias objects, such as matrix minor or block matrix, cannot be resized, thus must have the same dimensions as on the right hand side.

void std::swap(GenericMatrix&, GenericMatrix&);

Exchange the contents of two matrices in a most efficient way. If at least one non-persistent object is involved, the operands must have equal dimensions.

void Matrix::resize(int m, int n);

Extend or truncate to new dimensions (m rows, n columns.) Surviving elements keep their values, new elements are initialized with 0 (in sparse cases implicitly.)

SparseMatrix deploys an adaptive reallocation strategy similar to std::vector, reserving additional stock memory by every reallocation. A special case, looking at the first glance like a "no operation": M.resize(M.rows(), M.cols()), gets rid of this extra allocated storage.

void Matrix::clear();

Truncate to 0x0 matrix.

template <typename Scalar>

void GenericMatrix::fill(const Scalar& value);

Fill with given value without changing the dimensions. value can be of arbitrary type assignable to ElementType.

void GenericMatrix::remove0s();

Physically remove all zero elements from a sparse matrix that might have creeped in by some previous operation.

The following methods are defined for SparseMatrix objects only.

void squeeze();

template <typename row_index_consumer>

void squeeze(row_index_consumer rc);

template <typename row_index_consumer, typename col_index_consumer>

void squeeze(row_index_consumer rc, col_index_consumer cc);

Remove all empty (i.e., consisting entirely of implicit zeroes,) rows and columns, renumber the rest, and reduce the dimensions. If you need to know the exact mapping between the old and new row (column) indices, you can supply output iterators (e.g., back_inserter of a std::list.) They will get the old indices of non-empty rows (columns) assigned in the ascending order.

template <typename Iterator>

void permute_rows(Iterator perm);

template <typename Iterator>

void permute_cols(Iterator perm);

template <typename Iterator>

void permute_inv_rows(Iterator inv_perm);

template <typename Iterator>

void permute_inv_cols(Iterator inv_perm);

Permute the rows(columns) of the matrix without copying the elements. These operations are nevetherless expensive, as they need to visit each element and adjust its indices.

An iterator argument should run over a sequence of integer indices. In the straight variant (perm), it has the same meaning as in the select function. In the inverted variant (inv_perm), it is an inverse permutation: the first element specifies the new position of the 0-th row (column) etc.

If the index sequence does not build a proper permutation, the results are undefined. It is the responsibility of the application context to assure this condition.

The following methods are defined for ListMatrix objects only.

typename Rows<ListMatrix>::iterator insert_row (typename Rows<ListMatrix>::iterator pos, const GenericVector&);

Insert a new row at the given position: before the row pointed to by pos.

void ListMatrix::delete_row (typename Rows<ListMatrix>::iterator pos);

Delete the row pointed to by pos. Be aware that the iterator becomes void after the removal unless rescued in good time (with postincrement or postdecrement).

Element access

The concrete matrix type hidden behind the GenericMatrix.

int GenericMatrix::rows() const;

int GenericMatrix::cols() const;

Tell the current matrix dimensions.

Sequential access

The matrix classes as such do not implement any STL container interfaces: they are simply not defined for multi-dimensional containers. However, a matrix can be disguised in three different ways, yielding STL-conforming containers:

pseudo-container class	producing function	description
`Rows<Top>` `const Rows<Top>`	`rows(GenericMatrix&);` `rows(const GenericMatrix&);`	A sequence of row vectors; each element is a `GenericVector`.
`Cols<Top>` `const Cols<Top>`	`cols(GenericMatrix&);` `cols(const GenericMatrix&);`	A sequence of column vectors; each element is a `GenericVector`
`Serialized<Top>` `const Serialized<Top>`	`serialize(GenericMatrix&);` `serialize(const GenericMatrix&);`	A sequence of elements in rowwise order (the column index changes first). The `Serialized` object itself is a `GenericVector` of dimension `rows()*cols()`

`SparseMatrix` specialities

int Rows< SparseMatrix >::iterator::index();

int Cols< SparseMatrix >::iterator::index();

Tell the current row (column) number. The numbers start, as usual, with 0.

int SparseMatrix::row_type::iterator::index();

int SparseMatrix::col_type::iterator::index();

Tell the current column (row) number of the element pointed to by the given iterator. Remember that rows and columns of SparseMatrix are sparse containers, so the iterators over them visit only the (presumably non-zero) elements physically stored in the matrix.

SparseMatrix::col_type::iterator cross_direction(const SparseMatrix::row_type::iterator&);

SparseMatrix::row_type::iterator cross_direction(const SparseMatrix::col_type::iterator&);

For an element of a sparse matrix pointed to by the given iterator over the row, create an iterator over the column pointing to the same element (and vice versa.)

All these methods are, indeed, defined for the corresponding const_iterator classes too.

Random access

Random-row-access matrix: Matrix, SparseMatrix, block matrix built of random-row-access matrices, or a minor of a random-row-access matrix with a random-access row index set (sequence or series).

Random-column-access matrix: Matrix, SparseMatrix, block matrix built of random-column-access matrices, or a minor of a random-column-access matrix with a random-access column index set (sequence or series).

Random-access matrix: Matrix, block matrix built of random-access matrices, or a minor of a random-access matrix with random-access row and column index sets (sequence or series).

The real return type is a temporary proxy object, forwarding the read/write access operations to the sparse matrix. Prior to the destruction, which happens after the completion of the expression envolving the random access, it checks whether the element became zero; if so, it is removed from the matrix.

The real result type is a temporary alias object of type pm::VectorSlice. It inherits the const-ness from the original matrix.

ElementType& Matrix::operator() (int i, int j);

const ElementType& Matrix::operator() (int i, int j) const;

Get the element in i-th row, j-th column. The index range check can be activated.

ElementType& SparseMatrix::operator() (int i, int j);

const ElementType& SparseMatrix() (int i, int j) const;

The index range check can be activated. The const operator returns a reference to a static zero object if the specified element is an implicit 0.

GenericVector Matrix::row(int i);

GenericVector Matrix::operator[] (int i);

GenericVector Matrix::col(int i);

Get the i-th row or column respectively. The const variants of these methods create immutable slice objects. The index range check can be activated.

Rows and columns of a SparseMatrix offer the binary tree data access methods additionally to the GenericVector interface.

operator[] is a mere synonym for row(int), and is defined with the sole purpose to allow the C-style expressions M[i][j]. Please be aware that M(i,j) is a much more efficient way to access a single matrix element, as it doesn't involve creation of a temporary vector row object.

Operations

The real return type is a temporary object (expression template) implementing the lazy evaluation of the operation result. You should have it in mind when passing the result of a matrix expression as a GenericMatrix argument to a template function where it has a chance to be evaluated multiple times. See also prevent_lazy class.

The vector and matrix operators defined in the Polymake Template Library do take care of this, so this remark applies more likely to your own functions.

const GenericMatrix operator- (const GenericMatrix&);

const GenericMatrix operator+ (const GenericMatrix&, const GenericMatrix&);

const GenericMatrix operator- (const GenericMatrix&, const GenericMatrix&);

const GenericMatrix operator* (const Scalar&, const GenericMatrix&);

const GenericMatrix operator* (const GenericMatrix&, const Scalar&);

const GenericMatrix operator/ (const GenericMatrix&, const Scalar&);

const GenericMatrix operator* (const GenericMatrix&, const GenericMatrix&);

const GenericMatrix operator* (const GenericMatrix&, const GenericVector&);

const GenericMatrix operator* (const GenericVector&, const GenericMatrix&);

Do exactly what the mnemonics suggest. Every type combination is allowed, including mixing of sparse and dense vectors and matrices. ElementType of matrix and vector operands, and Scalar type can be different, provided the arithmetic operator with corresponding arguments exists. The dimension check can be activated.

Top& GenericMatrix::negate();

Top& GenericMatrix::operator+= (const GenericMatrix&);

Top& GenericMatrix::operator-= (const GenericMatrix&);

Top& GenericMatrix::operator*= (const Scalar&);

Top& GenericMatrix::operator/= (const Scalar&);

Corresponding assignment versions of the arithmetic operators.

bool operator== (const GenericMatrix&, const GenericMatrix&);

bool operator!= (const GenericMatrix&, const GenericMatrix&);

bool operator< (const GenericMatrix&, const GenericMatrix&);

bool operator> (const GenericMatrix&, const GenericMatrix&);

bool operator<= (const GenericMatrix&, const GenericMatrix&);

bool operator>= (const GenericMatrix&, const GenericMatrix&);

Compare two matrices lexicographically (iterating over elements in rowwise order).

Other constructions

The real result type is a masquerade reference Transposed&, referring to the originator matrix. It inherits the const-ness from the latter.

The real result type is a temporary alias object of type pm::MatrixMinor. It inherits the const-ness from the originator matrix.

The real result type is a temporary alias object of type pm::RowChain or pm::ColChain. It is mutable only if both operands are mutable too.

The real result type is a masquerade reference pm::SingleRow& or SingleCol&, pointing to the source vector.

The real result type is a temporary object of type pm::DiagMatrix.

The real result type is a temporary object of type pm::RepeatedRow or pm::RepeatedCol. It contains only a single copy of a vector, or just a reference to it, but makes it to appear as a sequence of identical row (column) vectors.

The real result type is a temporary object of type pm::SameElementSparseMatrix. It contains only a single copy of an element, or just a reference to an element, but makes it to appear as a matrix of identical elements.

Any combination of GenericSet, Complement, and a special identifier All (choosing all rows/columns) is allowed here.

GenericMatrix GenericMatrix::minor (const RowIndexSet& row_indices, const ColIndexSet& column_indices);

Select a matrix minor lying on the intersection of the given row and column subsets. The const variant of this method creates an immutable minor object. The index range check can be activated.

GenericVector GenericMatrix::diagonal (int i=0);

GenericVector GenericMatrix::anti_diagonal (int i=0);

Select the diagonal or anti-diagonal of a matrix:

i=0	the main diagonal (anti-diagonal)
i>0	the i-th diagonal below the main
i<0	the (-i)-th diagonal above the main

The const variants of these methods create immutable slice objects.

GenericMatrix T (GenericMatrix&);

Transpose the matrix. The roles of Rows and Cols get swapped.

GenericMatrix operator/ (const GenericMatrix&, const GenericMatrix&);

GenericMatrix operator/ (const GenericVector&, const GenericMatrix&);

GenericMatrix operator/ (const GenericMatrix&, const GenericVector&);

GenericMatrix operator/ (const GenericVector&, const GenericVector&);

Create a block matrix, virtually appending the rows of the second matrix after the last row of the first matrix. A vector argument is treated as a matrix with one row. Column dimensions of the operands must be equal. ElementType of the operands must be identical, otherwise you will get a rather cryptic message from the compiler.

Sorry for a rather weird mnemonics...

GenericMatrix operator| (const GenericMatrix&, const GenericMatrix&);

GenericMatrix operator| (const GenericVector&, const GenericMatrix&);

GenericMatrix operator| (const GenericMatrix&, const GenericVector&);

Create a block matrix, virtually appending the columns of the second matrix after the last column of the first matrix. A vector argument is treated as a matrix with one column. Row dimensions of the operands must be equal. ElementType of the operands must be identical, otherwise you will get a rather cryptic message from the compiler.

Matrix& Matrix::operator/= (const GenericMatrix&);

Matrix& Matrix::operator/= (const GenericVector&);

Matrix& Matrix::operator|= (const GenericMatrix&);

Matrix& Matrix::operator|= (const GenericVector&);

Append rows or columns to the matrix. All constraints for the non-assigning operators apply here too.

const GenericMatrix diag(const GenericMatrix&, const GenericMatrix&);

const GenericMatrix diag(const GenericMatrix&, const GenericVector&);

const GenericMatrix diag(const GenericVector&, const GenericMatrix&);

Create a block-diagonal matrix. Vector arguments are treated as square diagonal matrices.

GenericMatrix vector2row(const GenericVector& v);

GenericMatrix vector2col(const GenericVector& v);

Disguise v as a matrix with 1 row (column).

const GenericMatrix diag(const GenericVector& v);

Create a square diagonal matrix. The elements of v appear on the main diagonal.

const GenericMatrix repeat_row(const GenericVector& v, int n);

const GenericMatrix repeat_col(const GenericVector& v, int n);

Create a matrix with n rows (columns), each equal to v.

const GenericMatrix same_element_matrix(const ElementType& x, int m, int n);

const GenericMatrix ones_matrix<ElementType>(int m, int n);

const GenericMatrix zero_matrix<ElementType>(int m, int n);

Create a dense matrix with m rows and n columns, with all elements equal to x, 1, or 0, respectively. Note the obligatory explicit specification of the template parameter ElementType, where it cannot be deduced from the function arguments.

const GenericMatrix same_element_sparse_matrix<ElementType> (const GenericIncidenceMatrix& M);

const GenericMatrix same_element_sparse_matrix (const GenericIncidenceMatrix& M, const ElementType& x);

Create a sparse matrix with the same structure as M, with elements equal to x in true cells and 0 in false cells. If omitted, x is taken equal to 1; note the obligatory explicit template parameter specification in this case.

const GenericMatrix unit_matrix<ElementType> (int n);

Create a unit matrix n*n.

Input/output

std::ostream& operator<< (std::ostream&, const GenericMatrix&);

Print the matrix to the output stream rowwise. If std::ostream::width is not 0 before the operator is called, it is repeatedly restored for each row, thus yielding a column-aligned table.

explicit Matrix (Poly&);

Poly& operator>> (Poly&, Matrix&);

Read a property from polymake server. Change the matrix dimensions according to the obtained data. Parse errors are reported by raising the std::iostream::failure exception in the constructor, or via std::iostream::fail() flag in the input operator.

Restricted sparse matrix

Rows in the only_rows case, columns in the only_cols.

Columns in the only_rows case, rows in the only_cols.

enum sparse2d::restriction_kind { only_rows, only_cols };

template <typename ElementType, sparse2d::restriction_kind restriction=only_rows> class RestrictedSparseMatrix;

A special class allowing for efficient incremental building of sparse matrices.

If you have to fill a sparse matrix element for element, but do not know one of its dimensions in advance, you would have to resize it repeatedly, which would severely impact the performance. This special class maintains the internal data only for one direction, and thus can be filled cheaper than SparseMatrix.

On the other hand, almost all matrix operations described above need information about both directions (rows and columns.) Therefore, RestrictedSparseMatrix is deliberately not included in the GenericMatrix class family. Its only puspose is to gather elements and finally hand them over to a SparseMatrix object (without copying!)

RestrictedSparseMatrix();

Create a matrix with 0 rows and 0 columns.

explicit RestrictedSparseMatrix(int n);

Create a matrix with n rows and 0 columns (or vice versa in only_cols case), all elements are implicit zeroes.

template <typename Iterator>

RestrictedSparseMatrix(int n, Iterator src);

Create a matrix with n rows and 0 columns (or vice versa in only_cols case), initialize the rows (columns) from the input sequence. Each item supplied by src should be a GenericVector.

explicit RestrictedSparseMatrix(Poly&);

Read a property from the polymake server. Throw std::iostream::failure if the input data can't be parsed as items of ElementType. Currently this constructor is defined only for only_rows case.

void RestrictedSparseMatrix::clear();

Truncate to a 0x0 matrix.

void std::swap(RestrictedSparseMatrix&, RestrictedSparseMatrix&);

Swap the contents of two matrices in a most efficient way.

int RestrictedSparseMatrix::rows() const;

int RestrictedSparseMatrix::cols() const;

Tell the current matrix dimensions. For the unsupported direction, the value returned is the highest column (row) index ever seen in an element access or modification operation.

ElementType& RestrictedSparseMatrix::operator() (int i, int j);

const ElementType& RestrictedSparseMatrix() (int i, int j) const;

Random access to an element. The index range check for the main direction can be activated. The const operator returns a reference to a static zero object if the specified element is an implicit 0.

Rows< RestrictedSparseMatrix >& rows(RestrictedSparseMatrix&);

const Rows< RestrictedSparseMatrix >& rows(const RestrictedSparseMatrix&);

Cols< RestrictedSparseMatrix >& cols(RestrictedSparseMatrix&);

const Cols< RestrictedSparseMatrix >& cols(const RestrictedSparseMatrix&);

Sequential access to the rows and columns. Defined only for the main direction.

GenericVector RestrictedSparseMatrix::row (int i);

GenericVector RestrictedSparseMatrix::col (int i);

Random access to a row (column). Defined only for the main direction.

RestrictedSparseMatrix& RestrictedSparseMatrix::operator/= (const GenericMatrix&);

RestrictedSparseMatrix& RestrictedSparseMatrix::operator/= (const GenericVector&);

Append rows to the matrix. Allowed in both only_rows and only_cols modi.

RestrictedSparseMatrix& RestrictedSparseMatrix::operator|= (const GenericMatrix&);

RestrictedSparseMatrix& RestrictedSparseMatrix::operator|= (const GenericVector&);

Append columns to the matrix. Allowed in both only_rows and only_cols modi.

explicit SparseMatrix(RestrictedSparseMatrix<ElementType>&);

SparseMatrix& SparseMatrix::operator=(RestrictedSparseMatrix<ElementType>&);

The final and alone destination of a RestrictedSparseMatrix object is to hand over its internal data to a new created SparseMatrix. Any further action on the RestrictedSparseMatrix object (except destruction, obviously) is prohibited: it would lead to the immediate program crash.

SparseMatrix2x2

template <typename ElementType> struct SparseMatrix2x2 {
   int i,j;
   ElementType a_ii, a_ij, a_ji, a_jj;
};

A compact encoding of matrices of special structure, which occur as companion matrices in elimination, triangularization, etc. algorithms. Note that this class is not derived from GenericMatrix.

For two given integer indices i and j, it encodes a sparse square matrix A of arbitrary dimensions with 4 designated elements A_i,i, A_i,j, A_j,i, and A_j,j. Other elements are implicit zeroes except for the rest of the main diagonal, which is filled with ones.

The real return type is a masquerade reference const Transposed<SparseMatrix2x2>&.

SparseMatrix2x2(int _i, int _j,
                const ElementType& _a_ii, const ElementType& _a_ij,
                const ElementType& _a_ji, const ElementType& _a_jj);

The only constructor for this class.

SparseMatrix2x2 T (const SparseMatrix2x2&);

ElementType det(const SparseMatrix2x2&);

Elementary operations: transposition, determinant.

void GenericMatrix::multiply_from_left(const SparseMatrix2x2& U);

Multiply (*this) with U from the left, that is, apply U to the rows.

void GenericMatrix::multiply_from_right(const SparseMatrix2x2& U);

Multiply (*this) with U from the right, that is, apply U to the columns.

Prerequisits

Introduction

Generic type

Sequential access

SparseMatrix specialities

Random access

`SparseMatrix` specialities