MTL 4: Matrix Insertion
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.
// 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 argc, char* argv[]) { // 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.
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 argc, char* argv[]) { // 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).
// File: array_initialization.cpp #include <iostream> #include <boost/numeric/mtl/mtl.hpp> int main(int argc, char* argv[]) { 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
-- Generated on 24 Aug 2009 by Doxygen 1.5.9 -- Copyright 2008-09 by TU Dresden and the Trustees of Indiana University.