copy _changebasis_R implementation to _changebasis_I

src/_ChangeBasis.jl

@@ -314,18 +314,150 @@ function _changebasis_I(P::BSplineSpace{0,T,KnotVector{T}}, P′::BSplineSpace{p
     A⁰ = sparse(view(I,1:s-1), view(J,1:s-1), fill(one(U), s-1), n, n′)
     return A⁰
 end
-function _changebasis_I(P::BSplineSpace{p,T}, P′::BSplineSpace{p′,T′}) where {p,p′,T,T′}
+
+function _changebasis_I(P::BSplineSpace{p,T,KnotVector{T}}, P′::BSplineSpace{p′,T′,KnotVector{T′}}) where {p,p′,T,T′}
+    U = StaticArrays.arithmetic_closure(promote_type(T,T′))
     k = knotvector(P)
     k′ = knotvector(P′)
+    n = dim(P)
+    n′ = dim(P′)
+    K′ = [k′[i+p′] - k′[i] for i in 1:n′+1]
+    K = U[ifelse(k[i+p] ≠ k[i], U(1 / (k[i+p] - k[i])), zero(U)) for i in 1:n+1]
+    Aᵖ⁻¹ = _changebasis_R(_lower_R(P), _lower_R(P′))  # (n+1) × (n′+1) matrix
+    n_nonzero = exactdim_R(P′)*(p+1)  # This is a upper bound of the number of non-zero elements of Aᵖ (rough estimation).
+    I = Vector{Int32}(undef, n_nonzero)
+    J = Vector{Int32}(undef, n_nonzero)
+    V = Vector{U}(undef, n_nonzero)
+    s = 1
+    j_begin = 1
+    j_end = 1
+    for i in 1:n
+        # Skip for degenerated basis
+        isdegenerate_R(P,i) && continue
 
-    _P = BSplineSpace{p}(k[1+p:end-p] + p * KnotVector{T}([k[1+p], k[end-p]]))
-    _P′ = BSplineSpace{p′}(k′[1+p′:end-p′] + p′ * KnotVector{T′}([k′[1+p′], k′[end-p′]]))
-    _A = _changebasis_R(_P, _P′)
-    Asim = _changebasis_sim(P, _P)
-    Asim′ = _changebasis_sim(_P′, P′)
-    A = Asim * _A * Asim′
+        # The following indices `j*` have the following relationships.
+        #=
 
-    return A
+         1                                    n′
+         |--*-----------------------------*-->|
+            j_begin                       j_end
+            |----*----------------*------>|
+                 j_prev  j_mid    j_next
+                 |       |        |
+                  |--j₊->||<-j₋--|
+                 266666666777777773
+                 366666666777777773
+
+         1                                    n′
+         |--*-----------------------------*-->|
+            j_begin                       j_end
+           *|---------------*------------>|
+           j_prev  j_mid    j_next
+           |       |        |
+            |--j₊->||<-j₋--|
+           066666666777777773
+
+         1                                    n′
+         |--*-----------------------------*-->|
+            j_begin                       j_end
+            |-------------*-------------->|*
+                          j_prev  j_mid    j_next
+                          |       |        |
+                           |--j₊->||<-j₋--|
+                          266666666777777770
+                          366666666777777770
+
+         1                                    n′
+         *-----------------------------*----->|
+         j_begin                       j_end
+         *----------------*----------->|
+        j_prev  j_mid    j_next
+        |       |        |
+         |--j₊->||<-j₋--|
+        066666666777777773
+
+         1                                    n′
+         |------*---------------------------->*
+                j_begin                       j_end
+                |-------------*--------------->*
+                              j_prev  j_mid    j_next
+                              |       |        |
+                               |--j₊->||<-j₋--|
+                              266666666777777770
+                              366666666777777770
+
+        =#
+
+        # Precalculate the range of j
+        Pi = BSplineSpace{p}(view(k, i:i+p+1))
+        j_end::Int = findnext(j->Pi ⊆ BSplineSpace{p′}(view(k′, j_begin:j+p′+1)), 1:n′, j_end)
+        j_begin::Int = findprev(j->Pi ⊆ BSplineSpace{p′}(view(k′, j:j_end+p′+1)), 1:n′, j_end)
+        j_range = j_begin:j_end
+        # flag = 0
+
+        # Rule-0: outside of j_range
+        j_prev = j_begin-1
+        Aᵖᵢⱼ_prev = zero(U)
+        for j_next in j_range
+            # Rule-1: zero
+            if k′[j_next] == k′[j_next+p′+1]
+                # flag = 1
+                continue
+            # Rule-2: right limit
+            elseif k′[j_next] == k′[j_next+p′]
+                I[s] = i
+                J[s] = j_next
+                V[s] = Aᵖᵢⱼ_prev = bsplinebasis₊₀(P,i,k′[j_next+1])
+                s += 1
+                j_prev = j_next
+                # flag = 2
+            # Rule-3: left limit (or both limit)
+            elseif k′[j_next+1] == k′[j_next+p′]
+                j_mid = (j_prev + j_next) ÷ 2
+                # Rule-6: right recursion
+                for j₊ in (j_prev+1):j_mid
+                    I[s] = i
+                    J[s] = j₊
+                    V[s] = Aᵖᵢⱼ_prev = Aᵖᵢⱼ_prev + p * K′[j₊] * (K[i] * Aᵖ⁻¹[i, j₊] - K[i+1] * Aᵖ⁻¹[i+1, j₊]) / p′
+                    s += 1
+                end
+                I[s] = i
+                J[s] = j_next
+                V[s] = Aᵖᵢⱼ_prev = Aᵖᵢⱼ_next = bsplinebasis₋₀(P,i,k′[j_next+1])
+                s += 1
+                # Rule-7: left recursion
+                for j₋ in reverse((j_mid+1):(j_next-1))
+                    I[s] = i
+                    J[s] = j₋
+                    V[s] = Aᵖᵢⱼ_next = Aᵖᵢⱼ_next - p * K′[j₋+1] * (K[i] * Aᵖ⁻¹[i, j₋+1] - K[i+1] * Aᵖ⁻¹[i+1, j₋+1]) / p′
+                    s += 1
+                end
+                j_prev = j_next
+                # flag = 3
+            end
+        end
+        # Rule-0: outside of j_range
+        j_next = j_end + 1
+        Aᵖᵢⱼ_next = zero(U)
+        j_mid = (j_prev + j_next) ÷ 2
+        # Rule-6: right recursion
+        for j₊ in (j_prev+1):j_mid
+            I[s] = i
+            J[s] = j₊
+            V[s] = Aᵖᵢⱼ_prev = Aᵖᵢⱼ_prev + p * K′[j₊] * (K[i] * Aᵖ⁻¹[i, j₊] - K[i+1] * Aᵖ⁻¹[i+1, j₊]) / p′
+            s += 1
+        end
+        # Rule-7: left recursion
+        for j₋ in reverse((j_mid+1):(j_next-1))
+            I[s] = i
+            J[s] = j₋
+            V[s] = Aᵖᵢⱼ_next = Aᵖᵢⱼ_next - p * K′[j₋+1] * (K[i] * Aᵖ⁻¹[i, j₋+1] - K[i+1] * Aᵖ⁻¹[i+1, j₋+1]) / p′
+            s += 1
+        end
+    end
+
+    Aᵖ = sparse(view(I,1:s-1), view(J,1:s-1), view(V,1:s-1), n, n′)
+    return Aᵖ
 end
 
 ## Uniform B-spline space