Fix weighted FunctionProducts in restricted bases

Author: jagot

src/densities.jl

@@ -75,10 +75,14 @@ function FunctionProduct{Conjugated}(L::BasisOrRestricted, R::BasisOrRestricted,
         # Orthogonal case
         RV = RV[indices(R,2),:]
         uniform = all(isone, RV.data)
+        w = (w isa AbstractMatrix ? BandedMatrix(w)[indices(R,2),indices(R,2)] : w)
         C = if uniform
             w
         else
-            RV*(w isa AbstractMatrix ? w[indices(R,2),indices(R,2)] : w)
+            RV*w
+        end
+        if C isa AbstractMatrix && bandwidths(C) == (0,0)
+            C = Diagonal(C)
         end
         tmp = C === I ? nothing : tmpvec(T, size(RV, 1), nρ)
         FunctionProduct{Conjugated}(ρρ,

src/finite_differences.jl

@@ -198,6 +198,8 @@ end
 schafer_corner_fix(ρ, Z) = 1/ρ^2*Z*ρ/8*(1 + Z*ρ)
 
 function log_lin_grid!(r, ρₘᵢₙ, ρₘₐₓ, α, js)
+    0 < ρₘᵢₙ < Inf && 0 < ρₘₐₓ < Inf ||
+        throw(DomainError((ρₘᵢₙ, ρₘₐₓ), "Log–linear grid steps need to be finite and larger than zero"))
     δρ = ρₘₐₓ-ρₘᵢₙ
     for j = js
         r[j] = r[j-1] + ρₘᵢₙ + (1-exp(-α*r[j-1]))*δρ

test/densities.jl

@@ -1,46 +1,48 @@
 @testset "Densities" begin
-    f = x -> sin(2π*x)
-    g = x -> x*exp(-x)
-    h = x -> sin(2π*x)*x*exp(-x)
-
-    rmin = 0.01
-    rmax = 10.0
-    ρ = 0.05
-    N = ceil(Int, rmax/ρ)
-    k = 7
-    Nn = 71
-
-    ρmin=rmin
-    ρmax=0.5
-    α=0.01
-
-    @testset "$R" for (R,kind,rtol) in [
-        (FiniteDifferences(N, ρ), :orthogonal_uniform, 1e-14),
-        (StaggeredFiniteDifferences(ρmin, ρmax, α, rmax, 0.0), :orthogonal_non_uniform, 1e-14),
-        (FEDVR(range(0, stop=rmax, length=Nn), k), :orthogonal_non_uniform, 1e-14),
-        (BSpline(LinearKnotSet(k, 0, rmax, Nn)), :non_orthogonal, 1e-5),
-    ]
-        r = axes(R,1)
-
-        cf = R \ f.(r)
-        cg = R \ g.(r)
-        ch = R \ h.(r)
-
-        ρ = Density(applied(*,R,cf), applied(*,R,cg))
-        ch2 = ρ.ρ
-
-        @test eltype(ρ) == eltype(cf)
-
-        if kind == :orthogonal_uniform
-            @test ρ.LV == I
-            @test ρ.RV == I
-            @test ρ.C == I
-        elseif kind == :orthogonal_non_uniform
-            @test ρ.LV == I
-            @test ρ.RV == I
-        end
+    @testset "Single functions" begin
+        rmin = 0.01
+        rmax = 10.0
+        δr = 0.05
+        N = ceil(Int, rmax/δr)
+        k = 7
+        Nn = 71
+
+        δrmin=rmin
+        δrmax=0.5
+        α=0.01
+
+        @testset "$R" for (R,kind,rtol) in [
+            (FiniteDifferences(N, δr), :orthogonal_uniform, 1e-14),
+            (StaggeredFiniteDifferences(δrmin, δrmax, α, rmax, 0.0), :orthogonal_non_uniform, 1e-14),
+            (FEDVR(range(0, stop=rmax, length=Nn), k), :orthogonal_non_uniform, 1e-14),
+            (BSpline(LinearKnotSet(k, 0, rmax, Nn)), :non_orthogonal, 1e-5),
+        ]
+            f = x -> sin(2π*x)
+            g = x -> x*exp(-x)
+            h = x -> sin(2π*x)*x*exp(-x)
+
+            r = axes(R,1)
+
+            cf = R \ f.(r)
+            cg = R \ g.(r)
+            ch = R \ h.(r)
+
+            ρ = Density(applied(*,R,cf), applied(*,R,cg))
+            ch2 = ρ.ρ
+
+            @test eltype(ρ) == eltype(cf)
+
+            if kind == :orthogonal_uniform
+                @test ρ.LV == I
+                @test ρ.RV == I
+                @test ρ.C == I
+            elseif kind == :orthogonal_non_uniform
+                @test ρ.LV == I
+                @test ρ.RV == I
+            end
 
-        @test ch2 ≈ ch atol=1e-14 rtol=rtol
+            @test ch2 ≈ ch atol=1e-14 rtol=rtol
+        end
     end
 
     @testset "Multiple functions" begin
@@ -86,6 +88,43 @@
         end
     end
 
+    @testset "Restricted bases" begin
+        rmin = 0.0
+        rmax = 10.0
+        δr = 0.05
+        N = ceil(Int, rmax/δr)
+        k = 7
+        Nn = 71
+
+        δrmin=0.01
+        δrmax=0.5
+        α=0.01
+
+        f = x -> sinpi(x/rmax)
+        g = x -> cospi(x/rmax)
+        @testset "$(label)" for (w,label) in [(one,"Unweighted"),(exp,"Weighted")]
+            h = x -> f(x)*w(x)*g(x)
+
+            @testset "$R" for (R,rtol) in [
+                (FiniteDifferences(N, δr), 1e-15),
+                (StaggeredFiniteDifferences(δrmin, δrmax, α, rmax, 0.0), 1e-15),
+                (FEDVR(range(0, stop=rmax, length=Nn), k), 1e-15),
+                (BSpline(LinearKnotSet(k, 0, rmax, Nn)), 1e-14),
+            ]
+                display(R)
+                R̃ = R[:,1:ceil(Int, size(R,2))]
+
+                r = axes(R̃,1)
+
+                u = R̃ \ f.(r)
+                v = R̃ \ g.(r)
+
+                ρ = Density(applied(*, R̃, u), applied(*, R̃, v), w=w)
+                @test ρ.ρ ≈ (R̃ \ h.(r)) rtol=rtol
+            end
+        end
+    end
+
     @testset "Pretty-printing" begin
         R = FiniteDifferences(99, 0.1)
         u = R*rand(size(R,2))