JuliaApproximation / CompactBases.jl

Julia library for function approximation with compact basis functions
MIT License
16 stars 2 forks source link

transpose(::Basis) as opposed to adjoint(::Basis) #36

Open jagot opened 3 years ago

jagot commented 3 years ago

For discretization of non-Hermitian operators in non-uniform finite-differences, it can be advantageous to work with separate left- and right-vectors, where the inner product is computed as transpose(left)*right (i.e. the first argument is not conjugated), instead of u'*S*u. The reason is that the metric S is dense, and it's more efficient to work with a biorthogonal set of vectors.

How should this be represented in the ContinuumArrays framework? I'm thinking something like this:

R = Basis(...)

u = ...
# The current way of computing the norm (inner product) would be
u'*R'*R*u
# or equivalently
dot(u, R'R, u) # R'R is dense

uc = (R'R)*u # Left vector
# The new way of computing norms/inner products
transpose(uc)*transpose(R)*R*u
# or equivalently
dotu(uc, u)

I am not entirely satisfied by this, ideally one would always write an inner product as dot(u, S, v), and this would figure out if u is a left-vector that should be transposed in the biorthogonal case or adjointed in the normal case.

dlfivefifty commented 3 years ago

Why not just add dotu(u, A, v)?

jagot commented 3 years ago

Basically, I envision being able to run the same calculation in two "modes":

V = Matrix{ComplexF64}(...) # Right-vectors
# Right vectors and metric
U,S = if mode == :hermitian
    adjoint(V), R'R
else
    (R'R)*V, Diagonal(...)
end

for i = 1:steps
    propagate_right_vectors!(V)
    mode == :non_hermitian && propagate_left_vectors!(U)

    u = transpose(view(U, m, :)) # Some left-vector
    v = view(V, :, n) # Some right-vector
    uv = dot(u, S, v) # Correct inner product, regardless of propagation mode
end
jagot commented 3 years ago

Maybe this already works, as-is?