Skip to content
This repository was archived by the owner on Jul 19, 2023. It is now read-only.

Commit 331f135

Browse files
jagotjagot
committed
Concretization cleanup
1 parent 795a54a commit 331f135

File tree

2 files changed

+28
-90
lines changed

2 files changed

+28
-90
lines changed

‎Project.toml‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,9 +24,10 @@ julia = "1"
2424
[extras]
2525
FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
2626
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
27+
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
2728
SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
2829
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
2930
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
3031

3132
[targets]
32-
test = ["FillArrays", "OrdinaryDiffEq", "SafeTestsets", "SpecialFunctions", "Test"]
33+
test = ["FillArrays", "OrdinaryDiffEq", "Random", "SafeTestsets", "SpecialFunctions", "Test"]

‎src/derivative_operators/concretization.jl‎

Lines changed: 26 additions & 89 deletions
Original file line numberDiff line numberDiff line change
@@ -1,123 +1,60 @@
1-
function LinearAlgebra.Array(A::DerivativeOperator{T}, N::Int=A.len) where T
2-
L = zeros(T, N, N+2)
1+
function Base.copyto!(L::AbstractMatrix{T}, A::DerivativeOperator{T}, N::Int) where T
32
bl = A.boundary_point_count
43
stencil_length = A.stencil_length
4+
stencil_pivot = use_winding(A) ? (1 + stencil_length%2) : div(stencil_length,2)
55
bstl = A.boundary_stencil_length
6+
67
coeff = A.coefficients
7-
8-
if use_winding(A)
9-
stencil_pivot = 1 + A.stencil_length%2
8+
get_coeff = if coeff isa AbstractVector
9+
i -> coeff[i]
10+
elseif coeff isa Number
11+
i -> coeff
1012
else
11-
stencil_pivot = div(stencil_length,2)
13+
i -> true
1214
end
1315

1416
for i in 1:bl
15-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
17+
cur_coeff = get_coeff(i)
1618
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.low_boundary_coefs[i]) : A.low_boundary_coefs[i]
1719
L[i,1:bstl] = cur_coeff * cur_stencil
1820
end
1921

2022
for i in bl+1:N-bl
21-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
23+
cur_coeff = get_coeff(i)
2224
stencil = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i] : A.stencil_coefs
2325
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(stencil) : stencil
2426
L[i,i+1-stencil_pivot:i-stencil_pivot+stencil_length] = cur_coeff * cur_stencil
2527
end
2628

2729
for i in N-bl+1:N
28-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
30+
cur_coeff = get_coeff(i)
2931
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.high_boundary_coefs[i-N+bl]) : A.high_boundary_coefs[i-N+bl]
3032
L[i,N-bstl+3:N+2] = cur_coeff * cur_stencil
3133
end
32-
return L
34+
35+
L
3336
end
3437

35-
function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}, N::Int=A.len) where T
36-
L = spzeros(T, N, N+2)
37-
bl = A.boundary_point_count
38-
stencil_length = A.stencil_length
39-
stencil_pivot = div(stencil_length,2)
40-
bstl = A.boundary_stencil_length
41-
coeff = A.coefficients
38+
LinearAlgebra.Array(A::DerivativeOperator{T}, N::Int=A.len) where T =
39+
copyto!(zeros(T, N, N+2), A, N)
4240

43-
if use_winding(A)
44-
stencil_pivot = 1 + A.stencil_length%2
45-
else
46-
stencil_pivot = div(stencil_length,2)
47-
end
41+
SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}, N::Int=A.len) where T =
42+
copyto!(spzeros(T, N, N+2), A, N)
4843

49-
for i in 1:A.boundary_point_count
50-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
51-
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.low_boundary_coefs[i]) : A.low_boundary_coefs[i]
52-
L[i,1:bstl] = cur_coeff * cur_stencil
53-
end
54-
55-
for i in bl+1:N-bl
56-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
57-
stencil = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i] : A.stencil_coefs
58-
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(stencil) : stencil
59-
L[i,i+1-stencil_pivot:i-stencil_pivot+stencil_length] = cur_coeff * cur_stencil
60-
end
61-
62-
for i in N-bl+1:N
63-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
64-
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.high_boundary_coefs[i-N+bl]) : A.high_boundary_coefs[i-N+bl]
65-
L[i,N-bstl+3:N+2] = cur_coeff * cur_stencil
66-
end
67-
return L
68-
end
69-
70-
function SparseArrays.sparse(A::DerivativeOperator{T}, N::Int=A.len) where T
71-
SparseMatrixCSC(A,N)
72-
end
44+
SparseArrays.sparse(A::DerivativeOperator{T}, N::Int=A.len) where T = SparseMatrixCSC(A,N)
7345

7446
function BandedMatrices.BandedMatrix(A::DerivativeOperator{T}, N::Int=A.len) where T
75-
bl = A.boundary_point_count
7647
stencil_length = A.stencil_length
7748
bstl = A.boundary_stencil_length
78-
coeff = A.coefficients
79-
if use_winding(A)
80-
stencil_pivot = 1 + A.stencil_length%2
81-
else
82-
stencil_pivot = div(stencil_length,2)
83-
end
8449
L = BandedMatrix{T}(Zeros(N, N+2), (max(stencil_length-3,0,bstl),max(stencil_length-1,0,bstl)))
85-
86-
for i in 1:A.boundary_point_count
87-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
88-
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.low_boundary_coefs[i]) : A.low_boundary_coefs[i]
89-
L[i,1:bstl] = cur_coeff * cur_stencil
90-
end
91-
for i in bl+1:N-bl
92-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
93-
stencil = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i] : A.stencil_coefs
94-
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(stencil) : stencil
95-
L[i,i+1-stencil_pivot:i-stencil_pivot+stencil_length] = cur_coeff * cur_stencil
96-
end
97-
for i in N-bl+1:N
98-
cur_coeff = typeof(coeff) <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
99-
cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.high_boundary_coefs[i-N+bl]) : A.high_boundary_coefs[i-N+bl]
100-
L[i,N-bstl+3:N+2] = cur_coeff * cur_stencil
101-
end
102-
return L
50+
copyto!(L, A, N)
10351
end
10452

105-
function Base.convert(::Type{Array},A::DerivativeOperator{T}) where T
106-
Array(A)
107-
end
53+
Base.convert(::Type{Mat}, A::DerivativeOperator) where {Mat<:Union{Array,SparseMatrixCSC,BandedMatrix}} =
54+
Mat(A)
10855

109-
function Base.convert(::Type{SparseMatrixCSC},A::DerivativeOperator{T}) where T
110-
SparseMatrixCSC(A)
111-
end
112-
113-
function Base.convert(::Type{BandedMatrix},A::DerivativeOperator{T}) where T
56+
Base.convert(::Type{AbstractMatrix},A::DerivativeOperator) =
11457
BandedMatrix(A)
115-
end
116-
117-
function Base.convert(::Type{AbstractMatrix},A::DerivativeOperator{T}) where T
118-
BandedMatrix(A)
119-
end
120-
12158

12259
################################################################################
12360
# Boundary Padded Array concretizations
@@ -179,10 +116,10 @@ end
179116

180117
function BandedMatrices.BandedMatrix(Q::AffineBC{T}, N::Int) where {T}
181118
Q_l = BandedMatrix{T}(Eye(N), (length(Q.a_r)-1, length(Q.a_l)-1))
182-
inbands_setindex!(Q_L, Q.a_l, 1, 1:length(Q.a_l))
183-
inbands_setindex!(Q_L, Q.a_r, N, (N-length(Q.a_r)+1):N)
119+
BandedMatrices.inbands_setindex!(Q_l, Q.a_l, 1, 1:length(Q.a_l))
120+
BandedMatrices.inbands_setindex!(Q_l, Q.a_r, N, (N-length(Q.a_r)+1):N)
184121
Q_b = [Q.b_l; zeros(T,N); Q.b_r]
185-
return (Q_L, Q_b)
122+
return (Q_l, Q_b)
186123
end
187124

188125
function SparseArrays.sparse(Q::AffineBC{T}, N::Int) where {T}
@@ -383,4 +320,4 @@ function BlockBandedMatrices.BandedBlockBandedMatrix(A::DerivativeOperator{T,N},
383320
B = Kron(BandedMatrix(Eye(n)), BandedMatrix(A))
384321
end
385322
return BandedBlockBandedMatrix(B)
386-
end
323+
end

0 commit comments

Comments
 (0)