Fix indexing for PartialInverseOperator (#188)

* domain promotion for cached operator

* Fix PartialInverseOperator indexing

Author: jishnub

src/Operators/general/PartialInverseOperator.jl

@@ -1,5 +1,3 @@
-
-
 export PartialInverseOperator
 
 
@@ -23,18 +21,39 @@ rangespace(P::PartialInverseOperator)=domainspace(P.cache)
 domain(P::PartialInverseOperator)=domain(domainspace(P))
 bandwidths(P::PartialInverseOperator) = P.bandwidths
 
+# Compute the value at the (k,j)th index of inv(C), assumming that C is upper triangular
+function _getindexinv(C, k::Integer, j::Integer, ::Type{UpperTriangular})
+    j >= k || return zero(inv(one(eltype(C))))
+    j == k && return inv(C[k,k])
+    t = inv(C[k,k])
+    # k-th row of the inverse, starting from the diagonal
+    # but leaving out the j-th column
+    ret = zeros(eltype(t), j-k)
+    ret[1] = t
+    for m in k+1:j-1 # populate the k-th row of the inverse
+        t = zero(eltype(ret))
+        for i in k:j-1
+            t -= ret[i-k+1] * C[i,m]
+        end
+        ret[m - k + 1] = t/C[m,m]
+    end
+    t = zero(eltype(ret))
+    for (rind, i) in enumerate(k:j-1)
+        t -= ret[rind] * C[i,j]
+    end
+    t/C[j,j]
+end
+
 function getindex(P::PartialInverseOperator,k::Integer,j::Integer)
-    b = bandwidth(P.cache, 2)
+    b = bandwidth(P, 2)
     if k == j
-        return inv(P.cache[k,k])
+        inv(P.cache[k,k])
+    elseif j > k + b + 1
+        zero(eltype(P))
     elseif j > k
-        t = zero(T)
-        for i = max(k,j-b-1):j-1
-            t += ret[k,i]*P.cache[i,j]
-        end
-        return -t/P.cache[j,j]
+        _getindexinv(P.cache, k, j, UpperTriangular)
     else
-        return zero(eltype(P))
+        zero(eltype(P))
     end
 end
 

test/SpacesTest.jl

@@ -409,6 +409,27 @@ using LinearAlgebra
             end
         end
 
+        @testset "PartialInverseOperator" begin
+            @testset "sanity check" begin
+                A = UpperTriangular(rand(10, 10))
+                B = inv(A)
+                for I in CartesianIndices(B)
+                    @test B[I] ≈ ApproxFunBase._getindexinv(A, Tuple(I)..., UpperTriangular)
+                end
+            end
+            C = Conversion(Chebyshev(), Ultraspherical(1))
+            P = PartialInverseOperator(C, (0, 6))
+            Iapprox = (P * C)[1:10, 1:10]
+            @test all(isone, diag(Iapprox))
+            for k in axes(Iapprox,1), j in k + 1:min(k + bandwidths(P,2), size(Iapprox, 2))
+                @test Iapprox[k,j] ≈ 0 atol=eps(eltype(Iapprox))
+            end
+            B = AbstractMatrix(P[1:10, 1:10])
+            @testset for I in CartesianIndices(B)
+                @test B[I] ≈ P[Tuple(I)...] rtol=1e-8 atol=eps(eltype(B))
+            end
+        end
+
         @testset "istriu/istril" begin
             for D in Any[Derivative(Chebyshev()),
                     Conversion(Chebyshev(), Legendre()),