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

Commit f058ec5

Browse files
ajozefiakajozefiak
committed
fixed concretization types
1 parent efb1233 commit f058ec5

File tree

2 files changed

+63
-36
lines changed

2 files changed

+63
-36
lines changed

‎src/derivative_operators/concretization.jl‎

Lines changed: 29 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -275,7 +275,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N}, Mshape) where {T,N}
275275
end
276276
B = Kron(Eye(n), Array(A))
277277
end
278-
return BandedBlockBandedMatrix(B)
278+
return Array(B)
279279
end
280280

281281
function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N}, Mshape) where {T,N}
@@ -302,7 +302,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N}, Mshape) where
302302
end
303303
B = Kron(sparse(I,n,n), sparse(A))
304304
end
305-
return BandedBlockBandedMatrix(B)
305+
return sparse(B)
306306
end
307307

308308
function SparseArrays.sparse(A::DerivativeOperator{T,N}, Mshape) where {T,N}
@@ -325,6 +325,33 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N}, Mshape) where {
325325
B = Kron(Eye(n), B)
326326
end
327327

328+
# Case where A is differentiating along hte first dimension
329+
else
330+
n = 1
331+
for M_i in Mshape[2:end]
332+
n *= M_i
333+
end
334+
B = Kron(BandedMatrix(Eye(n)), BandedMatrix(A))
335+
end
336+
return BandedMatrix(B)
337+
end
338+
339+
function BlockBandedMatrices.BandedBlockBandedMatrix(A::DerivativeOperator{T,N}, Mshape) where {T,N}
340+
# Case where A is not differentiating along the first dimension
341+
if N != 1
342+
n = 1
343+
for M_i in Mshape[1:N-1]
344+
n *= M_i
345+
end
346+
B = Kron(BandedMatrix(A), Eye(n))
347+
if N != length(Mshape)
348+
n = 1
349+
for M_i in Mshape[N+1:end]
350+
n *= M_i
351+
end
352+
B = Kron(Eye(n), B)
353+
end
354+
328355
# Case where A is differentiating along hte first dimension
329356
else
330357
n = 1

‎test/concretization.jl‎

Lines changed: 34 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
using SparseArrays, DiffEqOperators, LinearAlgebra, Random,
2-
Test, BandedMatrices, FillArrays, LazyArrays
2+
Test, BandedMatrices, FillArrays, LazyArrays, BlockBandedMatrices
33

44
# This test file tests for the correctness of higher dimensional concretization.
55
# The tests verify that multiplication in the concretized case agrees with the matrix-free
@@ -15,15 +15,15 @@ using SparseArrays, DiffEqOperators, LinearAlgebra, Random,
1515
L2 = CenteredDifference(1,2,1.0,20)
1616
L3 = CenteredDifference(3,3,1.0,20)
1717

18-
@test L1*M ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*M ≈ BandedMatrix(L1, size(M))*M
19-
@test L2*M ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*M ≈ BandedMatrix(L2, size(M))*M
20-
@test L3*M ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*M ≈ BandedMatrix(L3, size(M))*M
18+
@test L1*M ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*M ≈ BandedMatrix(L1, size(M))*M ≈ BandedBlockBandedMatrix(L1,size(M))*M
19+
@test L2*M ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*M ≈ BandedMatrix(L2, size(M))*M ≈ BandedBlockBandedMatrix(L2,size(M))*M
20+
@test L3*M ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*M ≈ BandedMatrix(L3, size(M))*M ≈ BandedBlockBandedMatrix(L3,size(M))*M
2121

2222
M = rand(22,2,2,2,2)
2323

24-
@test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M)
25-
@test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M)
26-
@test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M)
24+
@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)
25+
@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)
26+
@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)
2727

2828
end
2929

@@ -35,21 +35,21 @@ end
3535
L2 = CenteredDifference{2}(1,2,1.0,20)
3636
L3 = CenteredDifference{2}(3,3,1.0,20)
3737

38-
@test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M)
39-
@test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M)
40-
@test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M)
38+
@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)
39+
@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)
40+
@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)
4141

4242
M = rand(2,22,2,2)
4343

44-
@test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M)
45-
@test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M)
46-
@test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M)
44+
@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)
45+
@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)
46+
@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)
4747

4848
M = rand(2,22,2,2,3)
4949

50-
@test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M)
51-
@test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M)
52-
@test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M)
50+
@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)
51+
@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)
52+
@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)
5353

5454
end
5555

@@ -62,21 +62,21 @@ end
6262
L3 = CenteredDifference{3}(3,3,1.0,20)
6363

6464

65-
@test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M)
66-
@test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M)
67-
@test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M)
65+
@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)
66+
@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)
67+
@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)
6868

6969
M = rand(3,2,22,2)
7070

71-
@test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M)
72-
@test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M)
73-
@test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M)
71+
@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)
72+
@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)
73+
@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)
7474

7575
M = rand(3,2,22,2,3)
7676

77-
@test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M)
78-
@test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M)
79-
@test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M)
77+
@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)
78+
@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)
79+
@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)
8080

8181
end
8282

@@ -88,20 +88,20 @@ end
8888
L2 = CenteredDifference{5}(1,2,1.0,20)
8989
L3 = CenteredDifference{5}(3,3,1.0,20)
9090

91-
@test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M)
92-
@test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M)
93-
@test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M)
91+
@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)
92+
@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)
93+
@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)
9494

9595
M = rand(3,2,3,2,22,3)
9696

97-
@test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M)
98-
@test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M)
99-
@test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M)
97+
@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)
98+
@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)
99+
@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)
100100

101101
M = rand(2,2,2,2,22,3,2)
102102

103-
@test vec(L1*M) ≈ Array(L1, size(M))*vec(M) ≈ sparse(L1,size(M))*vec(M) ≈ BandedMatrix(L1, size(M))*vec(M)
104-
@test vec(L2*M) ≈ Array(L2, size(M))*vec(M) ≈ sparse(L2,size(M))*vec(M) ≈ BandedMatrix(L2, size(M))*vec(M)
105-
@test vec(L3*M) ≈ Array(L3, size(M))*vec(M) ≈ sparse(L3,size(M))*vec(M) ≈ BandedMatrix(L3, size(M))*vec(M)
103+
@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)
104+
@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)
105+
@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)
106106

107107
end

0 commit comments

Comments
 (0)