Support for gemm and gemm!

[ayadlin]

Hi-
I am solving a complex system of DAE which calls on gemm. My parameters and initial conditions are all of type Measurement{Float64}. When solving toy systems, everyhitng works fine. issues are risen. Here is a toy example I run using Measurements and gemm:

A = rand(5,5)
B = rand(5,5)
# C = Array{Float64,2}
C=zeros(5,5)
D = BLAS.gemm('N', 'T', 1.0, A, B)
BLAS.gemm!('N','T',1.0, A, B,0.0,C)
D==C

Works and return true

when using measurement type

function addError(x;e=0.2)
    if typeof(x)!= Measurement{Float64}
        return x±e*x
    else
        return x
    end
A1=addError.(A;e=0.2)
B1=addError.(B;e=0.2)
C1=addError.(C;e=0.2)
D1 = BLAS.gemm('N', 'T', 1.0, A1, B1)
BLAS.gemm!('N','T',1.0, A1, B1,0.0,C1)
D1==C1

Either gemm or gemm! break because of the measurement type (errors below)
Is this something that can be addressed?
Making the library compatible with LinearAlgebra.BLAS would be very useful for complex problems that reuqire error propagations.

Thanks for your time.

GEMM:
MethodError: no method matching gemm(::Char, ::Char, ::Float64, ::Array{Measurement{Float64},2}, ::Array{Measurement{Float64},2})
Closest candidates are:
gemm(::AbstractChar, ::AbstractChar, ::Float64, !Matched::AbstractArray{Float64,2}, !Matched::AbstractArray{Float64,2}) at /opt/julia-1.0.0/share/julia/stdlib/v1.0/LinearAlgebra/src/blas.jl:1132
gemm(::AbstractChar, ::AbstractChar, !Matched::AbstractArray{Float64,2}, !Matched::AbstractArray{Float64,2}) at /opt/julia-1.0.0/share/julia/stdlib/v1.0/LinearAlgebra/src/blas.jl:1135
gemm(::AbstractChar, ::AbstractChar, !Matched::Float32, !Matched::AbstractArray{Float32,2}, !Matched::AbstractArray{Float32,2}) at /opt/julia-1.0.0/share/julia/stdlib/v1.0/LinearAlgebra/src/blas.jl:1132

GEMM!
MethodError: no method matching gemm!(::Char, ::Char, ::Float64, ::Array{Measurement{Float64},2}, ::Array{Measurement{Float64},2}, ::Float64, ::Array{Measurement{Float64},2})
Closest candidates are:
gemm!(::AbstractChar, ::AbstractChar, ::Float64, !Matched::Union{AbstractArray{Float64,1}, AbstractArray{Float64,2}}, !Matched::Union{AbstractArray{Float64,1}, AbstractArray{Float64,2}}, ::Float64, !Matched::Union{AbstractArray{Float64,1}, AbstractArray{Float64,2}}) at /opt/julia-1.0.0/share/julia/stdlib/v1.0/LinearAlgebra/src/blas.jl:1109
gemm!(::AbstractChar, ::AbstractChar, !Matched::Float32, !Matched::Union{AbstractArray{Float32,1}, AbstractArray{Float32,2}}, !Matched::Union{AbstractArray{Float32,1}, AbstractArray{Float32,2}}, !Matched::Float32, !Matched::Union{AbstractArray{Float32,1}, AbstractArray{Float32,2}}) at /opt/julia-1.0.0/share/julia/stdlib/v1.0/LinearAlgebra/src/blas.jl:1109
gemm!(::AbstractChar, ::AbstractChar, !Matched::Complex{Float64}, !Matched::Union{AbstractArray{Complex{Float64},1}, AbstractArray{Complex{Float64},2}}, !Matched::Union{AbstractArray{Complex{Float64},1}, AbstractArray{Complex{Float64},2}}, !Matched::Complex{Float64}, !Matched::Union{AbstractArray{Complex{Float64},1}, AbstractArray{Complex{Float64},2}}) at /opt/julia-1.0.0/share/julia/stdlib/v1.0/LinearAlgebra/src/blas.jl:1109

[giordano]

I think that BLAS will work only with the standard types. Try looking into equivalent functions in https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl, they should work out-of-the-box with Measurement objects (however don't expect super performance ☹️)

[ayadlin]

Thanks!
I'll look at it-

A

[giordano]

I think this can be closed, right? Reference: https://discourse.julialang.org/t/dae-with-measurement-type/25281

[ayadlin]

Yes, but i found a patchy solution if you want me to post it before you close it. I’m traveling, but I can post on Monday. A

[giordano]

I think you refer to (source)

import LinearAlgebra.BLAS: gemm!
function gemm!(tA, tB, alpha, A, B, beta, C)
#implement transpose logic here
#tX = 'N', no modifications
#tX = 'T', transpose(X)
#tX = 'C', adjoint(X)
    if tA=='T'
        A=transpose(A)
    elseif tB=='C'
        A=adjoint(A)
    else
        A=A
    end
    if tB=='T'
        B=transpose(B)
    elseif tB=='C'
        B=adjoint(B)
    else
        B=B
    end

C*=beta
C+=alpha*A*B
end

or is it another one?

[ayadlin]

Yes very patchy, but it worked :)

[giordano]

Thanks!