# DAE with measurement type

**URL:** <https://discourse.julialang.org/t/dae-with-measurement-type/25281>\
**Category:** Numerics\
**Tags:** diffeq\
**Created:** [June 13, 2019, 10:47pm UTC](https://discourse.julialang.org/t/dae-with-measurement-type/25281 "2019-06-13T22:47:27Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![alewolf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alewolf/32/9000_2.png) [@alewolf](https://discourse.julialang.org/u/alewolf)\
**Post date:** [June 13, 2019, 10:47pm UTC](https://discourse.julialang.org/t/dae-with-measurement-type/25281/1 "2019-06-13T22:47:28Z")

</div>

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](https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl), they should work out-of-the-box with Measurement objects (however don’t expect super performance ☹)”

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!

```julia
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

```julia
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  
…

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [June 14, 2019, 6:57am UTC](https://discourse.julialang.org/t/dae-with-measurement-type/25281/2 "2019-06-14T06:57:04Z")

</div>

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

```julia
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

```

---

<div class="post-metadata">

**Author:** ![alewolf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alewolf/32/9000_2.png) [@alewolf](https://discourse.julialang.org/u/alewolf)\
**Post date:** [June 14, 2019, 2:09pm UTC](https://discourse.julialang.org/t/dae-with-measurement-type/25281/3 "2019-06-14T14:09:51Z")

</div>

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 🙂 or add it to the measurements library. What would you recommend?

Thanks!

A

---

<div class="post-metadata">

**Author:** ![alewolf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alewolf/32/9000_2.png) [@alewolf](https://discourse.julialang.org/u/alewolf)\
**Post date:** [June 14, 2019, 4:06pm UTC](https://discourse.julialang.org/t/dae-with-measurement-type/25281/4 "2019-06-14T16:06:46Z")

</div>

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

A

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [June 14, 2019, 5:14pm UTC](https://discourse.julialang.org/t/dae-with-measurement-type/25281/5 "2019-06-14T17:14:55Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![alewolf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alewolf/32/9000_2.png) [@alewolf](https://discourse.julialang.org/u/alewolf)\
**Post date:** [June 14, 2019, 5:40pm UTC](https://discourse.julialang.org/t/dae-with-measurement-type/25281/6 "2019-06-14T17:40:29Z")

</div>

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

```julia

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

```