Matrix Programmers Guide Contents | Index |  Search

Category:containers Component type:concept
Description

The central concept in the Matrix Template Library is of course the Matrix. An MTL Matrix can be thought of as a Container of Containers, or a two-dimensional container. As an STL Container an MTL Matrix has begin() and end() functions which return iterators for traversing over the 2D container. These iterators dereferece to give a 1D container, which also has begin() and end() functions, since it to is an STL Container. The MTL implements a large variety of matrix formats, all with very different data representations and implementation details. However, the same MTL Matrix interface is provided for all of the matrices. This Matrix interface is described here.

1D and 2D Iterators

The following is a simplified example of how the Matrix iterators can be used to write a generic matrix-vector product.

``` template < class Matrix, class VecX, class VecY >
void matvec_mult(const Matrix& A, VecX x, VecY y) {
typename Matrix::const_iterator i;
typename Matrix::OneD::const_iterator j;
for (i = A.begin(); i != A.end(); i)
for (j = (*i).begin(); j != (*i).end(); j)
y[j.row()] += *j * x[j.column()];
}
```

There are two iterators used in the above algorithm, i and j. We refer to i as a 2D iterator since it iterates through the 2D container. We refer to j as the 1D iterator. *i gives a 1D Container, and the (*i).begin() and (*i).end() expressions define the iteration through the 1D part of the Matrix, which could be a Row, Column, or Diagonal depending on the Matrix type.

MTL matrices can also be sparse or dense, so the traversal behaviour of the 1D iterator j varies accordingly. j only iterates over the non-zero elements in the sparse case. In addition the row() and column() functions on the 1D iterator return the row and column index corresponding to the position of the iterator. This hides the differences in indexing between a large class of sparse and dense matrices.

Compare the MTL-style code to that of the BLAS-style and the SPARSKIT matrix-vector products. The BLAS and SPARSKIT algorithms include a lot of details that are specific to the matrix format being used. For instance, the sparse matrix algorithm must explicitly access the index and pointer arrays ia and ja.

``` // SPARSKIT-style sparse matrix-vector multiply
void matvec_mult(double* a, int n, int* ia, int* ja,
double* y, double* x) {
for (int i = 0; i < n; i)
for (int k = ia[i]; k < ia[i+1]; k)
y[i] +=  a[k] * x[ja[k]];
}
}
```

The dense matrix algorithm must use the leading dimension lda of the matrix to map to the appropriate place in memory.

``` // BLAS-style dense matrix-vector multiply
void matvec_mult(double* a, int m, int n, int lda,
double* y, double* x) {
for (int i = 0; i < m; i)
for (int j = 0; j < n; j)
y[i] += a[i*lda+j] * x[j];
}
}
```

The MTL-style algorithm has none of these format specific expressions because the implementation details are hidden under the MTL Matrix interface.

If one uses a dense MTL Matrix with the generic algorithm, the resulting assembly code looks similar to the BLAS-style algorithm. If one uses a sparse MTL Matrix with the generic algorithm, the resulting assembly code looks similar to the SPARSKIT-style algorithm.

operator()(i,j)

One may ask, what about the A(i,j) operator? MTL matrices do have an A(i,j) operator, but using it typically is not the most efficient or ``generic'' way to traverse and access the matrix, especially if the matrix is sparse. Furthermore, if the Matrix is banded, one would have to consider which indices of the Matrix are valid to access. With the MTL Matrix iterators, one does not have to worry about such considerations. Traversing with the MTL Matrix iterators gives you access to all matrix elements that are stored in any given storage format.

operator[](i)

This operator gives access to the OneD containers within a Matrix. For instance, if you have a column_major matrix A, and wish to find the index of the maximum element in the first column, one would do the following:
``` Matrix::Column first_column = A[0];
max_elt_index = max_index(first_column);
```

Accessing Rows vs. Columns

Most matrix types only provide access to either rows or columns, but not both. However, if one has a dense rectangle matrix, then the matrix can be easily converted back and forth from row_major to column_major using the rows and columns helper functions. The following finds the maximum element in the first row of the matrix.
``` max_elt_index = max_index(rows(A)[0]);
```

Submatrices

The sub_matrix() function returns a new matrix object that is a view into a particular portion of the matrix. The indexing within the new matrix is reset so that the first element is at (0,0). Here is an example of creating a submatrix:

```      [  1  2  3  4 ]
A =  [  5  6  7  8 ]
[  9 10 11 12 ]
[ 13 14 15 16 ]

A_00 = [  1  2 ]  A_01 = [ 3 4 ]
[  5  6 ]         [ 7 8 ]

A_10 = [  9 10 ]  A_11 = [ 11 12 ]
[ 13 14 ]         [ 15 16 ]

A_00 = A.sub_matrix(0,2,0,2);
A_01 = A.sub_matrix(2,4,0,2);
A_10 = A.sub_matrix(0,2,2,4);
A_11 = A.sub_matrix(2,4,2,4);
```

If one wants to create submatrices for the whole matrix, as above, MTL provides a short cut in the partition fuction. The function returns a Matrix of submatrices. The input is the row and column numbers that split up the matrices. The following code gives an example.

``` int splitrows[] = { 2 };
int splitcols[] = { 2 };
Matrix::partitioned Ap(2,2);
partition(array_to_vec(splitrows),
array_to_vec(splitcols), A, Ap);
```

Now Ap(0,0) is equivalent to A_00, Ap(0,1) is equivalent to A_01, etc.

Matrix Type Selection

To create an MTL Matrix, one uses the matrix type constructor to choose the particular storage format, element type, etc. This performs a shallow copy, since MTL objects are really just reference counted handles to the actually data.
Refinement of
Container
Associated types
Concept Type name Description
Tag X::shape See matrix_traits
Tag X::orientation Either row_tag or column_tag
Tag X::sparsity Either dense_tag or sparse_tag
Tag X::dimension Either oned_tag or twod_tag
Matrix X::transpose_type Used by trans helper function
Matrix X::strided_type Used by rows and columns helper functions
Matrix X::scaled_type Used by scaled helper function
Vector X::OneD The type for a OneD slice of the Matrix
Vector X::OneDRef The type for a reference to a OneD slice of the Matrix
Matrix X::submatrix_type The type for a submatrix of this Matrix
Matrix X::partition_type The type for a partitioned version of this Matrix
TrivialConcept X::value_type The element type
TrivialConcept& X::reference Reference to the element type
const TrivialConcept& X::const_reference Const reference to the element type
TrivialConcept* X::pointer Pointer to the element type
BidirectionalIterator X::iterator Iterator type, dereference gives a OneD
const BidirectionalIterator X::const_iterator Const iterator type
BidirectionalIterator X::reverse_iterator Reverse iterator type
const BidirectionalIterator X::const_reverse_iterator Const reverse iterator type
NonNegativeIntegralType X::size_type Size type
IntegralType X::difference_type Difference type
size_type X::M, N for static sized matrices
Notations

 X The type of a model of the Matrix concept A,B An object of type X stream A matrix stream for file IO
Definitions
Expression semantics
Description Expression Semantics
Normal Constructor X(m, n) or X A(m, n);
Constructor for banded matrices X(m, n, sub, super) or X A(m, n, sub, super);
Copy Constructor X(B) or X A(B);
Contstruct for external matrices X A(data, m, n);
Construct for external matrices, with different leading dimension (ld) X A(data, m, n, ld);
Constructor for external banded matrices X A(data, m, n, ld, sub, super);
Static sized dimension constructor X A(data);
Sparse external matrix constructor X A(m, n, nnz, val, ptrs, inds);
Banded view constructor X A(B, sub, super);
Matrix market stream constructors X A(stream);
Harwell Boeing stream constructor X A(stream);
Matrix market stream banded constructors X A(stream, sub, super);
Harwell Boeing stream banded constructor X A(stream, sub, super);
Iterate over OneD slices A.begin()
Const Iterator Begin A.begin()
Iterator End A.end()
Const Iterator End A.end()
Reverse Iterator Begin A.rbegin()
Const Reverse Iterator Begin A.rbegin()
Reverse Iterator End A.rend()
Const Reverse Iterator End A.rend()
Element Access A(i,j)
Const Element Access A(i,j)
OneD Access A[i]
Const OneD Access A[i]
Submatrix Access A.sub_matrix(row_start, row_finish, col_start, col_finish)
Number of rows A.nrows()
Number of columns A.ncols()
Number of non-zero elements (really number of stored elements) A.nnz()
Sub part of bandwidth A.sub()
Super part of bandwidth A.super()
Whether the main diagonal is all ones, and hence not stored A.is_unit()
Whether the matrix has elements stored only in the upper triangle A.is_upper()
Whether the matrix has elements stored only in the upper triangle A.is_lower()
Function specification
Name Function Complexity
Iterate over OneD slices iterator begin() constant time
Const Iterator Begin const_iterator begin() constant time
Iterator End iterator end() constant time
Const Iterator End const_iterator end() constant time
Reverse Iterator Begin reverse_iterator rbegin() constant time
Const Reverse Iterator Begin const_reverse_iterator rbegin() constant time
Reverse Iterator End reverse_iterator rend() constant time
Const Reverse Iterator End const_reverse_iterator rend() constant time
Element Access reference operator()(size_type i, size_type j) constant for dense, linear for sparse
Const Element Access const_reference operator()(size_type i, size_type j) constant for dense, linear for sparse
OneD Access OneDRef operator[](size_type i) constant
Const OneD Access const OneDRef operator[](size_type i) constant
Number of rows size_type nrows() constant time
Number of columns size_type ncols() constant time
Number of non-zero elements (really number of stored elements) size_type nnz() constant time
Invariants
Models
Notes
See also