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

using LowerTriangularArrays

L = rand(LowerTriangularMatrix{Float32}, 5, 5)
L2 = rand(LowerTriangularArray{Float32}, 5, 5, 5)

15×5 (LM+) LowerTriangularArray{Float32, 2, Array{...}}
 0.283031f0    0.6728255f0   0.6441826f0    0.23426145f0   0.45558852f0
 0.94413096f0  0.9830875f0   0.17703843f0   0.94946164f0   0.042922854f0
 0.26126927f0  0.86842453f0  0.6696464f0    0.9562146f0    0.30922252f0
 0.78367954f0  0.28755093f0  0.36937386f0   0.62999856f0   0.26224893f0
 0.8669858f0   0.71707296f0  0.4359895f0    0.98127156f0   0.5811741f0
 0.1672731f0   0.51947474f0  0.80713683f0   0.93548536f0   0.7179144f0
 0.56332207f0  0.15056598f0  0.06841022f0   0.097279966f0  0.54481596f0
 0.5962327f0   0.04488355f0  0.56018907f0   0.47993505f0   0.059660673f0
 0.23946643f0  0.5420443f0   0.69994724f0   0.486861f0     0.11177969f0
 0.8738056f0   0.38744575f0  0.5824845f0    0.7159746f0    0.8407996f0
 0.6532888f0   0.1558103f0   0.03211534f0   0.3215642f0    0.45011133f0
 0.63041574f0  0.70037895f0  0.22633296f0   0.9830787f0    0.05284667f0
 0.68313336f0  0.14320767f0  0.42466933f0   0.012978554f0  0.39536858f0
 0.27424031f0  0.3917879f0   0.012389481f0  0.37790138f0   0.23180753f0
 0.30372494f0  0.41088134f0  0.83429253f0   0.2553547f0    0.9878441f0

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)

M = rand(Float16, 3, 3)
L = LowerTriangularMatrix(M)

M2 = rand(Float16, 3, 3, 2)
L2 = LowerTriangularArray(M2)

6×2 (LM+) LowerTriangularArray{Float16, 2, Array{...}}
 Float16(0.1323)  Float16(0.5083)
 Float16(0.781)   Float16(0.797)
 Float16(0.641)   Float16(0.458)
 Float16(0.8545)  Float16(0.291)
 Float16(0.2944)  Float16(0.1284)
 Float16(0.647)   Float16(0.905)

Size of LowerTriangularArray

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

L = rand(LowerTriangularMatrix, 5, 5)

15-element, 5x5 LowerTriangularMatrix{Float32, Array{...}}
 0.570951  0.0        0.0        0.0       0.0
 0.630013  0.735788   0.0        0.0       0.0
 0.800062  0.0840679  0.0666391  0.0       0.0
 0.135972  0.284654   0.313994   0.819968  0.0
 0.121656  0.065399   0.859362   0.703236  0.885542

we have (additional dimensions follow naturally thereafter)

1-based vector indexing (default)

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

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.

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

L

L[2, 2]

0.73578805f0

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

L[6]

0.73578805f0

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

Matrix(L)[6]

0.0f0

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, see Iterators below.

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{Float32, Array{...}} at index [1, 2]

But reading from it will just return a zero

L[2, 3]     # in the upper triangle

0.0f0

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

L = rand(LowerTriangularArray{Float32}, 3, 3, 5)

L[2, 1] # second lower triangle element of the first lower triangle matrix

L[2, 1, 1] # (2,1) element of the first lower triangle matrix

0.07891822f0

The setindex! functionality follows accordingly.

Iterators

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

L = rand(LowerTriangularArray, 5, 5, 5)
for lm in eachindex(L)
    # do something
end

eachindex(L)

Base.OneTo(75)

In order to only loop over the harmonics (essentially the horizontal, ignoring other dimensions) use eachharmonic which will create iterators for both l, m via zip to be used like

L = zeros(LowerTriangularMatrix{Float32}, 3, 3)
for (l, m) in eachharmonic(L)
    L[l, m] = l+m
end
L

6-element, 3x3 LowerTriangularMatrix{Float32, Array{...}}
 2.0  0.0  0.0
 3.0  4.0  0.0
 4.0  5.0  6.0

But this needs to recalculate the running index lm for the non-zero harmonics, faster but otherwise identical is

for (lm, (l, m)) in enumerate(eachharmonic(L))
    L[lm] = l+m
end
L

6-element, 3x3 LowerTriangularMatrix{Float32, Array{...}}
 2.0  0.0  0.0
 3.0  4.0  0.0
 4.0  5.0  6.0

However note that this does not work on the GPU as scalar indexing generally does not work. Hence internally, SpeedyWeather will write such an operation either via broadcasting (see Broadcasting with LowerTriangularArray) if simple or with a custom kernel in general, see GPU and Architectures.

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

L = zeros(LowerTriangularArray, 3, 3, 3)
eachmatrix(L)

CartesianIndices((3,))

together they can be used as

for k in eachmatrix(L)
    for (lm, (l, m)) in enumerate(eachharmonic(L))
        L[lm, k]
    end
end

Note that k is a CartesianIndices 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.403491   0.0       0.0
 0.0367705  0.167762  0.0
 0.48649    0.452925  0.931146

julia> L * L
ERROR: MethodError: no method matching *(::LowerTriangularMatrix{Float32, Vector{Float32}, Spectrum{CPU{KernelAbstractions.CPU}, Vector{UnitRange{Int64}}, Vector{Int64}, Int64}, SpeedyWeatherInternals.ArrayDimensions.LM}, ::LowerTriangularMatrix{Float32, Vector{Float32}, Spectrum{CPU{KernelAbstractions.CPU}, Vector{UnitRange{Int64}}, Vector{Int64}, Int64}, SpeedyWeatherInternals.ArrayDimensions.LM})
The function `*` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  *(::Any, ::Any, ::Any, ::Oceananigans.Grids.AbstractGrid, ::Any, ::Any, ::Oceananigans.AbstractOperations.BinaryOperation, ::Oceananigans.AbstractOperations.BinaryOperation)
   @ Oceananigans ~/.julia/packages/Oceananigans/p6AAO/src/AbstractOperations/binary_operations.jl:75
  *(::Any, ::Any, ::Any, ::Oceananigans.Grids.AbstractGrid, ::Any, ::Any, ::Oceananigans.AbstractOperations.BinaryOperation, ::Oceananigans.Fields.AbstractField)
   @ Oceananigans ~/.julia/packages/Oceananigans/p6AAO/src/AbstractOperations/binary_operations.jl:81
  *(::Any, ::Any, ::Any, ::Oceananigans.Grids.AbstractGrid, ::Any, ::Any, ::Oceananigans.Fields.AbstractField, ::Oceananigans.AbstractOperations.BinaryOperation)
   @ Oceananigans ~/.julia/packages/Oceananigans/p6AAO/src/AbstractOperations/binary_operations.jl:86
  ...

julia> inv(L)
ERROR: MethodError: no method matching inv(::LowerTriangularMatrix{Float32, Vector{Float32}, Spectrum{CPU{KernelAbstractions.CPU}, Vector{UnitRange{Int64}}, Vector{Int64}, Int64}, SpeedyWeatherInternals.ArrayDimensions.LM})
The function `inv` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  inv(::CoordinateTransformations.PolarFromCartesian)
   @ CoordinateTransformations ~/.julia/packages/CoordinateTransformations/q6Edc/src/coordinatesystems.jl:78
  inv(::Rotations.CayleyMap)
   @ Rotations ~/.julia/packages/Rotations/wHQDQ/src/error_maps.jl:60
  inv(::Rotations.IdentityMap)
   @ Rotations ~/.julia/packages/Rotations/wHQDQ/src/error_maps.jl:64
  ...

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

L' * L

1.5011475f0

Summation with sum follows the flat, single index logic

L = rand(LowerTriangularArray{Float32}, 3, 3, 5)
sum(L, dims=2)

6×1 Matrix{Float32}:
 2.3928268
 2.6844249
 1.4734949
 2.169911
 1.743578
 2.7211347

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 front 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. In contrast to the grid-space rotate! (see Rotate and reverse Fields) which is restricted to multiples of 90˚, any rotation angle is possible in spectral space. For example start with

M = rand(LowerTriangularMatrix{ComplexF32}, 3, 3)

6-element, 3x3 LowerTriangularMatrix{ComplexF32, Array{...}}
 0.00618243+0.811319im       0.0+0.0im            0.0+0.0im
   0.866094+0.414133im  0.559887+0.557914im       0.0+0.0im
   0.803086+0.354882im  0.203311+0.508718im  0.825409+0.0491611im

Now rotate!(::LowerTriangularArray, degree)

rotate!(M, 45)

6-element, 3x3 LowerTriangularMatrix{ComplexF32, Array{...}}
 0.00618243+0.811319im       0.0+0.0im              0.0+0.0im
   0.866094+0.414133im  0.790404-0.0013946im        0.0+0.0im
   0.803086+0.354882im   0.50348+0.215955im   0.0491611-0.825409im

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

rotate!(M, 315)

6-element, 3x3 LowerTriangularMatrix{ComplexF32, Array{...}}
 0.00618243+0.811319im       0.0+0.0im            0.0+0.0im
   0.866094+0.414133im  0.559887+0.557914im       0.0+0.0im
   0.803086+0.354882im  0.203311+0.508718im  0.825409+0.0491611im

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 LowerTriangularArray 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. Both are consistent with reversing the corresponding field in grid space, see Rotate and reverse Fields.

reverse(M, dims=:lat)

6-element, 3x3 LowerTriangularMatrix{ComplexF32, Array{...}}
 0.00618243+0.811319im        0.0+0.0im            0.0+0.0im
  -0.866094-0.414133im   0.559887+0.557914im       0.0+0.0im
   0.803086+0.354882im  -0.203311-0.508718im  0.825409+0.0491611im

and in longitude

reverse(M, dims=:lon)

6-element, 3x3 LowerTriangularMatrix{ComplexF32, Array{...}}
 0.00618243-0.811319im       0.0+0.0im            0.0+0.0im
   0.866094-0.414133im  0.559887-0.557914im       0.0+0.0im
   0.803086-0.354882im  0.203311-0.508718im  0.825409-0.0491611im

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

L + L
2L

@. L + 2L - 1.1*L / L^2

6×5 (LM+) LowerTriangularArray{Float64, 2, Array{...}}
 -0.5655048710577639     0.7356071458526432  …   1.1460142345637232
  1.088774065928067      1.562425433713144      -0.11582119935829471
  1.5642812275241587   -18.2369814333732        -7.437883956271128
  0.07035814588005529  -30.144572318270193       0.691247146855072
 -4.434333998631767     -1.9750274948984377      0.2920995888056983
 -4.589838770172374     -3.6657477965690575  …   1.6815635257698263

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)

Array dimensions for LowerTriangularArray

Like Field (see Array dimensions of a Field in RingGrids), a LowerTriangularArray carries a dims::AbstractArrayDimensions field that records what the dimensions beyond the spherical harmonics (l, m) represent, see Array dimensions for a general overview of these dimension tags. The only difference is the name of the 2D (horizontal) dimension: LM instead of XY, since a LowerTriangularArray stores spherical harmonic coefficients rather than grid-point values. Correspondingly LMZ, LMT, and LMZT add a vertical and/or time dimension, exactly as XYZ, XYT, XYZT do for Field. LM is the default if no dims is given, e.g. for the L, L2 created above. Everything else – ArrayDimensions.hasvertical, ArrayDimensions.hastime, preservation through similar/zero/views/indexing – works the same way as described for Field.

L3 = zeros(LowerTriangularArray{Float32}, 5, 5, ArrayDimensions.LMZ(), 3)
ArrayDimensions.hasvertical(L3), ArrayDimensions.hastime(L3)

(true, false)

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:

spectrum = Spectrum(5, 5) # initailize

T4 Spectrum{...}
├ lmax = 5 (degrees)
├ mmax = 5 (orders)
└ architecture = CPU{KernelAbstractions.CPU}

L = rand(Float32, spectrum)

15-element, 5x5 LowerTriangularMatrix{Float32, Array{...}}
 0.336988   0.0       0.0        0.0       0.0
 0.647808   0.424832  0.0        0.0       0.0
 0.0796916  0.718515  0.879735   0.0       0.0
 0.606042   0.426354  0.0162833  0.109667  0.0
 0.872909   0.593394  0.967964   0.681627  0.533552

L = rand(ComplexF32, spectrum, 5)

15×5 (LM+) LowerTriangularArray{ComplexF32, 2, Array{...}}
  0.35410845f0 + 0.3869573f0im    …  0.10251653f0 + 0.97568387f0im
 0.041801453f0 + 0.8886602f0im       0.28862238f0 + 0.898713f0im
   0.4057095f0 + 0.8792975f0im       0.65794617f0 + 0.9296489f0im
  0.17064643f0 + 0.19452655f0im      0.37579894f0 + 0.8346743f0im
  0.46631438f0 + 0.6656572f0im        0.2517563f0 + 0.5091994f0im
 0.054786444f0 + 0.35033792f0im   …  0.47839397f0 + 0.19315451f0im
    0.535221f0 + 0.92554766f0im       0.9206643f0 + 0.41539896f0im
   0.3432979f0 + 0.2247439f0im         0.991665f0 + 0.13374346f0im
   0.6956607f0 + 0.8232098f0im       0.43988764f0 + 0.24385548f0im
   0.9882784f0 + 0.718436f0im        0.63473815f0 + 0.46474677f0im
   0.5945986f0 + 0.064747274f0im  …   0.0878042f0 + 0.63793325f0im
 0.122810066f0 + 0.54581964f0im      0.42790103f0 + 0.44703364f0im
   0.6199145f0 + 0.0498994f0im        0.7005154f0 + 0.4540599f0im
  0.38975435f0 + 0.01271379f0im        0.554838f0 + 0.99837154f0im
   0.3097093f0 + 0.66531825f0im      0.42191964f0 + 0.40844876f0im

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:

using SpeedyWeather    # SpectralGrid is not defined in LowerTriangularArrays
SG = SpectralGrid(trunc=5)
L = rand(Float32, SG.spectrum)

27-element, 7x6 LowerTriangularMatrix{Float32, Array{...}}
 0.141993  0.0        0.0        0.0       0.0        0.0
 0.698794  0.913593   0.0        0.0       0.0        0.0
 0.39406   0.0779604  0.887574   0.0       0.0        0.0
 0.297644  0.396918   0.112285   0.484569  0.0        0.0
 0.547624  0.925827   0.0594003  0.396436  0.0265943  0.0
 0.213309  0.527069   0.73134    0.836655  0.870764   0.600561
 0.691914  0.538687   0.76141    0.211872  0.300015   0.942151

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.LowerTriangularArrayWithTime Type

Type alias for all LowerTriangularArrays with a time dimension

source

LowerTriangularArrays.LowerTriangularArrayWithTimeAndVertical Type

Type alias for all LowerTriangularArrays with both time and vertical dimensions

source

LowerTriangularArrays.LowerTriangularArrayWithVertical Type

Type alias for all LowerTriangularArrays with a vertical dimension

source

LowerTriangularArrays.LowerTriangularMatrix Type

2D LowerTriangularArray of type T

source

LowerTriangularArrays.LowerTriangularMatrix Method

LowerTriangularMatrix(
    M::Array{T, 2}
) -> LowerTriangularArray{_A, 1, Vector{_A}, Spectrum{CPU{KernelAbstractions.CPU}, Vector{UnitRange{Int64}}, Vector{Int64}, Int64}, SpeedyWeatherInternals.ArrayDimensions.LM} where _A

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
) -> Union{LowerTriangularArray{_A, _B, var"#s179", Spectrum{A, O, L, IntType}, SpeedyWeatherInternals.ArrayDimensions.LM} where {_B, var"#s179"<:AbstractArray{_A, _B}, A, O, L, IntType}, LowerTriangularArray{_A, N, var"#s179", Spectrum{A, O, L, IntType}, SpeedyWeatherInternals.ArrayDimensions.LM} where {N, var"#s179"<:AbstractArray{_A, N}, A, O, L, IntType}} where _A

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}, _A, _B, <:Integer} where {_A, _B}

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}, _A, _B, <:Integer} where {_A, _B}

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 (mirror at the equator) when the elements in L represent the coefficients of the spherical harmonics. Reversal in latitude direction is obtained by flipping the sign of the harmonics with odd degree + order l + m. Reverses L in place.

source

Base._reverse! Method

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

Reverse the field represented by L in longitude (mirror at the 0˚ meridian) when the elements in L represent the coefficients of the spherical harmonics. Reversal in longitude direction is obtained by taking the complex conjugate of the harmonics, consistent with the grid-space reverse!(field, dims=:lon). 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 eastward by degree by multiplication of the spherical harmonics by exp(-i_m_2π*degree/360), with m the order of the spherical harmonic. Any degree is allowed, in contrast to the grid-space rotate! which is restricted to multiples of 90˚.

source

LowerTriangularArrays.eachharmonic Method

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

Iterator over all spherical harmonics of the LowerTriangularArrays provided as arguments, yielding (l, m) tuples of degree l and order m (both 1-based). Checks first that all arrays match in the horizontal, other dimensions may differ. Only loops over the horizontal dimension; combine with eachmatrix for the others.

source

LowerTriangularArrays.eachharmonic Method

eachharmonic(L::LowerTriangularArray) -> Any

Iterator over all spherical harmonics in L, yielding (l, m) tuples of degree l and order m (both 1-based) for every harmonic in the lower triangle. Only loops over the horizontal dimension; combine with eachmatrix for the other dimensions.

source

LowerTriangularArrays.eachharmonic Method

eachharmonic(S::Spectrum) -> Any

Iterator over all spherical harmonics in S, yielding (l, m) tuples of degree l and order m (both 1-based) for every harmonic in the lower triangle. To be used like

for (l, m) in eachharmonic(S)
    L[l, m]
end

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

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