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 built-in 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 two-dimensional 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 compile-time.

To dynamically adjust the size of the matrix, we can create a matrix using pointers. We can use:

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 University2, and the Matrix Expression Template implemented by Masakatsu Ito3.

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 time-consuming. 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 non-numerical types. When storing non-numerical 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.

STL-like 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 user-defined 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:

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.

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.

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="code-keyword"><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 Alexandrescu4. 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 non-numerical types, mathematical operations should not be available. This problem can be solved by using a mathematical policy. By implementing a mathematical policy, matrices containing non-numerical types cannot perform mathematical operations and will give compiler-error 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 2-dimensional. Multi-dimensional 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 STL-like 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 STL-like 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 hand-written 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 non-numerical 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 Veldhuizen5. 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:

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:

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:

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 Template3 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

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 STL-like container, it has to define various STL-like 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.

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 1-dimensional 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 performance-hungry 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 inter-convertible 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).

License

This article has no explicit license attached to it but may contain usage terms in the article text or the download files themselves. If in doubt please contact the author via the discussion board below.

I wanted a simple vector/matrix class to use for running a Monte Carlo simulation so I followed Rajesh's comments to modify this to work with g++. I was unhappy with how slow my simulation was compared to previous simulations I had written so I swapped this out for Boost::ublas and the difference was an immediate 30x speed-up in my code's execution.

Perhaps this matrix class served a purpose back when it was written, but given that Boost is free and easy to use, this project seems to be thoroughly obsolete.

Hunting for some template-use examples to learn from, stumble upon this set of tests in YMatrix_Demo.I find it puzzling why, when intializing an array using this system ("prepare"), the vector a1[9] = {1, 2, 3, 4, 5, 6, 7, 8, 9}, although arbitrary value-wise, is so essentially not your standard "double".

wherever you find error coming from undefined cols_ or rows_: change it to this->cols_ and this->rows_ respectively.

similarly any "undefined" instances should be attatched with "this->"for example swap(tempMatrix) should be changed to this->swap(tempMatrix),Update() should be this->Update()

Another problem comes with friend declaration. Instead of using typedef'd name such as EnclosingClass, replace it with actual name. For example at line number 733: friend class EnclosingClass; should be changed to friend class MET<exprt>; This trick should be used for almost all friend class declarations which is using typedef'd names.

Having done that, it is possible to compile and run the test code given on this webpage. But correctness needs to be checked. The changes I have suggested here is ad-hoc. Therefore the author who knows the philosophy of this code should comment on whether these changes will have any side-effect or not. I am sure visual studio compiler will not agree with such changes.

This implementation compiles with sunstudio-12 compiler toolset, but fails with g++ (>4.2) and intel c++ compiler 11.0 toolset (icc). However the kind of errors are different with different compilers. That means it might not be a good choice for a matrix class for general purpose use, because one does not know where the compiler-dependent problem can occur. However I find this effort very useful, at least conceptwise. If author can make it compatible with most of the compilers (especially g++) then it will be useful for a wide range of audience.

Promising! Not too extensive, and uses modern C++. I'm going to take a closer look.

How about getting rid of these many warnings in VS2005 (signed/unsigned mismatch, using this in initializer lists, deprecated functions)? Are you planning an update to fix it, or should I rather #pragma them away myself?

i basically wasted 40 min of my time trying it. try something else if you are quick and result driven.never gonna work for ya, there are other staff on the web that just workscheers-dm

BECAUSEi did not read that before i tried the class, wonderful motivation

Since there are already excellent matrix classes available, why is there a need to create another one? ---- 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.

[dmi@moscow geometry_lib]$ gcc MatrixTest.cppIn file included from MatrixTest.cpp:1:DataMatrix.h:76: declaration of `class ValueT'DataMatrix.h:58: shadows template parm `class ValueT'DataMatrix.h:81: syntax error before `>' tokenDataMatrix.h:82: syntax error before `>' tokenDataMatrix.h:107: syntax error before `&' tokenDataMatrix.h:109: missing ';' before right braceDataMatrix.h:118: syntax error before `operator'DataMatrix.h:67: warning: all member functions in class ` DataMatrix::MatrixStorageTemplate<valuet>' are privateDataMatrix.h: In constructor `DataMatrix::MatrixStorageTemplate<valuet>::RowIterator<ValueT, IteratorTraits>::RowIterator(ValueT*)':DataMatrix.h:97: class ` DataMatrix::MatrixStorageTemplate<valuet>::RowIterator<ValueT, IteratorTraits>' does not have any field named `current_'DataMatrix.h: At global scope:DataMatrix.h:122: ISO C++ forbids defining types within return typeDataMatrix.h:122: syntax error before `&' tokenDataMatrix.h:123: ISO C++ forbids declaration of `operator++' with no typeDataMatrix.h:123: `int& DataMatrix::operator++()' must have an argument of class or enumerated typeDataMatrix.h:123: `int& DataMatrix::operator++()' must take either one or two argumentsDataMatrix.h: In function `int& DataMatrix::operator++()':DataMatrix.h:124: invalid use of `this' in non-member functionDataMatrix.h:125: invalid use of `this' in non-member functionDataMatrix.h: At global scope:

I have been cruising around for different implementations of a matrix class, both my own implementations and others...

This implementation is one of the best I have seen, when also the compactness and simplicity is taken into account (THANKS!!!)

Though, my problem was to make this work without the Intel Compiler and without the STLport. Unfortunately, I haven't succeeded in installing STLport under Visual Studio .NET (v7.1). After having fiddled with STLport - withouth success - I decided to find an alternative...

In order for me to make this work simply implemented (read: added) the _Nonconst_traits and _Const_traits structs (taken from the STLport) to the header (hope it's ok concerning rights and so). Furthermore, the template definitions has to be added typename. Also, the bind2nd is part of the functional class.

All in all it did not take much more than half an hour (maybe one) to make the fixes.

It may be that this is an obvious way or that other exist. But for me at least, this was the only way to make it work!

Many years ago, in times of the first version of C language (and IBM compiler), I solved problem of transmitting library routines from FORTRAN and PL/I to C (not C++). Those languages are treated as lower level languages, but some their features has been more extended than those of C.

E.g. in Fortran IV, when allocating array, its size was naturally fixed. But there was possible to pass matrix dimension as argument to procedure (SUBROUTINE, how it was named in FORTRAN) and use it as file dimension in the style:

SUBROUTINE MATINV(A,N) INTEGER N REAL A(N,N) ....

Such approach solved many problems - in many cases, the dimension of matrix was definedby problem solved, so fixed dimension of matrix in caller was no problem,but it was possible to write universal function for any dimension.

When translating into C language, I passed naturally pointer and called it "a",but defined macro A(i,j) (a[n*i+j]) as acessor - so fortran program was simplyrewritten using A instead of a. The same approach can be used for simple matrixoperation. Finaly, that time format a(i,j) has been for me more natural than A[i][j], and it is for me more natural from some points of view up to now (or, say, A[i,j]).

More comfortable was situation in PL/I, where one was allowed allocate dynamically two or three dimensional arrays with variable lower and upper bounds with no problem - naturally using the only available internal array "class", with no chance of defining your own functions or operators. It was possible to use operators as +, - there, including * (multiplication member by member, no matrix multiplication). To say nothing about full structural feature of the language - freedom in definition of local functions inside other functions.

For this reason, I personally do not consider even C++ as so "high level" language as it is generally accepted. Naturally, many problems can be solved using "objects" - specific programmer class suited for their specific usage, but such basic concepts of structure programming as Dijkstra's coroutins has been practically forgetten.

There are also some good matrix routines available through Direct 3D. They might be worth accessing those for some applications.see D3DXMATRIX *D3DXMatrixMultiply (D3DXMATRIX *pOut, CONST D3DXMATRIX *pM1, CONST D3DXMATRIX *pM2);etc in the DirectX help

if you're looking for a small matrix class without too much functionability but will handle itself, search for apmatrix.h and apmatrix.cpp, they were made by the collegeboard, they work, and they're simple..... the only thing, this relies on their vector class (apvector.h and apvector.cpp)

Article is not wrong, neither is the project. Yes, this class cannot be compiled using Visual C++ compiler. But in Visual Studio, you can add-in the Intel C++ compiler and it is the Intel C++ compiler that was used to compiled the project.

It is not the authors problem that MSVC6 sucks errr- performs substandard when it comes to templates.With MSVC7.1 the class seems to work, using the tips from a forum entry here.GCC3.3 is actually quite good with templates, and STLPort is the native STL on most Linuxes. SO this class should work there also.

<hr size=1 /hr>"We trained hard, but it seemed that every time we were beginning to form up into teams we would be reorganised. I was to learn later in life that we tend to meet any new situation by reorganising: and a wonderful method it can be for creating the illusion of progress, while producing confusion, inefficiency and demoralisation."

LAPACK[^]: (Linear Algebra Package) is THE package for linear algebra built upon the BLAS to solve linear algebra problems: LU, Cholesky, QR, SVD, Schur, generalized Schur, well all you have ever dreamed off.

So here I come with the suggestion: could it be possible to build a template like DenseMatrix that has an internal structure compatible with LAPACK and could it be possible to build a policy MathLAPACK that wraps up the LAPACK methods ?

Thanks for your suggestion. I have heard about BLAS and LAPACK. From what I know, they are mostly fortran77 routines. LAPACK has been ported to C or C++ recently if I am not wrong. Both should be very good packages. However, to use these packages, you may have to interface fortran programs with C++, which is not an easy thing. There may also be performance penalty. So although these packages are good, they are of limited use to the C++ community.

It may be possible to build a policy that is compatible with LAPACK. Since the class Matrix itself is a template which depends on what policies you give it, it is possible to design a Storage policy that is compatible with LAPACK and possible to design a Math policy that wraps up LAPACK methods. It will not be easy though.

I have checked out MTL before. In fact, I thought it was a very good class and wanted to use it. In the end, the one thing that is sorely lacking in MTL killed my interest in it. That is the use of operators like +, -, and *. I believed that to simplify code and reduce coding errors and subsequent maintenance, the input of mathematical formulas should be intuitive and natural. MTL doesn't implement the mathematical operators because they believe it is bad for performance. Yes, performance is slower but with template expressions, the performance can be improved. Also, not every applications are that performance hungry.

Lastly, this is not only a matrix for performing mathematical calculations. It can also be used to store non-numerical data in a matrix form, something which other matrix class can't do without having odd behaviours like being able to perform mathematical operations

You really do want to use BLAS for the low-level calculation if you're doing compute-intensive work. If you have a build of ATLAS for your machine this is very very fast, comparable in many cases to hand-written (i.e. assembler) code. There is no performance penalty for calling Fortran from C/C++ -- the only difference is the calling convention which is pretty easy to deal with (all function names have _ appended and every argument is a pointer, although I'm not 100% sure this is right for all compilers) and the layout of the data in the matrix (column order).

I have written a template-based C++ matrix class which uses BLAS and also has some calls to LAPACK that we need (e.g. for Cholesky decomposition). Last I heard this code was going to be released under the GPL (as part of a larger system) but that is still some time off.