(untagged)

A generic, reusable, and extendable matrix class

Yap Chun Wei

0.00/5 (No votes)

19 Apr 2006

An article describing the use of a generic matrix class and how it is implemented.

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:

double* matrix = new double[100]

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++¹, the Matrix Template Library (MTL) developed by Open Systems Laboratory at Indiana University², and the Matrix Expression Template implemented by Masakatsu Ito³.

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:

STL-like 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 inter-convertible `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 inter-convertible `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 <string>

#include "Matrix.hpp"


using namespace std;
using namespace YMatrix;

int main(int argc, char* argv[])
{
     // typedef for non-mathematical matrix with string as datatype.

     typedef Matrix<string,DenseMatrix,NonMathMatrix> MatrixNM;
     // typedef for normal mathematical matrix with double as datatype.

     typedef Matrix<double,DenseMatrix,MathMatrix> MatrixM;
     // typedef for mathematical matrix implementing expression template.

     typedef Matrix<double,DenseMatrix,MathMETMatrix> MatrixEM;
     // typedef for normal mathematical row vector with double as datatype.

     typedef RowVector<double,DenseMatrix,MathMatrix> RVectorM;
     // typedef for normal mathematical col vector with double as datatype.

     typedef ColVector<double,DenseMatrix,MathMatrix> CVectorM;

     int i, j;

     // construct a 3 x 4 matrix and fill it with 'Hello' for all elements.

     MatrixNM stringMatrix(3, 4, "Hello");
     // resize matrix to 2,5 and fill new elements with 'World'

     stringMatrix.resize(2, 5, "World");
     // print matrix.

     cout << "First way to print matrix\n";
     for (i=0; i<stringMatrix.row.size(); ++i)
     {
          for (j=0; j<stringMatrix.col.size(); ++j)
          {
               cout << stringMatrix(i,j) << "\t";
          }
          cout << "\n";
     }
     cout << endl;

     vector<string> s(5, "New");
     stringMatrix.row.insert(1, s.begin());
     stringMatrix.col.insert(2, s.begin());
     // another way to print matrix

     cout << "Second way to print matrix\n";
     for (i=0; i<stringMatrix.row.size(); ++i)
     {
          copy (stringMatrix.row(i).begin(), 
                stringMatrix.row(i).end(), 
                ostream_iterator<string>(cout, "\t"));
          cout << "\n";
     }
     cout << endl;

     stringMatrix.col.erase(0);
     // another way to print matrix

     cout << "Third way to print matrix\n";
     cout << stringMatrix;

     // erase matrix.

     cout << "Deleting contents of matrix\n";
     stringMatrix.clear();
     // check if matrix is empty

     cout << (stringMatrix.empty() ? 
          "Matrix is empty\n" : "Matrix is not empty\n");

     // prepare two 3 x 3 arrays

     double a1[9] = {1,2,3, 4,5,6, 7,8,9};
     double** a2 = new double*[3];
     for (i=0; i<3; ++i) a2[i] = new double[3];
     for (i=0; i<3; ++i)
     {
          for (j=0; j<3; ++j)
          {
               a2[i][j] = 9 - (i*3+j);
          }
     }

     // construct three 3 x 3 matrix with no initial value

     MatrixM mathMatrix1(3,3), mathMatrix2(3,3), mathMatrix3;
     // fill matrices with array

     mathMatrix1 = a1;
     mathMatrix2 = a2;
     cout << "First matrix:\n" << mathMatrix1;
     cout << "Second matrix:\n" << mathMatrix2;
     cout << "Element(1,2) of first matrix: " << mathMatrix1.at(1,2) << endl;

     // perform mathematical operations on matrix and print the results

     cout << "Some mathematical operations\n";
     mathMatrix3 = mathMatrix1 + mathMatrix2; cout << mathMatrix3;
     mathMatrix3 = mathMatrix1 - mathMatrix2; cout << mathMatrix3;
     mathMatrix3 = mathMatrix1 * mathMatrix2; cout << mathMatrix3;
     mathMatrix3 = mathMatrix1 * 2; cout << mathMatrix3;
     mathMatrix3 = mathMatrix1 + mathMatrix2 * mathMatrix2;
     cout << mathMatrix3;
     mathMatrix3 += mathMatrix1; cout << mathMatrix3;
     mathMatrix3 -= mathMatrix1; cout << mathMatrix3;
     mathMatrix3 *= mathMatrix1; cout << mathMatrix3;
     mathMatrix3 *= mathMatrix1 + mathMatrix2; cout << mathMatrix3;
     mathMatrix3 *= 2; cout << mathMatrix3;

     // swap mathMatrix1 and mathMatrix2

     mathMatrix1.swap(mathMatrix2);
     cout << "Swapping first and second matrix:\n";
     cout << "First matrix:\n" << mathMatrix1;
     cout << "Second matrix:\n" << mathMatrix2;

     // swap row 0 and 1 of mathMatrix1

     cout << "Swapping first and second row of first matrix:\n";
     mathMatrix1.swaprows(0,1); cout << mathMatrix1;

     // swap column 0 and 2 of mathMatrix1

     cout << "Swapping first and third column of first matrix:\n";
     mathMatrix1.swapcols(0,2); cout << mathMatrix1;

     // print transpose of mathMatrix1

     cout << "Transpose\n" << Transpose<MatrixM>()(mathMatrix1);

     // print diagonal of mathMatrix1

     cout << "Diagonal\n" << Diagonal<MatrixM>()(mathMatrix1);

     // print covariance of mathMatrix1

     cout << "Covariance\n" << Covariance<MatrixM>()(mathMatrix1);

     // print power 2 of mathMatrix1

     cout << "Power 2\n" << Power<MatrixM>()(mathMatrix1, 2);

     // print mean of mathMatrix1

     cout << "Mean\n" << Mean<MatrixM>()(mathMatrix1);

     // print median of mathMatrix1

     cout << "Median\n" << Median<MatrixM>()(mathMatrix1);

     // print sum of mathMatrix1

     cout << "Sum\n" << Sum<MatrixM>()(mathMatrix1);

     // print cumulative sum of mathMatrix1

     cout << "CumulativeSum\n" << CumulativeSum<MatrixM>()(mathMatrix1);

     // construct three 3 x 3 matrix using the two arrays and mathMatrix1

     MatrixEM mathEmatrix1(3, 3, a1);
     MatrixEM mathEmatrix2(3, 3, a2);
     MatrixEM mathEmatrix3(mathMatrix1);
     cout << "First matrix:\n" << mathEmatrix1;
     cout << "Second matrix:\n" << mathEmatrix2;
     // perform mathematical operations on matrix and print the results

     cout << "Some mathematical operations\n";
     mathEmatrix3 = mathEmatrix1 + mathEmatrix2; cout << mathEmatrix3;
     mathEmatrix3 = mathEmatrix1 - mathEmatrix2; cout << mathEmatrix3;
     mathEmatrix3 = mathEmatrix1 * mathEmatrix2; cout << mathEmatrix3;
     mathEmatrix3 = mathEmatrix1 * 2; cout << mathEmatrix3;
     mathEmatrix3 = mathEmatrix1 + mathEmatrix2 * mathEmatrix2;
     cout << mathEmatrix3;
     mathEmatrix3 += mathEmatrix1; cout << mathEmatrix3;
     mathEmatrix3 -= mathEmatrix1; cout << mathEmatrix3;
     mathEmatrix3 *= mathEmatrix1; cout << mathEmatrix3;
     mathEmatrix3 *= 2; cout << mathEmatrix3;
     // invalid mathematical operation (compare with above)

     // mathEmatrix3 *= mathEmatrix1 + mathEmatrix2; cout << mathMatrix3;

     // demonstrate speed of expression template

     // (can't tell unless timing is done)

     mathEmatrix3 = mathEmatrix1 + mathEmatrix2 - 
                    mathEmatrix1 + mathEmatrix2;
     cout << mathMatrix3;
     // demonstrate interchangably of matrices

     mathEmatrix3 = mathMatrix1 + mathEmatrix2; cout << mathEmatrix3;
     // get transpose of mathEmatrix3

     MatrixEM mathEmatrix3Transposed;
     Transpose<MatrixEM>()(mathEmatrix3, mathEmatrix3Transposed);
     cout << "Transposed:\n" << mathEmatrix3Transposed;

     RVectorM rv(5, 3.0);
     cout << rv;
     rv.push_back(8); cout << "push_back 8\n" << rv;
     rv.push_front(1); cout << "push_front 1\n" << rv;
     rv.pop_back(); cout << "pop_back\n" << rv;
     rv.pop_front(); cout << "pop_front\n" << rv;
     cout << endl;

     CVectorM cv(5, 2.6);
     cout << cv;
     cv.push_back(3.5); cout << "push_back 3.5\n" << cv;
     cv.push_front(8.5); cout << "push_front 8.5\n" << cv;
     cv.pop_back(); cout << "pop_back\n" << cv;
     cv.pop_front(); cout << "pop_front\n" << cv;
     cout << endl;
}

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⁴. 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 Veldhuizen⁵. 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³ 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 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.

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 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).
Matrix Template Library.
Matrix Expression Template.
Andrei Alexandrescu, Modern C++ Design (Addison-Wesley, 2001).
Todd Veldhuizen, Expression Templates: C++ Report, Vol 7, No. 5, June 1995, pg 26-31.

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.

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.

A list of licenses authors might use can be found here