DAE with measurement type

Specific Domains Numerics
1

Hi there I am running a complex set of differential equations and at some point it fails because of gemm!.
The short story is that gemm! does not work with the meassurememnt type, and because the solver is calling gemm! at some point the equation can not be solved.

When I raise the issue within the Measurement repository they explained as much -
“I think that BLAS will work only with the standard types. Try looking into equivalent functions in GitHub - JuliaLinearAlgebra/GenericLinearAlgebra.jl: Generic numerical linear algebra in Julia, they should work out-of-the-box with Measurement objects (however don’t expect super performance :frowning:)”

Does anyone have some advice how to proceed.
I can solve small systems of equations using meassurment type, but when I extended to my rela problem it crashed because of gemm!

Thanks,

A

PS1 - Error form solving DAE system

error:

MethodError: no method matching gemm!(::Char, ::Char, ::Measurement{Float64}, ::SubArray{Measurement{Float64},2,Array{Measurement{Float64},2},Tuple{UnitRange{Int64},UnitRange{Int64}},false}, ::SubArray{Measurement{Float64},2,Array{Measurement{Float64},2},Tuple{UnitRange{Int64},UnitRange{Int64}},false}, ::Measurement{Float64}, ::SubArray{Measurement{Float64},2,Array{Measurement{Float64},2},Tuple{UnitRange{Int64},UnitRange{Int64}},false})
Closest candidates are:
gemm!(::AbstractChar, ::AbstractChar, !Matched::Float64, !Matched::Union{AbstractArray{Float64,1}, AbstractArray{Float64,2}}, !Matched::Union{AbstractArray{Float64,1}, AbstractArray{Float64,2}}, !Matched::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
…

Stacktrace:
[1] reckernel!(::Array{Measurement{Float64},2}, ::Val{true}, ::Int64, ::Int64, ::Array{Int64,1}, ::Base.RefValue{Int64}, ::Int64) at /opt/julia/packages/RecursiveFactorization/ow69g/src/lu.jl:98
[2] #lu!#3(::Bool, ::Int64, ::Function, ::Array{Measurement{Float64},2}, ::Array{Int64,1}, ::Val{true}) at /opt/julia/packages/RecursiveFactorization/ow69g/src/lu.jl:18
[3] #lu! at ./none:0 [inlined]
[4] #lu!#2 at /opt/julia/packages/RecursiveFactorization/ow69g/src/lu.jl:8 [inlined]
[5] lu! at /opt/julia/packages/RecursiveFactorization/ow69g/src/lu.jl:8 [inlined] (repeats 2 times)
[6] #call#436(::Nothing, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::DiffEqBase.DefaultLinSolve, ::Array{Measurement{Float64},1}, ::Array{Measurement{Float64},2}, ::Array{Measurement{Float64},1}, ::Bool) at /opt/julia/packages/DiffEqBase/6ewdP/src/linear_nonlinear.jl:40

PS2- An error and simple example code illustarting issue with 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 returns true

function addError(x;e=0.2)
    if typeof(x)!= Measurement{Float64}
        return x±e*x
    else
        return x
    end
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

Does not work:
Error
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
…

1 Like
2

what is the size of your problem? maybe you can get away with direct multiplication,defining the most generic code for gemm!:

import 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)
C*=beta
C+=alpha*A*B
end
3

That’s a great ideas. The problem has only 42 equations do this should work.
I guess i need to check where to modify the call in the differential equations library to substitute the gemm! call for this function :slight_smile: or add it to the measurements library. What would you recommend?

Thanks!

A

1 Like
4

Thanks this worked-
I just realized what the first line allowed me to do -
Thanks so much!

A

1 Like
5

It’s a hacky patch hahah, good to know, and thanks to Julia and his ability to patch even the core functions :sweat_smile:

6

It’s super simple, but here is the function I used in case it helps anyone else.
Should add some checks and balances for dimensions


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