Library

Domains

ApproxFunFourier.CircleType
Circle(c,r,o)

represents the circle centred at c with radius r which is positively (o=true) or negatively (o=false) oriented.

ApproxFun.CurveType

Curve Represents a domain defined by the image of a Fun. Example usage would be

x=Fun(1..2)
Curve(exp(im*x))  # represents an arc
ApproxFunBase.SegmentType
Segment(a,b)

represents a line segment from a to b. In the case where a and b are real and a < b, then this is is equivalent to an Interval(a,b).

IntervalSets.IntervalType

An Interval{L,R}(left, right) where L,R are :open or :closed is an interval set containg x such that

  1. left ≤ x ≤ right if L == R == :closed
  2. left < x ≤ right if L == :open and R == :closed
  3. left ≤ x < right if L == :closed and R == :open, or
  4. left < x < right if L == R == :open
ApproxFunOrthogonalPolynomials.RayType
Ray{a}(c,L,o)

represents a scaled ray (with scale factor L) at angle a starting at c, with orientation out to infinity (o = true) or back from infinity (o = false).

DomainSets.∂Function

Return the boundary of the given domain as a domain.

Accessing information about a spaces

ApproxFunBase.canonicalspaceFunction
canonicalspace(s::Space)

Return a space that is used as a default to implement missing functionality, e.g., evaluation. Implement a Conversion operator or override coefficients to support this.

Examples

julia> ApproxFunBase.canonicalspace(NormalizedChebyshev())
Chebyshev()
ApproxFunBase.itransformFunction
itransform(s::Space,coefficients::AbstractVector)

Transform coefficients back to values. Defaults to using canonicalspace as in transform.

Examples

julia> F = Fun(x->x^2, Chebyshev())
Fun(Chebyshev(), [0.5, 0.0, 0.5])

julia> itransform(Chebyshev(), coefficients(F)) ≈ values(F)
true

julia> itransform(Chebyshev(), [0.5, 0, 0.5])
3-element Vector{Float64}:
 0.75
 0.0
 0.75
ApproxFunBase.transformFunction
transform(s::Space, vals)

Transform values on the grid specified by points(s,length(vals)) to coefficients in the space s. Defaults to coefficients(transform(canonicalspace(space),values),canonicalspace(space),space)

Examples

julia> F = Fun(x -> x^2, Chebyshev());

julia> coefficients(F)
3-element Vector{Float64}:
 0.5
 0.0
 0.5

julia> transform(Chebyshev(), values(F)) ≈ coefficients(F)
true

julia> v = map(F, points(Chebyshev(), 4)); # custom grid

julia> transform(Chebyshev(), v)
4-element Vector{Float64}:
 0.5
 0.0
 0.5
 0.0
ApproxFunBase.evaluateFunction
evaluate(coefficients::AbstractVector, sp::Space, x)

Evaluate the expansion at a point x that lies in domain(sp). If x is not in the domain, the returned value will depend on the space, and should not be relied upon. See extrapolate to evaluate a function at a value outside the domain.

DomainSets.dimensionMethod
dimension(s::Space)

Return the dimension of s, which is the maximum number of coefficients.

Inbuilt spaces

ApproxFunOrthogonalPolynomials.ChebyshevType

Chebyshev() is the space spanned by the Chebyshev polynomials

    T_0(x),T_1(x),T_2(x),…

where T_k(x) = cos(k*acos(x)). This is the default space as there exists a fast transform and general smooth functions on [-1,1] can be easily resolved.

ApproxFunOrthogonalPolynomials.JacobiType

Jacobi(b,a) represents the space spanned by Jacobi polynomials P_k^{(a,b)}, which are orthogonal with respect to the weight (1+x)^β*(1-x)^α

ApproxFunOrthogonalPolynomials.LaguerreType

Laguerre(α) is a space spanned by generalized Laguerre polynomials Lₙᵅ(x) 's on (0, Inf), which satisfy the differential equations

    xy'' + (α + 1 - x)y' + ny = 0

Laguerre() is equivalent to Laguerre(0) by default.

ApproxFunOrthogonalPolynomials.UltrasphericalType

Ultraspherical(λ) is the space spanned by the ultraspherical polynomials

    C_0^{(λ)}(x),C_1^{(λ)}(x),C_2^{(λ)}(x),…

Note that λ=1 this reduces to Chebyshev polynomials of the second kind: C_k^{(1)}(x) = U_k(x). For λ=1/2 this also reduces to Legendre polynomials: C_k^{(1/2)}(x) = P_k(x).

ApproxFunFourier.TaylorType

Taylor() is the space spanned by [1,z,z^2,...]. This is a type alias for Hardy{true}.

ApproxFunFourier.HardyType

Hardy{false}() is the space spanned by [1/z,1/z^2,...]. Hardy{true}() is the space spanned by [1,z,z^2,...].

ApproxFunFourier.FourierType
Fourier()

The space spanned by the trigonemtric polynomials

    1, sin(θ), cos(θ), sin(2θ), cos(2θ), …

See also Laurent.


Fourier(d::Domain)

The space spanned by the trigonemtric polynomials

    1, sin(2pi/L*θ), cos(2pi/L*θ), sin(2pi/L*2θ), cos(2pi/L*2θ), …

for a domain with a period L.

ApproxFunFourier.LaurentType
Laurent()

The space spanned by the complex exponentials

    1,exp(-im*θ),exp(im*θ),exp(-2im*θ),…

See also Fourier.


Laurent(d::Domain)

The space spanned by the complex exponentials

    1, exp(-im * (2pi/L*θ)), exp(im * (2pi/L*θ)), exp(-2im * (2pi/L*θ)), …

for a domain with a period L.

ApproxFunFourier.CosSpaceType
CosSpace()

The space spanned by [1, cos θ, cos 2θ, ...]


CosSpace(d::Domain)

The space spanned by [1,cos(2pi/L*θ), cos(2pi/L*2θ), ...] for a domain with a period L

ApproxFunFourier.SinSpaceType
SinSpace()

The space spanned by [sin θ, sin 2θ, ...]


SinSpace(d::Domain)

The space spanned by [1, sin(2pi/L*θ), sin(2pi/L*2θ), ...] for a domain with a period L

ApproxFunSingularities.JacobiWeightType
JacobiWeight(β,α,s::Space)

weights a space s by a Jacobi weight, which on -1..1 is (1+x)^β*(1-x)^α. For other domains, the weight is inferred by mapping to -1..1.

ApproxFunSingularities.LogWeightType
LogWeight(β,α,s::Space)

represents a function on -1..1 weighted by log((1+x)^β*(1-x)^α). For other domains, the weight is inferred by mapping to -1..1.

ApproxFunBase.ArraySpaceType
ArraySpace(s::Space,dims...)

is used to represent array-valued expansions in a space s. The coefficients are of each entry are interlaced.

For example,

f = Fun(x->[exp(x),sin(x)],-1..1)
space(f) == ArraySpace(Chebyshev(),2)
ApproxFunBase.TensorSpaceType
TensorSpace(a::Space,b::Space)

represents a tensor product of two 1D spaces a and b. The coefficients are interlaced in lexigraphical order.

For example, consider

Fourier()*Chebyshev()  # returns TensorSpace(Fourier(),Chebyshev())

This represents functions on [-π,π) x [-1,1], using the Fourier basis for the first argument and Chebyshev basis for the second argument, that is, φ_k(x)T_j(y), where

φ_0(x) = 1,
φ_1(x) = sin x,
φ_2(x) = cos x,
φ_3(x) = sin 2x,
φ_4(x) = cos 2x
…

By Choosing (k,j) appropriately, we obtain a single basis:

φ_0(x)T_0(y) (= 1),
φ_0(x)T_1(y) (= y),
φ_1(x)T_0(y) (= sin x),
φ_0(x)T_2(y), …

Constructing a Fun

ApproxFunBase.FunType
Fun(s::Space, coefficients::AbstractVector)

Return a Fun with the specified coefficients in the space s

Examples

julia> f = Fun(Fourier(), [1,1]);

julia> f(0.1) == 1 + sin(0.1)
true

julia> f = Fun(Chebyshev(), [1,1]);

julia> f(0.1) == 1 + 0.1
true
Fun()

Return Fun(identity, Chebyshev()), which represents the identity function in -1..1.

Examples

julia> f = Fun(Chebyshev())
Fun(Chebyshev(), [0.0, 1.0])

julia> f(0.1)
0.1
Fun(s::Space)

Return Fun(identity, s)

Examples

julia> x = Fun(Chebyshev())
Fun(Chebyshev(), [0.0, 1.0])

julia> x(0.1)
0.1
Fun(f, d::Domain)

Return Fun(f, Space(d)), that is, it uses the default space for the specified domain.

Examples

julia> f = Fun(x->x^2, 0..1);

julia> f(0.1) ≈ (0.1)^2
true
Fun(f, s::Space)

Return a Fun representing the function, number, or vector f in the space s. If f is vector-valued, it Return a vector-valued analogue of s.

Examples

julia> f = Fun(x->x^2, Chebyshev())
Fun(Chebyshev(), [0.5, 0.0, 0.5])

julia> f(0.1) == (0.1)^2
true
Fun(f)

Return Fun(f, space) by choosing an appropriate space for the function. For univariate functions, space is chosen to be Chebyshev(), whereas for multivariate functions, it is a tensor product of Chebyshev() spaces.

Examples

julia> f = Fun(x -> x^2)
Fun(Chebyshev(), [0.5, 0.0, 0.5])

julia> f(0.1) == (0.1)^2
true

julia> f = Fun((x,y) -> x + y);

julia> f(0.1, 0.2) ≈ 0.3
true
Base.onesMethod
ones(d::Space)

Return the Fun that represents the function one on the specified space.

Examples

julia> ones(Chebyshev())
Fun(Chebyshev(), [1.0])
Base.zerosMethod
zeros(d::Space)

Return the Fun that represents the function one on the specified space.

Examples

julia> zeros(Chebyshev())
Fun(Chebyshev(), [0.0])

Accessing information about a Fun

ApproxFunBase.domainFunction
domain(f::Fun)

Return the domain that f is defined on.

Examples

julia> f = Fun(x->x^2);

julia> domain(f) == ChebyshevInterval()
true

julia> f = Fun(x->x^2, 0..1);

julia> domain(f) == 0..1
true
ApproxFunBase.coefficientsFunction
coefficients(cfs::AbstractVector, fromspace::Space, tospace::Space) -> Vector

Convert coefficients in fromspace to coefficients in tospace

Examples

julia> f = Fun(x->(3x^2-1)/2);

julia> coefficients(f, Chebyshev(), Legendre()) ≈ [0,0,1]
true

julia> g = Fun(x->(3x^2-1)/2, Legendre());

julia> coefficients(f, Chebyshev(), Legendre()) ≈ coefficients(g)
true
coefficients(f::Fun, s::Space) -> Vector

Return the coefficients of f in the space s, which may not be the same as space(f).

Examples

julia> f = Fun(x->(3x^2-1)/2);

julia> coefficients(f, Legendre()) ≈ [0, 0, 1]
true
coefficients(f::Fun) -> Vector

Return the coefficients of f, corresponding to the space space(f).

Examples

julia> f = Fun(x->x^2)
Fun(Chebyshev(), [0.5, 0.0, 0.5])

julia> coefficients(f)
3-element Vector{Float64}:
 0.5
 0.0
 0.5
ApproxFunBase.extrapolateFunction
extrapolate(f::Fun,x)

Return an extrapolation of f from its domain to x.

Examples

julia> f = Fun(x->x^2)
Fun(Chebyshev(), [0.5, 0.0, 0.5])

julia> extrapolate(f, 2) # 2 lies outside the domain -1..1
4.0
ApproxFunBase.ncoefficientsFunction
ncoefficients(f::Fun) -> Integer

Return the number of coefficients of a fun

Examples

julia> f = Fun(x->x^2)
Fun(Chebyshev(), [0.5, 0.0, 0.5])

julia> ncoefficients(f)
3
ApproxFunBase.pointsFunction
points(s::Space,n::Integer)

Return a grid of approximately n points, for which a transform exists from values at the grid to coefficients in the space s.

Examples

julia> chebypts(n) = [cos((2i+1)pi/2n) for i in 0:n-1];

julia> points(Chebyshev(), 4) ≈ chebypts(4)
true
points(f::Fun)

Return a grid of points that f can be transformed into values and back.

Examples

julia> f = Fun(x->x^2);

julia> chebypts(n) = [cos((2i+1)pi/2n) for i in 0:n-1];

julia> points(f) ≈ chebypts(ncoefficients(f))
true
ApproxFunBase.spaceFunction
space(f::Fun)

Return the space of f.

Examples

julia> f = Fun(x->x^2)
Fun(Chebyshev(), [0.5, 0.0, 0.5])

julia> space(f)
Chebyshev()
Base.valuesMethod
values(iterator)

For an iterator or collection that has keys and values, return an iterator over the values. This function simply returns its argument by default, since the elements of a general iterator are normally considered its "values".

Examples

julia> d = Dict("a"=>1, "b"=>2);

julia> values(d)
ValueIterator for a Dict{String, Int64} with 2 entries. Values:
  2
  1

julia> values([2])
1-element Vector{Int64}:
 2
values(f::Fun)

Return f evaluated at points(f).

Examples

julia> f = Fun(x->x^2)
Fun(Chebyshev(), [0.5, 0.0, 0.5])

julia> values(f)
3-element Vector{Float64}:
 0.75
 0.0
 0.75

julia> map(x->x^2, points(f)) ≈ values(f)
true
Base.strideMethod
stride(f::Fun)

Return the stride of the coefficients, checked numerically

Modify a Fun

ApproxFunBase.setdomainFunction
setdomain(f::Fun, d::Domain)

Return f projected onto domain.

Note

The new function may differ from the original one, as the coefficients are left unchanged.

Examples

julia> f = Fun(x->x^2);

julia> domain(f) == ChebyshevInterval()
true

julia> g = setdomain(f, 0..1);

julia> domain(g) == 0..1
true

julia> coefficients(f) == coefficients(g)
true
Base.chopMethod
chop(f::Fun[, tol = 10eps()]) -> Fun

Reduce the number of coefficients by dropping the tail that is below the specified tolerance.

Examples

julia> f = Fun(Chebyshev(), [1,2,3,0,0,0])
Fun(Chebyshev(), [1, 2, 3, 0, 0, 0])

julia> chop(f)
Fun(Chebyshev(), [1, 2, 3])

Bivariate Fun

ApproxFunBase.LowRankFunType
LowRankFun(f, space::TensorSpace)

Return an approximation to a bivariate function in a low-rank form

\[f(x,y) = \sum_i \sigma_i \phi_i(x) \psi_i(y)\]

where $\sigma_i$ represent the highest singular values, and $\phi_i(x)$ and $\psi_i(y)$ are orthogonal bases. The summation is truncated after an acceptable tolerance is reached.

Examples

julia> f = (x,y) -> x^2 * y^3;

julia> L = LowRankFun(f, Chebyshev() ⊗ Chebyshev());

julia> L(0.1, 0.2) ≈ f(0.1, 0.2)
true
ApproxFunBase.ProductFunType
ProductFun(coeffs::AbstractMatrix{T}, sp::AbstractProductSpace; [tol=100eps(T)], [chopping=false]) where {T<:Number}

Represent a bivariate function f in terms of the coefficient matrix coeffs, where the coefficients are obtained using a bivariate transform of the function f in the basis sp.

Examples

julia> P = ProductFun([0 0; 0 1], Chebyshev() ⊗ Chebyshev()) # corresponds to (x,y) -> x*y
ProductFun on Chebyshev() ⊗ Chebyshev()

julia> P(0.1, 0.2) ≈ 0.1 * 0.2
true
ProductFun(M::AbstractVector{<:Fun{<:UnivariateSpace}}, sp::UnivariateSpace)

Represent a bivariate function f(x,y) in terms of the univariate coefficient functions from M. The function f may be reconstructed as

\[f\left(x,y\right)=\sum_{i}M_{i}\left(x\right)b_{i}\left(y\right),\]

where $b_{i}\left(y\right)$ represents the $i$-th basis function for the space sp.

Examples

julia> P = ProductFun([zeros(Chebyshev()), Fun(Chebyshev())], Chebyshev()); # corresponds to (x,y)->x*y

julia> P(0.1, 0.2) ≈ 0.1 * 0.2
true

Operators

BandedMatrices.bandwidthsMethod
bandwidths(op::Operator)

Return the bandwidth of op in the form (l,u), where l ≥ 0 represents the number of subdiagonals and u ≥ 0 represents the number of superdiagonals.

ApproxFunBase.domainspaceFunction
domainspace(op::Operator)

Return the domain space of op. That is, op*f will first convert f to a Fun in the space domainspace(op) before applying the operator.

ApproxFunBase.rangespaceFunction
rangespace(op::Operator)

Return the range space of op. That is, op*f will return a Fun in the space rangespace(op), provided f can be converted to a Fun in domainspace(op).

Base.getindexMethod
(op::Operator)[k,j]

Return the kth coefficient of op*Fun([zeros(j-1);1],domainspace(op)).

Base.getindexMethod
(op::Operator)[f::Fun]

Construct the operator op * Multiplication(f), that is, it multiplies on the right by f first. Note that op * f is different: it applies op to f.

Examples

julia> x = Fun()
Fun(Chebyshev(), [0.0, 1.0])

julia> D = Derivative()
ConcreteDerivative : ApproxFunBase.UnsetSpace() → ApproxFunBase.UnsetSpace()

julia> Dx = D[x] # construct the operator y -> d/dx * (x * y)
TimesOperator : ApproxFunBase.UnsetSpace() → ApproxFunBase.UnsetSpace()

julia> twox = Dx * x # Evaluate d/dx * (x * x)
Fun(Ultraspherical(1), [0.0, 1.0])

julia> twox(0.1) ≈ 2 * 0.1
true
Base.:\Method
\(A::Operator,b;tolerance=tol,maxlength=n)

solves a linear equation, usually differential equation, where A is an operator or array of operators and b is a Fun or array of funs. The result u will approximately satisfy A*u = b.

LinearAlgebra.qrMethod
qr(A::Operator)

returns a cached QR factorization of the Operator A. The result QR enables solving of linear equations: if u=QR, then u approximately satisfies A*u = b.

LazyArrays.cacheMethod
cache(op::Operator)

Caches the entries of an operator, to speed up multiplying a Fun by the operator.

Inbuilt operators

ApproxFunBase.ConversionType
Conversion(fromspace::Space, tospace::Space)

Represent a conversion operator between fromspace and tospace, when available.

See also PartialInverseOperator that might be able to represent the inverse, even if this isn't banded.

ApproxFunBase.DerivativeType
Derivative(sp::Space, k::Int)

Return the k-th derivative operator on the space sp.

Examples

julia> Derivative(Chebyshev(), 2) * Fun(x->x^4) ≈ Fun(x->12x^2)
true
Derivative(sp::Space, k::AbstractVector{Int})

Return a partial derivative operator on a multivariate space. For example,

Dx = Derivative(Chebyshev()^2,[1,0]) # ∂/∂x
Dy = Derivative(Chebyshev()^2,[0,1]) # ∂/∂y
Tip

Using a static vector as the second argument would help with type-stability.

Examples

julia> ∂y = Derivative(Chebyshev()^2, [0,1]);

julia> ∂y * Fun((x,y)->x^2 + y^2) ≈ Fun((x,y)->2y)
true
Derivative(sp::Space)

Return the first derivative operator, equivalent to Derivative(sp,1).

Examples

julia> Derivative(Chebyshev()) * Fun(x->x^2) ≈ Fun(x->2x)
true
Derivative(k)

Return the k-th derivative, acting on an unset space. Spaces will be inferred when applying or manipulating the operator. If k is an Int, this returns a derivative in an univariate space. If k is an AbstractVector{Int}, this returns a partial derivative in a multivariate space.

Examples

julia> Derivative(1) * Fun(x->x^2) ≈ Fun(x->2x)
true

julia> Derivative([0,1]) * Fun((x,y)->x^2+y^2) ≈ Fun((x,y)->2y)
true
Derivative()

Return the first derivative on an unset space. Spaces will be inferred when applying or manipulating the operator.

Examples

julia> Derivative() * Fun(x->x^2) ≈ Fun(x->2x)
true
ApproxFunBase.DirichletType

Dirichlet(sp,k) is the operator associated with restricting the k-th derivative on the boundary for the space sp.

Dirichlet(sp) is the operator associated with restricting the the boundary for the space sp.

Dirichlet() is the operator associated with restricting on the the boundary.

ApproxFunBase.EvaluationType

Evaluation(sp,x,k) is the functional associated with evaluating the k-th derivative at a point x for the space sp.

Evaluation(sp,x) is the functional associated with evaluating at a point x for the space sp.

Evaluation(x) is the functional associated with evaluating at a point x.

ApproxFunBase.IntegralType
Integral(sp::Space, k::Int)

Return the k-th integral operator on sp. There is no guarantee on normalization.

Integral(sp::Space)

Return the first integral operator, equivalent to Integral(sp,1).

Integral(k::Int)

Return the k-th integral operator, acting on an unset space. Spaces will be inferred when applying or manipulating the operator.

Intergral()

Return the first integral operator on an unset space. Spaces will be inferred when applying or manipulating the operator.

ApproxFunBase.LaplacianType
Laplacian(sp::Space)

Return the laplacian operator on space sp.

Laplacian()

Return the laplacian operator on an unset space. Spaces will be inferred when applying or manipulating the operator.

ApproxFunBase.MultiplicationType

Multiplication(f::Fun,sp::Space) is the operator representing multiplication by f on functions in the space sp.

Multiplication(f::Fun) is the operator representing multiplication by f on an unset space of functions. Spaces will be inferred when applying or manipulating the operator.

ApproxFunBase.NeumannFunction

Neumann(sp) is the operator associated with restricting the normal derivative on the boundary for the space sp. At the moment it is implemented as Dirichlet(sp,1).

Neumann() is the operator associated with restricting the normal derivative on the boundary.

ApproxFunBase.PartialInverseOperatorType
PartialInverseOperator(O::Operator, bandwidths = (0, Infinities.ℵ₀))

Return an approximate estimate for inv(O), such that PartialInverseOperator(O) * O is banded, and is approximately I up to a bandwidth that is one less than the sum of the bandwidths of O and PartialInverseOperator(O).

Note

Only upper-triangular operators are supported as of now.

Examples

julia> C = Conversion(Chebyshev(), Ultraspherical(1));

julia> P = PartialInverseOperator(C); # default bandwidth

julia> P * C
TimesOperator : Chebyshev() → Chebyshev()
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋯
  ⋅   1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
  ⋅    ⋅   1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
  ⋅    ⋅    ⋅   1.0  0.0  0.0  0.0  0.0  0.0  0.0  ⋱
  ⋅    ⋅    ⋅    ⋅   1.0  0.0  0.0  0.0  0.0  0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅   1.0  0.0  0.0  0.0  0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  0.0  0.0  0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  0.0  0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  ⋱
  ⋮    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱    ⋱   ⋱

julia> P = PartialInverseOperator(C, (0, 4)); # specify an upper bandwidth

julia> P * C
TimesOperator : Chebyshev() → Chebyshev()
 1.0  0.0  0.0  0.0  0.0  0.0  -0.5    ⋅     ⋅     ⋅   ⋅
  ⋅   1.0  0.0  0.0  0.0  0.0   0.0  -1.0    ⋅     ⋅   ⋅
  ⋅    ⋅   1.0  0.0  0.0  0.0   0.0   0.0  -1.0    ⋅   ⋅
  ⋅    ⋅    ⋅   1.0  0.0  0.0   0.0   0.0   0.0  -1.0  ⋅
  ⋅    ⋅    ⋅    ⋅   1.0  0.0   0.0   0.0   0.0   0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅   1.0   0.0   0.0   0.0   0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    1.0   0.0   0.0   0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅     ⋅    1.0   0.0   0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅     ⋅     ⋅    1.0   0.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅     ⋅     ⋅     ⋅    1.0  ⋱
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅     ⋅     ⋅     ⋅     ⋅   ⋱