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

 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)

operator[](i)

 Matrix::Column first_column = A[0];
 max_elt_index = max_index(first_column);
 

Accessing Rows vs. Columns

 max_elt_index = max_index(rows(A)[0]);
 

Submatrices

      [  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

• row_matrix
• column_matrix
• diagonal_matrix
• triangle_matrix
• symmetric_matrix