Merge pull request #138 from JuliaDiffEq/kronecker

Higher Dimension Concretization

Author: ChrisRackauckas

Project.toml

@@ -5,8 +5,10 @@ version = "4.1.0"
 
 [deps]
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+LazyArrays =  "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

src/DiffEqOperators.jl

@@ -4,7 +4,7 @@ import Base: +, -, *, /, \, size, getindex, setindex!, Matrix, convert
 using DiffEqBase, StaticArrays, LinearAlgebra
 import LinearAlgebra: mul!, ldiv!, lmul!, rmul!, axpy!, opnorm, factorize, I
 import DiffEqBase: AbstractDiffEqLinearOperator, update_coefficients!, is_constant
-using SparseArrays, ForwardDiff, BandedMatrices, NNlib
+using SparseArrays, ForwardDiff, BandedMatrices, NNlib, LazyArrays, BlockBandedMatrices
 
 abstract type AbstractDerivativeOperator{T} <: AbstractDiffEqLinearOperator{T} end
 abstract type AbstractDiffEqCompositeOperator{T} <: AbstractDiffEqLinearOperator{T} end
@@ -15,7 +15,7 @@ include("matrixfree_operators.jl")
 include("jacvec_operators.jl")
 
 ### Utilities
-include("utils.jl") 
+include("utils.jl")
 
 ### Boundary Padded Arrays
 include("boundary_padded_arrays.jl")

src/derivative_operators/concretization.jl

@@ -50,7 +50,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}, N::Int=A.len) wh
     return L
 end
 
-function SparseArrays.sparse(A::AbstractDerivativeOperator{T}, N::Int=A.len) where T
+function SparseArrays.sparse(A::DerivativeOperator{T}, N::Int=A.len) where T
     SparseMatrixCSC(A,N)
 end
 
@@ -247,3 +247,118 @@ LinearAlgebra.Array(Q::ComposedMultiDimBC, Ns) = Tuple(Array.(Q.BCs, Ns))
 SparseArrays.SparseMatrixCSC(Q::ComposedMultiDimBC, Ns...) = Tuple(sparse.(Q.BCs, Ns))
 SparseArrays.sparse(Q::ComposedMultiDimBC, Ns) = SparseMatrixCSC(Q, Ns)
 BandedMatrices.BandedMatrix(Q::ComposedMultiDimBC, Ns) = Tuple(BandedMatrix.(Q.BCs, Ns))
+
+# HIgher Dimensional Concretizations. The following concretizations return two dimensional arrays
+# which operate on flattened vectors. Mshape is the size of the unflattened array on which A is operating on.
+
+function LinearAlgebra.Array(A::DerivativeOperator{T,N}, Mshape) where {T,N}
+    # Case where A is not differentiating along the first dimension
+    if N != 1
+        n = 1
+        for M_i in Mshape[1:N-1]
+            n *= M_i
+        end
+        B = Kron(Array(A), Eye(n))
+        if N != length(Mshape)
+            n = 1
+            for M_i in Mshape[N+1:end]
+                n *= M_i
+            end
+            B = Kron(Eye(n), B)
+        end
+
+    # Case where A is differentiating along hte first dimension
+    else
+        n = 1
+        for M_i in Mshape[2:end]
+            n *= M_i
+        end
+        B = Kron(Eye(n), Array(A))
+    end
+    return Array(B)
+end
+
+function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N}, Mshape) where {T,N}
+    # Case where A is not differentiating along the first dimension
+    if N != 1
+        n = 1
+        for M_i in Mshape[1:N-1]
+            n *= M_i
+        end
+        B = Kron(sparse(A), sparse(I,n,n))
+        if N != length(Mshape)
+            n = 1
+            for M_i in Mshape[N+1:end]
+                n *= M_i
+            end
+            B = Kron(sparse(I,n,n), B)
+        end
+
+    # Case where A is differentiating along hte first dimension
+    else
+        n = 1
+        for M_i in Mshape[2:end]
+            n *= M_i
+        end
+        B = Kron(sparse(I,n,n), sparse(A))
+    end
+    return sparse(B)
+end
+
+function SparseArrays.sparse(A::DerivativeOperator{T,N}, Mshape) where {T,N}
+    return SparseMatrixCSC(A,Mshape)
+end
+
+function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N}, Mshape) where {T,N}
+    # Case where A is not differentiating along the first dimension
+    if N != 1
+        n = 1
+        for M_i in Mshape[1:N-1]
+            n *= M_i
+        end
+        B = Kron(BandedMatrix(A), Eye(n))
+        if N != length(Mshape)
+            n = 1
+            for M_i in Mshape[N+1:end]
+                n *= M_i
+            end
+            B = Kron(Eye(n), B)
+        end
+
+    # Case where A is differentiating along hte first dimension
+    else
+        n = 1
+        for M_i in Mshape[2:end]
+            n *= M_i
+        end
+        B = Kron(BandedMatrix(Eye(n)), BandedMatrix(A))
+    end
+    return BandedMatrix(B)
+end
+
+function BlockBandedMatrices.BandedBlockBandedMatrix(A::DerivativeOperator{T,N}, Mshape) where {T,N}
+    # Case where A is not differentiating along the first dimension
+    if N != 1
+        n = 1
+        for M_i in Mshape[1:N-1]
+            n *= M_i
+        end
+        B = Kron(BandedMatrix(A), Eye(n))
+        if N != length(Mshape)
+            n = 1
+            for M_i in Mshape[N+1:end]
+                n *= M_i
+            end
+            B = Kron(Eye(n), B)
+        end
+
+    # Case where A is differentiating along hte first dimension
+    else
+        n = 1
+        for M_i in Mshape[2:end]
+            n *= M_i
+        end
+        B = Kron(BandedMatrix(Eye(n)), BandedMatrix(A))
+    end
+    return BandedBlockBandedMatrix(B)
+end

test/concretization.jl

@@ -0,0 +1,107 @@
+using SparseArrays, DiffEqOperators, LinearAlgebra, Random,
+      Test, BandedMatrices, FillArrays, LazyArrays, BlockBandedMatrices
+
+# This test file tests for the correctness of higher dimensional concretization.
+# The tests verify that multiplication in the concretized case agrees with the matrix-free
+# multiplication
+
+@testset "First Dimension" begin
+
+      # Test that even when we have a vector, the concretizations using the higher dimension dispatch still function
+      # correctly
+      M = rand(22)
+
+      L1 = CenteredDifference(1,2,1.0,20)
+      L2 = CenteredDifference(1,2,1.0,20)
+      L3 = CenteredDifference(3,3,1.0,20)
+
+      @test L1*M ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*M ≈ BandedMatrix(L1, size(M))*M ≈ BandedBlockBandedMatrix(L1,size(M))*M
+      @test L2*M ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*M ≈ BandedMatrix(L2, size(M))*M ≈ BandedBlockBandedMatrix(L2,size(M))*M
+      @test L3*M ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*M ≈ BandedMatrix(L3, size(M))*M ≈ BandedBlockBandedMatrix(L3,size(M))*M
+
+      M = rand(22,2,2,2,2)
+
+      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
+      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
+      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)
+
+end
+
+@testset "Second Dimension" begin
+
+      M = rand(2,22,2)
+
+      L1 = CenteredDifference{2}(1,2,1.0,20)
+      L2 = CenteredDifference{2}(1,2,1.0,20)
+      L3 = CenteredDifference{2}(3,3,1.0,20)
+
+      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
+      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
+      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)
+
+      M = rand(2,22,2,2)
+
+      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
+      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
+      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)
+
+      M = rand(2,22,2,2,3)
+
+      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
+      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
+      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)
+
+end
+
+@testset "Third Dimension" begin
+
+      M = rand(3,2,22)
+
+      L1 = CenteredDifference{3}(1,2,1.0,20)
+      L2 = CenteredDifference{3}(1,2,1.0,20)
+      L3 = CenteredDifference{3}(3,3,1.0,20)
+
+
+      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
+      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
+      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)
+
+      M = rand(3,2,22,2)
+
+      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
+      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
+      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)
+
+      M = rand(3,2,22,2,3)
+
+      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
+      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
+      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)
+
+end
+
+@testset "Fifth Dimension" begin
+
+      M = rand(3,2,3,2,22)
+
+      L1 = CenteredDifference{5}(1,2,1.0,20)
+      L2 = CenteredDifference{5}(1,2,1.0,20)
+      L3 = CenteredDifference{5}(3,3,1.0,20)
+
+      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
+      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
+      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)
+
+      M = rand(3,2,3,2,22,3)
+
+      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
+      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
+      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)
+
+      M = rand(2,2,2,2,22,3,2)
+
+      @test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L1,size(M))*vec(M)
+      @test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L2,size(M))*vec(M)
+      @test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M) ≈ BandedBlockBandedMatrix(L3,size(M))*vec(M)
+
+end

test/runtests.jl

@@ -16,3 +16,4 @@ import Base: isapprox
 @time @safetestset "Convolutions" begin include("convolutions.jl") end
 @time @safetestset "Differentiation Dimension" begin include("differentiation_dimension.jl") end
 @time @safetestset "2D and 3D fast multiplication" begin include("2D_3D_fast_multiplication.jl") end
+@time @safetestset "Higher Dimensional Concretization" begin include("concretization.jl") end