MTL 4: Matrix Types

Right now, MTL4 provides five matrix types:

- matrix::dense2D;
- matrix::morton_dense;
- matrix::compressed2D; and
- multi_vector, see Type Multivector
- element_structure.

The type matrix::dense2D defines regular row-major and column-major matrices:

// File: dense2D.cpp #include <iostream> #include <boost/numeric/mtl/mtl.hpp> int main(int, char**) { using namespace mtl; // A is a row-major matrix dense2D<double> A(10, 10); // Matrices are not initialized by default A= 0.0; // Assign a value to a matrix element A(2, 3)= 7.0; // You can also use a more C-like notation A[2][4]= 3.0; std::cout << "A is \n" << A << "\n"; // B is a column-major matrix dense2D<float, matrix::parameters<tag::col_major> > B(10, 10); // Assign the identity matrix times 3 to B B= 3; std::cout << "B is \n" << B << "\n"; return 0; }

If no matrix parameters are defined, dense matrices are by default row-major. There are more matrix parameters besides the orientation. As they are not yet fully supported we refrain from discussing them.

Matrix elements can be accessed by a(i, j) or in the more familiar form a[i][j]. The second form is internally transformed into the first one at compile time so that the run-time performance is not affected (unless the compiler does not inline completely which we never observed so far). Also, the compile time is not conceivably increased by this transformation.

Please notice that overwriting single matrix elements is only defined for dense matrix types. For a generic way to modify matrices see Matrix Insertion.

Assigning a scalar value to a matrix stores a multiple of the identity matrix, i.e. the scalar is assigned to all diagonal elements and all off-diagonal elements are 0. Diagonal elements are matrix entries with identical row and column index. Therefore scalars can also be assigned to non-square matrices. This operation is generic (i.e. applicable to all matrix types including sparse).

Just in case you wonder why the scalar value is only assigned to the diagonal elements of the matrix not to all entries, this becomes quite clear when you think of a matrix as a linear operator (from one vector space to another one). For instance, consider the multiplication of vector x with the scalar alpha:

y= alpha * x;

where y is a vector, too. This operation is equivalent to assigning alpha to the matrix A and multiplying x with A:

A= alpha; y= A * x;

In other words, the matrix A has the same impact on x as the scalar alpha itself.

If the matrix is not square, i.e. the linear operator's domain and image have different dimensions, the equivalence with the scalar multiplication applies accordingly. In case that the image has a lower dimension, say m, then only the first m entries of the vector from the domain are scaled with alpha and the others are ignored. If the image has an higher dimension then the last m-n entries are zero with n the dimension of the domain. When you rely on this behavior please check the revision of your MTL4 library: old versions, i.e. before revision 6843, considered it erroneous to store a non-zero scalar to a non-square matrix.

From a more pragmatic prospective:

A= 0;

can be used to clear any matrix, square or rectangular, sparse and dense.

Dense matrices with a recursively designed memory layout can be defined with the type matrix::morton_dense :

// File: morton_dense.cpp #include <iostream> #include <boost/numeric/mtl/mtl.hpp> int main(int, char**) { using namespace mtl; // Z-order matrix morton_dense<double, recursion::morton_z_mask> A(10, 10); A= 0; A(2, 3)= 7.0; A[2][4]= 3.0; std::cout << "A is \n" << A << "\n"; // B is an N-order matrix with column-major 4x4 blocks, see paper morton_dense<float, recursion::doppled_4_col_mask> B(10, 10); // Assign the identity matrix times 3 to B B= 3; std::cout << "B is \n" << B << "\n"; return 0; }

In the pure Morton order format 2 by 2 sub-matrices are stored contiguously in memory. 4 by 4 matrices constitute of 4 2-by-2-matrices and use consecutive memory. The continuation of this recursive scheme provides square sub-matrices with power of two sizes that are in contiguous memory and allow for cache-efficient recursive algorithms. On the other hand, algorithms that are implemented fully recursively create considerable overhead for function calls. We therefore recommend using mixed schemes of recursion and iteration. Particularly efficient are algorithms that operate on mid-size blocks, e.g. 64 by 64, with regular row-major or column-major layout. MTL4 provides a virtually infinite number of memory layouts for dense matrices that are specified by a bitmask. A detailed description and discussion of recursive matrices and algorithm is provided in this conference paper.

Sparse matrices are defined with the type matrix::compressed2D :

// File: compressed2D.cpp #include <iostream> #include <boost/numeric/mtl/mtl.hpp> int main(int, char**) { using namespace mtl; // CRS matrix compressed2D<double> A(12, 12); // Laplace operator discretized on a 3x4 grid matrix::laplacian_setup(A, 3, 4); std::cout << "A is \n" << A; // Element access is allowed for reading std::cout << "A[3][2] is " << A[3][2] << "\n\n"; // CCS matrix compressed2D<float, matrix::parameters<tag::col_major> > B(10, 10); // Assign the identity matrix times 3 to B B= 3; std::cout << "B is \n" << B << "\n"; return 0; }

Matrix a is stored as compressed row storage (CRS). Its assigned values correspond to a discretized Laplace operator. To change or insert single elements of a compressed matrix is not supported. Especially for very large matrices, this would result in an unbearable performance burden.

However, it is allowed to assign a scalar value to the entire matrix given it is square as in the example. Matrix b is stored in compressed column storage (CCS).

Which orientation is favorable dependents on the performed operations and might require some experimentation. All operations are provided in the same way for both formats

All matrices have free functions for the number of rows and columns and the matrix size, which is understood as the product of the former and not the number of non-zeros.

To find out the number of rows use

unsigned r= num_rows(A);

It returns an unsigned integer (more precisely the size_type of the matrix type). Likewise the number of columns is given

unsigned c= num_cols(A);

The matrix size is given by

unsigned s= size(A); assert (s == r * c);

and is defined as product of the numbers of rows and columns. These definitions are consistent with the according functions for vectors (Vector Types).

The matrices can take a second (in case of matrix::morton_dense a third) template argument that allows for certain specialization. The argument should be an instance of the template class matrix::parameters. There are the following arguments:

- Orientation: whether the matrix is row_major or column-major (col_major), diagonal orientation might be added later, default is row_major;
- Index: was intended for easier handling of 1-based Fortran-like indexing but turned out to be too error-prone and unreadable in a generic context; might be removed in the future;
- Dimensions: non_fixed::dimensions for giving the dimension at run time or fixed::dimensions to specify during compilation, default is non_fixed::dimensions;
- OnStack: whether data should be stored on stack or on heap; storage on stack requires fixed::dimensions (although some compilers like g++ also support arrays with run-time sizes which is not standard-compliant); default is true for static dimensions and false otherwise;
- SizeType: the type for storing indices; default is std::size_t.

Some arguments in certain matrix types or have little impact, for instance:

- Orientation is ignored in matrix::morton_dense since the layout is determined by the mask;
- OnStack is ignored in matrix::compressed2D and fixed::dimensions reduce the overall memory need only marginally;
- Reversely, the choice of SizeType has no effect on the performance of dense matrices are very little on their memory requirements.

Using only 32 bit integers instead of 64 bit can accelerate sparse matrix operations significantly because twice as much indices can be loaded from memory at the same time (and as we all know, memory bandwidth is the limiting factor in sparse algebra), see Reducing the Size Type. Multiple operations are specialized for dense matrices with fixed dimensions, see Using Fixed-size Matrices and Vectors.

How to fill sparse matrices is shown on page Matrix Insertion.

As a new matrix type we have an element structure. The use of this new structure is to be explained with the following 2 by 3 grid.

#include <iostream> #include <boost/numeric/mtl/mtl.hpp> using namespace std; int main(int, char**) { typedef double value_type; typedef int size_type; const int nb_elements= 6, nb_nodes= 12; value_type array[][4]= {{2, 3, 4, 5}, {4, 10, 13, 16}, {6, 25, 38, 46}, {8, 32, 77, 100}}; mtl::dense2D<value_type> E_mat(array); std::cout<<"E_mat=\n"<< E_mat <<"\n"; mtl::matrix::element_structure<value_type> A; typedef mtl::matrix::element<value_type> element_type; element_type* elements = new element_type[nb_elements]; mtl::dense_vector<size_type> index_a(4, 0), index_b(4, 0), index_c(4, 0), index_d(4, 0), index_e(4, 0), index_f(4, 0); // construct nodes for every element index_a[0]= 0; index_a[1]= 1; index_a[2]= 4; index_a[3]= 5; index_b= index_a + 1; index_c= index_a + 2; index_d= index_a + 4; index_e= index_a + 5; index_f= index_a + 6; //construct the 6 elements from the example grid element_type a(0, index_a, E_mat); element_type b(1, index_b, E_mat); element_type c(2, index_c, E_mat); element_type d(3, index_d, E_mat); element_type e(4, index_e, E_mat); element_type f(5, index_f, E_mat); //construct neighborhood information for each element a.add_neighbors(&b, &d, &e); b.add_neighbors(&a, &c, &d, &e, &f); c.add_neighbors(&a, &b, &e); d.add_neighbors(&a, &b, &e); e.add_neighbors(&a, &b, &c, &d, &f); f.add_neighbors(&b, &c, &e); std::cout<< "a=" << a << "\n"; //construct array of elements elements[0]=a; elements[1]=b; elements[2]=c; elements[3]=d; elements[4]=e; elements[5]=f; //construct element_structure from the 6 single elements A.consume(nb_elements, nb_nodes, elements); mtl::dense_vector<value_type> x(nb_nodes, 1.0), test(A * x); std::cout<< "test="<< test << "\n"; return 0; }

Return to Vector Types Table of Content Proceed to Type Multivector

*
*