Dear julia users,
When looking at e.g. BLAS functions, they have a general format of adding a
scaled result array to an existing array, which itself can also be scaled;
e.g. AB could be the result of a matrix multiplication of matrices A and B,
and the BLAS gemm! function (using Julia's name) allows to store
alpha*AB+beta*C into C. Looking at the high-performant microkernels of
BLAS, they specialize on the different possibilities for alpha and beta
using if constructions:
if beta == 0
# do not actually compute C*0 since that could produce NaN if C was not
initialized and contains NaN or infinity
# instead set C to zero
elseif beta != 1
# only scale C if beta is not 1
end
if alpha == 1
# immediately add AB, do not scale with alpha
else
# C+=alpha * AB
end
In writing something related in Julia, I thought this would be a perfect
case for dispatch. The high performant low level function could just
contain the general statement C=beta*C+alpha*AB and it would be called in a
high level function as (note that I am not actually trying to reimplement
gemm!, this is just as an example)
gemm!((beta == 1 ? _one : (beta == 0 ?_zero : beta)) , C, (alpha == 1 ?_one
: alpha), A,B)
in combination with the following definitions
immutable Zero <: Integer
end
immutable One <: Integer
end
const _zero = Zero()
const _one = One()
Base.promote_rule{T<:Number}(::Type{Zero}, ::Type{T}) = T
Base.promote_rule{T<:Number}(::Type{One}, ::Type{T}) = T
Base.convert{T<:Number}(::Type{T}, ::Zero) = zero(T)
Base.convert{T<:Number}(::Type{T}, ::One) = one(T)
# add special rules, the most essential of which are:
+(::Zero, a::Number) = a
+(::Number, a::Zero) = a
*(::Zero, a::Number) = _zero
*(a::Number, ::Zero) = _zero
*(::One, a::Number) = a
*(a::Number, ::One) = a
Unfortunately this produces a lot of ambiguity warnings. Of course I need
to add the rules for +(::Zero,::Zero) etc, i.e. the interaction of my newly
defined types with themselves. But this is not sufficient. The remaining
addition and multiplication ambiguities can be avoided by defining them
separately for a::Real and a::Complex. But the most difficult one is the
convert definition, which is ambiguous with almost any other number type
out there.
I am aware of the many ambiguity related issues and the proposal to make
those runtime errors instead of warnings, but I have the feeling that there
might be a more specific solution for what I am trying to accomplish and
that I am missing something trivial. I know there are several other number
types defined in packages, and I assume the corresponding authors must have
faced similar problems? So I appreciate any tips or input.
