expression class for dgemv

[Peter]

Currently the expression classes are minimalistic that each class represents one operator and two argument classes. For instance, Matrix*Matrix is represented by a MatrixBinary<Matrix, Matrix, Multiplies> and then for each operation in the expression another branch is added in the template class tree.

Typically (ignoring plus and minus for a second) this means that for every operation a GSL vector/matrix is allocated and computed via appropriate blas function. This is a bit wasteful as it allocates memory for structs that are just needed in intermediate steps and calls. If we for example have:

alpha * A * x + beta * y

it creates ScaledMatrix representing alpha * A, then MatrixVector created representing ScaledMatrix*x, then a ScaledVector representing beta*y, and then lastly the resulting gsl_vector is allocated and computed as MatrixVector+ScaledVector element by element.

The idea is to instead recognize that we can compute the result with one BLAS call to dgemv, and to accomplish that we need a class that represents this expression.

The way to accomplish consist of two things. We need specializations of the operators. As an example, double*ScaledMatrix<Matrix> can return ScaledMatrix<Matrix> rather than ScaledMatrix<ScaledMatrix<Matrix> > (using the fact that alpha*A*beta equals (alpha*beta)*A. The second thing is that we need to bridge from class from taking two objects to class that takes five objects via classes that take 3 or 4 objects. To start with we need a ScaledMatrixVector to represent alpha*A*x and then similarly operator+ or operator- should can return DGEMV class when argument classes are appropriate (remember cases when scalars are absent, i.e., implicitly being 1.0.

related to ticket #892

[Peter]

For this to work, constructors of these expression classes must be lazy and not compute anything. MatrixBinary<T,T,Multiplies> is not lazy at the moment.

[Peter]

Replying to peter:

For this to work, constructors of these expression classes must be lazy and not compute anything. MatrixBinary<T,T,Multiplies> is not lazy at the moment.

This is not true since r3653.

[Peter]

Note that y appears on both LHS and RHS in

y = alpha * A * x + beta * y`

which limits the complexity of this class. A finite beta seems mostly useful for

VectorMutable::operator+=(expression&&)
VectorMutable::operator-=(expression&&)

which is off-pist for this ticket.

[Peter]

In 3908:

#closes #881

Lift out expression classes to their own files, so they can be
included separately (and avoiding cyclic inclusion).

Replace the two traits classes for matriox expression with one
class. Use this class in MatrixVector? expression, which makes
transposition a no-op and we can therefore implement Vector * Matrix
as transpose(Matrix) * Vector (since there is no distinction between
row and column vectors in yat).