LowerTriangularArrays

LowerTriangularArrays is a package that has been developed for SpeedyWeather.jl but can also be used standalone.

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 LowerTriangularArrays
julia> L = rand(LowerTriangularMatrix{Float32}, 5, 5)15-element, 5x5 LowerTriangularMatrix{Float32}
 0.431945  0.0       0.0       0.0       0.0
 0.373974  0.751113  0.0       0.0       0.0
 0.567648  0.884494  0.846708  0.0       0.0
 0.686385  0.269477  0.741976  0.412334  0.0
 0.927618  0.984433  0.683471  0.77768   0.805296
julia> L2 = rand(LowerTriangularArray{Float32}, 5, 5, 5)15×5 LowerTriangularArray{Float32, 2, Matrix{Float32}, Spectrum{CPU{KernelAbstractions.CPU}, Vector{UnitRange{Int64}}, Vector{Int64}}}
 0.9665051f0    0.09185016f0  0.65675944f0   0.6918121f0    0.13463825f0
 0.49905884f0   0.44240206f0  0.44376028f0   0.11803311f0   0.5451692f0
 0.18484777f0   0.5911841f0   0.5160752f0    0.5680251f0    0.18442726f0
 0.056351483f0  0.23988068f0  0.03368616f0   0.7145026f0    0.5064971f0
 0.043976188f0  0.11256951f0  0.46404934f0   0.8937672f0    0.22286618f0
 0.19179171f0   0.5260416f0   0.049563527f0  0.8212688f0    0.22976404f0
 0.3599372f0    0.7569856f0   0.7384545f0    0.13504118f0   0.7110862f0
 0.26300055f0   0.21832329f0  0.44815207f0   0.9149946f0    0.27490604f0
 0.52976906f0   0.16888243f0  0.77007663f0   0.60937274f0   0.4511652f0
 0.20589817f0   0.7690395f0   0.67822605f0   0.5953352f0    0.4827767f0
 0.49433088f0   0.27963996f0  0.61867696f0   0.6387465f0    0.98527884f0
 0.094744444f0  0.2357468f0   0.65639466f0   0.06884819f0   0.04322648f0
 0.11909282f0   0.423909f0    0.23065567f0   0.10755545f0   0.7532134f0
 0.24254924f0   0.5394788f0   0.28993505f0   0.9864153f0    0.31207544f0
 0.57483786f0   0.08467996f0  0.8466997f0    0.035477042f0  0.12895006f0

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

Note how for a matrix both the upper triangle and the lower triangle are shown in the terminal. The zeros are evident. However, for higher dimensional LowerTriangularArray we fall back to show the unravelled first two dimensions. Hence, here, the first column is the first matrix with 15 elements forming a 5x5 matrix, but the zeros are not shown.

Alternatively, it can be created through conversion from Array, which drops the upper triangle entries and sets them to zero (which are not stored however)

julia> M = rand(Float16, 3, 3)3×3 Matrix{Float16}:
 0.083   0.2485  0.871
 0.5195  0.413   0.888
 0.849   0.5127  0.1694
julia> L = LowerTriangularMatrix(M)6-element, 3x3 LowerTriangularMatrix{Float16}
 0.083   0.0     0.0
 0.5195  0.413   0.0
 0.849   0.5127  0.1694
julia> M2 = rand(Float16, 3, 3, 2)3×3×2 Array{Float16, 3}:
[:, :, 1] =
 0.5684   0.3345  0.6064
 0.9443   0.268   0.1855
 0.03076  0.6274  0.2998

[:, :, 2] =
 0.4185  0.9956  0.8374
 0.9023  0.2461  0.396
 0.896   0.6045  0.887
julia> L2 = LowerTriangularArray(M2)6×2 LowerTriangularArray{Float16, 2, Matrix{Float16}, Spectrum{CPU{KernelAbstractions.CPU}, Vector{UnitRange{Int64}}, Vector{Int64}}}
 Float16(0.5684)   Float16(0.4185)
 Float16(0.9443)   Float16(0.9023)
 Float16(0.03076)  Float16(0.896)
 Float16(0.268)    Float16(0.2461)
 Float16(0.6274)   Float16(0.6045)
 Float16(0.2998)   Float16(0.887)

Size of LowerTriangularArray

There are three different ways to describe the size of a LowerTriangularArray. For example with L

julia> L = rand(LowerTriangularMatrix, 5, 5)15-element, 5x5 LowerTriangularMatrix{Float64}
 0.303165  0.0        0.0        0.0        0.0
 0.571405  0.531802   0.0        0.0        0.0
 0.435465  0.546216   0.0375098  0.0        0.0
 0.736736  0.0669808  0.390149   0.992316   0.0
 0.223564  0.282781   0.686346   0.0873902  0.510838

we have (additional dimensions follow naturally thereafter)

1-based vector indexing (default)

julia> size(L)    # equivalently size(L, OneBased, as=Vector)(15,)

The lower triangle is unravelled hence the number of elements in the lower triangle is returned.

1-based matrix indexing

julia> size(L, as=Matrix)    # equivalently size(L, OneBased, as=Matrix)(5, 5)

If you think of a LowerTriangularMatrix as a matrix this is the most intuitive size of L, which, however, does not agree with the size of the underlying data array (hence it is not the default).

0-based matrix indexing

Because LowerTriangularArrays are used to represent the coefficients of spherical harmonics which are commonly indexed based on zero (i.e. starting with the zero mode representing the mean), we also add ZeroBased to get the corresponding size.

julia> size(L, ZeroBased, as=Matrix)(4, 4)

which is convenient if you want to know the maximum degree and order of the spherical harmonics in L. 0-based vector indexing is not implemented.

Indexing LowerTriangularArray

We illustrate the two types of indexing LowerTriangularArray supports.

• Matrix indexing, by denoting two indices, column and row [l, m, ..]
• Vector/flat indexing, by denoting a single index [lm, ..].

The matrix index works as expected

julia> L15-element, 5x5 LowerTriangularMatrix{Float64}
 0.303165  0.0        0.0        0.0        0.0
 0.571405  0.531802   0.0        0.0        0.0
 0.435465  0.546216   0.0375098  0.0        0.0
 0.736736  0.0669808  0.390149   0.992316   0.0
 0.223564  0.282781   0.686346   0.0873902  0.510838
julia> L[2, 2]0.5318024871866116

But the single index skips the zero entries in the upper triangle, i.e. a 2, 2 index points to the same element as the index 6

julia> L[6]0.5318024871866116

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

julia> Matrix(L)[6]0.0

which would point to the first element in the upper triangle (hence zero).

In performance-critical code a single index should be used, as this directly maps to the index of the underlying data vector. The matrix 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, as=Matrix)(5, 5)
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).

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 15-element, 5x5 LowerTriangularMatrix{Float64} at index [1, 2]

But reading from it will just return a zero

julia> L[2, 3]     # in the upper triangle0.0

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

julia> L = rand(LowerTriangularArray{Float32}, 3, 3, 5)6×5 LowerTriangularArray{Float32, 2, Matrix{Float32}, Spectrum{CPU{KernelAbstractions.CPU}, Vector{UnitRange{Int64}}, Vector{Int64}}}
 0.9869555f0   0.87066007f0  0.9969142f0   0.67151076f0  0.6004017f0
 0.970635f0    0.30506748f0  0.56721866f0  0.4967081f0   0.609935f0
 0.3664909f0   0.8996582f0   0.59946454f0  0.8767897f0   0.33051252f0
 0.39431113f0  0.13089699f0  0.9859782f0   0.1016112f0   0.40536368f0
 0.15392935f0  0.58854145f0  0.8007999f0   0.5794253f0   0.82851845f0
 0.7566889f0   0.40112936f0  0.37226796f0  0.479097f0    0.17071092f0
julia> L[2, 1] # second lower triangle element of the first lower triangle matrix0.970635f0
julia> L[2, 1, 1] # (2,1) element of the first lower triangle matrix0.970635f0

The setindex! functionality follows accordingly.

Iterators

An iterator over all entries in the array can be created with eachindex

julia> L = rand(LowerTriangularArray, 5, 5, 5)15×5 LowerTriangularArray{Float64, 2, Matrix{Float64}, Spectrum{CPU{KernelAbstractions.CPU}, Vector{UnitRange{Int64}}, Vector{Int64}}}
 0.8280740936385513    0.5812411578752604   …  0.4250565446576843
 0.4646394506757542    0.44719200978056106     0.09362839137291812
 0.017556819845855998  0.957190345355617       0.42129958893841013
 0.12504404765566912   0.9972983087211111      0.5217362720681946
 0.22830381071667571   0.13210639266014435     0.5168103343626186
 0.49533808837207793   0.6627232170239564   …  0.07078460253188268
 0.8967254291326145    0.07187048084790848     0.2266686226271949
 0.21695038928841825   0.8425128451250508      0.856932979921012
 0.8441835797734459    0.7087684552303365      0.7947851194019955
 0.9515659377335692    0.3993838006688908      0.909236529205353
 0.35263448782421036   0.2632528613262278   …  0.8321311400192318
 0.8743126223549444    0.41967821254312987     0.9134503613638271
 0.09978044442496603   0.23277425063764912     0.8293968047129032
 0.41407219489941705   0.23542047213499184     0.8674395228328761
 0.7831943920541231    0.7625876752245078      0.7106761171209531
julia> for ij in eachindex(L)
           # do something
       end
julia> eachindex(L)Base.OneTo(75)

In order to only loop over the harmonics (essentially the horizontal, ignoring other dimensions) use eachharmonic

julia> eachharmonic(L)Base.OneTo(15)

If you only want to loop over the other dimensions use eachmatrix

julia> eachmatrix(L)CartesianIndices((5,))

together they can be used as

julia> for k in eachmatrix(L)
           for lm in eachharmonic(L)
               L[lm, k]
           end
       end

Note that k is a CartesianIndex that will loop over all other dimensions, whether there's only 1 (representing a 3D variable) or 5 (representing a 6D variable with the first two dimensions being a lower triangular matrix).

Linear algebra with LowerTriangularArray

The LowerTriangularArrays module's main purpose is not linear algebra, and typical matrix operations will not work with LowerTriangularMatrix because it's treated as a vector not as a matrix, meaning that the following will not work as expected

julia> L = rand(LowerTriangularMatrix{Float32}, 3, 3)6-element, 3x3 LowerTriangularMatrix{Float32}
 0.34962   0.0       0.0
 0.342062  0.157181  0.0
 0.187859  0.764076  0.276006
julia> L * LERROR: MethodError: no method matching *(::LowerTriangularMatrix{Float32, Vector{Float32}, Spectrum{CPU{KernelAbstractions.CPU}, Vector{UnitRange{Int64}}, Vector{Int64}}}, ::LowerTriangularMatrix{Float32, Vector{Float32}, Spectrum{CPU{KernelAbstractions.CPU}, Vector{UnitRange{Int64}}, Vector{Int64}}})
The function `*` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:596
  *(::ChainRulesCore.NotImplemented, ::Any)
   @ ChainRulesCore ~/.julia/packages/ChainRulesCore/Vsbj9/src/tangent_arithmetic.jl:37
  *(::Any, ::ChainRulesCore.NotImplemented)
   @ ChainRulesCore ~/.julia/packages/ChainRulesCore/Vsbj9/src/tangent_arithmetic.jl:38
  ...
julia> inv(L)ERROR: MethodError: no method matching inv(::LowerTriangularMatrix{Float32, Vector{Float32}, Spectrum{CPU{KernelAbstractions.CPU}, Vector{UnitRange{Int64}}, Vector{Int64}}})
The function `inv` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  inv(::Crayons.ANSIStyle)
   @ Crayons ~/.julia/packages/Crayons/u3AH8/src/crayon.jl:76
  inv(::CoordinateTransformations.CartesianFromCylindrical)
   @ CoordinateTransformations ~/.julia/packages/CoordinateTransformations/q6Edc/src/coordinatesystems.jl:300
  inv(::CoordinateTransformations.CylindricalFromSpherical)
   @ CoordinateTransformations ~/.julia/packages/CoordinateTransformations/q6Edc/src/coordinatesystems.jl:301
  ...

And many other operations that require L to be a AbstractMatrix which it isn't. In contrast, typical vector operations like a scalar product between two "LowerTriangularMatrix" vectors does work

julia> L' * L0.9592282f0

Summation with sum follows the flat, single index logic

julia> L = rand(LowerTriangularArray{Float32}, 3, 3, 5)ERROR: UndefVarError: `LowerTriangularArray` not defined in `Main`
Suggestion: check for spelling errors or missing imports.
Hint: a global variable of this name also exists in LowerTriangularArrays.
Hint: a global variable of this name also exists in SpeedyWeather.
julia> sum(L, dims=2)ERROR: UndefVarError: `L` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

sums along the second dimension of the underlying vector, not of the full matrix representation.

Rotation of LowerTriangularArray

LowerTriangularArrays are used to describe spherical harmonics. In that case each element represents the coefficient in fron of the respective harmonic describing a field on the sphere when transformed to grid space. We implement rotate! (and rotate for an allocating version) for LowerTriangularArray to rotate these coefficients in complex number space to represent a longitude rotation of the represented grid space field. For example start with

julia> M = rand(LowerTriangularMatrix{ComplexF32}, 3, 3)6-element, 3x3 LowerTriangularMatrix{ComplexF32}
 0.391699+0.689666im       0.0+0.0im            0.0+0.0im
 0.129261+0.713228im   0.83777+0.835835im       0.0+0.0im
 0.848531+0.136385im  0.741117+0.963409im  0.481724+0.788217im

Now rotate!(::LowerTriangularArray, degree)

julia> rotate!(M, 45)6-element, 3x3 LowerTriangularMatrix{ComplexF32}
 0.391699+0.689666im      0.0+0.0im              0.0+0.0im
 0.129261+0.713228im  1.18342-0.00136858im       0.0+0.0im
 0.848531+0.136385im  1.20528+0.157184im    0.788217-0.481724im

represents the same (up to rounding errors from the rotation when not rotating by $\pm 90, \pm 180, ...$) field in grid space but rotated by 45˚ eastward. Note how the zonal modes (the first column) aren't rotated because they are zonally constant anyway (in fact their imaginary part can be dropped, but isn't here as created by the rand) and for the other modes this amounts to a multiplication with

\[\exp(-i\frac{2π}{360}\phi m)\]

With $\phi$ the rotation angle in degrees (positive eastward) and $m$ the zonal wavenumber (order of the spherical harmonic, the zero-based column index). Rotating again by 315˚ yields the original array

julia> rotate!(M, 315)6-element, 3x3 LowerTriangularMatrix{ComplexF32}
 0.391699+0.689666im       0.0+0.0im            0.0+0.0im
 0.129261+0.713228im   0.83777+0.835835im       0.0+0.0im
 0.848531+0.136385im  0.741117+0.963409im  0.481724+0.788217im

Reverse of LowerTriangularArray

A LowerTriangularArray is an AbstractArray, as such reverse and reverse! (in-place) are defined, reversing all elements of these arrays in the way how they are indexed with a single index. For LowerTriangularArra representing the coefficients of the spherical harmonics this does not make much sense, however, we describe here the functionality to reverse these arrays as they represent spherical harmonics, adding methods for dims=:lat for reversal in latitude direction and dims=:lon in longitude direction. Spherical harmonics are reversed in latitude by flipping the sign of the odd harmonics, which are the ones that are not symmetric around the equator. Spherical harmonics are reversed in longitude by taking the complex conjugate of every element as this flips the sign of the imaginary parts, which effectively mirrors the rotation of that harmonic around 0˚E.

julia> reverse(M, dims=:lat)6-element, 3x3 LowerTriangularMatrix{ComplexF32}
  0.391699+0.689666im        0.0+0.0im            0.0+0.0im
 -0.129261-0.713228im    0.83777+0.835835im       0.0+0.0im
  0.848531+0.136385im  -0.741117-0.963409im  0.481724+0.788217im

and in longitude

julia> reverse(M, dims=:lon)6-element, 3x3 LowerTriangularMatrix{ComplexF32}
 0.391699-0.689666im       0.0+0.0im            0.0+0.0im
 0.129261-0.713228im   0.83777-0.835835im       0.0+0.0im
 0.848531-0.136385im  0.741117-0.963409im  0.481724-0.788217im

Broadcasting with LowerTriangularArray

In contrast to linear algebra, many element-wise operations work as expected thanks to broadcasting, so operations that can be written in . notation whether implicit +, 2*, ... or explicitly written .+, .^, ... or via the @. macro

julia> L + L6-element, 3x3 LowerTriangularMatrix{Float32}
 0.69924   0.0       0.0
 0.684124  0.314362  0.0
 0.375717  1.52815   0.552012
julia> 2L6-element, 3x3 LowerTriangularMatrix{Float32}
 0.69924   0.0       0.0
 0.684124  0.314362  0.0
 0.375717  1.52815   0.552012
julia> @. L + 2L - 1.1*L / L^26-element, 3x3 LowerTriangularMatrix{Float64}
 -2.09741   0.0        0.0
 -2.18961  -6.52675    0.0
 -5.29189   0.852579  -3.1574

GPU

LowerTriangularArray{T, N, ArrayType, S} 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 (work in progress). 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 existing LowerTriangularArray.

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

The Spectrum type

Internally, a LowerTriangularArray is represented by an array that holds all non-zero elements of the matrices and a Spectrum type that holds all spectral discretization information and the architecture the array is on. The Spectrum can also be used to create new LowerTriangularArrays with the same spectral discretization:

julia> spectrum = Spectrum(5, 5) # initailizeT4 Spectrum
├ lmax=5 (degrees)
├ mmax=5 (orders)
└ architecture: CPU{KernelAbstractions.CPU}

julia> L = rand(Float32, spectrum)15-element, 5x5 LowerTriangularMatrix{Float32}
 0.548299  0.0       0.0       0.0        0.0
 0.41656   0.451309  0.0       0.0        0.0
 0.579934  0.06595   0.544848  0.0        0.0
 0.612964  0.207314  0.305646  0.163982   0.0
 0.915852  0.3305    0.413271  0.0686005  0.174597

julia> L = rand(ComplexF32, spectrum, 5)15×5 LowerTriangularArray{ComplexF32, 2, Matrix{ComplexF32}, Spectrum{CPU{KernelAbstractions.CPU}, Vector{UnitRange{Int64}}, Vector{Int64}}}
 0.77169997f0 + 0.5891573f0im    …   0.3805048f0 + 0.9890936f0im
  0.1407786f0 + 0.8091592f0im        0.6230468f0 + 0.8239662f0im
  0.7794726f0 + 0.055342555f0im      0.7439987f0 + 0.9099983f0im
 0.97537977f0 + 0.9082579f0im       0.17868459f0 + 0.8457592f0im
  0.2458018f0 + 0.994555f0im         0.3127923f0 + 0.13422251f0im
 0.76875776f0 + 0.82194805f0im   …   0.7871114f0 + 0.48232245f0im
 0.31201243f0 + 0.6664442f0im       0.29432726f0 + 0.21469063f0im
 0.27333224f0 + 0.08920771f0im      0.34725034f0 + 0.49068838f0im
    0.48824f0 + 0.8937826f0im        0.6192591f0 + 0.44708544f0im
  0.5070303f0 + 0.663254f0im        0.92056805f0 + 0.60454315f0im
 0.14705628f0 + 0.8031362f0im    …  0.12131572f0 + 0.44583803f0im
  0.9858015f0 + 0.5756321f0im       0.32028204f0 + 0.6841664f0im
 0.12962598f0 + 0.4497428f0im       0.87851703f0 + 0.07772893f0im
  0.3411851f0 + 0.91684014f0im      0.91468585f0 + 0.87955165f0im
 0.13270545f0 + 0.09332347f0im      0.07185066f0 + 0.85488623f0im

In the SpeedyWeather.jl model, the Spectrum is stored just once in the SpectralGrid type, and all LowerTriangularArrays are created with the same Spectrum. Therefore, once you've initialized the SpectralGrid, you can create LowerTriangularArrays with the same spectral discretization as follows:

julia> SG = SpectralGrid(trunc=31)ERROR: UndefVarError: `SpectralGrid` not defined in `Main`
Suggestion: check for spelling errors or missing imports.
Hint: a global variable of this name also exists in SpeedyWeather.
julia> L = rand(Float32, SG.spectrum)ERROR: UndefVarError: `SG` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

Function and type index

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

LowerTriangularMatrix(
    M::Array{T, 2}
) -> LowerTriangularMatrix{T, ArrayType, Spectrum{CPU{KernelAbstractions.CPU}, Vector{UnitRange{Int64}}, Vector{Int64}}} where {T, ArrayType<:Vector{T}}

LowerTriangularMatrix(
    M::Array{T, 2},
    spectrum::AbstractSpectrum
) -> LowerTriangularMatrix{T, ArrayType} where {T, ArrayType<:Vector{T}}

Spectrum(
    lmax::Integer,
    mmax::Integer;
    architecture
) -> Spectrum{CPU{KernelAbstractions.CPU}}

Spectrum(
    trunc::Integer;
    one_degree_more,
    kwargs...
) -> Spectrum{CPU{KernelAbstractions.CPU}}

Spectrum(
    spectrum::Spectrum;
    architecture
) -> Spectrum{CPU{KernelAbstractions.CPU}}

_reverse!(
    L::LowerTriangularArray,
    _::Val{:lat}
) -> LowerTriangularArray

_reverse!(
    L::LowerTriangularArray,
    _::Val{:lon}
) -> LowerTriangularArray

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

length(L::LowerTriangularArray) -> Any

size(L::LowerTriangularArray; ...) -> Any
size(
    L::LowerTriangularArray,
    base::Type{<:LowerTriangularArrays.IndexBasis};
    as
) -> Any

rotate!(
    L::LowerTriangularArray,
    degree::Real
) -> LowerTriangularArray

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

eachharmonic(L::LowerTriangularArray) -> Any

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

eachmatrix(L::LowerTriangularArray) -> Any

for k in eachmatrix(L)
    L[1, k]

eachorder(L1::LowerTriangularArray) -> Any

for lms in eachorder(L)
    for lm in lms
        L[lm] 
    end
end

get_2lm_range(m, lmax) -> Any

get_lm_range(m, lmax) -> Any

i2lm(k::Integer, mmax::Integer) -> Tuple{Any, Any}

lm2i(l::Integer, m::Integer, lmax::Integer) -> Any

lowertriangular_match(
    L1::LowerTriangularArray,
    L2::LowerTriangularArray;
    horizontal_only
) -> Any

lowertriangular_match(
    L1::LowerTriangularArray,
    Ls::LowerTriangularArray...;
    kwargs...
) -> Any