3.7 Expression templates for matrix types
3.7.2 Participating classes and operators
A = B*C*D + E*F - G;
the expression tree shown in Figure 3.7 will be built by the operators-,+and *. It will be evaluated when the = operator (a member of the MatrixType class) is called. The compute()method of the MatrixExpressionobject on the right side of the =sign is called and the expression tree is traversed from its root to its leaves, descending first to the right and then to the left (this is a pure and arbitrary implementation decision, of course it could have been done also first to the left and then to the right).
If the algorithm descend to a leaf, nothing is done, if it descends to an operatorop, the result of the expression ’left childop right child’ is computed.
-+ G
*
* E F
B *
C D
Figure 3.7: Expression tree
3.7.2 Participating classes and operators The MatrixExpression class
First and foremost, we have the MatrixExpression class, that is intended to represent the matrix expressions.
Matrix expression class template <class Left, class Op, class Right>
class MatrixExpression {
public:
typedef typename Left::VerticalDimensionType VerticalDimensionType;
typedef typename Right::HorizontalDimensionType HorizontalDimensionType;
typedef RowWise Orientation;
typedef typename Left::index_type index_type;
typedef typename Right::value_type value_type;
typedef typename Right::WeightMemoryType WeightMemoryType;
private:
typedef typename vector<map<index_type,value_type> > TmpWeightValuesContainer;
typedef typename TmpWeightValuesContainer::iterator TmpWeightValuesIterator;
public:
typedef typename
two_dim_map_indexed_iterator_adapter<TmpWeightValuesIterator> RowIterator;
private:
Left& left;
Right& right;
bool computed;
bool providememory;
TmpWeightValuesContainer* tmpvalues;
public:
MatrixExpression(Left& l, Right& r)
: left(l), right(r), computed(false), providememory(false), tmpvalues(0) {}
};
Within this class, the public types were defined in order to ensure that MatrixExpression can behave like a normal matrix class. Note that MatrixExpression is yet only designed to hold results which are row wise oriented.
The class stores references to the two involved operands of type Left and Right. The member TmpWeightValuesContaineris a pointer to the temporary data structure mentioned in Section 3.7.1.
If the stored values in the matrix type are not pointers, then the normal procedure would be the following: the = assignment operator of the matrix class calls the compute() function of the MatrixExpressionclass. This in turn computes the expression from the leaves of the expression tree recursively to its root. Then the values are copied from the root expression object by using the member functions TmpWeightValues begin() and TmpWeightValues end(). Thereafter, the deleteTmp() member function is called and frees the allocated memory for the temporary expression(s).
The boolean variable providememory is only of interest when computing expressions with block matrices, or more precisely, if the WeightMemoryType is of type Pointer. If the member function useMemory()is called,providememoryis set to true. This has the consequence that (in theroot expres-sion object) the allocated memory to which the pointers in the vectors of TmpWeightValuesContainer point to, won’t be destroyed by the deleteTmp()function. The reason for this mechanism is that in this case, there is no need to physically copy the already allocated submatrices from the root expression object into the result matrix object. Instead, only the pointers are copied.
In addition to the publicly defined types, the member functionsRow begin(),Row end(), vertical-Dimension(), andhorizontalDimension(), were defined in order to provide a matrix interface to the involved operator given by the class Op. Table 3.1 gives a summary of the publicly provided member functions.
return type function name
void compute()
void deleteTmp()
void useMemory()
TmpWeightValuesIterator TmpWeightValues begin() const TmpWeightValuesIterator TmpWeightValues end() const Left::VerticalDimensionType verticalDimension() const Right::HorizontalDimensionType horizontalDimension() const
RowIterator Row begin() const
RowIterator Row end() const
Table 3.1: Member functions of classMatrixExpression
The matrix product functor structure
The matrix productstructure acts as a functor that is given to theMatrixExpression class by the
’*’ operator.
Matrix product functor template <class Matrix1Type, class Matrix2Type,
class Orientation, class WeightMemoryType>
struct matrix_product {
template<class RowIterator>
static void computeResult(Matrix1Type& A, Matrix2Type& B, RowIterator it) {
matrix_product(A.Row_begin(), A.Row_end(),
B.Row_begin(), it, WeightMemoryType());
} };
The matrix product functor, as well as the matrix sumand the matrix differencefunctor solely has the static member functioncomputeResult. It calls the template functionmatrix product, that computes the sparse structure of the product and the new matrix entries on the fly by storing them in the according vector of STL maps (the auxiliary *tmpvaluesinMatrixExpression).
In this default version, the functor implements a sparse matrix-matrix multiplication in which the two operands are row wise oriented. If at least one operand is oriented column wise, we need a specialization. Any combination of column wise and row wise oriented matrices would require an own specialization with an own algorithm, since the direction of traversal through the matrix is different each time. Thus, not all combinations are yet implemented. Table 3.2 gives an overview of the implemented algorithms.
2. operand row wise 2. operand column wise 1. operand row wise implemented not implemented
1. operand column wise implemented not implemented
Table 3.2: Implemented combinations of sparse matrix-matrix multiplication
The specialization for the variant with the first operand being column wise oriented is declared as follows:
Matrix product functor
template <class Matrix1Type, class Matrix2Type, class WeightMemoryType>
struct matrix_product<Matrix1Type,Matrix2Type, ColumnWise,WeightMemoryType>
{
//...
};
Thetransposed matrix producttemplate function is used instead to calculate the result within this implemenation.
The matrix sum functor structure
Thematrix sum structure is a functor that is given to theMatrixExpressionclass by the+operator.
Matrix sum functor template <class Matrix1Type, class Matrix2Type,
class Orientation, class WeightMemoryType>
struct matrix_sum {
//...
}
The generateResultStructuremember function uses the generateMatrixSumPositionstemplate function to compute the sparsity pattern of the matrix sum A+B. The template function matrix sum then computes the actual values of the sum and is used by thecomputeResult member function.
Again, we have restricted the implemented combinations of the row wise/column wise algorithms, see Table 3.3.
2. operand row wise 2. operand column wise 1. operand row wise implemented not implemented
1. operand column wise not implemented partly implemented Table 3.3: Implemented combinations of sparse matrix-matrix addition
Note that the column wise/column wise matrix addition algorithm is the same as the row wise/row wise version. However, since the MatrixExpression class up to now only supports row wise result matrices, there would have to be at least a conversion routine from column wise to row wise somewhere in the functor structure or in the MatrixExpressionclass.
Thematrix sumtemplate function makes use of a template metaprogramming technique. Its dec-laration looks like the following.
Matrix sum function template
template <class M1Iterator, class M2Iterator, class M3Iterator, class WeightMemoryType, class Positive2ndOperand>
inline void matrix_sum(M1Iterator matrix1_begin, M1Iterator matrix1_end, M2Iterator matrix2_begin, M3Iterator matrix3_begin,
WeightMemoryType wt,
Positive2ndOperand positivesecond);
Depending now on the type of Positive2ndOperand respectively thepositivesecond parameter it chooses different helper functions inside the function body. If Positive2ndOperand is of type True, then an addition is performed (by calling functions that do an addition of two values or assign a positive value). The function call looks like this:
Call from matrix sum functor matrix_sum(A.Row_begin(), A.Row_end(), B.Row_begin(), cit,
WeightMemoryType(), True());
If Positive2ndOperand is of type False, then the 2nd operand is formally negated and a sub-traction is performed (by calling functions that do a subsub-traction of two values or assign a negative value). Now since the type of positivesecondis known at compile time, the according functions can be inlined and optimized by the compiler.
The matrix difference functor structure
Thematrix differencestructure is a functor that is given to the MatrixExpressionclass by the -operator.
Matrix difference functor template <class Matrix1Type, class Matrix2Type,
class Orientation, class WeightMemoryType>
struct matrix_difference {
//...
}
Here we can reuse the generateMatrixSumPositionstemplate function for computing the sparsity pattern of the matrix difference since it is the same as for the sum.
The computeResultmember function uses thematrix sumtemplate function as described above.
It is called by
Call from matrix difference functor matrix_sum(A.Row_begin(), A.Row_end(), B.Row_begin(), cit,
WeightMemoryType(), False());
The table of implemented row wise/column wise combinations is the same as for the matrix sum functor.
The operators +, - and *
The operators +, -and *are implemented as template functions. Their purpose is to instantiate an object of the class MatrixExpression with the according template arguments for the left and right operand type and the operator type.
At this point we have to take care of the possible combinations of types that can instantiate a binary expression. Let op be an arbitrary operator type, MatrixType1 and MatrixType2 some
arbitrary matrix types, andMatrixExpression1and MatrixExpression2some arbitrary expression types then there are the possibilities
• MatrixType1 op MatrixType2
• MatrixExpression1 op MatrixType2
• MatrixType1 op MatrixExpression2
• MatrixExpression1 op MatrixExpression2
For the *operator and the third possibility e.g., the implementation looks like this:
The operator* for MatrixType * MatrixExpression template <class Left, class Op, class Right, class RightRight>
MatrixExpression<Left,
matrix_product<Left,
MatrixExpression<Right,Op,RightRight>, typename Left::Orientation,
typename Left::WeightMemoryType>, MatrixExpression<Right,Op,RightRight> >
operator*(Left& A, MatrixExpression<Right,Op,RightRight> B) {
typedef typename Left::Orientation Orientation;
typedef typename Left::WeightMemoryType WeightMemoryType;
return MatrixExpression<Left,
matrix_product<Left,
MatrixExpression<Right,Op,RightRight>, Orientation,
WeightMemoryType>,
MatrixExpression<Right,Op,RightRight> >(A,B);
}
For the implementation of the other operators we refer to the library code for further details.
The assignment layer
To the sparse matrix expression template framework of course also belongs an assignment opera-tor = which has to be a member function of the according matrix type. It is implemented in the AssignmentLayerwhich is described in Section 3.3.7.