Introduction
There are times when it is necessary to store data in a matrix so as to make programming easier. This is especially true in scientific programming where matrices are often used in calculations. However, there is no builtin Matrix type in C++.
There are several ways to implement a matrix in C++. The easiest involve arrays or pointers. For example, to create a simple matrix using an array, just do double matrix[10][10];
to create a variable ‘matrix’ that has 10 rows and 10 columns. Programming books will tell you that the data isn't exactly stored in a twodimensional form, but that is of no interest to us. What is important is that we can easily access the elements in the matrix; for example, to access the element in row 2 and column 5, we do matrix[1][4]
taking into account that arrays start from 0. For ease of explaining, for now on, all matrices start from zero, thus if I refer to row 2 and column 5, it is matrix[2][5]. The disadvantage of using arrays as matrices is that the size of the rows and columns must be fixed at compiletime.
To dynamically adjust the size of the matrix, we can create a matrix using pointers. We can use:
double* matrix = new double[100]
or
double** matrix = new double*[10];
for(int i=0; i<10;++i)
matrix[i] = new double[10];
The first form is a special type of matrix in that, to access an element in row i
and column j
, you need to do matrix[i*10+j]
. This is obviously not as intuitive as matrix[i][j]
. However, it has its own merits in that the matrix is simpler to create. The second form is slightly harder to create a matrix, but accessing an element is the same as in the array, i.e. matrix[i][j]
. The disadvantage of using pointers is that you must remember to call delete[]
on the matrix after using it, otherwise there will be a memory leak.
A third way of creating a matrix is to use STL containers like vectors. We can nest two vectors together, like vector<vector<double> > matrix
. To access an element, we can still use matrix[i][j]
. In addition, we can get the size of the matrix relatively easily, through the size
function of the vector container. We can also easily increase the size of the matrix by using the resize
function. Another advantage of using STL containers is it is not prone to memory leaks unlike pointers.
Mathematical operations are harder to perform on matrices implemented using arrays, pointers, or STL containers. For example, to add two matrices together, you can't write "matrixA + matrixB
". Instead, you must write two loops to add the individual elements together. To enable ease of calculations and also to provide useful functions for manipulating matrices, we can create a matrix class. Many programmers have created their own matrix classes. Among the more notable are the simple matrix class implemented for use in Numerical Recipes in C++^{1}, the Matrix Template Library (MTL) developed by Open Systems Laboratory at Indiana University^{2}, and the Matrix Expression Template implemented by Masakatsu Ito^{3}.
Since there are already excellent matrix classes available, why is there a need to create another one? There are several motivations for creating a new matrix class which is the focus of this article. First, the available matrix classes doesn't have all the features that I need. To add extra functions to these matrix classes will require an understanding of the implementation of the classes, which is timeconsuming. Second, as an amateur programmer, creating a new matrix class will improve my programming skills as I will have to consider the implementation problems thoroughly. The third reason is that I will not have to spend time learning about how to use the other matrix classes. Although some are quite straightforward, some are not. Implementing my own matrix will eliminate this learning curve.
The matrix class that I have implemented is not perfect. It is not meant to be the best matrix class. Before proceeding further, the deficiencies of the matrix class are listed below:
 Only dense matrix storage type has been implemented. Other storage types like sparse matrix, triangular matrix, diagonal matrix etc., have not been implemented.
 The performance of the matrix is not the fastest among the available matrix classes. Although I have not done extensive comparison tests, I believe that my matrix class is not the fastest in terms of calculations.
 The number of functions available for manipulating matrices is less than that in some available matrix classes.
 The most important deficiency is that this class can only be compiled by compilers that are compliant to the C++ standard. I have only tested the present version on Visual C++ 2005.
Other than the motivations for creating this matrix class, my matrix class hase advantages over other matrix classes. Otherwise, I would not have published this class. The advantages are:
 The matrix can be used to store numerical types as well as nonnumerical types. When storing nonnumerical types, the mathematical functions are not available for use, and will give compiler errors if accidentally used, instead of causing errors during runtime.
 The functions available for performing matrix calculations are intuitive to use, unlike those in some matrix classes.
 Although other storage types are not implemented, they can be created relatively easily by the user by following the
MatrixStorageTemplate
provided.
 STLlike functions and typedefs are available for users familiar with STL containers. Certain STL algorithms can also be performed on the matrix.
For those who have decided that this class may be suitable for them or who wishes to learn certain programming techniques like using policies with templates and expression templates, the rest of the article is divided into the following sections. The first section gives an overview of the functions available. The second section gives an example of how to use the matrix class. The third section discusses the implementation of the matrix class, giving details of using policies for generic programming and the use of expression templates to speed up certain matrix calculations.
Available Matrix Functions
The only class that the user needs to use is the template Matrix<ValueT,MatrixTypeP,MathP>
class. ValueT
is a type name for the type of data that will be stored in the matrix. It can be numerical types like int, double, or complex, or numerical types like char*
, or even userdefined classes. MatrixTypeP
is a template class policy for the data storage type. At the present moment, only DenseMatrix
is available. MathP
is a template class policy for mathematical functions. Three choices are available: MathMatrix
, MathMETMatrix
, and NonMathMatrix
.
No matter what type of matrix is created, the following are common for all types of matrixes:
STLlike typedefs:
value_type

type of data stored in matrix. 
reference

reference to type of data stored in matrix. 
const_reference

const reference to type of data stored in matrix. 
pointer

pointer to type of data stored in matrix. 
const_pointer

const pointer to type of data stored in matrix. 
difference_type

difference types for pointer calculations. 
size_type

size type. 
Constructors, destructor, copy constructors, and assignment:
Matrix<ValueT, MatrixTypeP, MathP>

Default constructor 
Matrix<ValueT, MatrixTypeP, MathP>(size_type sizeRow, size_type sizeCol, value_type x=value_type)

Constructor requiring number of rows, columns, and initial value for elements. 
Matrix<ValueT, MatrixTypeP, MathP>(size_type sizeRow, size_type sizeCol, value_type* am)

Constructor requiring number of rows, columns, and a 1D array containing initial values for elements. 
Matrix<ValueT, MatrixTypeP, MathP>(size_type sizeRow, size_type sizeCol, value_type** am)

Constructor requiring number of rows, and a 2D array containing initial values for elements. 
~Matrix<ValueT, MatrixTypeP, MathP>

Destructor. 
Matrix<ValueT, MatrixTypeP, MathP>(const Matrix<ValueT, MatrixTypeP, MathP>& m)

Copy constructor. 
template<typename ValueTT, template<typename> class MatrixTypePT, template<typename, typename> class MathPT>
Matrix(const Matrix<ValueTT, MatrixTypePT, MathPT>& m)

Copy constructor taking any matrix of interconvertible ValueT . 
Matrix<ValueT, MatrixTypeP, MathP>& operator=(Matrix<ValueT, MatrixTypeP, MathP> m)

Assignment operator taking the same type of matrix. 
template<typename ValueTT,template<typename> class MatrixTypePT,template<typename,typename> class MathPT>
Matrix<ValueT,MatrixTypeP,MathP>& operator=(Matrix<ValueTT,MatrixTypePT,MathPT> m)

Assignment operator taking any matrix of interconvertible ValueT . 
Matrix<ValueT, MatrixTypeP,MathP>& operator=(const value_type* am)

Copy a 1D array into the matrix. The number of array elements must be equal or larger than the number of matrix elements, and the sequence of array elements must be a concatenation of rows to be put into the matrix. 
Matrix<ValueT, MatrixTypeP, MathP>& operator=(const value_type** am)

Copy a 2D array into the matrix. The size of the array must be the same as the matrix. 
Matrix member functions
void clear

Erase matrix. 
bool empty const

Check if matrix is empty. 
void resize(size_type sizeRow, size_type sizeCol, value_type x=value_type)

Resize matrix, filling new elements with supplied or default values and removing extra elements. 
void swap(Matrix<ValueT, MatrixTypeP,MathP>& m)

Swap the contents of another matrix with the present matrix. The matrix to be swapped must be of the same type. 
void swaprows(size_type i1, size_type i2)

Swap two rows in the matrix. 
void swapcols(size_type j1, size_type j2)

Swap two columns in the matrix. 
bool operator==(const Matrix<ValueT,MatrixTypeP,MathP>& m) const
bool operator!=(const Matrix<ValueT,MatrixTypeP,MathP>& m) const

Comparison operators. 
void Update

Update matrix. 
const_reference operator(size_type posRow, size_type posCol)
const reference operator(size_type posRow, size_type posCol)
const_reference at(size_type posRow, size_type posCol)
const reference at(size_type posRow, size_type posCol)

Member accessor functions 
Row and Column iterator functions
size_type size() const 
Return the total number of rows/columns. 
template<typename ForwardIterator> void insert(size_type rowNo, const ForwardIterator& first) 
Insert row/column at specified row/column. 
void erase(size_type rowNo) 
Erase row/column at specific row/column. 
template<typename ForwardIterator> void push_back(const ForwardIterator& first) 
Add row/column at end of matrix. 
template<typename ForwardIterator> void push_front(const ForwardIterator& first) 
Add row/column at beginning of matrix. 
void pop_back() 
Remove row/column at end of matrix. 
void pop_front() 
Remove row/column at beginning of matrix. 
iterator begin()
const_iterator begin() const
iterator end()
const_iterator end() const
reverse_iterator rbegin()
const_reverse_iterator rbegin() const
reverse_iterator rend()
const_reverse_iterator rend() const

Functions to obtain iterators to data in matrix. 
Functors operating on Matrix
MatrixT Transpose<MatrixT>()(const MatrixT& matrix)
void Transpose<MatrixT>()(const MatrixT& matrix, MatrixT& transposedMatrix)

Get the transpose of a matrix, i.e., swap rows with columns and vice versa. The first form returns the transposed matrix. The second saves the transposed matrix into a variable which is passed in. 
MatrixT Diagonal<MatrixT>()(const MatrixT& matrix)
void Diagonal<MatrixT>()(const MatrixT& matrix, MatrixT& diagonalMatrix)

Get the main diagonal of a matrix, or put a vector into the main diagonal of a matrix. The first form returns the diagonal matrix. The second saves the diagonal matrix into a variable which is passed in. 
MatrixT ovariance<MatrixT>()(const MatrixT& matrix)
void Covariance<MatrixT>()(const MatrixT& matrix, MatrixT& covarianceMatrix)

Get the covariance of a matrix. The first form returns the covariance matrix. The second saves the covariance matrix into a variable which is passed in. 
MatrixT Power<MatrixT>()(const MatrixT& matrix, const value_type& power)
void Power<MatrixT>()(const MatrixT& matrix, MatrixT& powerMatrix, const value_type& power)

Get a matrix with all its elements raised to a defined power. The first form returns the powered matrix. The second saves the powered matrix into a variable which is passed in. 
MatrixT Mean<MatrixT>()(const MatrixT& matrix)
void Mean<MatrixT>()(const MatrixT& matrix, MatrixT& meanMatrix)

Get a row vector with the mean of each column of a matrix. The first form returns the row vector. The second saves the row vector into a variable which is passed in. 
MatrixT Median<MatrixT>()(const MatrixT& matrix)
void edian<MatrixT>()(const MatrixT& matrix, MatrixT& medianMatrix)

Get a row vector with the median of each column of a matrix. The first form returns the row vector. The second saves the row vector into a variable which is passed in. 
MatrixT Sum<MatrixT>()(const MatrixT& matrix)
void Sum<MatrixT>()(const MatrixT& matrix, MatrixT& sumMatrix)

Get a row vector with the sum over each column of a matrix. The first form returns the row vector. The second saves the row vector into a variable which is passed in. 
MatrixT CumulativeSum<MatrixT>()(const MatrixT& matrix)
void CumulativeSum<MatrixT>()(const MatrixT& matrix, MatrixT& cumulativeSumMatrix)

Get the cumulative sum of elements of a matrix. The first form returns the cumulative sum matrix. The second saves the cumulative sum matrix into a variable which is passed in. 
Matrix Class Example
##include <span class="codekeyword"><string>
</span>
Implementation Details
Introduction
This matrix class was designed to be generic and with reusability in mind. The inspiration for this design comes from Andrei Alexandrescu^{4}. Since the matrix was to be generic, templates are essential for the implementation. The matrix must be able to store different types of data, not just numerical data. Thus, a typename ValueT
must be present to indicate the type of data that is stored.
There are different types of matrices: dense matrix, sparse matrix, triangular matrix, diagonal matrix, etc. Each differs in how the matrix is organized, and how the data is stored. Thus, there should be a policy which determines how the matrix is stored and how the data can be accessed. Hence, a separate class for each kind of matrix should be created. This class will handle the storage and access of the data. The second element in the matrix class template parameter is thus template<typename> class MatrixTypeP
. At the present moment, MatrixTypeP
can only be DenseMatrix
.
The third thing that must be considered when designing the matrix class is that for the matrix to be useful in scientific programming, the matrix should have mathematical operations. However, if the matrix contains nonnumerical types, mathematical operations should not be available. This problem can be solved by using a mathematical policy. By implementing a mathematical policy, matrices containing nonnumerical types cannot perform mathematical operations and will give compilererror if this is attempted, whereas mathematical operations will be available for numerical types. Initially, two classes are created to achieve this policy: MathMatrix
and NonMathMatrix
, and the third element in the matrix class template parameter is the template<typename,typename>
class MathP
. In the initial version of the matrix class, MathP
was derived from MatrixTypeP
, and then the matrix class was derived from either of the math classes. This was later changed to the present version where the matrix class was derived from both MatrixTypeP
and MathP
, and MathP
no longer derived from MatrixTypeP
. In addition, a third math class, MathMETMatrix
, was created which implements expression templates to increase the performance when performing certain types of matrix calculations. The use of expression templates will be discussed later.
Thus, the skeleton of the matrix class is formed. The matrix class will take in three template arguments which denote the type of data which will be stored, the policy for storing and accessing the data, and the policy for determining whether mathematical operations are available for the matrix. The only thing which was not considered was the dimension of the matrix. The matrix was assumed to be 2dimensional. Multidimensional matrixes were not considered as this would increase the difficulty of the implementation.
The operations of the matrix class was designed to be similar to a STL container. The matrix class implements the various STLlike typedefs like value_type
, size_type
, etc. Iterators for iterating through the elements in the matrix were also implemented. Certain common STL container functions like clear
, empty
, resize
, and swap
were implemented. Others like assign
, erase
, insert
, front
, back
, push_back
, push_front
, pop_back
, pop_front
, begin
, and end
were implemented in the row and column iterators.
Implementation of MatrixTypeP
There are a few basic public functions that any MatrixTypeP
matrix must have. There must be a default constructor which takes in no argument. There should be constructors which take in the initial size of the matrix as well as the initial values. The initial values can be a default value for all the elements, or from a matrix implemented using arrays or pointers. A copy constructor is essential because of the presence of a raw pointer as a member variable. There should be an assignment operator for copying another similar MatrixTypeP
matrix, and an assignment operator for copying any matrix class. The MatrixTypeP
matrix must also implement the resize, swap, and comparison operators. In addition, an Update
function was added to allow the user to inform the matrix to update itself whenever data is changed. This Update
function is not used by DenseMatrix
but it may be necessary for other MatrixTypeP
.
There are other basic public functions which can theoretically be placed in the final matrix class instead of the MatrixTypeP
matrix. They are clear
, empty
, the various member accessor functions, operator()
, and at
, etc. They can theoretically be placed in the final matrix class because their code doesn’t need to be changed for all MatrixTypeP
matrices. Their implementation is based on common member variables and functions that will be present in all MatrixTypeP
matrices. However, they are still placed in the MatrixTypeP
matrix at the present moment, because they provide functions which may be useful for the MatrixTypeP
matrix to use in the implementation of other functions.
There are two more private functions which should be present in a MatrixTypeP
matrix. They are ResizeIterators
and UpdateIterators
. ResizeIterators
is used to resize the four private member variables: rowIteratorStart
, rowIteratorFinish
, colIteratorStart
, and colIteratorFinish
. UpdateIterators
is used to update the data stored by these four member variables. Initially, there is only UpdateIterators
. The resizing of the four private member variables is done inside UpdateIterators
. Later on, the resizing is separate from UpdateIterators
to form ResizeIterators
in order to increase the performance of the matrix. The four private member variables are raw pointers which keep information on the iterators which point to the start and end of each row and column. By storing these iterators, the member accessor functions and the various iterator functions can be speeded up and made generic for all MatrixTypeP
. Initially, the four private member variables were implemented using the STL vector. Then it was changed to use raw pointers, to increase performance.
Finally, each MatrixTypeP
matrix must implement a row iterator class and a column iterator class. The row iterator is responsible for iterating through a row in the matrix, and the column iterator for iterating through a column. These two iterator classes are the most important aspects of the implementation of a MatrixTypeP
matrix. I will try to give a detailed explanation on how to implement the iterator classes.
Since the iterator classes are to be similar to an STL container, it should have the various STLlike typedefs. Thus, iterator_category
, value_type
, reference, pointer, size_type
, and difference_type
are defined.
The iterator should have appropriate constructors. It must have a default constructor which takes in no arguments. It should have a constructor which takes in an iterator (VectorTypeIterator
) which iterates through a row or column of the type which stores the matrix data. For example, in DenseMatrix
, the matrix data is stored by a raw pointer of type value_type*
. Thus, the VectorTypeIterator
should be of type value_type*
. If the matrix data is stored by an STL vector, then the VectorTypeIterator
is of type std::vector::iterator
. If the matrix data is stored by an STL deque<deque>
, then the VectorTypeIterator
is of type std::deque::iterator
. This VectorTypeIterator
is stored in a member variable for use in various functions. In addition to the VectorTypeIterator
, the constructor may take in any variable which may help in the operation of the iterator. For example, the constructor in the column iterator of the DenseMatrix
takes in the total number of columns in the matrix. It uses this information in its Subtract
, Increment
, Decrement
, and Advance
functions.
Since the default constructor may be used to construct the iterator, there should be an assignment operator which takes in VectorTypeIterator
. This is to allow the update of the member variable which stores the VectorTypeIterator
. Other functions to allow updates of any variables which is necessary for the iterator operation should also be available. For example, there is a SetCols
function to update the total number of column variables in the column iterator of DenseMatrix
.
The type of functions which should be present in the iterator will depend on the type of the iterator, i.e., input, output, bidirectional, or random. All types of iterators should implement the operator*
to allow dereferencing, and operator>
to allow accessing of matrix elements’ member variables or functions. Friend comparison functions must also be implemented for all types of iterators.
For a random access iterator, it should have pre and post increment and decrement operators. It should also have operator
to determine the distance between two similar iterators. In addition, it should implement operator+=
, operator=
, operator+
, and operator
to allow the iterator to increase or decrease by a certain distance. These functions can be made generic by delegating the jobs to other functions like Subtract
, Increment
, Decrement
, and Advance
. A random access iterator should also implement the operator[]
to allow it to act like an array.
template<typename ValueT> class DenseMatrix
The implementation of dense matrix has gone through several stages. Initially, deque<deque<ValueT> >
was used to store the matrix elements. Then, deque<ValueT>
was used. In the end, ValueT*
was used to increase the performance of the matrix as the use of STL containers slow matrix calculations greatly. However, the usage of raw pointers in the storage increases the difficulty of implementing the DenseMatrix
class, slightly. Care must be taken during the implementation to avoid memory leaks. Also, memory leaks may occur when exceptions occur. Since performance is crucial in scientific programming, this disadvantage of using raw pointers is overlooked in this case.
The implementation of iterators is slightly different from that of the template described above. The column iterator uses the template above for implementation, but the row iterator uses a raw pointer ValueT*
. In the initial version, a class following the template format was used to create the row iterator. Later, it was found that in matrix calculations, repeated use of functions in the row iterator class actually slows down the calculations. Since the storage of data allows the use of raw pointers for the row iterators, raw pointers are used to increase the speed of matrix calculations. The same cannot be done for the column iterator because iterating through a column requires knowledge of the total number of columns in the matrix, and this information needs to be stored and used. Thus, a proper class is needed for the column iterator. Since using the row iterator is faster than using the column iterator, the row iterator should be used when using the DenseMatrix
storage policy unless it is absolutely necessary to use the column iterator so as to increase the efficiency of the matrix.
Approximate analysis of the speed of matrix calculations using the DenseMatrix
storage policy indicates that it is nowhere as fast as handwritten loops involving raw pointers or arrays. This is to be expected since the underlying matrix elements are stored under many layers of code, which add a significant overhead. In order to make this matrix class appealing to scientific programmers who require the fast speed of raw pointers, three member variables and a function, which normally should be private or protected, is made public. This is against the principle of encapsulation, but the end justifies the means. The three variables are rows_
and cols_
, which contain the total number of rows and columns in the matrix, respectively, and matrix_
which is a raw pointer ValueT*
containing the matrix elements. The function which is made public is the ResizeIterators
. Scientific programmers who desire speed, yet convenience, can make use of the various functions of the matrix class to help maintain the matrix, and when doing calculations, can directly access the matrix storage variable, matrix_
, to perform fast lookup, assignment, or calculations, by using it as an ordinary raw pointer. The ResizeIterators
is exposed so that if the programmer ever resizes the matrix_
or points the matrix_
to another matrix manually, the internal iterators can be updated so that the structure of the matrix is maintained.
Implementation of MathP
The functions which are to be implemented in a MathP
class must support addition, subtraction, dot multiplication, and scalar multiplication of the matrix when the data is of numerical type. When the data is of nonnumerical type, all these functions should not be available. There are two MathP
classes which provide mathematical operations, and one which does not. The two that does, implement the operator+=
, operator+
, operator=
, operator
, operator*=
, and operator*
.
The implementation of a MathP
class presents a unique problem. The MathP
class is a base class from which the final matrix class inherits from. It is separated from the storage class, thus it has no idea of how the matrix elements are stored. Since the storage class is also in charge of providing accessor functions, the MathP
class effectively has no way of accessing the matrix elements. But to perform matrix calculations, it has to have access to the matrix elements. The way to overcome this problem is to pass in a pointer representing the final matrix class to the MathP
class so that it can have access to the matrix elements. The MathP
class must maintain a pointer variable that points to its child matrix. Thus, the final declaration of a MathP
class is template<typename ValueT, typename MatrixT>
, where MatrixT
is the type of the final matrix class. Since the MathP
class cannot exist without a pointer to its child matrix, there is no default constructor for a MathP
class which takes in no arguments. There is only a constructor which requires a pointer to the child matrix.
template<typename ValueT, typename MatrixT> class MathMatrix
The MathMatrix
class is a relatively simple class. The main function of this class is to provide mathematical operators to the final matrix class so that matrix calculations can be done intuitively using mathematical operators.
The mathematical operators are implemented in such a way that the operator+
and operator
calls the functions operator+=
and operator=
respectively to perform their functions. This technique is also used to implement the operator*=
and operator*
for scalar multiplication. However, for dot multiplication, the situation is reversed. The operator*=
calls the functions operator*
to perform its function. The reason for this reversal is because it is more efficient. During matrix dot multiplication, a temporary matrix is created to store the results of the multiplication. If this code resides in operator*=
instead of in operator*
, when using operator*
, three temporary matrices will be created, instead of the usual two, in operator+
and operator
. This reduces the efficiency of the matrix dot multiplication. By putting the code in operator*
, only two temporary matrixes are created, so it is more efficient. Approximate analysis reveals that this does not degrade the performance of the operator*=
, so there is no need to duplicate the multiplication code in operator*=
.
The MathMatrix
class implements the checking of the matrix sizes before performing the various matrix calculations. It will return an exception if the two matrixes cannot be added, subtracted, or dot multiplied together because of wrong sizes. This checking does not slow down the matrix calculations appreciably, but can be removed if really necessary for performance.
template<typename ValueT, typename MatrixT> class MathMETMatrix
MathMETMatrix
uses a technique called expression templates to increase the performance of certain types of matrix addition and subtraction. Expression templates was introduced by Todd Veldhuizen^{5}. To explain expression templates, consider the problem below.
A, B, C, D, E, F are matrices. Consider the code A = B + C + D + E + F;
. To perform matrix calculations, the compiler will first add B
and C
together, then add the results to D
, then the results to E
, and the results to F
, before putting the final result into A
. Thus, several temporay matrices are created during the addition, and this decreases the performance greatly. What expression templates hope to achieve is to eliminate all these unnecessary temporary matrices by creating code equivalent to the following:
for (size_t i=0; i<rows; ++i)
{
for (size j=0; j<cols; ++j)
{
A(i,j) = B(i,j) + C(i,j) + D(i,j) +E(i,j) + F(i,j);
}
}
To achieve the above, templates are used, and several supporting classes are needed. I will not go into the theory of expression templates. Rather, I will explain the creation of the various supporting template classes.
The first template class that is needed is a wrapping class for the Operations type. It is named MET. In the MathMETMatrix
class, there is only one type of Operations, that is METElementBinaryOp
, which perform binary operations on the elements in the matrix. The MET main job is to wrap around METElementBinaryOp
, providing access to the functions present in METElementBinaryOp
. To do that, it maintains a variable for the template class METElementBinaryOp
. MET is necessary because it helps to bind different METElementBinaryOp
template classes together. Without MET, the whole concept of expression templates cannot work.
The next template class is the Operations types. Presently, there is only METElementBinaryOp
. METElementBinaryOp
takes in two template variables in its constructor, and saves them in pointers that point to the data in those variables. The two template variables can either be a matrix class or a MET template class. METElementBinaryOp
template argument also takes in an operator type class. This operator type class determines the kind of operation to perform on the two template variables.
There are three operator types present in MathMETMatrix
: METAdd
, METSubtract
, and METMultiply
. METMultiply
is not used yet because a good algorithm implementing expression templates for performing dot multiplication has not been found. METAdd
and METSubtract
are structures which provide a static function Evaluate
. Evaluate
takes in two matrix elements and returns the results of their addition or subtraction.
To implement expression templates, the operator+
and operator
are overloaded so that it will return an instance of a MET object. There is a operator+
and operator
which perform operations between a MathMETMatrix
derived matrix and any other matrix. There are also friend operator+
and operator
which perform operations between any matrix and any MET object, and between two MET objects. Each operator overload will return a MET object which contains information on the Operations types to do (only METElementBinaryOp
at the present moment), the operator types to perform (METAdd
for operator+
and METSubtract
for operator
), and the variables to perform the operator types on (a matrix or an expression).
Thus, when the compiler meets the code above: A = B + C + D + E + F;
, it will generate the following code:
A = MET< METElementBinaryOp< METAdd<value_type>, Matrix,
Matrix> >(METElementBinaryOp< METAdd<value_type>,
MET< METElementBinaryOp< METAdd<value_type>,
Matrix, Matrix> >(METElementBinaryOp< METAdd<value_type>,
MET< METElementBinaryOp< METAdd<value_type>, Matrix,
Matrix> >(METElementBinaryOp< METAdd<value_type>, MET<
METElementBinaryOp< METAdd<value_type>, Matrix, Matrix> >,
Matrix>, Matrix>(& MET< METElementBinaryOp<
METAdd<value_type>, Matrix, Matrix> >(METElementBinaryOp<
METAdd<value_type>, MET< METElementBinaryOp<
METAdd<value_type>, Matrix, Matrix> >, Matrix>,
Matrix>(& MET< METElementBinaryOp<
METAdd<value_type>, Matrix, Matrix>
>(METElementBinaryOp< METAdd<value_type>, MET<
METElementBinaryOp< METAdd<value_type>, Matrix, Matrix>
>(METElementBinaryOp< METAdd<value_type>, MET<
METElementBinaryOp< METAdd<value_type>,
Matrix, Matrix> >, Matrix>,
Matrix>(& MET< METElementBinaryOp<
METAdd<value_type>, Matrix, Matrix>
>(METElementBinaryOp< METAdd<value_type>,
MET< METElementBinaryOp< METAdd<value_type>,
Matrix, Matrix> >, Matrix>(& MET<
METElementBinaryOp< METAdd<value_type>,
Matrix, Matrix> >(METElementBinaryOp<
METAdd<value_type>, Matrix, Matrix>(&B,
&C)), &D)), &E)), &F)); &C)), &D)), &E)), &F))
Don’t worry about the above code. What it basically does is that it generates a MET object which encapsulates all the information required for doing B + C + D + E + F
. All that is required now is to implement the calculation of the results and put it in matrix A
.
To implement the calculation, the assignment operator (operator=
) was overloaded in the final matrix class. The assignment operator will take in a MET object as its argument. Since there are many different types of MET objects, the operator=
has to be a template function. The resulting code is as follows:
template<typename ExprT> Self& operator=(MET<ExprT> expr)
{
size_t sizeRow = expr.rowsize;
size_t sizeCol = expr.colsize;
Self tempMatrix(sizeRow,sizeCol);
for (size_t i=0; i<sizeRow; ++i)
{
for (size_t j=0; j<sizeCol; ++j)
{
tempMatrix(i,j) = expr(i,j);
}
}
swap(tempMatrix);
Update;
return *this;
}
The code is relatively simple to understand. What it does is that it gets the required total number of rows and columns from the MET object and creates a temporary matrix to hold the results of the matrix calculations. Then it loops through each individual element, and finds the results for the calculation. Using the example above, the results of the calculation is the result of B(i,j) + C(i,j) + D(i,j) +E(i,j) + F(i,j)
. The temporary matrix is then swapped with the current matrix (which is matrix A
, in the example) so that the results of the calculation is now assigned to the current matrix.
To digress slightly, this function is slower than that of the Matrix Expression Template^{3} because of the creation of the temporary matrix. However, this temporary matrix is necessary to make the operation more intuitive, i.e., it doesn’t matter what the size of matrix A
is before the calculation, because it will be resized to the correct size required by the matrix calculation of B + C + D + E + F
. This is not the case for the Matrix Expression Template class. In that class, matrix A
is required to be the same size as B, C, D, E, and F, otherwise the operation will fail. The present matrix class can achieve better efficiency than the Matrix Expression Template if the same stance was adopted. The temporary matrix just needs to be removed, and substitute tempMatrix(i,j)
with (*this)(i,j)
. However, this was not done because it is more essential for the operation to be intuitive in this case than for it to be fast.
Now, let’s go back to our discussion on expression templates. From the operator=
code, it can be seen that there is a need to implement some functions in the various expression template supporting classes. Two functions, row.size
and col.size
need to be implemented to find out the total number of rows and columns needed for the matrix receiving the results. There is also a need to implement an operator which takes in the matrix element position and outputs the results of the mathematical operation.
For the MET class, it is simple to implement these three functions. These three functions are just passed on to the METElementBinaryOp
for processing.
METElementBinaryOp
implements row.size
and col.size
by passing it on to one of its two template variables. The variable left_
is chosen. There is no reason why the variable right_
cannot be chosen. If left_
is a matrix, it will have the functions row.size
and col.size
, and thus will return the number of rows and columns respectively. If left_
is a MET object, it will pass it on to METElementBinaryOp
, which will pass in on to left_
recursively and so forth, until it reaches a matrix.
METElementBinaryOp
implements an operator by applying it on the left_
and right_
variables. If left_
and right_
are matrices, they will return the element at position (i,j). If they are a MET object, they will pass it on to METElementBinaryOp
until it finally reaches a matrix and the element at position (i,j) is returned. The returned elements are then sent to the Evaluate
function of the operator type which will then add or subtract them and return the results. Thus, the innocently looking expr(i,j)
actually performs all the necessary calculations on the individual elements (e.g., B(i,j) + C(i,j) + D(i,j) + E(i,j) + F(i,j)
).
That concludes the discussion on the implementation of expression templates for the matrix class. One problem which can be seen is that there is no proper way to ensure that all the matrices are of the same size, before the matrix calculation is done. Thus, it is up to the user to ensure that the matrices are of the same size, before performing the calculations. Otherwise, unexpected errors may arise. Another problem is that expression templates may not speed up all matrix addition or subtraction. For short additions like A = B + C;
, using expression templates may actually slow down the calculation. Thus, the user should experiment with both MathMatrix
and MathMETMatrix
to determine the best mathematical matrix to use for each situation. Both contains code to interconvert one form to another, and both types of matrices can be used in the same calculation, so there is no penalty incurred in having both types available for use in the final matrix class. There is currently no good algorithm for implementing dot multiplication using expression templates. Thus, at the present moment, dot multiplication for MathMETMatrix
is using the same code as MathMatrix
.
template<typename ValueT, typename MatrixT> NonMathMatrix
This is the simplest of the three MathP
classes. It contains the bare essential functions so that it presents the same interface as MathMatrix
and MathMETMatrix
. Other than that, it does not contain any mathematical operators so that a matrix deriving from it cannot perform any mathematical operations and will give a compiler error if tried.
Implementation of the Matrix
template<typename ValueT, template<typename> class MatrixTypeP, template<typename, typename> class MathP> class Matrix
We have reached the most important class, the Matrix
class. It is the class that the user interacts with. The final functions that are present in this class is dependent on which storage policy and which mathematical policy that it inherits from.
Even though Matrix
inherits most of its functions from MatrixTypeP
and MathP
, constructors, the destructor, copy constructors, and assignment operators are not inherited. Thus, it needs to implement these. Since it is to be a STLlike container, it has to define various STLlike typedefs. These typedefs are obtained from MatrixTypeP
by redefining them. The MathP
policy does not have any constructor or assignment operator which must be mirrored by Matrix
, thus Matrix
only has to define constructors, copy constructors, and assignment operators which are present in MatrixTypeP
and pass the operations to it. Hence, there must be a default constructor which takes in no arguments. There should be constructors which take in the initial size of the matrix as well as the initial values. The initial values can be a default value for all elements, or from a matrix implemented using arrays or pointers. A copy constructor should also be present. There should be an assignment operator for copying another similar matrix, and an assignment operator for copying any matrix class.
In addition, Matrix
implements one additional copy constructor which is not found in MatrixTypeP
. This copy constructor takes in any matrix class. The code is implemented by passing the matrix to be copied to the template operator=
in MatrixTypeP
instead of a similar copy constructor in MatrixTypeP
. The reason for this is that if there is a copy constructor in MatrixTypeP
which can take in all matrix classes, it will prevent the normal copy constructor in Matrix
from working properly because all the matrixes passed into MatrixTypeP
will be processed by the template version of the copy constructor instead of the normal version. The reason for this is simple. The construct of the normal copy constructor in Matrix
is Matrix(const Self& m) : MathBase(this), StorageBase(m)
. As can be seen, StorageBase
is passed a const reference to m
. If there is no template version of the copy constructor in StorageBase
, m
will be processed by the normal copy constructor in StorageBase
. However, if a template version is present in StorageBase
, the StorageBase
will use the template version to process m
since m
is of type Matrix
, so it does not fit into the normal copy constructor which requires the type StorageBase
. Hence, a template version of the copy constructor should not be present in StorageBase
, otherwise all copies may result in slow operations or error. However, a template version of the copy constructor will be useful since there could be different Matrix
types each differing only in their MatrixTypeP
or MathP
, and it should be possible to copy one type to another. Thus, a template version of the copy constructor is created only in the Matrix
class, and it passes the operation to a template operator=
in the MatrixTypeP
which normally should perform a copy of the matrix using the common functions available to all matrices.
Matrix
also implements an additional assignment operator (operator=
) which takes in a MET object. The reason is discussed above in MathMETMatrix
.
Two other functions which are useful for all matrices but are not essential for MatrixTypeP
operations are implemented in Matrix
. They are swaprows
and swapcols
. These two functions use only iterators for their operations so that they can be applicable to all types of MatrixTypeP
matrices.
Implementation of RowVector and ColVector
template<typename ValueT,template<typename> class MatrixTypeP,template<typename,typename> class MathP> class RowVector
template<typename ValueT,template<typename> class MatrixTypeP,template<typename,typename> class MathP> class ColVector
Both RowVector
and ColVector
are derived from the Matrix
class. Their existence is almost unnecessary as the Matrix
class can be a row vector or a column vector, by setting the number of rows or columns to 1, respectively. Also, STL already implements three useful and powerful 1dimensional containers. So, why is there a need for RowVector
and ColVector
? The initial idea for creating them is for performance and ease of use.
RowVector
and ColVector
can be faster than STL containers when they use the DenseMatrix
class as a storage policy. As mentioned earlier, DenseMatrix
exposes its storage mechanism to performancehungry users. Thus, users can directly access the elements in the matrix using the exposed variable instead of passing through accessor functions. STL containers do not expose their storage mechanism, thus the only way to access their elements is through accessor functions. This adds overhead and slows down performance.
For variables destined to be a vector throughout their lifetime, RowVector
and ColVector
are easier to use than the Matrix
class. For these variables, either their number of rows or columns is always 1. Thus, the flexibility offered by the Matrix
class in accessing the elements or obtaining iterators is unnecessary. To enable ease of use, RowVector
and ColVector
implement several functions present in STL containers:operator[]
, at
, push_back
, resize
, size
, begin
, end
, rbegin
, and rend
. RowVector
and ColVector
also implement constructors requiring only a single size variable instead of the two in the Matrix
class (one for the number of rows and one for the number of columns).
Although RowVector
and ColVector
have simplified functions for accessing and modifying its elements and iterators, they retain all the functions available in the Matrix
class in order to be compatible with the Matrix
class. They are also interconvertible with the Matrix
class. However, the user must ensure that when copying or assigning a Matrix
class to a RowVector
or ColVector
class, the matrix class must be either a row vector or column vector (i.e. either its number of rows or columns must be 1) respectively. Otherwise, the RowVector
or ColVector
may end up being a matrix instead. This will not cause any errors when using the class, and the user may not even notice it, but it is logically wrong as a row vector or column vector should not store a matrix.
Functions Operating On Matrix
template<typename MatrixT> struct Transpose
Transpose
is a functor which finds the transpose of a matrix (a transposed matrix is a matrix where the rows and columns are swapped). It has two overloaded operators, the first taking in a matrix to be transposed and returns the transposed matrix, the second taking in both a matrix to be transposed and a matrix to hold the transposed matrix. The first operator is suitable to use in matrix calculations where it can be used intuitively, or in any appropriate STL algorithms. The second operator is faster, and should be used when the transpose of a matrix needs to be found but need not be used immediately in calculations.
template<typename MatrixT> struct Diagonal
Diagonal
is a functor which finds the main diagonal of a matrix or puts a vector into the main diagonal of a matrix. It has two overloaded operators, the first taking in a matrix and returns the diagonal matrix, the second taking in both a matrix and a matrix to hold the diagonal matrix. The first operator is suitable to use in matrix calculations where it can be used intuitively or in any appropriate STL algorithms. The second operator is faster, and should be used when the diagonal of a matrix needs to be found but need not be used immediately in calculations.
template<typename MatrixT> struct Covariance
Covariance
is a functor which finds the covariance of a matrix. It has two overloaded operators, the first taking in a matrix and returns the covariance matrix, the second taking in both a matrix and a matrix to hold the covariance matrix. The first operator is suitable to use in matrix calculations where it can be used intuitively, or in any appropriate STL algorithms. The second operator is faster, and should be used when the covariance of a matrix needs to be found but need not be used immediately in calculations.
template<typename MatrixT> struct Power
Power
is a functor which calculates the matrix with all its elements raised to a defined power. It has two overloaded operators, the first taking in a matrix and returns the powered matrix, the second taking in both a matrix and a matrix to hold the powered matrix. The first operator is suitable to use in matrix calculations where it can be used intuitively, or in any appropriate STL algorithms. The second operator is faster, and should be used when the powered matrix needs to be found but need not be used immediately in calculations.
template<typename MatrixT> struct Mean
Mean
is a functor which finds a row vector with the mean of each column of a matrix. It has two overloaded operators, the first taking in a matrix and returns the row vector, the second taking in both a matrix and a matrix to hold the row vector. The first operator is suitable to use in matrix calculations where it can be used intuitively, or in any appropriate STL algorithms. The second operator is faster, and should be used when the row vector needs to be found but need not be used immediately in calculations.
template<typename MatrixT> struct Median
Median
is a functor which finds a row vector with the median of each column of a matrix. It has two overloaded operators, the first taking in a matrix and returns the row vector, the second taking in both a matrix and a matrix to hold the row vector. The first operator is suitable to use in matrix calculations where it can be used intuitively, or in any appropriate STL algorithms. The second operator is faster, and should be used when the row vector needs to be found but need not be used immediately in calculations.
template<typename MatrixT> struct Sum
Sum
is a functor which finds a row vector with the sum over each column of a matrix. It has two overloaded operators, the first taking in a matrix and returns the row vector, the second taking in both a matrix and a matrix to hold the row vector. The first operator is suitable to use in matrix calculations where it can be used intuitively, or in any appropriate STL algorithms. The second operator is faster, and should be used when the row vector needs to be found but need not be used immediately in calculations.
template<typename MatrixT> struct CumulativeSum
CumulativeSum
is a functor which finds a matrix with the cumulative sum of the elements of a matrix. It has two overloaded operators, the first taking in a matrix and returns the cumulative sum matrix, the second taking in both a matrix and a matrix to hold the cumulative sum matrix. The first operator is suitable to use in matrix calculations where it can be used intuitively, or in any appropriate STL algorithms. The second operator is faster, and should be used when the cumulative sum matrix needs to be found but need not be used immediately in calculations.
Conclusions
This concludes the discussion on the Matrix
class. The Matrix
class makes use of concepts like policies and expression templates in its implementation. It is designed to be generic, reusable, and easily extendable. Though some of the concepts may not be implemented correctly or in the best possible ways, it is an attempt to consolidate various C++ techniques together. This article serves as an example of how to use the Matrix
class, a documentation for all the various classes and functions, and also a documentation of all the various phases and thoughts in the design and implementation of the classes. It is hoped that this article will provide sufficient information for all who wish to use, modify, or extend this class. It is also hoped that the brief explanations on policies, iterators, and expression templates will spur those who are interested to read more about these techniques and use them in their own work.
Any comments about this class are greatly appreciated, and any requests for new functions will be considered and added as soon as possible, if appropriate.
References
 William H. Press, Numerical recipes in C++: the art of scientific computing (Cambridge University Press, 2002).
 Matrix Template Library.
 Matrix Expression Template.
 Andrei Alexandrescu, Modern C++ Design (AddisonWesley, 2001).
 Todd Veldhuizen, Expression Templates: C++ Report, Vol 7, No. 5, June 1995, pg 2631.
Change logs
 Version 1.1
 Removed the identity matrix implementation because it was not well implemented.
 Version 1.2
 Added
RowVector
and ColVector
classes.
 Fixed errors involving reverse iterators.
 Changed function
empty()
.
 Version 1.3
 Added member spaces
row
and col
to StoragePolicy and Expression Templates.
 Replaced
rowsize
, colsize
, rowbegin
, rowend
, colbegin
, and colend
with row::size
, col::size
, row(i).begin
, row(i).end
, col(i).begin
, and col(i).end
.
 Added
row::insert
, row::erase
, row::push_back
, row::pop_back
, col::insert
, col::erase
, col::push_back
, and col::pop_back
.
 Added
pop_back
for RowVector
and ColVector
.
 Version 1.4
 Changed member spaces
row
and col
to Row
and Col
to prevent conflicts causing row(i)
not to work when RowVector
or ColVector
is declared.
 Version 1.5 (16 April 2006)
 Fixed code so that it is no longer dependent on STLport and will work on Visual Studio 2005, and GCC 3.3.2.
 Because of the fix, the constancy of the row and col iterators may not really work.
 Added
row::push_front
, row::pop_front
, col::push_front
, and col::pop_front
.
 Added
push_front
and pop_front
for RowVector
and ColVector
.