MTL 4: Recursion
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
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, char**) { 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, char**) { 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
-- Gen. with
rev. 7542
on 7 Apr 2011 by doxygen 1.5.9 -- © 2010 by SimuNova UG.