LowerTriangularMatrices

LowerTriangularMatrices is a submodule that has been developed for SpeedyWeather.jl which is technically independent (SpeedyWeather.jl however imports it and so does SpeedyTransforms) and can also be used without running simulations. It is just not put into its own respective repository.

This module defines LowerTriangularArray, a lower triangular matrix format, which in contrast to LinearAlgebra.LowerTriangular does not store the entries above the diagonal. SpeedyWeather.jl uses LowerTriangularArray which is defined as a subtype of AbstractArray to store the spherical harmonic coefficients (see Spectral packing). For 2D LowerTriangularArray the alias LowerTriangularMatrix exists. Higher dimensional LowerTriangularArray are 'batches' of 2D LowerTriangularMatrix. So, for example a $(10\times 10\times 10)$ LowerTriangularArray holds 10 LowerTriangularMatrix of size $(10\times 10)$ in one array.

Creation of `LowerTriangularArray`

A LowerTriangularMatrix and LowerTriangularArray can be created using zeros, ones, rand, or randn

julia> using SpeedyWeather.LowerTriangularMatrices
julia> L = rand(LowerTriangularMatrix{Float32}, 5, 5)5×5 LowerTriangularMatrix{Float32}:
 0.499474   0.0       0.0         0.0       0.0
 0.518371   0.413868  0.0         0.0       0.0
 0.827691   0.117217  0.355117    0.0       0.0
 0.998721   0.272952  0.452217    0.394454  0.0
 0.0650921  0.171766  0.00929779  0.852978  0.5402
julia> L2 = rand(LowerTriangularArray{Float32}, 5, 5, 5)5×5×5 LowerTriangularArray{Float32, 3, Matrix{Float32}}:
[:, :, 1] =
 0.290082  0.242957  0.803701  0.922044  0.900542
 0.982923  0.408152  0.063689  0.900542  0.176469
 0.76264   0.375047  0.922044  0.176469  0.410569
 0.401429  0.803701  0.900542  0.410569  0.898092
 0.242957  0.063689  0.176469  0.898092  0.311436

[:, :, 2] =
 0.177292  0.915217  0.742044  0.173701  0.550479
 0.275552  0.955204  0.749312  0.550479  0.480088
 0.419899  0.371084  0.173701  0.480088  0.435227
 0.281232  0.742044  0.550479  0.435227  0.909262
 0.915217  0.749312  0.480088  0.909262  0.702046

[:, :, 3] =
 0.779754   0.0204619  0.859188  0.831609  0.093742
 0.743511   0.437347   0.18614   0.093742  0.900595
 0.386988   0.246032   0.831609  0.900595  0.631347
 0.880061   0.859188   0.093742  0.631347  0.581569
 0.0204619  0.18614    0.900595  0.581569  0.48943

[:, :, 4] =
 0.850312   0.0372618  0.47023    0.427034   0.0507299
 0.325737   0.413884   0.0924453  0.0507299  0.573355
 0.886221   0.953995   0.427034   0.573355   0.621187
 0.145165   0.47023    0.0507299  0.621187   0.468975
 0.0372618  0.0924453  0.573355   0.468975   0.268479

[:, :, 5] =
 0.803995   0.683818   0.112888  0.62992   0.326685
 0.499729   0.766977   0.395839  0.326685  0.384049
 0.247414   0.0221771  0.62992   0.384049  0.616339
 0.0842742  0.112888   0.326685  0.616339  0.465186
 0.683818   0.395839   0.384049  0.465186  0.728698

or the undef initializor LowerTriangularMatrix{Float32}(undef, 3, 3). The element type is arbitrary though, you can use any type T too.

Alternatively, it can be created through conversion from Array, which drops the upper triangle entries and sets them to zero

julia> M = rand(Float16, 3, 3)3×3 Matrix{Float16}:
 0.8574  0.4556  0.1406
 0.4111  0.4336  0.0713
 0.621   0.995   0.971
julia> L = LowerTriangularMatrix(M)3×3 LowerTriangularMatrix{Float16}:
 0.8574  0.0     0.0
 0.4111  0.4336  0.0
 0.621   0.995   0.971
julia> M2 = rand(Float16, 3, 3, 5)3×3×5 Array{Float16, 3}:
[:, :, 1] =
 0.2866  0.547     0.651
 0.879   0.003906  0.3818
 0.8857  0.1675    0.2866

[:, :, 2] =
 0.332   0.1577  0.7764
 0.4263  0.1411  0.874
 0.5796  0.1113  0.1714

[:, :, 3] =
 0.1523  0.7075  0.001953
 0.703   0.1919  0.7095
 0.3496  0.671   0.7563

[:, :, 4] =
 0.417   0.4946  0.1582
 0.2134  0.4888  0.442
 0.2231  0.4116  0.4658

[:, :, 5] =
 0.6787  0.5815   0.774
 0.4385  0.02734  0.0635
 0.3267  0.3271   0.2798
julia> L2 = LowerTriangularArray(M)3×3 LowerTriangularMatrix{Float16}:
 0.8574  0.0     0.0
 0.4111  0.4336  0.0
 0.621   0.995   0.971

Indexing `LowerTriangularArray`

LowerTriangularArray supports two types of indexing: 1) by denoting two indices, column and row [l, m, ..] or 2) by denoting a single index [lm, ..]. The double index works as expected

julia> L3×3 LowerTriangularMatrix{Float16}:
 0.8574  0.0     0.0
 0.4111  0.4336  0.0
 0.621   0.995   0.971
julia> L[2, 2]Float16(0.4336)

But the single index skips the zero entries in the upper triangle, i.e.

julia> L[4]Float16(0.4336)

which, important, is different from single indices of an AbstractMatrix

julia> Matrix(L)[4]Float16(0.0)

In performance-critical code a single index should be used, as this directly maps to the index of the underlying data vector. The double index is somewhat slower as it first has to be converted to the corresponding single index.

Consequently, many loops in SpeedyWeather.jl are build with the following structure

julia> n, m = size(L)(3, 3)
julia> ij = 00
julia> for j in 1:m, i in j:n
           ij += 1
           L[ij] = i+j
       end

which loops over all lower triangle entries of L::LowerTriangularArray and the single index ij is simply counted up. However, one could also use [i, j] as indices in the loop body or to perform any calculation (i+j here). An iterator over all entries in the lower triangle can be created by

julia> for ij in eachindex(L)
           # do something
       end

The setindex! functionality of matrixes will throw a BoundsError when trying to write into the upper triangle of a LowerTriangularArray, for example

julia> L[2, 1] = 0    # valid index0
julia> L[1, 2] = 0    # invalid index in the upper triangleERROR: BoundsError: attempt to access 3×3 LowerTriangularMatrix{Float16} at index [1, 2]

Higher dimensional LowerTriangularArray can be indexed with multidimensional array indices like most other arrays types. Both the single index and the double index for the lower triangle work as shown here

julia> L = rand(LowerTriangularMatrix{Float32}, 3, 3, 5)ERROR: MethodError: no method matching Vector{Float32}(::Matrix{Float32})

Closest candidates are:
  Array{T, N}(::AbstractArray{S, N}) where {T, N, S}
   @ Base array.jl:673
  Vector{T}(::UndefInitializer, ::Tuple{Int64}) where T
   @ Core boot.jl:486
  Vector{T}(::UndefInitializer, ::Int64) where T
   @ Core boot.jl:477
  ...
julia> L[2, 1] # second lower triangle element of the first lower triangle matrixFloat16(0.0)
julia> L[2, 1, 1] # (2,1) element of the first lower triangle matrixERROR: BoundsError: attempt to access 3×3 LowerTriangularMatrix{Float16} at index [2, 1, 1]

The setindex! functionality follows accordingly.

Linear algebra with `LowerTriangularArray`

The LowerTriangularMatrices module's main purpose is not linear algebra, and it's implementation may not be efficient, however, many operations work as expected

julia> L = rand(LowerTriangularMatrix{Float32}, 3, 3)3×3 LowerTriangularMatrix{Float32}:
 0.740298  0.0        0.0
 0.410585  0.0190914  0.0
 0.233839  0.294438   0.992309
julia> L + L3×3 LowerTriangularMatrix{Float32}:
 1.4806    0.0        0.0
 0.82117   0.0381827  0.0
 0.467677  0.588875   1.98462
julia> L * L3×3 Matrix{Float32}:
 0.648516  0.183196  0.101929
 0.380645  0.183069  0.305633
 0.526042  0.352475  1.07583

Note, however, that the latter includes a conversion to Matrix, which is true for many operations, including inv or \. Hence when trying to do more sophisticated linear algebra with LowerTriangularMatrix we quickly leave lower triangular-land and go back to normal matrix-land.

GPU

LowerTriangularArray{T,N,ArrayType} wraps around an array of type ArrayType. If this array is a GPU array (e.g. CuArray), all operations are performed on GPU as well. The implementation was written so that scalar indexing is avoided in almost all cases, so that GPU operation should be performant. To use LowerTriangularArray on GPU you can e.g. just adapt an exisiting LowerTriangularArray.

using Adapt
L = rand(LowerTriangularArray{Float32}, 5, 5, 5)
L_gpu = adapt(CuArray, L)

Function and type index

SpeedyWeather.LowerTriangularMatrices.LowerTriangularArray — Type

A lower triangular array implementation that only stores the non-zero entries explicitly. L<:AbstractMatrix although in general we have L[i] != Matrix(L)[i], the former skips zero entries, tha latter includes them.

Supports n-dimensional lower triangular arrays, so that for all trailing dimensions L[:, :, ..] is a matrix in lower triangular form, e.g. a (5x5x3)-LowerTriangularArray would hold 3 lower triangular matrices.

source

SpeedyWeather.LowerTriangularMatrices.LowerTriangularArray — Method

LowerTriangularArray(
    M::AbstractArray{T, N}
) -> LowerTriangularArray

Create a LowerTriangularArray L from Array M by copying over the non-zero elements in M.

source

SpeedyWeather.LowerTriangularMatrices.LowerTriangularMatrix — Type

2-dimensional LowerTriangularArray of type Twith its non-zero entries unravelled into aVector{T}`

source

SpeedyWeather.LowerTriangularMatrices.LowerTriangularMatrix — Method

LowerTriangularMatrix(M::Array{T, 2}) -> Any

Create a LowerTriangularArray L from Matrix M by copying over the non-zero elements in M.

source

Base.fill! — Method

fill!(L::LowerTriangularArray, x) -> LowerTriangularArray

Fills the elements of L with x. Faster than fill!(::AbstractArray, x) as only the non-zero elements in L are assigned with x.

source

Base.length — Method

length(L::LowerTriangularArray) -> Any

Length of a LowerTriangularArray defined as number of non-zero elements.

source

Base.size — Method

size(
    L::LowerTriangularArray
) -> Tuple{Int64, Int64, Vararg{Any}}

Size of a LowerTriangularArray defined as matrix size including zero upper triangle.

source

SpeedyWeather.LowerTriangularMatrices.add! — Method

add!(
    L::LowerTriangularArray,
    L1::LowerTriangularArray,
    L2::LowerTriangularArray
) -> LowerTriangularArray

Adds L1 and L2 to L, so that L .+= L1 .+ L2

source

SpeedyWeather.LowerTriangularMatrices.add! — Method

add!(
    L::LowerTriangularArray,
    L1::LowerTriangularArray
) -> LowerTriangularArray

Adds L1 to L, so that L .+= L1

source

SpeedyWeather.LowerTriangularMatrices.eachharmonic — Method

eachharmonic(
    L1::LowerTriangularArray,
    Ls::LowerTriangularArray...
) -> Any

creates unit_range::UnitRange to loop over all non-zeros in the LowerTriangularMatrices provided as arguments. Checks bounds first. All LowerTriangularMatrix's need to be of the same size. Like eachindex but skips the upper triangle with zeros in L.

source

SpeedyWeather.LowerTriangularMatrices.eachharmonic — Method

eachharmonic(L::LowerTriangularArray) -> Any

creates unit_range::UnitRange to loop over all non-zeros/spherical harmonics numbers in a LowerTriangularArray L. Like eachindex but skips the upper triangle with zeros in L.

source

SpeedyWeather.LowerTriangularMatrices.ij2k — Method

ij2k(i::Integer, j::Integer, m::Integer) -> Any

Converts the index pair i, j of an mxn LowerTriangularMatrix L to a single index k that indexes the same element in the corresponding vector that stores only the lower triangle (the non-zero entries) of L.

source

SpeedyWeather.LowerTriangularMatrices.k2ij — Method

k2ij(k::Integer, m::Integer) -> Tuple{Any, Any}

Converts the linear index k in the lower triangle into a pair (i, j) of indices of the matrix in column-major form. (Formula taken from Angeletti et al, 2019, https://hal.science/hal-02047514/document)

source

LowerTriangularMatrices

Creation of LowerTriangularArray

Indexing LowerTriangularArray

Linear algebra with LowerTriangularArray

GPU

Function and type index

Creation of `LowerTriangularArray`

Indexing `LowerTriangularArray`

Linear algebra with `LowerTriangularArray`