Complex Linear Algebra

[mipals]

Hello Boris,

First of all thank you for making this package available. HSS-matrices are something that I for a long time have wanted to play around with. This package finally makes that possible :smile:

My main application are acoustical problems. Here we mostly deal with complex matrices. Do you believe it would be possible to use this package for complex linear algebra?

Even some simple linear algebra such as multiplication seem to initially cause some problems

julia> n = 200
julia> A = rand(ComplexF64,n,n)
julia> hssA = hss(A)
julia> hssA*ones(n)
ERROR: MethodError: no method matching _matmatdown!(::Matrix{Float64}, ::HssMatrix{ComplexF64}, ::Matrix{Float64}, ::BinaryNode{Matrix{ComplexF64}}, ::Nothing, ::Float64, ::Float64)
Closest candidates are:
  _matmatdown!(::AbstractMatrix{T}, ::HssMatrix{T}, ::AbstractMatrix{T}, ::BinaryNode{AT}, ::Union{Nothing, AbstractMatrix{T}}, ::Number, ::Number) where {T, AT<:(AbstractArray{T, N} where N)} at /home/mpasc/.julia/packages/HssMatrices/1ZvgD/src/matmul.jl:44
  _matmatdown!(::HssMatrix{T}, ::HssMatrix{T}, ::BinaryNode{Matrix{T}}, ::Matrix{T}) where T at /home/mpasc/.julia/packages/HssMatrices/1ZvgD/src/matmul.jl:96
Stacktrace:
 [1] mul!(C::Matrix{Float64}, hssA::HssMatrix{ComplexF64}, B::Matrix{Float64}, α::Float64, β::Float64)
   @ HssMatrices ~/.julia/packages/HssMatrices/1ZvgD/src/matmul.jl:26
 [2] *
   @ ~/.julia/packages/HssMatrices/1ZvgD/src/matmul.jl:13 [inlined]
 [3] *(hssA::HssMatrix{ComplexF64}, x::Vector{Float64})
   @ HssMatrices ~/.julia/packages/HssMatrices/1ZvgD/src/matmul.jl:15
 [4] top-level scope
   @ REPL[128]:1

From the error it is obvious that this can be solved by casting the vector elements to have the same time as the matrix elements

julia> x = hssA*ones(eltype(hssA),n)

Unfortunately things go wrong when I try to solve the system

julia> hssA\x
ERROR: MethodError: no method matching lu!(::HssMatrix{ComplexF64}, ::Val{true}; check=true)
Closest candidates are:
  lu!(::StridedMatrix{T}, ::Union{Val{false}, Val{true}}; check) where T<:Union{Float32, Float64, ComplexF32, ComplexF64} at /builddir/build/BUILD/julia-1.6.1/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:79
  lu!(::Union{Hermitian{T, S}, Symmetric{T, S}} where {T, S}, ::Union{Val{false}, Val{true}}; check) at /builddir/build/BUILD/julia-1.6.1/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:88
  lu!(::StridedMatrix{T} where T, ::Union{Val{false}, Val{true}}; check) at /builddir/build/BUILD/julia-1.6.1/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:130
  ...
Stacktrace:
 [1] lu(A::HssMatrix{ComplexF64}, pivot::Val{true}; check::Bool)
   @ LinearAlgebra /builddir/build/BUILD/julia-1.6.1/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:273
 [2] lu(A::HssMatrix{ComplexF64}, pivot::Val{true}) (repeats 2 times)
   @ LinearAlgebra /builddir/build/BUILD/julia-1.6.1/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lu.jl:272
 [3] \(A::HssMatrix{ComplexF64}, B::Vector{ComplexF64})
   @ LinearAlgebra /builddir/build/BUILD/julia-1.6.1/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/generic.jl:1136
 [4] top-level scope
   @ REPL[131]:1

I have not been able to find a easy work-around (not that I have spend any reasonable time on it ). Is this something do you believe could be easily fixed?

Cheers, Mikkel

[bonevbs]

Hi Mikkel,

thanks for the interest. You're right - I have not really played around with complex linear algebra, I merely kept the door open so to say. I would be interested to hear your experience with this. I believe that HssMatrices.jl can be easily adapted to accommodate complex linear algebra.

On second inspection, yes this should be easily adapted... I believe that there was a type stability issue and that's why I disabled the Complex option for the ULV algorithm. I will have a look at it!

Cheers Boris

[mipals]

Hi again Boris,

I found that the reason the complex-solve did not work was simply due to the way the memory of the output was allocated. It was easily fixed as you mention :smile: (see my PR)

I also fixed so that you can solve with a Vector as a righ-hand side. The fix here is not pretty, but it works: Its a simple reshape, however it does mean that the output will be a Matrix and not a Vector. This can probably also be fixed, however I do not have time to look into it before a week or two. For the same reason I did not have time to write actual tests, but I added the PR anyways.

[bonevbs]

Hi Mikkel, cool! glad to hear that! I will have a look at your PR and the other one in the coming days :)

[mipals]

Ahh, I did not see that the existing PR dealt with the exact same issue 👍 I guess the only new thing about my PR then is the allowance of having a vector as a RHS 😄

[BaptisteLamic]

Hi, My PR should solve most of the issues for working with complex valued matrices. I am currently using it to solve quantum transport problem, where everything is complex valued.

But, there is still one test that does not pass with ComplexF32 ( have a look to the test file) at leat when using OpenBlas. Hence I never push my last version.

By the way there is typo in the code to get an individual element from a complex valued matrix ( extra conjugaison take place). I can push a fix to my PR for this last problem by the end of the week, I am a bit far from my computer for the moment. Best, Baptiste

Le ven. 13 août 2021 à 14:15, Mikkel Paltorp Schmitt < @.***> a écrit :

Ahh, I did not see that the existing PR dealt with the exact same issue 👍 I guess the only new thing about my PR then is the allowance of having a vector as a RHS 😄

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/bonevbs/HssMatrices.jl/issues/16#issuecomment-898417673, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEVAE64YHRTL3TTUQLF57CDT4UEGBANCNFSM5B6LISAQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email .

[mipals]

Hello Baptiste,

Sounds good. Sorry for overlooking your PR! No need to stress, enjoy your vacation (or whatever thing that is keeping you away from your computer 😉 )

[bonevbs]

Hi guys,

sorry for the silence - I am preparing to defend my thesis and had very little time to work on code. I should have some time soon to check the commits and incorporate them into main code.

Cheers, Boris

[bonevbs]

Closing as addressed by #17