Some operations do not work although they (probably) should

[andreasvarga]

The following code illustrates this issue:

M = ones(Rational,3,3)
v = ones(Int,3)
M*v  # OK, the result is a rational vector
LinearOperator(M)*v # failure

The execution of this code produces:

ERROR: TypeError: in typeassert, expected Array{Real,1}, got a value of type Array{Rational,1}
Stacktrace:
 [1] *(::LinearOperator{Rational}, ::Array{Int64,1}) at C:\Users\Andreas\.julia\packages\LinearOperators\z6AM5\src\operations.jl:7
 [2] top-level scope at REPL[126]:1
 [3] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1088

I wonder what is the right approach in determining the type of the elements of the result of an operation M*v. For example, in my concrete case, M is an inverse operator to solve a matrix equation, with the underlying matrices of a certain element type T (say T = Float32). The vector v can have any element type T1 (say T1 = Float64) such that it can be still converted to the element type T, without error. In this case, the result M*v would have element type T, which would be, for me, the desirable type. Still note that for M a matrix, the result of M*v has element type T (and not T1), which involves a conversion of M (generally not desirable in my case).

[abelsiqueira]

The issue here is that Rational does not provide enough information. It should be Rational{Int} instead:

M = ones(Rational{Int},3,3)
v = ones(Int,3)
LinearOperator(M)*v # works

The issue is not due to LinearOperators but due to the promotion rules:

julia> promote_type(Rational, Int)
Real

julia> promote_type(Rational{Int}, Int)
Rational{Int64}

[abelsiqueira]

Regarding the more general issue

...underlying matrices of a certain element type T (say T = Float32). The vector v can have any element type T1 (say T1 = Float64) such that it can be still converted to the element type T, without error. In this case, the result M*v would have element type T...

The default Julia promotion rules should apply here:

julia> M = ones(Float32, 3, 3)
3×3 Array{Float32,2}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

julia> v = ones(Float64, 3)
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

julia> M * v
3-element Array{Float64,1}:
 3.0
 3.0
 3.0

Is there any additional operation involved?

[andreasvarga]

Many thanks for your explanation. Abel Siqueira <notifications@github.com> schrieb am Mo., 5. Okt. 2020, 17:00:

…

The issue here is that Rational does not provide enough information. It should be Rational{Int} instead: M = ones(Rational{Int},3,3) v = ones(Int,3)LinearOperator(M)*v # works The issue is not due to LinearOperators but due to the promotion rules: julia> promote_type(Rational, Int) Real julia> promote_type(Rational{Int}, Int) Rational{Int64} — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#153 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ALJDHECQ7TFIOVCHU764Y3LSJHNRJANCNFSM4SEZDWCA> .

[andreasvarga]

In MatrixEquations.jl the operation exits with error for a Lyapunov operator of type Float32 and vector of type Float64. Abel Siqueira <notifications@github.com> schrieb am Mo., 5. Okt. 2020, 17:07:

…

Regarding the more general issue ...underlying matrices of a certain element type T (say T = Float32). The vector v can have any element type T1 (say T1 = Float64) such that it can be still converted to the element type T, without error. In this case, the result M*v would have element type T... The default Julia promotion rules should apply here: julia> M = ones(Float32, 3, 3) 3×3 Array{Float32,2}: 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 julia> v = ones(Float64, 3) 3-element Array{Float64,1}: 1.0 1.0 1.0 julia> M * v 3-element Array{Float64,1}: 3.0 3.0 3.0 Is there any additional operation involved? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#153 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ALJDHEBZGED4FNBSBHY7UUDSJHOL3ANCNFSM4SEZDWCA> .