CompactBases.jl icon indicating copy to clipboard operation
CompactBases.jl copied to clipboard

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

Open jagot opened this issue 4 years ago • 3 comments

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.

jagot avatar Feb 19 '21 10:02 jagot

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

dlfivefifty avatar Feb 19 '21 10:02 dlfivefifty

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

  • Hermitian mode where one pretends that left-vectors are simply the adjoints of right-vectors and the metric is dense (or some diagonal approximation; this will lower order of convergence), or
  • non-Hermitian mode where left-vectors are computed explicitly and the metric is the identity matrix by construction.
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 avatar Feb 19 '21 11:02 jagot

Maybe this already works, as-is?

jagot avatar Feb 19 '21 11:02 jagot