MTL 4: Matrix Types
The type dense2D defines regular row-major and column-major matrices:
// File: dense2D.cpp #include <iostream> #include <boost/numeric/mtl/mtl.hpp> int main(int, char**) { using namespace mtl; // A is a row-major matrix dense2D<double> A(10, 10); // Matrices are not initialized by default A= 0.0; // Assign a value to a matrix element A(2, 3)= 7.0; // You can also use a more C-like notation A[2][4]= 3.0; std::cout << "A is \n" << A << "\n"; // B is a column-major matrix dense2D<float, matrix::parameters<tag::col_major> > B(10, 10); // Assign the identity matrix times 3 to B B= 3; std::cout << "B is \n" << B << "\n"; return 0; }
If no matrix parameters are defined, dense matrices are by default row-major. There are more matrix parameters besides the orientation. As they are not yet fully supported we refrain from discussing them.
Matrix elements can be accessed by a(i, j) or in the more familiar form a[i][j]. The second form is internally transformed into the first one at compile time so that the run-time performance is not affected (unless the compiler does not inline completely which we never observed so far). Also, the compile time is not conceivably increased by this transformation.
Please notice that overwriting single matrix elements is only defined for dense matrix types. For a generic way to modify matrices see Matrix Insertion.
Assigning a scalar value to a matrix stores a multiple of the identity matrix, i.e. the scalar is assigned to all diagonal elements and all off-diagonal elements are 0. Diagonal elements are matrix entries with identical row and column index. Therefore scalars can also be assigned to non-square matrices. This operation is generic (i.e. applicable to all matrix types including sparse).
Just in case you wonder why the scalar value is only assigned to the diagonal elements of the matrix not to all entries, this becomes quite clear when you think of a matrix as a linear operator (from one vector space to another one). For instance, consider the multiplication of vector x with the scalar alpha:
y= alpha * x;
A= alpha; y= A * x;
If the matrix is not square, i.e. the linear operator's domain and image have different dimensions, the equivalence with the scalar multiplication applies accordingly. In case that the image has a lower dimension, say m, then only the first m entries of the vector from the domain are scaled with alpha and the others are ignored. If the image has an higher dimension then the last m-n entries are zero with n the dimension of the domain. When you rely on this behavior please check the revision of your MTL4 library: old versions, i.e. before revision 6843, considered it erroneous to store a non-zero scalar to a non-square matrix.
From a more pragmatic prospective:
A= 0;
Dense matrices with a recursively designed memory layout can be defined with the type morton_dense:
// File: morton_dense.cpp #include <iostream> #include <boost/numeric/mtl/mtl.hpp> int main(int, char**) { using namespace mtl; // Z-order matrix morton_dense<double, recursion::morton_z_mask> A(10, 10); A= 0; A(2, 3)= 7.0; A[2][4]= 3.0; std::cout << "A is \n" << A << "\n"; // B is an N-order matrix with column-major 4x4 blocks, see paper morton_dense<float, recursion::doppled_4_col_mask> B(10, 10); // Assign the identity matrix times 3 to B B= 3; std::cout << "B is \n" << B << "\n"; return 0; }
In the pure Morton order format 2 by 2 sub-matrices are stored contiguously in memory. 4 by 4 matrices constitute of 4 2-by-2-matrices and use consecutive memory. The continuation of this recursive scheme provides square sub-matrices with power of two sizes that are in contiguous memory and allow for cache-efficient recursive algorithms. On the other hand, algorithms that are implemented fully recursively create considerable overhead for function calls. We therefore recommend using mixed schemes of recursion and iteration. Particularly efficient are algorithms that operate on mid-size blocks, e.g. 64 by 64, with regular row-major or column-major layout. MTL4 provides a virtually infinite number of memory layouts for dense matrices that are specified by a bitmask. A detailed description and discussion of recursive matrices and algorithm is provided in this conference paper.
Sparse matrices are defined with the type compressed2D:
// File: compressed2D.cpp #include <iostream> #include <boost/numeric/mtl/mtl.hpp> int main(int, char**) { using namespace mtl; // CRS matrix compressed2D<double> A(12, 12); // Laplace operator discretized on a 3x4 grid matrix::laplacian_setup(A, 3, 4); std::cout << "A is \n" << A; // Element access is allowed for reading std::cout << "A[3][2] is " << A[3][2] << "\n\n"; // CCS matrix compressed2D<float, matrix::parameters<tag::col_major> > B(10, 10); // Assign the identity matrix times 3 to B B= 3; std::cout << "B is \n" << B << "\n"; return 0; }
Matrix a is stored as compressed row storage (CRS). Its assigned values correspond to a discretized Laplace operator. To change or insert single elements of a compressed matrix is not supported. Especially for very large matrices, this would result in an unbearable performance burden.
However, it is allowed to assign a scalar value to the entire matrix given it is square as in the example. Matrix b is stored in compressed column storage (CCS).
Which orientation is favorable dependents on the performed operations and might require some experimentation. All operations are provided in the same way for both formats
All matrices have free functions for the number of rows and columns and the matrix size, which is understood as the product of the former and not the number of non-zeros.
To find out the number of rows use
unsigned r= num_rows(A);
unsigned c= num_cols(A);
unsigned s= size(A); assert (s == r * c);
How to fill sparse matrices is shown in section Matrix Insertion.
Return to Vector Types Table of Content Proceed to Type Multivector
Matrix Types -- MTL 4 -- Peter Gottschling and Andrew Lumsdaine
-- Gen. with
rev. 7542
on 7 Apr 2011 by doxygen 1.5.9 -- © 2010 by SimuNova UG.