Conjugate transpose

[mipals]

Hi there,

I am currently using this nice package for various matrices in relation to the acoustical boundary element method. An example of such matrix is

$$ B = \begin{bmatrix} G(k,t_1,y_1) & G(k,t_1,y_2) & \dots & G(k,t1,y{NQ}) \ G(k,t_2,y_1) & G(k,t_2,y_2) & \dots & G(k,t2,y{NQ}) \ \vdots & \vdots & \ddots & \vdots\ G(k,t_n,y_1) & G(k,t_n,y_2) & \dots & G(k,tn,y{NQ}) \ \end{bmatrix}, $$

where G is the Green's function for the scalar Helmholtz equation. Notice that the system is not square, meaning that i need to solve the resulting linear system using e.g. lsqr. Doing so requires products with the conjugate transpose, which fails for the HMatrix format (due to getindex not being defined). Was there another reason than time why the conjugate transpose was left out? Would it not be "easily" by implemented by conjugating each of the blocks?

My workaround so far is to simply define another HMatrix as

$$ B^\text{H} = \begin{bmatrix} G(-k, t_1, y_1) & G(-k, t_2, y_1) & \dots & G(-k, tn, y{1}) \ G(-k, t_1, y_2) & G(-k, t_2, y2) & \dots & G(-k, t{n}, y_{2}) \ \vdots & \vdots & \ddots & \vdots\ G(-k, t1, y{NQ}) & G(-k, t2, y{NQ}) & \dots & G(-k, tn, y{NQ}) \ \end{bmatrix}, $$

but this of course doubles both assembly time and storage, so I would rather not have to do that 😄

Cheers, Mikkel

[maltezfaria]

Calling adjoint on an HMatrix should be supported: if that is not the case, we should add that 👍.

p.s. In your formula above, don't you also need to invert t and y for the adjoint?

[maltezfaria]

Ok, I just checked, and indeed while adjoint(H) creates a (lazy) adjoint of the HMatrix, the multiplication is not implemented. I will patch it up soon. And I am glad you got the error about getindex: otherwise the code falls back to some generic multiplication algorithm which would have been terribly slow for your problem 😉

[maltezfaria]

@mipals can you test with the branch in #50 see if it resolves your issue?

[mipals]

Hi @maltezfaria ,

Regarding the notion of t and y, then I am not so sure. From a purely an algebraic point of view the position of t and y do not matter as only the distance between the two is used. But I also get what you are saying. When I use FMM3D to do the above I interchange targets and sources.

That was quick. I unfortunately still run into an error of the adjoint not having a .partition. Which makes sense since the instead the partition is accessed as .parent.partition for the adjoint?

julia> Fh.H'*pI
ERROR: type Adjoint has no field partition
Stacktrace:
 [1] getproperty
   @ ./Base.jl:37 [inlined]
 [2] mul!(y::Vector{…}, A::Adjoint{…}, x::Vector{…}, a::Int64, b::Int64; global_index::Bool, threads::Bool)
   @ HMatrices ~/.julia/packages/HMatrices/NfIiG/src/multiplication.jl:404
 [3] mul! (repeats 2 times)
   @ ~/.julia/packages/HMatrices/NfIiG/src/multiplication.jl:365 [inlined]
 [4] *(A::Adjoint{ComplexF64, HMatrices.HMatrix{HMatrices.ClusterTree{3, Float64}, ComplexF64}}, x::Vector{ComplexF64})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:57
 [5] top-level scope
   @ REPL[21]:1

Cheers, Mikkel

[maltezfaria]

Regarding the notion of t and y, then I am not so sure. From a purely an algebraic point of view the position of t and y do not matter as only the distance between the two is used. But I also get what you are saying. When I use FMM3D to do the above I interchange targets and sources.

AFAIU, purely based on the expected size of $B^H$, you need to swap targets and sources in the definition above.

That was quick. I unfortunately still run into an error of the adjoint not having a .partition. Which makes sense since the instead the partition is accessed as .parent.partition for the adjoint?

julia> Fh.H'*pI
ERROR: type Adjoint has no field partition
Stacktrace:
 [1] getproperty
   @ ./Base.jl:37 [inlined]
 [2] mul!(y::Vector{…}, A::Adjoint{…}, x::Vector{…}, a::Int64, b::Int64; global_index::Bool, threads::Bool)
   @ HMatrices ~/.julia/packages/HMatrices/NfIiG/src/multiplication.jl:404
 [3] mul! (repeats 2 times)
   @ ~/.julia/packages/HMatrices/NfIiG/src/multiplication.jl:365 [inlined]
 [4] *(A::Adjoint{ComplexF64, HMatrices.HMatrix{HMatrices.ClusterTree{3, Float64}, ComplexF64}}, x::Vector{ComplexF64})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:57
 [5] top-level scope
   @ REPL[21]:1

Cheers, Mikkel

My bad... I only tested the serial case, which does not run into that part of the code. I will fix it now. In the meantime, would you mind testing the serial version (run e.g. HMatrices.use_threads() = false to redefine the global option)

[mipals]

Hmm, are they not already swapped? I wrote it quite quickly so there might be a mistake. At least I think we understand each other 😄

The serial version does run! Comparing it to a dense calculation it seems like it also computes the correct thing 🥳

[maltezfaria]

Hmm, are they not already swapped? I wrote it quite quickly so there might be a mistake. At least I think we understand each other 😄

Oops, you are correct, they are already swapped! I got thrown off by their position as arguments not being swapped, but as you mentioned that does not matter since the function depends only on the distance. Sorry.

The serial version does run! Comparing it to a dense calculation it seems like it also computes the correct thing 🥳

I just pushed some changes that should make the parallel version work. I don't necessarily like the implementation, but at least we can test it. FWIW, it may be cleaner to implement $x^B H$, which uses only $H$, and then say that $H^B x = (x^B H)^B$. The problem with the current implementation is that I need to keep track of whether I am dealing with an adjoint or not throughout the multiplication if the partition is used. If I have time I will try to come up with a better implementation later.

[mipals]

Oops, you are correct, they are already swapped! I got thrown off by their position as arguments not being swapped, but as you mentioned that does not matter since the function depends only on the distance. Sorry.

No worries. I guess the interchanging of t and y will have an impact when I go for the normal derivative of the Green's function?

I just pushed some changes that should make the parallel version work. I don't necessarily like the implementation, but at least we can test it. FWIW, it may be cleaner to implement $x^BH$, which uses only $H$, and then say that $H^Bx=(x^BH)B$. The problem with the current implementation is that I need to keep track of whether I am dealing with an adjoint or not throughout the multiplication if the partition is used. If I have time I will try to come up with a better implementation later.

The parallel version does indeed seem to work 👍 Hmm, yeah, that might be cleaner. But no need to stress on getting it implemented from my end. I am just playing a little around 😄