BandedMatrices.jl Documentation

Creating banded matrices

BandedMatrices.BandedMatrix — Type

BandedMatrix{T}(undef, (n, m), (l, u))

returns an uninitialized n×m banded matrix of type T with bandwidths (l,u).

source

BandedMatrix{T}(kv::Pair, (m,n), (l,u))

Construct a m × n BandedMatrix with bandwidths (l,u) from Pairs of diagonals and vectors. Vector kv.second will be placed on the kv.first diagonal.

source

BandedMatrix(kv::Pair{<:Integer,<:AbstractVector}...)

Construct a square matrix from Pairs of diagonals and vectors. Vector kv.second will be placed on the kv.first diagonal.

source

BandedMatrices.brand — Function

brand(T,n,m,l,u)

Creates an n×m banded matrix with random numbers in the bandwidth of type T with bandwidths (l,u)

source

BandedMatrices.brandn — Function

brandn(T,n,m,l,u)

Creates an n×m banded matrix with random normals in the bandwidth of type T with bandwidths (l,u)

source

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 — Function

bandwidths(A)

Returns a tuple containing the lower and upper bandwidth of A, in order.

source

BandedMatrices.bandwidth — Function

bandwidth(A,i)

Returns the lower bandwidth (i==1) or the upper bandwidth (i==2).

source

BandedMatrices.bandrange — Function

bandrange(A)

Returns the range -bandwidth(A,1):bandwidth(A,2).

source

BandedMatrices.band — Function

band(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

source

BandedMatrices.BandRange — Constant

BandRange

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

source

BandedMatrices.isbanded — Function

isbanded(A)

returns true if a matrix implements the banded interface.

source

BandedMatrices.BandSlice — Type

BandSlice(band::Band, indices)

Represent a range of indices corresponding to a band.

Upon calling to_indices, Bands 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

source

BandedMatrices.colstart — Function

colstart(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

source

BandedMatrices.colstop — Function

colstop(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

source

BandedMatrices.colrange — Function

colrange(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

source

BandedMatrices.collength — Function

collength(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

source

BandedMatrices.rowstart — Function

rowstart(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

source

BandedMatrices.rowstop — Function

rowstop(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

source

BandedMatrices.rowrange — Function

rowrange(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

source

BandedMatrices.rowlength — Function

rowlength(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

source

BandedMatrices.bandeddata — Function

bandeddata(A)

returns a matrix containing the data of a banded matrix, in the BLAS format.

This is required for gbmv! support

source

BandedMatrices.BandedMatrixBand — Type

BandedMatrixBand

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

source

BandedMatrices.dataview — Function

dataview(V::BandedMatrices.BandedMatrixBand)

Forward a view of a band of a BandedMatrix to the parent's data matrix.

Warn

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

source

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.

Bandwidth minimization

The reverse Cuthill-McKee algorithm permutes a symmetric matrix into a band matrix with small bandwidth. This algorithm is implemented by the function symrcm.

julia> using BandedMatrices, CliqueTrees

julia> matrix = [
           1 0 1 0 0 0 0 1 1
           0 1 1 0 0 1 1 1 0
           1 1 1 0 0 0 0 0 0
           0 0 0 1 0 0 0 1 1
           0 0 0 0 1 0 1 1 0
           0 1 0 0 0 1 1 0 0
           0 1 0 0 1 1 1 0 0
           1 1 0 1 1 0 0 1 0
           1 0 0 1 0 0 0 0 1
       ];

julia> bandedmatrix, perm, invp = symrcm(matrix);

julia> bandedmatrix
9×9 LinearAlgebra.Symmetric{Int64, BandedMatrix{Int64, Matrix{Int64}, Base.OneTo{Int64}}}:
 1  1  1  0  ⋅  ⋅  ⋅  ⋅  ⋅
 1  1  1  1  0  ⋅  ⋅  ⋅  ⋅
 1  1  1  0  1  1  ⋅  ⋅  ⋅
 0  1  0  1  0  1  0  ⋅  ⋅
 ⋅  0  1  0  1  0  1  0  ⋅
 ⋅  ⋅  1  1  0  1  1  1  0
 ⋅  ⋅  ⋅  0  1  1  1  0  1
 ⋅  ⋅  ⋅  ⋅  0  1  0  1  1
 ⋅  ⋅  ⋅  ⋅  ⋅  0  1  1  1

julia> bandedmatrix == matrix[perm, perm]
true

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`
`BandedMatrices.bandedrowsdata(A)`	Override to return a matrix of the entries reshaped to a row-major format. Required if `MemoryLayout(A)` returns `BandedRowsMajor`

Note that certain SubArrays of BandedMatrix are also banded matrices. The banded matrix interface is implemented for such SubArrays 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.