Concretization of L operators

Project.toml

@@ -10,6 +10,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 
 [compat]
 julia = "1"

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
+using SparseArrays, ForwardDiff, BandedMatrices
 
 abstract type AbstractDerivativeOperator{T} <: AbstractDiffEqLinearOperator{T} end
 abstract type AbstractDiffEqCompositeOperator{T} <: AbstractDiffEqLinearOperator{T} end
@@ -30,6 +30,7 @@ include("derivative_operators/derivative_irreg_operator.jl")
 include("derivative_operators/derivative_operator.jl")
 include("derivative_operators/abstract_operator_functions.jl")
 include("derivative_operators/convolutions.jl")
+include("derivative_operators/concretization.jl")
 
 ### Composite Operators
 include("composite_operators.jl")

src/derivative_operators/concretization.jl

@@ -1,36 +1,57 @@
-function Base.convert(::Type{Array}, A::AbstractDerivativeOperator{T}, N::Int=A.dimension) where T
-    @assert N >= A.stencil_length # stencil must be able to fit in the matrix
-    mat = zeros(T, (N, N+2))
-    v = zeros(T, N+2)
-    for i=1:N+2
-        v[i] = one(T)
-        #=
-            calculating the effect on a unit vector to get the matrix of transformation
-            to get the vector in the new vector space.
-        =#
-        mul!(view(mat,:,i), A, v)
-        v[i] = zero(T)
-    end
-    return mat
+function LinearAlgebra.Array(A::DerivativeOperator{T}) where T
+    N = A.dimension
+    L = zeros(T, N, N+2)
+    bl = A.boundary_length
+    stl = A.stencil_length
+    stl_2 = div(stl,2)
+    for i in 1:A.boundary_length
+        L[i,1:stl] = A.low_boundary_coefs[i]
+    end
+    for i in bl+1:N-bl
+        L[i,i+1-stl_2:i+1+stl_2] = A.stencil_coefs
+    end
+    for i in N-bl+1:N
+        L[i,N-stl+3:N+2] = A.high_boundary_coefs[i-N+bl]
+    end
+    return L / A.dx^A.derivative_order
 end
 
-function SparseArrays.sparse(A::AbstractDerivativeOperator{T}) where T
+function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}) where T
     N = A.dimension
-    mat = spzeros(T, N, N)
-    v = zeros(T, N)
-    row = zeros(T, N)
-    for i=1:N
-        v[i] = one(T)
-        #=
-            calculating the effect on a unit vector to get the matrix of transformation
-            to get the vector in the new vector space.
-        =#
-        mul!(row, A, v)
-        copyto!(view(mat,:,i), row)
-        @. row = 0 * row;
-        v[i] = zero(T)
-    end
-    return mat
+    L = spzeros(T, N, N+2)
+    bl = A.boundary_length
+    stl = A.stencil_length
+    stl_2 = div(stl,2)
+    for i in 1:A.boundary_length
+        L[i,1:stl] = A.low_boundary_coefs[i]
+    end
+    for i in bl+1:N-bl
+        L[i,i+1-stl_2:i+1+stl_2] = A.stencil_coefs
+    end
+    for i in N-bl+1:N
+        L[i,N-stl+3:N+2] = A.high_boundary_coefs[i-N+bl]
+    end
+    return L / A.dx^A.derivative_order
 end
 
-# BandedMatrix{A,B,C,D}(A::DerivativeOperator{A,B,C,D}) = BandedMatrix(convert(Array, A, A.stencil_length), A.stencil_length, div(A.stencil_length,2), div(A.stencil_length,2))
+function SparseArrays.sparse(A::AbstractDerivativeOperator{T}) where T
+    SparseMatrixCSC(A)
+end
+
+function BandedMatrices.BandedMatrix(A::DerivativeOperator{T}) where T
+    N = A.dimension
+    bl = A.boundary_length
+    stl = A.stencil_length
+    stl_2 = div(stl,2)
+    L = BandedMatrix{T}(Zeros(N, N+2), (max(stl-3,0),max(stl-1,0)))
+    for i in 1:A.boundary_length
+        L[i,1:stl] = A.low_boundary_coefs[i]
+    end
+    for i in bl+1:N-bl
+        L[i,i+1-stl_2:i+1+stl_2] = A.stencil_coefs
+    end
+    for i in N-bl+1:N
+        L[i,N-stl+3:N+2] = A.high_boundary_coefs[i-N+bl]
+    end
+    return L / A.dx^A.derivative_order
+end

src/derivative_operators/convolutions.jl

@@ -37,7 +37,7 @@ function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::
     coeffs = A.high_boundary_coefs
     Threads.@threads for i in 1 : A.boundary_length
         xtempi = coeffs[i][end]*x[end]
-        @inbounds for idx in A.stencil_length:-1:1
+        @inbounds for idx in A.stencil_length-1:-1:1
             xtempi += coeffs[i][end-idx] * x[end-idx]
         end
         x_temp[end-A.boundary_length+i] = xtempi

test/derivative_operators_interface.jl

@@ -1,4 +1,4 @@
-using SparseArrays, DiffEqOperators, LinearAlgebra, Random, Test
+using SparseArrays, DiffEqOperators, LinearAlgebra, Random, Test, BandedMatrices
 
 function second_derivative_stencil(N)
   A = zeros(N,N+2)
@@ -30,15 +30,37 @@ end
     N = 10
     d_order = 2
     approx_order = 2
-    x = collect(1:1.0:N).^2
     correct = second_derivative_stencil(N)
     A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N)
 
     @test convert_by_multiplication(Array,A,N) == correct
-    @test_broken convert(Array, A, N) == second_derivative_stencil(N)
-    @test_broken sparse(A) == second_derivative_stencil(N)
+    @test Array(A) == second_derivative_stencil(N)
+    @test sparse(A) == second_derivative_stencil(N)
+    @test BandedMatrix(A) == second_derivative_stencil(N)
     @test_broken opnorm(A, Inf) == opnorm(correct, Inf)
 
+
+    # testing higher derivative and approximation concretization
+    N = 20
+    d_order = 4
+    approx_order = 4
+    A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N)
+    correct = convert_by_multiplication(Array,A,N)
+
+    @test Array(A) == correct
+    @test sparse(A) == correct
+    @test BandedMatrix(A) == correct
+
+    N = 40
+    d_order = 8
+    approx_order = 8
+    A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N)
+    correct = convert_by_multiplication(Array,A,N)
+
+    @test Array(A) == correct
+    @test sparse(A) == correct
+    @test BandedMatrix(A) == correct
+
     # testing correctness of multiplication
     N = 1000
     d_order = 4