Opened 3 years ago
Closed 3 weeks ago
#881 closed discussion (fixed)
expression class for dgemv
Reported by: | Peter | Owned by: | Peter |
---|---|---|---|
Priority: | major | Milestone: | yat 0.18 |
Component: | utility | Version: | trunk |
Keywords: | Cc: |
Description (last modified by )
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
Change History (7)
comment:1 follow-up: 2 Changed 3 years ago by
comment:2 Changed 3 years ago by
comment:3 Changed 3 years ago by
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.
comment:4 Changed 3 years ago by
Description: | modified (diff) |
---|
comment:5 Changed 4 weeks ago by
Milestone: | yat 0.x+ → yat 0.18 |
---|
comment:6 Changed 4 weeks ago by
Owner: | changed from Jari Häkkinen to Peter |
---|---|
Status: | new → accepted |
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.