Logo MTL4

Matrix Insertion

Setting the values of a dense matrix is an easy task since each element has its dedicated location in memory. Setting sparse matrices, esp. compressed ones, is a little more complicated. There exist two extreme approaches:

The former approach has the advantage that it is handier and that the set-up of sparse matrices can be handled like dense matrices (which eases the development of generic code). However, when matrices grow larger, the insertion becomes more and more expensive, up to the point of being unusable. Most high-performance libraries use therefore the second approach. In practice, a sparse matrix is usually the result of discretization (FEM, FDM, ...) that is set up once and then used many times in linear or non-linear solvers. Many libraries even establish a two-phase set-up: first building the sparsity pattern and then populating the non-zero elements with values.

The MTL4 approach lies somewhere between. Sparse matrices can be either written (inserted) or read. However, there can be multiple insertion phases.

Element-wise Insertion

Before giving more details, we want to show you a short example:

// File: insert.cpp

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

using namespace mtl;

template <typename Matrix>
void fill(Matrix& m)
{
    // Matrices are not initialized by default
    m= 0.0;

    // Create inserter for matrix m
    matrix::inserter<Matrix> ins(m);
    
    // Insert value in m[0][0]
    ins[0][0] << 2.0;
    ins[1][2] << 0.5;
    ins[2][1] << 3.0;

    // Destructor of ins sets final state of m
}

template <typename Matrix>
void modify(Matrix& m)
{
    // Type of m's elements
    typedef typename Collection<Matrix>::value_type value_type;

    // Create inserter for matrix m
    // Existing values are not overwritten but inserted
    matrix::inserter<Matrix, update_plus<value_type> > ins(m, 3);
    
    // Increment value in m[0][0]
    ins[0][0] << 1.0;

    // Elements that doesn't exist (in sparse matrices) are inserted
    ins[1][1] << 2.5;
    ins[2][1] << 1.0;
    ins[2][2] << 4.0;

    // Destructor of ins sets final state of m
}

int main(int, char**)
{
    // Matrices of different types
    compressed2D<double>              A(3, 3);
    dense2D<double>                   B(3, 3);
    morton_dense<float, morton_mask>  C(3, 3);

    // Fill the matrices generically
    fill(A); fill(B); fill(C);
    std::cout << "A is \n" << A << "\nB is \n" << B << "\nC is \n" << C;

    // Modify the matrices generically
    modify(A); modify(B); modify(C);
    std::cout << "\n\nAfter modification:\nA is \n" << A 
              << "\nB is \n" << B << "\nC is \n" << C;

    return 0;
}

The first aspect worth pointing at is that sparse and dense matrices are treated the same way. If we cannot handle sparse matrices like dense (at least not efficiently), we can treat dense matrices like sparse ones. For performance reasons, matrices are not initialized by default. Therefore, the first operation in the function fill is to set the matrix to zero.

Internally the inserters for dense and sparse matrices are implemented completely differently but the interface is the same. Dense inserters insert the value directly and there is not much to say about.

Sparse inserters are more complicated. The constructor stretches the matrix so that the first five elements in a row (in a CCS matrix likewise the first 5 elements in a column) are inserted directly. During the live time of the inserter, new elements are written directly into empty slots. If all slots of a row (or column) are filled, new elements are written into an std::map. During the entire insertion process, no data is shifted.

If an element is inserted twice then the existing element is overwritten, regardless if the element is stored in the matrix itself or in the overflow container. Overwriting is only the default. The function modify() illustrates how to use the inserter incrementally. Existing elements are incremented by the new value. We hope that this ability facilitates the development of FEM code.

For performance reasons it is advisable to customize the number of elements per row (or column), e.g., ins(m, 13). Reason being, the overflow container consumes more memory per element then the regular matrix container. In most applications, an upper limit can be easily given. However, the limit is not that strict: if some rows need more memory than the slot size it only results in slightly higher memory need for the overflow container. If the number of elements per row is very irregular we recommend a slot size over the average (and maybe under the maximum). Since only a small part of the data is copied during the compression, sparse matrices can be created that fill almost the entire memory.

Nota bene: inserters for dense matrices are not much more than facades for the matrices themselves in order to provide the same interface as for sparse ones. However, dense inserters can be also very useful in the future for extending the library to parallel computations. Then the inserter can be used to write values into remote matrix elements.

INSERTERS MUST BE DESTROYED

A mistake that many people did with inserters was using the matrix before the inserter was destroyed, e.g.:
using namespace mtl;
typedef compressed2D<double> matrix_type;

matrix_type A(5, 5);
matrix::inserter<matrix_type> ins(A);
ins[0][0] << 7.3; // .... more insertions

do_something_with(A);  // TROUBLE!!!
Then the matrix A is not ready and an exception is thrown (or an assertion fails depending on compile flags).

The issue was apparently not sufficiently discussed in the tutorial and we have to blamed not the users for doing this wrong.

The insertion problem is circumvented by defining the inserter in a separate function as we did in the previous section. If we accessed the matrix within the fill-in function we would experience the same problem.

If no separate function shall be defined for brevity, one can define the inserter in an extra block. The following program implements the function "fill" of the example insert.cpp with a compressed matrix:

// Filename: insert_scope.cpp

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

using namespace mtl;

int main(int, char**)
{
    compressed2D<double>              A(3, 3);
    {
        matrix::inserter<compressed2D<double> > ins(A);  
        ins[0][0] << 2.0;
        ins[1][2] << 0.5;
        ins[2][1] << 3.0;
    } // ins is destroyed here

    std::cout << "A is \n" << A;  // we can access A now

    return 0;
}

Alternatively, one can handle the insertion destruction explicitly with pointer as will be explained later.

Block-wise Insertion

A more powerful method to fill sparse (and dense) matrices provide the two functions element_matrix() and element_array().

The following program illustrates how to use them:

// File: element_matrix.cpp

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

using namespace mtl;

template <typename Matrix>
void fill(Matrix& m)
{
    // Matrices are not initialized by default
    m= 0.0;

    // Type of m's elements
    typedef typename Collection<Matrix>::value_type value_type;

    // Create inserter for matrix m
    // Existing values are not overwritten but inserted
    matrix::inserter<Matrix, update_plus<value_type> > ins(m, 3);
    
    // Define element matrix (array)
    double m1[2][2]= {{1.0, -.4}, {-0.5, 2.0}}; 

    // Corresponding indices of the elements
    std::vector<int> v1(2);
    v1[0]= 1; v1[1]= 3;

    // Insert element matrix
    ins << element_array(m1, v1);

    // Insert same array with different indices
    v1[0]= 0; v1[1]= 2;
    ins << element_array(m1, v1);

    // Use element matrix type with dynamic size
    dense2D<double> m2(2, 3);
    m2[0][0]= 1; m2[0][1]= 0.2; m2[0][2]= 0.1; 
    m2[1][0]= 2; m2[1][1]= 1.2; m2[1][2]= 1.1;

    // Vector for column indices 
    dense_vector<int> v2(3);
    // Indices can be out of order
    v2[0]= 4; v2[1]= 1; v2[2]= 3;

    // Use element_matrix and separate vectors for row and column indices
    ins << element_matrix(m2, v1, v2);
}

int main(int, char**)
{
    // Matrices of different types
    compressed2D<double>              A(5, 5);
    dense2D<double>                   B(5, 5);
    morton_dense<float, morton_mask>  C(5, 5);

    // Fill the matrices generically
    fill(A); fill(B); fill(C);
    std::cout << "A is \n" << with_format(A, 4, 3) 
              << "\nB is \n" << with_format(B, 4, 3)
              << "\nC is \n" << with_format(C, 4, 3);

    return 0;
}

The function element_array is designed for element matrices that are stored as a 2D C/C++ array. The entries of such an element matrix are accessed by A[i][j], while the entries are accessed by A(i, j) if the function element_matrix is used. Element matrices stored in MTL4 types can be accessed both ways and either element_array or element_matrix can be used.

Both functions can be called with two or three arguments. In the former case the first argument is the element matrix and the second argument a vector containing the indices that correspond to the rows and columns of the assembled matrix. With three arguments, the second one is a vector of row indices and the third one a vector with column indices. Evidently, the size of the vector with the row/column indices should be equal to the number of rows/columns of the element matrix.

The vector type must provide a member function size and a bracket operator. Thus, mtl::dense_vector and std::vector can used (are models).

Insertion in Multiple Function calls

If a matrix is set up by means of multiple function calls as it happens often in finite element assembly. Say we have a class world_matrix that contains a compressed matrix which is set by calling add_entry several times. We can define an inserter in the function add_entry and the inserter will be destroyed after leaving the function:

// Filename: insert_class_expensive.cpp

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

using namespace mtl;

class world_matrix
{
  public:
    world_matrix(unsigned nrows, unsigned ncols) : A(nrows, ncols) {}

    void add_entry(unsigned row, unsigned col, double value) 
    {
        // Extremely expensive -> must not be done 
        matrix::inserter<compressed2D<double>, update_plus<double> > ins(A); 
        ins[row][col] << value;
    }

    friend inline std::ostream& operator<<(std::ostream& os, const world_matrix& w)
    { return os << w.A; }

  private:
    compressed2D<double>              A;
};


int main(int, char**)
{
    world_matrix              A(3, 3);

    A.add_entry(0, 0, 2.0);
    A.add_entry(1, 2, 0.5);
    A.add_entry(2, 1, 3.0);

    std::cout << "A is \n" << A;  // we can access A now

    return 0;
}

This works correctly but it is horribly slow because every value inserted need the creation of an inserter which is extremely expensive.

Defining an inserter as member of the class does not work at all because the inserter will live as long as the containing object and the matrix cannot be accessed during its life time.

The solution is to define a pointer to an inserter as member. The pointer lives as long as the object and the life time of the referred inserter can be controlled manually with new and delete. This said, we need a function that allocates the pointer that must be called before starting the insertion. Accordingly, a function is required that deallocates the pointer so that the inserter is destroyed. This function must be called after terminating the insertion phase and before accessing the matrix:

// Filename: insert_scope.cpp

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

using namespace mtl;

class world_matrix
{
    typedef matrix::inserter<compressed2D<double>, update_plus<double> > inserter_type;
  public:
    world_matrix(unsigned nrows, unsigned ncols) : A(nrows, ncols), ins(0) {}

    void add_entry(unsigned row, unsigned col, double value) 
    {  (*ins)[row][col] << value;    }

    void start_insertion() // must be called before first insertion
    { if (!ins) ins= new inserter_type(A); }

    void finish_insertion() // must be called before first usage
    {
        if (ins) {
            delete ins;
            ins= 0;
        }
    }

    friend inline std::ostream& operator<<(std::ostream& os, const world_matrix& w)
    { return os << w.A; }

  private:
    compressed2D<double>      A;
    inserter_type*            ins;
};


int main(int, char**)
{
    world_matrix              A(3, 3);

    A.start_insertion();
    A.add_entry(0, 0, 2.0);
    A.add_entry(1, 2, 0.5);
    A.add_entry(2, 1, 3.0);
    A.finish_insertion();

    std::cout << "A is \n" << A;  // we can access A now

    return 0;
}

Note that "ins" must be dereferred in add_entry. We find this approach more error-prone than defining the inserter in an extra scope but im some situations this is the only feasible way.

This approach is used in several software projects such as AMDiS and FEniCS.

Initializing Matrices with Arrays

For small matrices in examples it is more convenient to initialize the matrix from a 2D C/C++ array instead of filling it element-wise:

// File: array_initialization.cpp

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

int main(int, char**)
{
    using namespace mtl;

    double array[][4]= {{3, 7.2,   0, 6}, 
                        {2, 4.444, 5, 3},
                        {1, 7,     9, 2}};

    dense2D<double> A(array);

    std::cout << "A = \n" << A << "\n";

    return 0;
}

C/C++ arrays can be initialized be nested lists. All MTL4 matrices provide construction from arrays. Unfortunately, it is not (yet) possible to initialize user-defined types with lists. This is proposed for the next C++ standard and we will incorporate this feature as soon as it is generally available.

Return to Vector Insertion                                Table of Content                                Proceed to Vector Assignment


Matrix Insertion -- MTL 4 -- Peter Gottschling and Andrew Lumsdaine -- Gen. with rev. 7542 on 7 Apr 2011 by doxygen 1.5.9 -- © 2010 by SimuNova UG.