Logo MTL4

Matrix Types

Right now, MTL4 provides five matrix types:

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 Scalar Values to Matrices

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.

Recursive Memory Layout

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 :

Compressed Sparse Matrices

// 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).

Matrix Parameters

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:

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

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.

Matrix Insertion

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

Elementstructure

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.

2by3grid.png
#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


Matrix Types -- MTL 4 -- Peter Gottschling and Andrew Lumsdaine -- Gen. with rev. 9518 on 24 Apr 2014 by doxygen 1.6.3 -- © 2010-2013 by SimuNova UG.