GitHub

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 LowerTriangularMatrix, a lower triangular matrix, which in contrast to LinearAlgebra.LowerTriangular does not store the entries above the diagonal. SpeedyWeather.jl uses LowerTriangularMatrix which is defined as a subtype of AbstractMatrix to store the spherical harmonic coefficients (see Spectral packing).

Creation of LowerTriangularMatrix

A LowerTriangularMatrix 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.155149   0.0       0.0       0.0       0.0
 0.127645   0.749576  0.0       0.0       0.0
 0.862919   0.550036  0.610428  0.0       0.0
 0.0363234  0.505492  0.905685  0.716618  0.0
 0.532488   0.990047  0.070392  0.756007  0.820987

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 Matrix, which drops the upper triangle entries and sets them to zero

julia> M = rand(Float16, 3, 3)3×3 Matrix{Float16}:
 0.9146   0.6567  0.58
 0.5957   0.802   0.4072
 0.02832  0.0928  0.5625
julia> L = LowerTriangularMatrix(M)3×3 LowerTriangularMatrix{Float16}:
 0.9146   0.0     0.0
 0.5957   0.802   0.0
 0.02832  0.0928  0.5625

Indexing LowerTriangularMatrix

LowerTriangularMatrix 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.9146   0.0     0.0
 0.5957   0.802   0.0
 0.02832  0.0928  0.5625
julia> L[2, 2]Float16(0.802)

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

julia> L[4]Float16(0.802)

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::LowerTriangularMatrix 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 LowerTriangularMatrix, 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]

Linear algebra with LowerTriangularMatrix

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.840596  0.0       0.0
 0.499203  0.60787   0.0
 0.366153  0.169193  0.241796
julia> L + L3×3 LowerTriangularMatrix{Float32}:
 1.68119   0.0       0.0
 0.998405  1.21574   0.0
 0.732305  0.338386  0.483592
julia> L * L3×3 Matrix{Float32}:
 0.706602  0.0       0.0
 0.723078  0.369506  0.0
 0.480782  0.143758  0.0584654

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.

Function and type index

SpeedyWeather.LowerTriangularMatrices.LowerTriangularMatrix — Type
L = LowerTriangularMatrix{T}(v::Vector{T}, m::Int, n::Int)

A lower triangular matrix 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.

source
Base.fill! — Method
fill!(L::LowerTriangularMatrix, x)

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
SpeedyWeather.LowerTriangularMatrices.eachharmonic — Method
unit_range = eachharmonic(Ls::LowerTriangularMatrix...)

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.ij2k — Method
k = ij2k(   i::Integer,     # row index of matrix
            j::Integer,     # column index of matrix
            m::Integer)     # number of rows in matrix

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