Logo MTL4
Matrix Expressions

The following program illustrates how to add matrices, including scaled matrices:

#include <iostream>
#include <boost/numeric/mtl/mtl.hpp>

int main(int, char**)
{
    using namespace mtl;
    
    const unsigned n= 10;
    compressed2D<double>                         A(n, n);
    dense2D<int, matrix::parameters<col_major> > B(n, n);
    morton_dense<double, 0x555555f0>             C(n, n), D(n, n);

    matrix::laplacian_setup(A, 2, 5);
    matrix::hessian_setup(B, 1); matrix::hessian_setup(C, 2.0); matrix::hessian_setup(D, 3.0);

    D+= A - 2 * B + C;

    std::cout << "The matrices are: A=\n" << A << "B=\n" << B << "C=\n" << C << "D=\n" << D;

    return 0;
}

The example shows that arbitrary combinations of matrices can be added, regardless their orientation, recursive or non-recursive memory layout, and sparseness.

Matrix multiplication can be implemented as elegantly:

#include <boost/numeric/mtl/mtl.hpp>

int main(int, char**)
{
    using namespace mtl; using namespace mtl::matrix;
    
    const unsigned n= 20;
    dense2D<double>                            A(n, n), B(n, n);
    morton_dense<double, doppled_64_row_mask>  C(n, n);

    hessian_setup(A, 3.0); hessian_setup(B, 1.0); 
    hessian_setup(C, 2.0);

    // Corresponds to A= B * B;
    mult(B, B, A);

    A= B * B;   // use BLAS
    A= B * C;   // use recursion + tiling from MTL4

    A+= B * C;  // Increment A by the product of B and C
    A-= B * C;  // Likewise with decrement

    return 0;
}

Arbitrary matrix types can be multiplied in MTL4. Let's start with the operation that is the holy grail in high-performance computing: dense matrix multiplication. This is also the operation shown in the example above. The multiplication can be executed with the function mult where the first two arguments are the operands and the third the result. Exactly the same is performed with the operator notation below.

Warning: the arguments and the result must be different! Expressions like A= A*B will throw an exception. More subtle aliasing, e.g., partial overlap of the matrices might not be detected and result in undefined mathematical behavior.

Products of three matrices are supported now. Internally they are realized by binary products creating temporaries (thus, sequences of two-term products should provide better performance). Moreover, products can be arbitrarily added and subtracted:

#include <boost/numeric/mtl/mtl.hpp>

int main(int, char**)
{
    using namespace mtl; using namespace mtl::matrix;
    
    const unsigned n= 40;
    dense2D<double>                            A(n, n), B(n, n);
    morton_dense<double, doppled_64_row_mask>  C(n, n), D(n, n);

    hessian_setup(A, 3.0); hessian_setup(B, 1.0); 
    hessian_setup(C, 2.0); hessian_setup(D, 11.0);

    A+= B * B + C * B - B * B * C * D;   

    return 0;
}

Return to Rich Vector Expressions                                Table of Content                                Proceed to Matrix-Vector Expressions


Matrix Expressions -- MTL 4 -- Peter Gottschling and Andrew Lumsdaine -- Gen. with rev. 7542 on Sat Aug 11 2012 by doxygen 1.7.6.1 -- © 2010 by SimuNova UG.