Inner products between functions on (un)restricted bases

Assume a linear operator Ĥ with a discrete and continuous
spectrum. Then any function may be expanded over its eigenfunctions:

When calculating inner products and matrix elements between such
expansions, one encounters inner products of the kinds

Usually, the eigenfunctions of the discrete spectrum have compact
support, or are at least very localized (they may have an
exponentially decaying tail, which is negligible). It would be nice to
be able to calculate such inner products, without having to represent
the discrete eigenfunctions on the entire interval. Maybe something
like this?:

B <: AbstractQuasiMatrix
B̃ = B[1:100]
i = B̃*rand(100)
k = B*rand(size(B,2))

# This is the syntax I'd like to automagically work:
i'k

As a follow-up to this, I'll need to solve the Poisson problem on the
entire interval, with source terms of the form

where the source term will be non-zero only on the interval where
is.

Start with the AbstractQuasiMatrix example and the desired i'k operation described in the issue. Investigate how restricted and unrestricted bases currently interact, then determine the scope for supporting these inner products and the related Poisson source terms. Done should include a decided API and validated behavior for both use cases.