barycentricweights(r)

Create an array analogous to r of the barycentric weights associated with the array of points r. See also spectralderivative and spectralinterpolation.

Examples

Create barycentric weights for equally spaced points

julia> Bennu.barycentricweights(collect(-1.0:1.0))
3-element Vector{Float64}:
  0.5
 -1.0
  0.5

or create barycentric weights for Legendre–Gauss–Lobatto points

julia> x, _ = legendregausslobatto(4);
julia> Bennu.barycentricweights(x)
4-element Vector{Float64}:
 -0.6250000000000002
  1.3975424859373686
 -1.3975424859373689
  0.6250000000000006

References

Jean-Paul Berrut & Lloyd N. Trefethen, "Barycentric Lagrange Interpolation", SIAM Review 46 (2004), pp. 501-517. https://doi.org/10.1137/S0036144502417715

batchedbandedlu!(A::AbstractArray{T, 4} [, kl = div(size(mat, 2) - 1, 2)])

Compute the banded LU factors of the batched matrix A. Uses the A as storage for the factors. The optional argument kl is the lower bandwidth of the matrix; upper bandwidth is calculated as ku = size(mat, 2) - 1 - kl.

The matrix A should have indexing A[i, b, v, h] where the indices i and h are the batch indices with the size of the first index defining the kernel workgroup size the last index the number of workgroups. The index b is the band index for column v; storage uses the LAPACK format.

For example the matrix A would be stored in banded storage B with lower bandwidth kl = 2.

julia> A = [11 12  0  0  0  0  0
            21 22 23  0  0  0  0
            31 32 33 34  0  0  0
            41 42 43 44 45  0  0
             0 52 53 54 55 56  0
             0  0 63 64 65 66 67
             0  0  0 74 75 76 77]
7×7 Matrix{Int64}:
 11  12   0   0   0   0   0
 21  22  23   0   0   0   0
 31  32  33  34   0   0   0
 41  42  43  44  45   0   0
  0  52  53  54  55  56   0
  0   0  63  64  65  66  67
  0   0   0  74  75  76  77

julia> B = [ 0 12 23 34 45 56 66 67
            11 22 33 44 55 65 77  0
            21 32 43 54 64 76  0  0
            31 42 53 63 75  0  0  0
            41 52  0 74  0  0  0  0]
5×8 Matrix{Int64}:
  0  12  23  34  45  56  66  67
 11  22  33  44  55  65  77   0
 21  32  43  54  64  76   0   0
 31  42  53  63  75   0   0   0
 41  52   0  74   0   0   0   0

batchedbandedmatrix!(matvec!, A::BatchedBandedMatrix, y, x[, nthreads=1024])

Update the banded matrix A using the matvec! function. See also batchedbandedmatrix

batchedbandedmatrix(matvec!, y::AbstractArray{T, 3}, x::AbstractArray{T, 3},
                    kl, ku[, nthreads=1024])

Return the banded matrix with lower and upper bandwidths kl and ku defined by the matrix-vector multiplication matvec! where matvec!(y, x) should set y = A * x where x and y are batched vectors of with size(x) == size(y) == [Nqh, n, Neh].

The optional argument nthreads defines the number of total threads to use per workgroup in the kernel launch.

x, w = gaussrule(lo, hi, a, b, endpoint, maxiterations=100)

Generates the points x and weights w for a Gauss rule with weight function w(x) on the interval lo < x < hi.

The arrays a and b hold the coefficients (as given, for instance, by legendrecoefficients) in the three-term recurrence relation for the monic orthogonal polynomials p(0,x), p(1,x), p(2,x), ... , that is,

p(k, x) = (x-a[k]) p(k-1, x) - b[k]² p(k-2, x),    k ≥ 1,

where p(0, x) = 1 and, by convention, p(-1, x) = 0 with

          hi
b[1]^2 = ∫  w(x) dx.
         lo

Thus, p(k, x) = xᵏ + lower degree terms and

 hi
∫  p(j, x) p(k, x) w(x) dx = 0 if j ≠ k.
lo

get_batched_array(x::StructArray, Nqh, Neh)

Get the array that backs fieldarray x with horiztonal workgroup size Nqh and horizonal number workgroups Neh in the format used by the batched routines.

The main functionality of this routine is to permuate the vertical dof and field dimensions for more efficient banded solvers.

hilbertcode(Y::AbstractArray{T}; bits=8sizeof(T)) where T

Given an array of axes coordinates Y stored as bits-bit integers the function returns the associated Hilbert integer H.

The encoding of the Hilbert integer is best described by example. If 5-bits are used from each of 3 coordinates then the function performs

 X[2]|                       H[0] = A B C D E
     | /X[1]       ------->  H[1] = F G H I J
axes |/                      H[2] = K L M N O
     0------ X[0]                   high low

where the 15-bit Hilbert integer = A B C D E F G H I J K L M N O is stored in H.

This function is based on public domain code from John Skilling which can be found in https://dx.doi.org/10.1063/1.1751381.

Examples

We can generate the two-dimensional Hilbert code for a 1-bit integers with

julia> hilbertcode([0,0], bits=1)
2-element Vector{Int64}:
 0
 0

julia> hilbertcode([0,1], bits=1)
2-element Vector{Int64}:
 0
 1

julia> hilbertcode([1,1], bits=1)
2-element Vector{Int64}:
 1
 0

julia> hilbertcode([1,0], bits=1)
2-element Vector{Int64}:
 1
 1

References

John Skilling, "Programming the Hilbert curve", AIP Conference Proceedings 707, 381 (2004). https://doi.org/10.1063/1.1751381

x, w = legendregauss(n, endpoint::EndPoint=Bennu.NeitherEndPoint())

Convenience function with type T = Float64:

x, w = legendregauss(T, n, endpoint::EndPoint=Bennu.NeitherEndPoint())

Returns points x and weights w for the n-point Gauss-Legendre rule for the interval -1 < x < 1 with weight function w(x) = 1.

Use endpoint=LeftEndPoint(), RightEndPoint() or BothEndPoints() for the left Radau, right Radau, or Lobatto rules, respectively.

x, w = legendregausslobatto(n)

Convenience function with type T = Float64:

x, w = legendregausslobatto(T, n)

Returns points x and weights w for the n-point Legendre-Gauss-Lobatto rule for the interval -1 ≤ x ≤ 1 with weight function w(x) = 1.

minmaxflip(x, y)

Returns x, y sorted lowest to highest and a bool that indicates if a swap was needed.

partition(r, p, P)

Partition r into P pieces and return the pth piece.

This will provide an equal partition when N = length(r) is divisible by P and otherwise the ranges will have lengths of either div(N, P) or cld(N, P).

Examples

We can get the second piece of 10:20 partitioned into three pieces with

julia> partition(10:20, 2, 3)
13:16

We can get the third piece of 10:2:20 partitioned into four pieces with

julia> partition(10:2:20, 3, 4)
16:2:16

partition(r, P)

Partition r into P pieces.

This will provide an equal partition when N = length(r) is divisible by P and otherwise the ranges will have lengths of either div(N, P) or cld(N, P).

Examples

We can partition 10:20 into three pieces with

julia> partition(10:20, 3)
3-element Vector{UnitRange{Int64}}:
 10:12
 13:16
 17:20

We can partition 10:2:20 into four pieces with

julia> partition(10:2:20, 4)
4-element Vector{StepRange{Int64, Int64}}:
 10:2:10
 12:2:14
 16:2:16
 18:2:20

quantize([T=UInt64], x)

Quantize a number x, between zero(x) and one(x), to an integer of type T between zero(T) and typemax(T). If x is an array or tuple each element is quantized.

Examples

julia> quantize(0.0)
0x0000000000000000

julia> quantize(UInt32, 0.5f0)
0x7fffffff

julia> quantize((1.0, 0.5))
(0xffffffffffffffff, 0x7fffffffffffffff)

spectralderivative(r, w=barycentricweights(r))

Create the differentiation matrix, type analogous to r, for polynomials defined on the points r with associated barycentric weights w. See also barycentricweights and spectralinterpolation.

Examples

Create the differentiation matrix for Legendre–Gauss–Lobatto points

julia> x, _ = legendregausslobatto(4);
julia> spectralderivative(x)
4×4 Matrix{Float64}:
 -3.0        4.04508      -1.54508       0.5
 -0.809017   7.77156e-16   1.11803      -0.309017
  0.309017  -1.11803      -1.55431e-15   0.809017
 -0.5        1.54508      -4.04508       3.0

References

Jean-Paul Berrut & Lloyd N. Trefethen, "Barycentric Lagrange Interpolation", SIAM Review 46 (2004), pp. 501-517. https://doi.org/10.1137/S0036144502417715

spectralinterpolation(r, s, w=barycentricweights(r))

Create the interpolation matrix, type analogous to r, for polynomials defined on the points r with associated barycentric weights w evaulated at the points s. See also barycentricweights and spectralderivative.

Examples

Create the interpolation matrix which evaluates polynomials on Legendre–Gauss–Lobatto points at an equally spaced grid.

julia> x, _ = legendregausslobatto(4);
julia> y = collect(-1:0.4:1);
julia> spectralinterpolation(x,y)
6×4 Matrix{Float64}:
  1.0           0.0           0.0           0.0
  0.16          0.936656     -0.136656      0.04
 -0.12          0.868328      0.331672     -0.08
 -0.08          0.331672      0.868328     -0.12
  0.04         -0.136656      0.936656      0.16
 -1.11022e-16   3.43078e-16  -8.98189e-16   1.0

References

Jean-Paul Berrut & Lloyd N. Trefethen, "Barycentric Lagrange Interpolation", SIAM Review 46 (2004), pp. 501-517. https://doi.org/10.1137/S0036144502417715

tuplesort(a)

Returns the tuple a as a sorted tuple.

tuplesortperm(a)

Returns the permutation to sort a.

tuplesortpermnum(a)

Returns the lexicographic permutation number p to sort the tuple a.