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    LowerTriangularArray holds 10 LowerTriangularMatrix of size   in one array.

LowerTriangularMatrix is actually a vector

LowerTriangularMatrix and LowerTriangularArray can in many ways be used very much like a Matrix or Array, however, because they unravel the lower triangle into a vector their dimensionality is one less than their Array counterparts. A LowerTriangularMatrix should therefore be treated as a vector rather than a matrix with some (limited) added functionality to allow for matrix-indexing (vector or flat-indexing is the default though). More details below.

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, Array{...}}
 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, Array{...}}
 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, Array{...}}
 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, Array{...}}
 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, Array{...}}
 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> L
15-element, 5x5 LowerTriangularMatrix{Float64, Array{...}}
 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 = 0
0

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).

end doesn't work for matrix indexing

Indexing LowerTriangularMatrix and LowerTriangularArray in matrix style ([i, j]) with end doesn't work. It either returns an error or wrong results as the end is lowered by Julia to the size of the underlying flat array dimension.

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

julia> L[2, 1] = 0    # valid index
0

julia> L[1, 2] = 0    # invalid index in the upper triangle
ERROR: BoundsError: attempt to access 15-element, 5x5 LowerTriangularMatrix{Float64, Array{...}} at index [1, 2]

But reading from it will just return a zero

julia> L[2, 3]     # in the upper triangle
0.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, Array{...}}
 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 matrix
0.970635f0

julia> L[2, 1, 1] # (2,1) element of the first lower triangle matrix
0.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, Array{...}}
 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, Array{...}}
 0.34962   0.0       0.0
 0.342062  0.157181  0.0
 0.187859  0.764076  0.276006

julia> L * L
ERROR: 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.NoTangent, ::Any)
   @ ChainRulesCore ~/.julia/packages/ChainRulesCore/Vsbj9/src/tangent_arithmetic.jl:64
  *(::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(::BigFloat)
   @ Base mpfr.jl:592
  inv(::CoordinateTransformations.CylindricalFromCartesian)
   @ CoordinateTransformations ~/.julia/packages/CoordinateTransformations/q6Edc/src/coordinatesystems.jl:299
  inv(::Proj.PJ_DIRECTION)
   @ Proj ~/.julia/packages/Proj/e2sVT/src/coord.jl:356
  ...

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' * L
0.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, Array{...}}
 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, Array{...}}
 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  ) 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

With the rotation angle in degrees (positive eastward) and 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, Array{...}}
 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, Array{...}}
  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, Array{...}}
 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 + L
6-element, 3x3 LowerTriangularMatrix{Float32, Array{...}}
 0.69924   0.0       0.0
 0.684124  0.314362  0.0
 0.375717  1.52815   0.552012

julia> 2L
6-element, 3x3 LowerTriangularMatrix{Float32, Array{...}}
 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^2
6-element, 3x3 LowerTriangularMatrix{Float64, Array{...}}
 -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) # initailize
T4 Spectrum
├ lmax=5 (degrees)
├ mmax=5 (orders)
└ architecture: CPU{KernelAbstractions.CPU}

julia> L = rand(Float32, spectrum)
15-element, 5x5 LowerTriangularMatrix{Float32, Array{...}}
 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, Array{...}}
 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

LowerTriangularArrays.LowerTriangularArray Type

A lower triangular array implementation that only stores the non-zero entries explicitly. L<:AbstractArray{T,N-1} although we do allow both "flat" N-1-dimensional indexing and additional N-dimensional or "matrix-style" indexing.

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

LowerTriangularArrays.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

LowerTriangularArrays.LowerTriangularMatrix Type

2D LowerTriangularArray of type T

source

LowerTriangularArrays.LowerTriangularMatrix Method

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

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

source

LowerTriangularArrays.LowerTriangularMatrix Method

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

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

source

LowerTriangularArrays.OneBased Type

Abstract type to dispatch for 1-based indexing of the spherical harmonic degree l and order m, i.e. l=m=1 is the mean, the zonal modes are m=1 etc. This indexing matches Julia's 1-based indexing for arrays.

source

LowerTriangularArrays.Spectrum Type

Encodes the spectral trunction, orders and degrees of the spherical harmonics. Is used by every LowerTriangularArray and also defines the architecture on which the data of the LowerTriangularArray is stored.

source

LowerTriangularArrays.Spectrum Method

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

Create a Spectrum from the spectral truncation lmax and mmax. Both are assumed to be one-based, i.e. lmax=5 and mmax=5 will create a spectrum with T4 truncation.

source

LowerTriangularArrays.Spectrum Method

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

Create a Spectrum for the spectral truncation trunc. trunc is assumed to be zero-based, i.e. trunc=4 will create a Spectrum with T4 truncation. With one_degree_more==true the Spectrum wil have an lmax increased by one, which is needed for spectral gradients.

source

LowerTriangularArrays.Spectrum Method

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

Create a Spectrum from another Spectrum but with a new architecture.

source

LowerTriangularArrays.ZeroBased Type

Abstract type to dispatch for 0-based indexing of the spherical harmonic degree l and order m, i.e. l=m=0 is the mean, the zonal modes are m=0 etc. This indexing is more common in mathematics.

source

Base._reverse! Method

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

Reverse the field represented by L in latitude when the element in L represent the coefficients of the spherical harmonics. Reversal in latitude direction is obtained by flipping the sign of the odd harmonics. Reverses L in place.

source

Base._reverse! Method

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

Reverse the field represented by L in longitude (mirror at 0˚) when the element in L represent the coefficients of the spherical harmonics. Reversal in longitude direction is obtained by the complex harmonics. Reverses L in place.

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 Function

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

Size of a LowerTriangularArray defined as size of the flattened array if as <: AbstractVector and as if it were a full matrix when as <: AbstractMatrix.

source

LinearAlgebra.rotate! Method

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

Rotate the field(s) represented by a LowerTriangularArray zonally by degree by multiplication of the spherical harmonics by exp(-i_m_2π*degree/360), with m the order of the spherical harmonic.

source

LowerTriangularArrays.eachharmonic Method

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

creates unit_range::UnitRange to loop over all non-zeros in the LowerTriangularArrays 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

LowerTriangularArrays.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

LowerTriangularArrays.eachmatrix Method

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

Iterator for the non-horizontal dimensions in LowerTriangularArrays. Checks that the LowerTriangularArrays match according to lowertriangular_match.

source

LowerTriangularArrays.eachmatrix Method

eachmatrix(L::LowerTriangularArray) -> Any

Iterator for the non-horizontal dimensions in LowerTriangularArrays. To be used like

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

to loop over every non-horizontal dimension of L.

source

LowerTriangularArrays.eachorder Method

eachorder(L1::LowerTriangularArray) -> Any

Iterator for the order m, for each m return all ls, therefore the columns in the lower triangular matrix.

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

to loop over every order of L.

source

LowerTriangularArrays.find_L Method

L = find_L(Ls) returns the first LowerTriangularArray among the arguments. Adapted from Julia documentation of Broadcast interface

source

LowerTriangularArrays.get_2lm_range Method

get_2lm_range(m, lmax) -> Any

range of the doubled running indices 2lm in a l-column (degrees of spherical harmonics) given the column index m (order of harmonics)

source

LowerTriangularArrays.get_lm_range Method

get_lm_range(m, lmax) -> Any

range of the running indices lm in a l-column (degrees of spherical harmonics) given the column index m (order of harmonics)

source

LowerTriangularArrays.i2lm Method

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

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

source

LowerTriangularArrays.interpolate Method

interpolate(
    _::Type{NF},
    alms::LowerTriangularArray{T, N, ArrayType, S},
    ltrunc::Integer,
    mtrunc::Integer
) -> Any

Returns a LowerTriangularArray that is interpolated from alms to the size (ltrunc+1) x (mtrunc+1), both inputs are 0-based, by padding zeros for higher wavenumbers. If ltrunc or mtrunc are smaller than the corresponding size of alms than truncate is automatically called instead, returning a smaller LowerTriangularArray.

source

LowerTriangularArrays.lm2i Method

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

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

source

LowerTriangularArrays.lowertriangular_match Method

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

True if both L1 and L2 are of the same size (as matrix), but ignores singleton dimensions, e.g. 5x5 and 5x5x1 would match. With horizontal_only=true (default false) ignore the non-horizontal dimensions, e.g. 5x5, 5x5x1, 5x5x2 would all match.

source

LowerTriangularArrays.lowertriangular_match Method

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

True if all lower triangular matrices provided as arguments match according to lowertriangular_match wrt to L1 (and therefore all).

source

LowerTriangularArrays.truncate! Method

truncate!(
    A::AbstractMatrix
) -> LowerTriangularArray{T, 2, ArrayType} where {T, ArrayType<:AbstractMatrix{T}}

Sets the upper triangle of A to zero.

source

LowerTriangularArrays.truncate! Method

truncate!(
    alms::LowerTriangularArray,
    ltrunc::Integer,
    mtrunc::Integer
) -> LowerTriangularMatrix{T} where T

Triangular truncation to degree ltrunc and order mtrunc (both 0-based). Truncate spectral coefficients alms in-place by setting all coefficients for which the degree l is larger than the truncation ltrunc or order m larger than the truncaction mtrunc.

source

LowerTriangularArrays.truncate! Method

truncate!(
    alms::LowerTriangularArray,
    trunc::Integer
) -> LowerTriangularArray

Triangular truncation of alms to degree and order trunc in-place.

source

LowerTriangularArrays.truncate! Method

truncate!(
    alms::LowerTriangularArray
) -> LowerTriangularArray

Triangular truncation of alms to the size of it, sets additional rows to zero.

source

LowerTriangularArrays.truncate Method

truncate(
    _::Type{NF},
    alms::LowerTriangularArray{T, N, ArrayType, S},
    ltrunc::Integer,
    mtrunc::Integer
) -> Any

Returns a LowerTriangularArray that is truncated from alms to the size (ltrunc+1) x (mtrunc+1), both inputs are 0-based. If ltrunc or mtrunc is larger than the corresponding size ofalms than truncate is automatically called instead, returning a LowerTriangularArray padded zero coefficients for higher wavenumbers.

source

LowerTriangularArrays.zero_last_degree! Method

zero_last_degree!(L::LowerTriangularArray)

Zeros the largest degree (last row, l = lmax) of a LowerTriangularArray L. This sets all elements where l = lmax to zero

source