Logo MTL4

Recursion

Recursion is an important theme in MTL4. Besides matrices with recursive recursive memory layout -- cf. Matrix Types and morton_dense -- recursion with regard to algorithms plays a decisive role.

To support the implementation of recursive algorithms we introduced -- in collaboration with David S. Wise -- the concept to Recursator, an analogon of Iterator. The class matrix::recursator enables recursive subdivision of all matrices with a sub_matrix function (e.g., dense2D and morton_dense). We refrained from providing the sub_matrix functionality to compressed2D; this would possible but very inefficient and therefor not particularly useful. Thus matrixrecursator of compressed2D cannot be declared. A recursator for vectors is planned for the future.

Generally spoken, the matrix::recursator (cf. recursion::matrix::recursator) consistently divides a matrix into four quadrants

with the self-evident cartographic meaning (from here on we abreviate matrix recursator to recursator). The quadrants itself can be sub-divided again providing the recursive sub-division of matrices into scalars (or blocks with user-defined maximal size).

The following program illustrates how to divide matrices via recursator:

// File: recursator.cpp

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

int main(int argc, char* argv[])
{
    using namespace mtl; using std::cout;

    // Z-order matrix
    typedef morton_dense<double, recursion::morton_z_mask>  matrix_type;
    matrix_type                                             A(10, 10);
    matrix::hessian_setup(A, 3.0);

    // Define a recursator over A
    matrix::recursator<matrix_type>                          rec(A);

    // Access a quadrant of the matrix
    cout << "Upper right quadrant (north_east) of A is \n" << *north_east(rec) << "\n";

    // Access a quadrant's quadrant of the matrix
    cout << "Lower left (south_west) of upper right quadrant (north_east) of A is \n" 
         << *south_west(north_east(rec)) << "\n";

    cout << "The virtual bound of 'rec' is " << rec.bound() << "\n";

    return 0;
}

The functions north_west(), north_east(), south_west(), and south_east() return recursators that refer to sub-matrices. The sub-matrices can be accessed by dereferring the recursator, i.e. *rec. Only then a sub-matrix is created.

As the example shows, the quadrant (represented by a recursator) can be sub-divided further (returning another recursator). Block-recursive algorithms can be implemented efficiently by sub-dividing large matrices recursively into blocks of decreasing size until a block size is reached that allows efficient iterative treatment. Sub-matrices are only created at the base case and not during the recursive descent because the creation of sub-matrix might be a relatively expensive operation (e.g., with morton_dense) while the creation of a new recursator requires only a few integer operations.

The recursator uses internally a virtual bound that is a power of 2 and at least as large as the number of rows and columns. In the example, the bound is 16 (as shown by the member function bound). When computing a quadrant the bound is halved and the starting row and column are potentially increased. For instance, the north_east quadrant is a virtual 8 by 8 matrix starting at row 0 and column 8. The sub-matrix referred by the north_east recursator is the intersection of this virtual quadrant with the original matrix A, i.e. an 8 by 2 matrix starting in row 0 and column 8.

More functionality of recursators is shown in the following example:

// File: recursator2.cpp

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

int main(int argc, char* argv[])
{
    using namespace mtl; using std::cout;

    typedef morton_dense<double, recursion::morton_z_mask>  matrix_type;
    matrix_type                                             A(10, 10);
    matrix::hessian_setup(A, 3.0);
    matrix::recursator<matrix_type>                          rec(A);

    // Create a recursator for the north_east quadrant of A
    matrix::recursator<matrix_type>                          ne(north_east(rec));

    cout << "Test if recursator 'ne' refers to an empty matrix (shouldn't): " << is_empty(ne) << "\n";
    cout << "Test if north_east of 'ne' refers to an empty matrix (it should): " << is_empty(north_east(ne)) << "\n";

    cout << "Number of rows and columns of north_east quadrant is: " << num_rows(ne)
         << " and " << num_cols(ne) << "\n";

    cout << "Test if 'ne' fills ils virtual quadrant (shouldn't): " << is_full(ne) << "\n";
    cout << "Test if north_west fills its virtual quadrant (it should): " << is_full(north_west(rec)) << "\n";

    return 0;
}

The function is_empty applied on a recursator computes whether the referred sub-matrix is empty, i.e. the intersection of the virtual quadrant and the original matrix A is empty. The sub-matrix itself is not generated since this test can be performed from size and index information. In the same way, number of rows and columns of the referred sub-matrix can be computed without its creation.

The function is_full() comes in handy in block-recursive algorithms. Assume we have a base case of 64 by 64, i.e. matrices with at most 64 rows and columns are treated iteratively. Then it is worthwile to write a blazingly fast iterative implementation for 64 by 64 matrices, in other words when the sub-matrix fills the entire virtual quadrant (when bound is 64). Thus, the function is_full() can be used to dispatch between this optimized code and the (hopefully not much slower) code for smaller matrices.

Return to Iteration                                Table of Content                                Proceed to Why and How we use Functors


Recursion -- MTL 4 -- Peter Gottschling and Andrew Lumsdaine -- Generated on 24 Aug 2009 by Doxygen 1.5.9 -- Copyright 2008-09 by TU Dresden and the Trustees of Indiana University.