Scaled Apply + GEMV
The proposal of scaled operators in !703 initiated this issue.
Having the standard BLAS level 2 routine gemv
(see e.g. here)
y \gets \alpha A x + \beta y
as part of the Operator
interface would be great for obvious reasons.
What we currently can do is
y \gets A x + \gamma y
where \gamma \in \{0, 1\}
(a.k.a., Replace
and Add
).
Implementing the full gemv
would be the optimal fix, however, since many of the matrix-free operator implementations work in an "additive" fashion (particularly all elementwise operators including the new generated versions) I do not see an obvious way around having to touch y
(i.e., the "dst" vector) twice to apply scaling with \beta
. If there is no obvious solution, we are forced to first scale y
(if \beta \neq 0
) and then perform a scaled matrix-vector multiplication after that. We can do this scaling of y
inside the operator, but I do not see a real reason to do that inside the operator and not just call assign
beforehand.
MR !703 proposes an implementation that enables scaled operators
y \gets \alpha A x + \gamma y
that is in my opinion not optimal since \alpha
becomes part of the operator's state.
To resolve the "state"-issue mentioned above, I propose to actually extend the apply
routine by a parameter \alpha
instead of changing the operator's state:
virtual void apply( const ValueType& alpha // new parameter
const SourceFunction& src,
const DestinationFunction& dst,
size_t level,
DoFType flag,
UpdateType updateType = Replace ) const
Refactoring this requires some work, but should hopefully be doable. I'd keep the current apply, and just make it call the new one with \alpha = 1
.
EDIT: We basically already do this:
We can do this scaling of
y
inside the operator, but I do not see a real reason to do that inside the operator and not just callassign
beforehand.
if UpdateType == Replace
since we interpolate 0
.
Therefore, for completeness I'll also propose an interface for the gemv
:
virtual void gemv( const ValueType& alpha // new parameter
const SourceFunction& src,
const ValueType& beta // new parameter
const DestinationFunction& dst,
size_t level,
DoFType flag ) const
with the perk that it gets rid of the UpdateType
parameter.