Added metric, metric_shape, matrix_element_metric

Author: jagot

src/CompactBases.jl

@@ -69,14 +69,47 @@ is simply [`locs`](@ref), i.e. the location of the quadrature nodes.
 """
 centers(B::BasisOrRestricted) = locs(B)
 
+"""
+    metric(B)
+
+Returns the metric or mass matrix of the basis `B`, equivalent to
+`S=B'B`.
+"""
+metric(B::BasisOrRestricted) = B'B
+
+"""
+    metric_shape(B)
+
+Returns the shape of the metric or mass matrix of the basis `B`,
+i.e. `Diagonal`, `BandedMatrix`, etc.
+"""
+metric_shape(B::BasisOrRestricted) = typeof(metric(B))
+
+function matrix_element_metric(B::BasisOrRestricted)
+    # For non-orthogonal bases, we need to apply the metric inverse,
+    # after computing the action of a matrix representation of a
+    # linear operator on a ket. However, when computing the inner
+    # product with the bra afterwards, we again use the metric. Here,
+    # we short-circuit this by combining them into them identity
+    # operator I, if the operator metric equals the metric. For
+    # orthogonal bases, with the integration weights not built into the
+    # coefficients, and an identity matrix operator metric, we /do/
+    # need to use the metric for the inner product.
+    S = metric(B)
+    oS = operator_metric(B)
+    oS == S ? I : (oS \ S)
+end
+
 export AbstractFiniteDifferences,
     FiniteDifferences, StaggeredFiniteDifferences, ImplicitFiniteDifferences,
     Derivative, dot, QuasiDiagonal, Inclusion, .., distribution,
     FEDVR, Derivative, @elem, dot,
     BSpline,
     ArbitraryKnotSet, LinearKnotSet, ExpKnotSet,
     order, numintervals, numfunctions, nonempty_intervals,
-    centers, vandermonde, FunctionProduct, Density,
+    centers, vandermonde,
+    metric, metric_shape, matrix_element_metric,
+    FunctionProduct, Density,
     LinearOperator, DiagonalOperator, ShiftAndInvert
 
 end # module

src/bsplines.jl

@@ -288,6 +288,8 @@ end
     end
 end
 
+metric(B::BSpline) = B.S
+
 # * Diagonal operators
 
 @materialize function *(Ac::AdjointBSplineOrRestricted,

src/fedvr.jl

@@ -201,6 +201,8 @@ inverse_weight(B::FEDVR, j) = B.n[j]
 inverse_weights(B::FEDVR) = B.n
 
 vandermonde(B::FEDVR) = BandedMatrix(0 => B.n)
+metric(B::FEDVROrRestricted{T}) where T = Diagonal(Ones{T}(size(B,2)))
+metric_shape(B::FEDVROrRestricted) = Diagonal
 
 # * Basis functions
 

test/inner_products.jl

@@ -4,7 +4,13 @@
     N = ceil(Int, rmax/ρ)
     k = 7
 
-    @testset "$(grid_type)" for grid_type = [:fd, :sfd, :implicit_fd, :fedvr, :bsplines]
+    @testset "$(grid_type)" for (grid_type,MS,MatElMetric) = [
+        (:fd, Diagonal, ρ*I),
+        (:sfd, Diagonal, ρ*I),
+        (:implicit_fd, Diagonal, ρ*I),
+        (:fedvr, Diagonal, I),
+        (:bsplines, BandedMatrix, I)
+    ]
         R = if grid_type == :fd
             FiniteDifferences(N, ρ)
         elseif grid_type == :sfd
@@ -51,6 +57,10 @@
         @test isapprox(f'g, 0.0, atol=1e-14)
         @test g'g isa Real
         @test isapprox(g'g, 0.5, rtol=1e-3)
+
+        @test metric(R) == S
+        @test metric_shape(R) <: MS
+        @test matrix_element_metric(R) == MatElMetric
     end
 
     @warn "Need to test inner products with restricted bases as well"