The Matrix class is declared by:
template <typename Real>
class Matrix;
Operator overloading is essential in being able to imitate the brief and coincise notation of mathematics. Using operator overloading, one gets to write A * B in infix notation rather than multiply(A, B) in functional notation. Pastel overloads the operators for the matrix class in a way that resembles mathematical notation.
There are several useful ways to view the data in matrices. We enumerate them as follows:
Matrix as a 2d-array of numbers. Thus an individual element (i, j) of the matrix can be referenced. Also, the underlying Array object can be accessed as a non-mutable reference when that is needed for a function interface.
Matrix as a 1d-array of numbers (in row-major order). Via iterator abstraction, elements can be accessed as a sequence of values. Specifically, standard library algorithms can be used to read and modify matrix elements.
Matrix as a collection of row vectors. One can handle a matrix row as a vector.
Matrix as a collection of column vectors. One can handle a matrix column as a vector.
Matrix as a matrix expression. This makes it automatically work with whatever new matrix expression are introduced in the future.
The main technique in making matrix computations efficient in both space and time is the use of expression templates. Consider a matrix expression such as
A := A * (B + D)
Using a naive implementation, the result of B + D is computed in its entirety, requiring a temporary matrix E to hold the result. Next the product A * E is computed whose result again requires a new temporary matrix F. Finally F is assigned to A. Because the matrices can be of considerable size in memory, such a strategy is wasteful in memory space. There is also a problem with efficiency. First, the dynamic memory allocations can be slow. Second, some computations could be computed much faster if they could be combined (inlined) into same place.
These problems can be solved elegantly by so-called expression templates. Rather than actually carrying out the computation, one constructs an expression tree describing the computation. Every expression then represents an abstract matrix which exists only as a computational description. Such a description allows, for each element (i, j), to compute the value of the expression. Specifically, an expression tree allows to evaluate an expression breadth-first rather than depth-first. This lazy evaluation approach to computation allows one to avoid the temporary matrix E. Moreover, the expression tree is a compile-time construct which after code generation is for the most part removed by the compiler’s optimizer.
Because both sides of our example assignment contains A, the expression on the right necessarily needs to be evaluated to a temporary matrix before being assigned to A: otherwise the assigned values could possibly affect the computation on the right. Pastel automatically takes care of this temporary evaluation problem by inspecting the expression tree of the assignment parameters for self-references. However, even if there is a self-reference in the assigned expression, creation of a temporary can be avoided if the reference is trivial (i.e. to compute value for (i, j), it only uses the element (i, j) of itself). For example, the following assignment does not generate temporaries at all:
A := A + B + (2 * A)
Matrix expressions can also be used other things than to implement fast operations between matrices. For example, the identity matrix is implemented as an expression which returns value 1 for the diagonal positions and 0 for the other positions. Assigning such an expression to a matrix allows to initialize the matrix to an identity matrix without actually creating a temporary one in memory.
The matrix data can be stored in either the row-major, or the column-major order. By default, the data is stored using the the row-major convention. This means that the rows are stored one after each other in ascending order, and the elements in each row are stored in ascending order in column position. For example, the elements of a 2x2 matrix are stored in the order (0, 0); (0, 1); (1, 0); (1, 1). This convention corresponds to the convention that programming language C (and C++) uses with its multi-dimensional tables. The column-major convention is useful when interfacing with Matlab; there the storage convention is column-major.
Two matrix classes can share the same data. This is called aliasing. Because we know the storage convention, an aliasing matrix can also reshape its view of the data area somewhat. Specifically, the aliasing matrix can choose to skip some rows from the start and from the end. However, the aliasing matrix can not choose to skip any columns. The aliasing matrix will not attempt to free the memory after it is destroyed. On the other hand, if the aliased matrix is destroyed before the aliasing matrix, there is a danger of accessing unallocated memory. One must be wary for this not to happen.