From 1196cd3d16218d73de9f54c46054e4e307657051 Mon Sep 17 00:00:00 2001 From: Adam Jozefiak Date: Mon, 10 Jun 2019 03:57:54 -0700 Subject: [PATCH 1/3] constructor overloads for Array, SparseMatricCSC, and BandedMatrix implemented for DerivativeOperator --- Project.toml | 1 + src/DiffEqOperators.jl | 2 +- .../derivative_operator.jl | 54 +++++++++++++++++++ 3 files changed, 56 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 938c6f503..9a3790e09 100644 --- a/Project.toml +++ b/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" diff --git a/src/DiffEqOperators.jl b/src/DiffEqOperators.jl index fe6b47444..11eb46542 100644 --- a/src/DiffEqOperators.jl +++ b/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 diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index c6a2a0e11..3f69f5a92 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -59,3 +59,57 @@ function LinearAlgebra.opnorm(A::DerivativeOperator{T,S}, p::Real=2) where {T,S} opnorm(convert(Array,A), p) end end + +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.SparseMatrixCSC(A::DerivativeOperator{T}) where T + N = A.dimension + 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 + +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}(undef, (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 From 9cc822005e75150807779c5602204f3cf48cb194 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 10 Jun 2019 08:02:29 -0400 Subject: [PATCH 2/3] some refactors --- src/derivative_operators/concretization.jl | 83 ++++++++++++------- .../derivative_operator.jl | 54 ------------ test/derivative_operators_interface.jl | 2 +- 3 files changed, 53 insertions(+), 86 deletions(-) diff --git a/src/derivative_operators/concretization.jl b/src/derivative_operators/concretization.jl index 2df10c8ae..0e3806859 100644 --- a/src/derivative_operators/concretization.jl +++ b/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}(undef, (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 diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index 3f69f5a92..c6a2a0e11 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -59,57 +59,3 @@ function LinearAlgebra.opnorm(A::DerivativeOperator{T,S}, p::Real=2) where {T,S} opnorm(convert(Array,A), p) end end - -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.SparseMatrixCSC(A::DerivativeOperator{T}) where T - N = A.dimension - 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 - -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}(undef, (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 diff --git a/test/derivative_operators_interface.jl b/test/derivative_operators_interface.jl index 98fe8f46e..344ca7391 100644 --- a/test/derivative_operators_interface.jl +++ b/test/derivative_operators_interface.jl @@ -35,7 +35,7 @@ end 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 Array(A) == second_derivative_stencil(N) @test_broken sparse(A) == second_derivative_stencil(N) @test_broken opnorm(A, Inf) == opnorm(correct, Inf) From 5c569d90acec190df308a2b31728c4f63c884352 Mon Sep 17 00:00:00 2001 From: Adam Jozefiak Date: Thu, 13 Jun 2019 12:12:50 -0700 Subject: [PATCH 3/3] concretization tests and fixed BandedMatrix initialization and convolve_right_BC --- src/DiffEqOperators.jl | 1 + src/derivative_operators/concretization.jl | 2 +- src/derivative_operators/convolutions.jl | 2 +- test/derivative_operators_interface.jl | 30 +++++++++++++++++++--- 4 files changed, 29 insertions(+), 6 deletions(-) diff --git a/src/DiffEqOperators.jl b/src/DiffEqOperators.jl index 11eb46542..3d484bd98 100644 --- a/src/DiffEqOperators.jl +++ b/src/DiffEqOperators.jl @@ -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") diff --git a/src/derivative_operators/concretization.jl b/src/derivative_operators/concretization.jl index 0e3806859..0843a45d9 100644 --- a/src/derivative_operators/concretization.jl +++ b/src/derivative_operators/concretization.jl @@ -43,7 +43,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T}) where T bl = A.boundary_length stl = A.stencil_length stl_2 = div(stl,2) - L = BandedMatrix{T}(undef, (N, N+2), (max(stl-3,0),max(stl-1,0))) + 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 diff --git a/src/derivative_operators/convolutions.jl b/src/derivative_operators/convolutions.jl index 3bf322b13..7285b3cc4 100644 --- a/src/derivative_operators/convolutions.jl +++ b/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 diff --git a/test/derivative_operators_interface.jl b/test/derivative_operators_interface.jl index 344ca7391..1a8e34785 100644 --- a/test/derivative_operators_interface.jl +++ b/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 Array(A) == 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