Docstrings

AbstractLinearRecurrence{T} <: AbstractRecurrence{T}

A recurrence where the results are linear w.r.t. the initial conditions. To support pseudostablization, the following methods are needed:

• rdiv!(::AbstractLinearRecurrence, x::Number): divides the system by x. It is used to prevent floating point overflow.
• LinearAlgebra.norm(step!(::AbstractLinearRecurrence), Inf): computes the ∞-norm of a step. It is used to get the amplification of the system.

AbstractLinearRecurrencePlan <: AbstractRecurrencePlan

A recurrence plan where the results are linear w.r.t. the initial conditions. Pseudostablization algorithm on linear recurrences is implemented.

AbstractRecurrence{T}

The abstract type where a recurrence instance runs. It needs to support the following methods:

• step!(::AbstractRecurrence): step forward the recurrence. Returns nothing if the recurrence terminates, otherwise returns a view of the stepping result.

AbstractRecurrencePlan

The abstract type of recurrence plans. A recurrence plan should include all the information that can perform a recurrence with a precision. It needs to support the following methods:

• init(::AbstractRecurrencePlan; T=Float64, init=:default): generates an AbstractRecurrence object where the actual recurrence runs on. Returns a vector of the view of the initial steps as well.
  • T is the underlying type where the precision is implied.
  • init tells how the initial values are generated. Must support :default for the forward recurrence to run. It is also suggested to support :rand for the pseudostablization algorithm.

BandedSylvesterRecurrence{T, TA<:AbstractMatrix{T}, TB<:AbstractMatrix{T}, TC<:AbstractMatrix{T}, TX<:AbstractMatrix{T}} <: AbstractLinearRecurrence{slicetype(TX)}

The recurrence generated from the infinite Sylvester equation $AX+XB+C=0$, assuming $X$ has infinite number of columns. The upper bandwidth of $A$ has to be finite, the lower bandwidth of $B$ has to be positive and finite and the lower bounding band of $B$ can't contain zero. The recurrence is basically a cross-shaped stencil recurrence, where the width of the stencil is determined by The total bandwidth of $B$ and the height by that of $A$. See also BandedSylvesterRecurrencePlan.

Properties

• A::TA, B::TB, C::TC: the matrices $A$, $B$ and $C$. It's recommended to use BandedMatrices.jl to boost performance. Their dimensions don't have to match. Just make sure that no BoundsError happens during the recurrence.
• buffer::TX: the buffer that stores temp results.
• sliceind::Int: the current column index.
• lastind::Int: the last column index to be computed.

BandedSylvesterRecurrencePlan{FA<:Function, FB<:Function, FC<:Function, INIT<:Function} <: AbstractLinearRecurrencePlan

See also BandedSylvesterRecurrence.

Properties

• fA::FA, fB::FB, fC::FC: the functions that generate A, B and C for the BandedSylvesterRecurrence. The functions should have the form f(eltype, size...).
• init::INIT: the function that generates the first few columns of $X$ in order to start the recurrence. It should have the form f(eltype, size...).
• size::Dims{2}: the size of $X$.
• bandB::NTuple{2,Int}: the bandwidths of B. Although B can have infinite upper bandwidth, that will cause the stencil to have infinite width and hence in practice, the upper bandwidth of B is always limited by the finite width of $X$.
• firstind::Int: the first column to be computed. The default value is bandB[2]+1.

StencilRecurrence{N,T,S,
    COEF<:NTuple{S,AbstractArray{T,N}},
    TB<:AbstractArray{T,N},}
    (stencil, coef, buffer, slicestart, sliceend, lastind)

Parameters

• N: system dimension
• T: eltype
• S: stencil size

Properties

For coef and slicesupport, tt's suggested to use lazy arrays for performance.

• stencil::NTuple{S, CartesianIndex{N}}: The relative index of the stencil. Can contain (0,0) (see coef)
• coef::COEF<:NTuple{S,AbstractArray{T,N}}: The coefficient associated with each relative index. The one associated with CartesianIndex(0,0) refers to a constant added to that entry. It's suggested to use lazy arrays for performance.
• buffer::CircularArray{T,N,TB<:AbstractArray{T,N}}: a buffer to store temp results.
• slicestart::MVector{N, Int} and sliceend::MVector{N, Int}: marks the current range of entries to be determined. Technically NTuple{N-1, Int} should work, but Julia doesn't support computed type parameters.
• lastslice::Int: marks the index of the slice where the recurrence terminates.

StencilRecurrencePlan{N, D, S, COEF<:NTuple{S,Function}, INIT<:Function} <: AbstractLinearRecurrencePlan

Parameters

• N: system dimension
• D: number domain, e.g. Real, Complex, etc.
• S: stencil size

Properties

• stencil::SVector{S, CartesianIndex{N}}: The relative index of the stencil. Can contain (0,0) (see coef)
• coef::COEF<:NTuple{S,Function}: The coefficient associated with each relative index. The one associated with CartesianIndex(0,0) refers to a constant added to that entry. The functions should be in the form f(T, I...) where I is the index of the stencil and T is the suggested return type. Coefficients should be at least as accurate as T. Exact-value types such as Irrational, Rational or Integer would do the job, and if that's not possible, BigFloat would work as well.
• init::INIT<:Function: the function used for initial values. The functions should be in the form f(I...) or f(T, I...) where I is the size of the array, excluding the last dimension, and T is the precisiontype. For the former case, f should return exact values, i.e. Integer, Rational or Irrational.
• size::Dims{N}: the size of the whole array.
• offset::NTuple{N,Int}: the very first index where the recurrence starts at.

ToPrecision{F}(f::F) <: Function

Create a function where the first argument specifies the returned precisiontype: (f::ToPrecision)(T, args...) = convert_precisiontype(precisiontype(T), f.f(args...))

Examples

julia> f() = π
f (generic function with 1 method)

julia> f(n) = 1//n + (n-1)//n*im
f (generic function with 2 methods)

julia> g = ToPrecision(f)
(::PseudostableRecurrences.ToPrecision{typeof(f)}) (generic function with 1 method)

julia> g(Float64)
3.141592653589793

julia> g(BigFloat)
3.141592653589793238462643383279502884197169399375105820974944592307816406286198

julia> g(Float16, 7)
Float16(0.1428) + Float16(0.857)im

julia> g(Float32, 3)
0.33333334f0 + 0.6666667f0im

_dotu(x, y) = mapreduce(*, +, x, y)

Not to be confused with LinearAlgebra.BLAS.dotu. See https://github.com/JuliaLang/julia/pull/27677

method_to_precision(f, argtypes)

If hasmethod(f, argtypes), return ToPrecision(f). Otherwise return f.

Examples

julia> f() = π

julia> Tf = methodtoprecision(BigFloat, f, Tuple{})

precision_shift(P::AbstractLinearRecurrencePlan)

Estimates log2 of the amplification of P by performing a full recurrence based on random initial conditions.

slicetype(T)

Get the type of a slice of T

Examples

julia> using PseudostableRecurrences: slicetype

julia> slicetype(Matrix{Float64})
Vector{Float64} (alias for Array{Float64, 1})

julia> slicetype(UnitRange{Int})
Int64

julia> slicetype(BitArray{3})
BitMatrix (alias for BitArray{2})