JuliaLinearAlgebra
Constructors from Kronecker Products?
Usually these types of matrices from from PDEs where they can be constructed from the Kronecker product of the actions along each dimension, like:
using SparseArrays
Iy = SparseMatrixCSC(I,N,N)
Ix = SparseMatrixCSC(I,N,N)
fJ = ones(3,3)
Dz = [1 0 0
0 0 0
0 0 0]
A = kron(Dz,Iy,sparse(Mx)) + kron(Dz,sparse(My),Ix) + kron(fJ,Iy,Ix)
Could there be a way like that last statement to directly form a BlockBandedMatrix or BandedBlockBandedMatrix?
What's Mx and My?
Note that there is support for Kronecker products via LazyBandedMatrices:
using LazyBandedMatrices, BandedMatrices, BlockBandedMatrices, FillArrays, LazyArrays
N = 2
Ix = Iy = BandedMatrix(Eye(N,N))
julia> BandedBlockBandedMatrix(Kron(Ix,Iy))
2×2-blocked 4×4 BandedBlockBandedMatrix{Float64,BlockArrays.PseudoBlockArray{Float64,2,Array{Float64,2},Tuple{BlockArrays.BlockedUnitRange{Array{Int64,1}},BlockArrays.BlockedUnitRange{Array{Int64,1}}}},BlockArrays.BlockedUnitRange{Array{Int64,1}}}:
1.0 ⋅ │ ⋅ ⋅
⋅ 1.0 │ ⋅ ⋅
──────────┼──────────
⋅ ⋅ │ 1.0 ⋅
⋅ ⋅ │ ⋅ 1.0
const Mx = Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1])
const My = copy(Mx)
Mx[2,1] = 2.0
Mx[end-1,end] = 2.0
My[1,2] = 2.0
My[end,end-1] = 2.0
Seems like there's some missing methods:
using LinearAlgebra, LazyBandedMatrices, BandedMatrices, BlockBandedMatrices, FillArrays, LazyArrays
N = 32
Mx = Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1])
My = copy(Mx)
Mx[2,1] = 2.0
Mx[end-1,end] = 2.0
My[1,2] = 2.0
My[end,end-1] = 2.0
Ix = BandedMatrix(Eye(N,N))
Iy = BandedMatrix(Eye(N,N))
Ix = SparseMatrixCSC(I,N,N)
fJ = ones(3,3)
Dz = [1 0 0
0 0 0
0 0 0]
A = BlockBandedMatrix(Kron(Dz,Iy,sparse(Mx))) +
BlockBandedMatrix(Kron(Dz,sparse(My),Ix)) +
BlockBandedMatrix(Kron(fJ,Iy,Ix))
MethodError: no method matching LazyBandedMatrices.BlockKron{Float64,A,B} where B where A(::Array{Int64,2}, ::BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}, ::SparseMatrixCSC{Float64,Int64})
Closest candidates are:
LazyBandedMatrices.BlockKron{Float64,A,B} where B where A(::AA, ::BB) where {T, AA, BB} at C:\Users\accou\.julia\packages\LazyBandedMatrices\l1EHv\src\LazyBandedMatrices.jl:331
LazyBandedMatrices.BlockKron(::ApplyArray{Float64,2,typeof(kron),Tuple{Array{Int64,2},BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}},SparseMatrixCSC{Float64,Int64}}}) at LazyBandedMatrices.jl:333
BlockSkylineMatrix{T,Array{T,1},BlockBandedMatrices.BlockSkylineSizes{Tuple{BlockArrays.BlockedUnitRange{Array{Int64,1}},BlockArrays.BlockedUnitRange{Array{Int64,1}}},Fill{Int64,1,Tuple{Base.OneTo{Int64}}},Fill{Int64,1,Tuple{Base.OneTo{Int64}}},BandedMatrix{Int64,Array{Int64,2},Base.OneTo{Int64}},Array{Int64,1}}} where T(::ApplyArray{Float64,2,typeof(kron),Tuple{Array{Int64,2},BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}},SparseMatrixCSC{Float64,Int64}}}) at LazyBandedMatrices.jl:359
top-level scope at test.jl:102
This example is trying to do https://github.com/JuliaDiffEq/DiffEqBenchmarks.jl/issues/91 but create a BandedBlockBandedMatrix, instead of doing a flat sparse matrix from SparsityDetection.jl
You don’t want dense arguments (e.g. Dz ) in the construction of a BandedBlockBandedMatrix
Also, is a triple kron correspondent to a 3D PDE? That would require a more sophisticated data structure...
It's essentially
u_t = u_xx + u_yy + f1(u,v,w)`
v_t = f2(u,v,w)
w_t = f3(u,v,w)
The dense matrix is then the sparsity pattern of f.
Yes that still requires a more sophisticated data structure: I believe if we interlace the rows one gets a banded-block-banded structure, which is what ApproxFun does but with a hacky implementation (InterlaceOperator). I had started a while ago to make an Interlace matrix but never wrapped it up
Got it. I thought it ended up as a block banded matrix in the end (but not banded block banded).
It should be blocks of the block banded matrix of the 2D Laplacian operator.