Logo MTL4
Iteration

Iterative traversal of collections is implemented in MTL4 in two ways:

The latter is more general and allows especially for sparse structures a cleaner abstraction. Initially MTL4 was implemented entirely with this paradigm but it has shown that algorithms running exclusively on dense structures are easier to implement in terms of iterators.

All cursors and iterators are handled by:

They are all templated by a tag that determines the form of traversal. The following tags are currently available for cursors:

For iterators:

And finally for constant iterators:

Let's consider cursors in more detail.

Cursors

The approach was proposed by David Abrahams in order to separate the form of traversal from the manner of access. A cursor is a tool that can be used to visit different objects of a collection. In an array it can be compared with a position rather than a pointer because it is not fixed how one accesses the values. The traversal is essential the same as with iterators, e.g.:

    for (Cursor cursor(begin(x)), cend(end(x)); cursor != cend; ++cursor)
       do_something(cursor);

We will come back to the type Cursor later (please be patient).

In order to have more flexibility we templatized the begin and end functions:

    for (Cursor cursor(begin<tag::all>(x)), cend(end<tag::all>(x)); cursor != cend; ++cursor)
       do_something(cursor);

This cursor for instance goes over all elements of a matrix or vector, including structural zeros.

Nested Cursors

Several cursors can be used to create other cursors. This is necessary to traverse multi-dimensional collections like matrices. In most cases you will use nested cursors via the tags tag::row and tag::col. The returned cursor can be a certain collection (e.g. a vector) or just a place-holder that only contains some index and reference to a collection but cannot be used directly in operations. If the type and orientation permits, one can access the elements with tag::all or tag::nz, e.g.:

    for (Cursor cursor(begin<tag::row>(x)), cend(end<tag::row>(x)); cursor != cend; ++cursor)
       for (ICursor icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor)
           do_something(icursor);

Often it is more efficient to adapt an algorithm to the orientation of a matrix. Then it is convenient to use tag::major instead of dispatching for row-major and column major matrices:

    for (Cursor cursor(begin<tag::major>(x)), cend(end<tag::major>(x)); cursor != cend; ++cursor)
       for (ICursor icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor)
           do_something(icursor);

Property Maps

The concept of property maps has not only the advantage to allow for different forms of accessibility of values but also to provide different views or details of this value. Matrices have four property maps:

They are all accessed by dereferenced cursors, e.g.

    for (Cursor cursor(begin<tag::nz>(x)), cend(end<tag::nz>(x)); cursor != cend; ++cursor)
        cout << "matrix[" << row(*cursor) << ", " << col(*cursor) << "] = " 
             << const_value(*cursor) << '\n';

Three of the property maps are constant (guess which). Obviously only value can be changed. The syntax is the following:

    value(*cursor, 7);

Range Generator

The type traits traits::range_generator<Tag, Collection> is used to determine the type of cursor:

    typedef typename traits::range_generator<tag::row, Matrix>::type c_type;
    typedef typename traits::range_generator<tag::nz, c_type>::type  ic_type;

    for (c_type cursor(begin<tag::row>(x)), cend(end<tag::row>(x)); cursor != cend; ++cursor)
       for (ic_type icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor)
           do_something(icursor);

As can be seen in the examples, cursors that represents sub-collections (e.g. rows) can be used as collection type.

Iterators

In some contexts, especially with dense data only, iterators are simpler to use. With the property map syntax, one cannot apply operators like += or a modifying function. Therefore we provide iterators for dense matrices and vectors. For sparse matrices there was no use case so far because iterators do not reveal which matrix element they are pointing at.

The usage of iterators is very similar to those of cursors:

    for (Iter iter(begin<tag::const_iter::nz>(x)), iend(end<tag::const_iter::nz>(x)); 
         iter != iend; ++iter)
        cout << "matrix value = " << *iter << '\n';

In contrast to the previous examples we can only output the value without the indices. The type of Iter can be determined with range_generator in the same way.

Nested Iterators

Nesting of iterators is also analog to cursors. However, iterators only exist to access elements not sub-collections. The nesting is therefore realized by mixing cursors and iterators.

    for (Cursor cursor(begin<tag::major>(x)), cend(end<tag::major>(x)); cursor != cend; ++cursor)
        for (Iter iter(begin<tag::const_iter::nz>(cursor)), iend(end<tag::const_iter::nz>(cursor)); 
             iter != iend; ++iter)
            cout << "matrix value = " << *iter << '\n';

In the example we iterate over the rows by a cursor and then iterate over the elements with an iterator.

Advanced topic: Choosing traversal by complexity

Range generators in MTL4 have a notion of complexity. That is for a given collection and a given form of traversal it can be said at compile time which complexity this traversal has.

Dense matrices are traversed with linear or cached_linear complexity. The latter is used for contiguous memory access over strided ones, which is also linear but considerably slower. This distinction is mathematically questionable but useful in practical contexts.

Sparse matrices have linear complexity when traversed along the orientation. Traversing compressed matrices perpendicular to the orientation (e.g. a CRS matrix column-wise) has infinite complexity because it is not implemented. Moreover, the default (non-spezialized) range_generator has infinite complexity so that it is per se defined for arbitrary collections and tags. Whether the range generator is actually really implemented can be tested by comparing the complexity with infinite (by using MPL functions).

The following example shows a simpler way to find out the best traversal:

// MTL4 example: minimize complexity of traversal

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

using namespace mtl;
    
template <typename Matrix>
void f(Matrix& m)
{
    using traits::range_generator; using traits::range::min;

    // Set values in diagonal
    m= 7.0;
    
    // Choose between row and column traversal
    typedef typename min<range_generator<tag::row, Matrix>, range_generator<tag::col, Matrix> >::type range_type;
    range_type                                                      my_range;

    // Type of outer cursor
    typedef typename range_type::type                               c_type;
    // Type of inner cursor
    typedef typename traits::range_generator<tag::nz, c_type>::type ic_type;

    // Define the property maps
    typename traits::row<Matrix>::type                              row(m); 
    typename traits::col<Matrix>::type                              col(m);
    typename traits::const_value<Matrix>::type                      value(m); 

    // Now iterate over the matrix    
    for (c_type cursor(my_range.begin(m)), cend(my_range.end(m)); cursor != cend; ++cursor)
       for (ic_type icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor)
           std::cout << "matrix[" << row(*icursor) << ", " << col(*icursor) << "] = " << value(*icursor) << '\n';    
}


int main(int, char**)
{
    // Define a CRS matrix
    compressed2D<double>                                  A(3, 3);
    // And a CCS matrix
    compressed2D<double, matrix::parameters<col_major> >  B(3, 3);

    f(A);
    f(B);
    
    return 0;
}

Please not that the example uses compressed sparse matrices and not all forms of traversion are supported. Obviously a linear complexity is lower than an infinite and the range generator without implemenation is never used. As the free functions begin() and end() are internally always implemented by member functions of range_generator (free template functions cannot be spezialized partially) we used directly the member functions in the example.

The range generator can also be minimized recursively between three and more alternatives:

    typedef typename min<range_generator<tag::row, Matrix>, 
                         typename min<range_generator<tag::col, Matrix>,
                                      range_generator<tag::major, Matrix> >::type 
                        >::type range_type;

In many cases there is no need for explicitly minimizing the complexity because tag::major usually will yield the same results (but this is not so cool).

Return to Using Predefined Linear Solvers                                Table of Content                                Proceed to Recursion

-----------------------------------------------------------

/*!


Iteration -- 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.