BandedMatrices.jl Documentation
Creating banded matrices
BandedMatrices.BandedMatrix
— TypeBandedMatrix{T}(undef, (n, m), (l, u))
returns an uninitialized n
×m
banded matrix of type T
with bandwidths (l,u)
.
BandedMatrix{T}(kv::Pair, (m,n), (l,u))
Construct a m × n BandedMatrix with bandwidths (l,u) from Pair
s of diagonals and vectors. Vector kv.second
will be placed on the kv.first
diagonal.
BandedMatrix(kv::Pair{<:Integer,<:AbstractVector}...)
Construct a square matrix from Pair
s of diagonals and vectors. Vector kv.second
will be placed on the kv.first
diagonal.
BandedMatrices.brand
— Functionbrand(T,n,m,l,u)
Creates an n×m
banded matrix with random numbers in the bandwidth of type T
with bandwidths (l,u)
BandedMatrices.brandn
— Functionbrandn(T,n,m,l,u)
Creates an n×m
banded matrix with random normals in the bandwidth of type T
with bandwidths (l,u)
To create a banded matrix of all zeros, identity matrix, or with a constant value use the following constructors:
julia> BandedMatrix(Zeros(5,5), (1,2))
5×5 BandedMatrix{Float64} with bandwidths (1, 2):
0.0 0.0 0.0 ⋅ ⋅
0.0 0.0 0.0 0.0 ⋅
⋅ 0.0 0.0 0.0 0.0
⋅ ⋅ 0.0 0.0 0.0
⋅ ⋅ ⋅ 0.0 0.0
julia> BandedMatrix(Eye(5), (1,2))
5×5 BandedMatrix{Float64} with bandwidths (1, 2):
1.0 0.0 0.0 ⋅ ⋅
0.0 1.0 0.0 0.0 ⋅
⋅ 0.0 1.0 0.0 0.0
⋅ ⋅ 0.0 1.0 0.0
⋅ ⋅ ⋅ 0.0 1.0
julia> BandedMatrix(Ones(5,5), (1,2))
5×5 BandedMatrix{Float64} with bandwidths (1, 2):
1.0 1.0 1.0 ⋅ ⋅
1.0 1.0 1.0 1.0 ⋅
⋅ 1.0 1.0 1.0 1.0
⋅ ⋅ 1.0 1.0 1.0
⋅ ⋅ ⋅ 1.0 1.0
To create a banded matrix of a given size with constant bands (such as the classical finite difference approximation of the one-dimensional Laplacian on the unit interval [0,1]), you can use the following:
n = 128
h = 1/n
A = BandedMatrix{Float64}(undef, (n,n), (1,1))
A[band(0)] .= -2/h^2
A[band(1)] .= A[band(-1)] .= 1/h^2
Creating a large banded matrix from a dense matrix should be avoided because that costs time and memory:
julia> @time BandedMatrix(ones(10000,10000),(0,0));
0.775120 seconds (10 allocations: 763.016 MiB, 28.04% gc time)
Try to use structured matrices to get around this:
julia> @time BandedMatrix(Ones(10000,10000),(0,0));
0.000074 seconds (9 allocations: 78.469 KiB)
Another example:
julia> @time BandedMatrix([Ones(10000,5000) Zeros(10000,5000)],(1,1));
0.346918 seconds (22 allocations: 763.169 MiB, 8.38% gc time)
julia> using LazyArrays
julia> @time BandedMatrix(Hcat(Ones(10000,5000),Zeros(10000,5000)),(1,1));
0.012627 seconds (30.01 k allocations: 1.374 MiB, 92.24% gc time)
See LazyArrays, FillArrays for other implemented structured matrices.
Accessing banded matrices
BandedMatrices.bandwidths
— Functionbandwidths(A)
Returns a tuple containing the lower and upper bandwidth of A
, in order.
BandedMatrices.bandwidth
— Functionbandwidth(A,i)
Returns the lower bandwidth (i==1
) or the upper bandwidth (i==2
).
BandedMatrices.bandrange
— Functionbandrange(A)
Returns the range -bandwidth(A,1):bandwidth(A,2)
.
BandedMatrices.band
— Functionband(i)
Represents the i
-th band of a banded matrix.
julia> using BandedMatrices
julia> A = BandedMatrix(0=>1:4, 1=>5:7, -1=>8:10)
4×4 BandedMatrix{Int64} with bandwidths (1, 1):
1 5 ⋅ ⋅
8 2 6 ⋅
⋅ 9 3 7
⋅ ⋅ 10 4
julia> A[band(1)]
3-element Vector{Int64}:
5
6
7
julia> A[band(0)]
4-element Vector{Int64}:
1
2
3
4
julia> A[band(-1)]
3-element Vector{Int64}:
8
9
10
BandedMatrices.BandRange
— ConstantBandRange
Represents the entries in a row/column inside the bands.
julia> using BandedMatrices
julia> A = BandedMatrix(0=>1:4, 1=>5:7, -1=>8:10)
4×4 BandedMatrix{Int64} with bandwidths (1, 1):
1 5 ⋅ ⋅
8 2 6 ⋅
⋅ 9 3 7
⋅ ⋅ 10 4
julia> A[2, BandRange]
3-element Vector{Int64}:
8
2
6
BandedMatrices.isbanded
— Functionisbanded(A)
returns true if a matrix implements the banded interface.
BandedMatrices.BandSlice
— TypeBandSlice(band::Band, indices)
Represent a range of indices corresponding to a band.
Upon calling to_indices
, Band
s are converted to BandSlice
objects to represent the indices over which the Band
spans.
This mimics the relationship between Colon
and Base.Slice
.
Example
julia> B = BandedMatrix(0 => 1:4, 1=>1:3);
julia> bs = to_indices(B, (Band(1),))[1];
julia> bs isa BandedMatrices.BandSlice
true
julia> using LinearAlgebra
julia> bs == diagind(B, 1)
true
BandedMatrices.colstart
— Functioncolstart(A, i::Integer)
Return the starting row index of the filled bands in the i
-th column, bounded by the actual matrix size.
Examples
julia> A = BandedMatrix(0=>1:4, 1=>5:7)
4×4 BandedMatrix{Int64} with bandwidths (0, 1):
1 5 ⋅ ⋅
⋅ 2 6 ⋅
⋅ ⋅ 3 7
⋅ ⋅ ⋅ 4
julia> BandedMatrices.colstart(A, 3)
2
julia> BandedMatrices.colstart(A, 4)
3
BandedMatrices.colstop
— Functioncolstop(A, i::Integer)
Return the stopping row index of the filled bands in the i
-th column, bounded by the actual matrix size.
Examples
julia> A = BandedMatrix(0=>1:4, 1=>5:7)
4×4 BandedMatrix{Int64} with bandwidths (0, 1):
1 5 ⋅ ⋅
⋅ 2 6 ⋅
⋅ ⋅ 3 7
⋅ ⋅ ⋅ 4
julia> BandedMatrices.colstop(A, 3)
3
julia> BandedMatrices.colstop(A, 4)
4
BandedMatrices.colrange
— Functioncolrange(A, i::Integer)
Return the range of rows in the i
-th column that correspond to filled bands.
Examples
julia> A = BandedMatrix(0=>1:4, 1=>5:7)
4×4 BandedMatrix{Int64} with bandwidths (0, 1):
1 5 ⋅ ⋅
⋅ 2 6 ⋅
⋅ ⋅ 3 7
⋅ ⋅ ⋅ 4
julia> colrange(A, 1)
1:1
julia> colrange(A, 3)
2:3
BandedMatrices.collength
— Functioncollength(A, i::Integer)
Return the number of filled bands in the i
-th column.
Examples
julia> A = BandedMatrix(0=>1:4, 1=>5:7)
4×4 BandedMatrix{Int64} with bandwidths (0, 1):
1 5 ⋅ ⋅
⋅ 2 6 ⋅
⋅ ⋅ 3 7
⋅ ⋅ ⋅ 4
julia> BandedMatrices.collength(A, 1)
1
julia> BandedMatrices.collength(A, 2)
2
BandedMatrices.rowstart
— Functionrowstart(A, i::Integer)
Return the starting column index of the filled bands in the i
-th row, bounded by the actual matrix size.
Examples
julia> A = BandedMatrix(0=>1:4, 1=>5:7)
4×4 BandedMatrix{Int64} with bandwidths (0, 1):
1 5 ⋅ ⋅
⋅ 2 6 ⋅
⋅ ⋅ 3 7
⋅ ⋅ ⋅ 4
julia> BandedMatrices.rowstart(A, 2)
2
julia> BandedMatrices.rowstart(A, 3)
3
BandedMatrices.rowstop
— Functionrowstop(A, i::Integer)
Return the stopping column index of the filled bands in the i
-th row, bounded by the actual matrix size.
Examples
julia> A = BandedMatrix(0=>1:4, 1=>5:7)
4×4 BandedMatrix{Int64} with bandwidths (0, 1):
1 5 ⋅ ⋅
⋅ 2 6 ⋅
⋅ ⋅ 3 7
⋅ ⋅ ⋅ 4
julia> BandedMatrices.rowstop(A, 2)
3
julia> BandedMatrices.rowstop(A, 4)
4
BandedMatrices.rowrange
— Functionrowrange(A, i::Integer)
Return the range of columns in the i
-th row that correspond to filled bands.
Examples
julia> A = BandedMatrix(0=>1:4, 1=>5:7)
4×4 BandedMatrix{Int64} with bandwidths (0, 1):
1 5 ⋅ ⋅
⋅ 2 6 ⋅
⋅ ⋅ 3 7
⋅ ⋅ ⋅ 4
julia> rowrange(A, 1)
1:2
julia> rowrange(A, 4)
4:4
BandedMatrices.rowlength
— Functionrowlength(A, i::Integer)
Return the number of filled bands in the i
-th row.
Examples
julia> A = BandedMatrix(0=>1:4, 1=>5:7)
4×4 BandedMatrix{Int64} with bandwidths (0, 1):
1 5 ⋅ ⋅
⋅ 2 6 ⋅
⋅ ⋅ 3 7
⋅ ⋅ ⋅ 4
julia> BandedMatrices.rowlength(A, 1)
2
julia> BandedMatrices.rowlength(A, 4)
1
BandedMatrices.bandeddata
— Functionbandeddata(A)
returns a matrix containing the data of a banded matrix, in the BLAS format.
This is required for gbmv! support
BandedMatrices.BandedMatrixBand
— TypeBandedMatrixBand
Type to represent a view of a band of a BandedMatrix
Examples
julia> B = BandedMatrix(0=>1:3);
julia> view(B, band(0)) isa BandedMatrices.BandedMatrixBand
true
BandedMatrices.dataview
— Functiondataview(V::BandedMatrices.BandedMatrixBand)
Forward a view of a band of a BandedMatrix
to the parent's data matrix.
This will error if the indexing is out-of-bounds for the data matrix, even if it is inbounds for the parent BandedMatrix
Examples
julia> A = BandedMatrix(0=>1:4, 1=>5:7, -1=>8:10)
4×4 BandedMatrix{Int64} with bandwidths (1, 1):
1 5 ⋅ ⋅
8 2 6 ⋅
⋅ 9 3 7
⋅ ⋅ 10 4
julia> v = view(A, band(1));
julia> BandedMatrices.dataview(v)
3-element view(::Matrix{Int64}, 1, 2:4) with eltype Int64:
5
6
7
To loop over the nonzero elements of a BandedMatrix, you can use colrange(A, c)
and rowrange(A, r)
.
Creating symmetric banded matrices
Use Symmetric(::BandedMatrix)
to work with symmetric banded matrices.
Banded matrix interface
Banded matrices go beyond the type BandedMatrix
: one can also create matrix types that conform to the banded matrix interface, in which case many of the utility functions in this package are available. The banded matrix interface consists of the following:
Required methods | Brief description |
---|---|
bandwidths(A) | Returns a tuple containing the sub-diagonal and super-diagonal bandwidth |
BandedMatrices.isbanded(A) | Override to return true |
Optional methods | Brief description |
---|---|
inbands_getindex(A, k, j) | Unsafe: return A[k,j] , without the need to check if we are inside the bands |
inbands_setindex!(A, v, k, j) | Unsafe: set A[k,j] = v , without the need to check if we are inside the bands |
BandedMatrices.MemoryLayout(A) | Override to get banded lazy linear algebra, e.g. y .= Mul(A,x) |
BandedMatrices.bandeddata(A) | Override to return a matrix of the entries in BLAS format. Required if MemoryLayout(A) returns BandedColumnMajor |
Note that certain SubArray
s of BandedMatrix
are also banded matrices. The banded matrix interface is implemented for such SubArray
s to take advantage of this.
Eigenvalues
To compute efficiently a selection of eigenvalues for a BandedMatrix
, you may use any Krylov method that relies on a sequence of matrix * vector operations. For instance, using the package KrylovKit:
using KrylovKit
A = BandedMatrix(Eye(5), (1, 1))
KrylovKit.eigsolve(A, 1, :LR)
Implementation
Currently, only column-major ordering is supported: a banded matrix B
[ a_11 a_12 ⋅ ⋅
a_21 a_22 a_23 ⋅
a_31 a_32 a_33 a_34
⋅ a_42 a_43 a_44 ]
is represented as a BandedMatrix
with a field B.data
representing the matrix as
[ ⋅ a_12 a_23 a_34
a_11 a_22 a_33 a_44
a_21 a_32 a_43 ⋅
a_31 a_42 ⋅ ⋅ ]
B.l
gives the number of subdiagonals (2) and B.u
gives the number of super-diagonals (1). Both B.l
and B.u
must be non-negative at the moment.