_last_index_below_degree(set::MultiindexSet, max_total_deg::Int) -> Int

Return the last column index i in the Grlex‑sorted set such that the total degree of set[i] is strictly less than max_total_deg. If no such index exists, return 0. O(1) via the precomputed degree_offsets field.

_last_index_below_degree(set::MultiindexSet, max_total_deg::Int,
						 allowed_indices::AbstractVector{Int}) -> Int

Same as the 2‑argument version, but restricted to the indices listed in allowed_indices (must be sorted in increasing order). Uses the O(1) degree_offsets boundary then searchsortedlast to locate the answer; no per-element degree sums are computed.

_lex_cmp_same_degree(v::SVector{N,Int}, e) -> Int

Three‑way descending‑lexicographic comparison of two exponents of equal total degree: -1 if v precedes e, +1 if e precedes v, 0 if they are equal. Since the degrees are known to match, no degree sums are computed.

_multinomial(e::Int, k::Vector{Int}) -> Int

Multinomial coefficient: e! / (k₁! k₂! … kₚ!) where sum(k) = e. Uses iterative multiplication of binomial coefficients to avoid overflow.

_total_degree(exp, ::Val{N}) -> Int

Total degree of an N-component exponent. Written as an explicit loop so that it works uniformly for Vector, SVector and NTuple, and returns 0 for N == 0, where sum on an empty tuple is not defined.

all_multiindices_in_box(bound::Vector{Int}) -> MultiindexSet

Generate all multi-indices v of length length(bound) such that 0 ≤ v[i] ≤ bound[i] for all i. The vectors are generated and then sorted according to graded lexicographic order.

all_multiindices_up_to(nvars::Int, max_degree::Int) -> MultiindexSet

Generate all exponent vectors with nvars variables whose total degree ≤ max_degree, sorted according to graded lexicographic order. Returns a MultiindexSet.

bounded_index_tuples(M, exp::SVector{N,Int})

Enumerate all index tuples of length M whose component counts are bounded by exp.

This function generates all tuples (i₁, ..., i_M) with entries in 1:N such that:

• Each index k appears at most exp[k] times

• The total length of the tuple is exactly M

Each tuple is represented in its canonical sorted form (non-decreasing order), so that it uniquely corresponds to a count vector (multiindex).

For each valid tuple, the function also returns:

• The associated multiindex alpha, where alpha[k] is the number of times k appears

• The number of distinct permutations of the tuple

Returns a vector of tuples: (indextuple, multiindex, permutationcount)

where:

• index_tuple::NTuple{M,Int} is the sorted tuple

• multiindex::SVector{N,Int} contains the counts

• permutation_count::Int is the number of distinct permutations

build_exponent_index_map(set::MultiindexSet) -> Dict{SVector{N,Int}, Int}

Build a dictionary mapping each exponent SVector to its 1-based column index in set. Useful for O(1) lookups without repeated binary searches.

delete_multiindices(set::MultiindexSet{N}, exps) -> MultiindexSet{N}

Return a new MultiindexSet holding every exponent of set except those listed in exps. set itself is never modified — there is no in-place variant.

exps may be a single exponent (SVector, Vector{Int} or NTuple{N,Int}), any iterable of such exponents, or another MultiindexSet{N}. Exponents that are not members of set are ignored, matching Base.setdiff; an exponent whose length differs from N throws an ArgumentError.

Deleting cannot reorder a Grlex-sorted list, so the result is built through the pre-sorted path: no re-sort is performed, and degree_offsets is rebuilt so the O(1) degree-boundary queries stay correct even when an entire degree block disappears.

Removing an exponent can break downward closure, which parametrise(model, spectral, mset) requires — check the result with is_downward_closed.

Examples

S = all_multiindices_up_to(2, 2)              # 6 monomials
T = delete_multiindices(S, [[2, 0], [0, 2]])  # drop the two pure squares
length(S), length(T)                          # (6, 4) — S is unchanged

delete_multiindices(pred, set::MultiindexSet{N}) -> MultiindexSet{N}

Return a new MultiindexSet holding the exponents of set for which pred(α) is false; set itself is never modified. pred receives each exponent as an SVector{N,Int}. This is the exact complement of filter(pred, set): together the two partition set.

Examples

S = all_multiindices_up_to(3, 4)
delete_multiindices(α -> sum(α) > 2, S) == all_multiindices_up_to(3, 2)   # true

divides(a::AbstractVector{Int}, b::AbstractVector{Int}) -> Bool

Check whether a divides b componentwise, i.e., a[i] ≤ b[i] for all i.

factorisations_asymmetric(set, exp, num_factors, candidate_indices) -> Vector{FactorisationEntry}

Return every ordered num_factors-tuple of indices from candidate_indices whose exponent vectors sum to exp. Each ordering is a separate entry with multiplier = 1.

factorisations_fully_symmetric(set, exp, num_factors, candidate_indices) -> Vector{FactorisationEntry}

Return all unordered (non-decreasing index) num_factors-tuples from candidate_indices whose exponent vectors sum to exp. multiplier is the multinomial coefficient num_factors! / (m₁! m₂! … mₖ!) counting the distinct ordered arrangements of each combination.

Duplicate entries in candidate_indices are removed internally.

factorisations_groupwise_symmetric(set, exp, group_sizes, candidate_indices) -> Vector{FactorisationEntry}

Return factorisations of exp into M = length(group_sizes) groups, where group i has group_sizes[i] factor slots and is internally symmetric.

factor_indices in each entry is the concatenation of per-group indices in group order (each group sorted non-decreasingly). multiplier is the product of the per-group permutation counts — the total number of ordered arrangements.

filter(pred, set::MultiindexSet{N}) -> MultiindexSet{N}

Return a new MultiindexSet holding the exponents of set for which pred(α) is true, leaving set untouched. pred receives each exponent as an SVector{N,Int}.

Filtering preserves Grlex order, so the result is built through the pre-sorted path. Use it to intersect independent conditions, e.g. a total-degree bound on the master coordinates together with a per-parameter box:

mset = filter(α -> α[1] + α[2] ≤ 4 && α[3] ≤ 2, all_multiindices_in_box([4, 4, 2]))

See delete_multiindices for the removal-side counterpart.

find_in_set(set::MultiindexSet, exp::AbstractVector{Int}) -> Union{Int, Nothing}

Return the column index of exp in set.exponents, or nothing if exp is not present (including when its length differs from the number of variables).

The total degree of exp brackets the search to a single degree block via degree_offsets, after which a binary search in Grlex (here plain descending lexicographic, the degrees being equal) locates the exponent.

find_in_set(set::MultiindexSet, exp::NTuple{N,Int}) where N -> Union{Int, Nothing}

Tuple version – avoids allocating a vector for the exponent.

grlex_precede(a::AbstractVector{<:Integer}, b::AbstractVector{<:Integer}) -> Bool

Graded lexicographic order: compare total degree first, then lexicographic descending. Returns true if a comes before b in this order.

indices_in_box_with_bounded_degree(set::MultiindexSet, box_upper::AbstractVector{Int},
								   degree_lower_bound::Int, total_deg_upper::Int) -> Vector{Int}

Return the column indices of all multiindices v in set such that

• v ≤ box_upper componentwise,

• degree_lower_bound ≤ sum(v[i]) < total_deg_upper.

Uses the degree bounds to limit the search to relevant columns.

indices_in_box_with_bounded_degree(set::MultiindexSet, box_upper::AbstractVector{Int},
								   degree_lower_bound::Int, total_deg_upper::Int,
								   allowed_indices::AbstractVector{Int}) -> Vector{Int}

Same as the 4‑argument version, but the search is restricted to the indices listed in allowed_indices (which must be sorted). Only those indices that also satisfy the componentwise bound are returned.

is_conjugate_closed(set::MultiindexSet{N}, perm::AbstractVector{Int}) -> Bool

Return true when the conjugate partner of every member of set is itself a member.

perm is an N-length permutation pairing conjugate coordinates (self-paired entries for real ones). It acts on the components of an exponent,

(P·α)[k] = α[perm[k]],

the same convention the cohomological solve uses when it fills a conjugate monomial from its partner. Monomials with P·α == α are self-paired and trivially closed.

parametrise(model, spectral, mset; conjugate_permutation = ...) requires this: a member whose partner is absent is solved directly instead of being filled by conjugation, which costs the pairing optimisation and returns a parametrisation without the conjugate structure the permutation asserts. It is enforced through validate_multiindex_set, which reports which partner is missing; use this predicate when a plain Bool is enough.

Throws an ArgumentError if length(perm) != N.

is_constant(exp::AbstractVector{Int}) -> Bool

Return true if the exponent vector is all zeros.

is_downward_closed(set::MultiindexSet{N}) -> Bool

Return true when every divisor of every member of set is itself a member.

This is the downward closure property that parametrise(model, spectral, mset) requires: the graded solve reads lower-order coefficients W[α - eᵢ] while working on α, so a missing divisor would be silently read as zero. Combinatorial truncations (all_multiindices_up_to, all_multiindices_in_box) are closed by construction, but a spectral criterion such as |⟨λ, α⟩| ≤ R generally is not.

parametrise and solve_parametrisation enforce this through validate_multiindex_set, which reports which divisor is missing; use this predicate when a plain Bool is enough.

The zero multiindex is exempt: DPIM sets are built with min_degree = 1, so a degree-1 exponent is not asked to have its constant divisor present.

Only the immediate predecessors α - eᵢ are probed; by induction over the total degree that is equivalent to checking every divisor. Cost is O(N · |set| · log) through find_in_set.

monomial_rank(exp::AbstractVector{Int}, nvars::Int, max_degree::Int) -> Int

Return the 1‑based rank of the exponent vector exp in the complete set of all multiindices of length nvars with total degree ≤ max_degree, sorted according to graded lexicographic order (Grlex).

The vector exp must satisfy length(exp) == nvars and sum(exp) ≤ max_degree. Uses combinatorial counting formulas for efficiency (O(nvars) time).

If the computed rank exceeds typemax(Int), an error is thrown.

multiindex(components::Int...) -> Vector{Int}

Convenience constructor for an exponent vector.

multiindices_with_total_degree(nvars::Int, deg::Int) -> MultiindexSet

Generate all exponent vectors of length nvars with total degree exactly deg, sorted in lexicographic order (the tie‑breaker for Grlex). Returns a MultiindexSet.

num_multiindices_up_to(nvars::Int, max_degree::Int) -> Integer

Return the number of exponent vectors of length nvars with total degree ≤ max_degree. For nvars = 0 the set contains one element if max_degree ≥ 0, otherwise zero.

zero_multiindex(n::Int) -> Vector{Int}

Return the zero exponent vector of length n (all components zero).

_monomial_vector!(m, poly, vals)

Fill the pre-allocated vector m (length L) with the value of each monomial in poly.multiindex_set evaluated at vals. This is the central loop; called once per evaluate.

_precompute_powers(T, vals, max_exps) -> NTuple of Vectors

For each variable j, store vals[j]^e for e = 0 … max_exps[j]. Avoids repeated ^ calls inside the inner loop.

coeff_shape(p) -> NTuple

Leading axes of the coefficient array. () for scalar, (K,) for K-vector, (m,n) for a matrix-valued polynomial, etc.

coefficient(p, exp) -> scalar or view

Return the coefficient of the monomial with exponent exp. For scalar polynomials this is a T; for vector-valued it is a Vector{T} view into the backing array — no allocation.

compose_linear(poly::DensePolynomial, M::Matrix{TA}) where TA -> DensePolynomial

Compose a multivariate polynomial with a linear map.

Arguments

• poly: polynomial in variables x₁, …, x_n. (The coefficient type can be numeric or array‑valued.)

• M: an n × p matrix. Composition means replacing x_i by ∑_{j=1}^p M[i,j] * y_j, where y₁, …, y_p are new variables.

Returns

A new polynomial in the variables y₁, …, y_p. The returned polynomial has the same coefficient type as the input poly.

each_term(poly) -> generator of (exponent, coefficient)

Yields (SVector{NVAR,Int}, coeff) for every nonzero monomial.

• Scalar (N=1): coeff is a T.

• Vector (N=2): coeff is a SubArray{T,1} view (no allocation).

evaluate(poly::DensePolynomial{T,NVAR,1}, vals) -> T

Scalar polynomial evaluation: sum(c_i * m_i) without conjugation.

evaluate(poly::DensePolynomial{T,NVAR,2}, vals) -> Vector{T}

Vector-valued polynomial evaluation via BLAS gemv:

result = coefficients  (K × L)  *  m  (L,)  →  (K,)

A single BLAS call, no per-term allocation.

evaluate(poly::DensePolynomial{T,NVAR,N}, vals) -> Array

General tensor-valued case: reshape to (prod(leading_dims), L), gemv, reshape back.

evaluate(poly::DensePolynomial{T,NVAR,2}, vals, component::Int) -> T

Evaluate a single component of a vector-valued polynomial without allocating a full output vector.

extract_component(poly::DensePolynomial{T,NVAR,2}, idx) -> DensePolynomial{T,NVAR,1}

Return the idx-th component as a scalar polynomial. The coefficient vector is a copy of row idx of the matrix.

linear_matrix_of_polynomial(poly::DensePolynomial{T,NVAR,2}) -> Matrix{T}

Return the (K × NVAR) matrix A such that the linear part equals A * x. Reads directly from the coefficient matrix — no intermediate arrays.

mmap_polynomial(path, coeff_size::NTuple, T::Type, mset; write=false)

Return a DensePolynomial whose coefficient array is backed by a memory-mapped file at path. coeff_size is the leading axes of the coefficient array (empty () for scalar, (K,) for K-vector, …).

The OS page-cache handles RAM pressure automatically: only accessed pages are loaded, so polynomials larger than available RAM work transparently.

Example

# 5-vector polynomial with 1 000 000 monomials, stored on disk
p = mmap_polynomial("coeffs.bin", (5,), Float64, mset; write=true)

restrict_polynomial_to_degree(poly::DensePolynomial, max_degree::Int) -> DensePolynomial

Returns a polynom that contains all the monomials of poly that are of degree lower or equal than max_degree.

similar_poly(dict::Dict{SVector{NVAR,Int}, T}) where {NVAR, T<:Number}

Scalar polynomial from a dictionary mapping exponent vectors to scalar coefficients.

similar_poly(dict::Dict{SVector{NVAR,Int}, SVector{K,T}}) where {NVAR, K, T<:Number}

Fixed-size-vector polynomial from a dictionary mapping exponent vectors to SVector{K,T} coefficients. The output coefficient matrix has size K × L.

similar_poly(dict::Dict{SVector{NVAR,Int}, Vector{T}}) where {NVAR, T<:Number}

Vector-valued polynomial from a dictionary mapping exponents to Vector{T} coefficients. All vectors must have the same length K.

zero(DensePolynomial{T}, mset)              # scalar
zero(DensePolynomial{T}, coeff_shape, mset) # tensor-valued

_arity_description(fixed, va) -> String

Human-readable summary of the argument counts a callable accepts, e.g. "3, ≥ 1".

_build_multilinear_map(f!, multiindex, multiplicity_external, fully_asymmetric)

Validate a fully resolved term and build it. The single place where the invariants of MultilinearMap are enforced; every public constructor ends up here.

Checks run in the order negatives → degree → arity so that sum(multiindex) is meaningful in every message.

_call_signature(multiindex, multiplicity_external) -> String

Render the call f! receives, e.g. "f!(res, x^(0), x^(0), x^(1), r)" for multiindex = (2, 1) and multiplicity_external = 1.

_check_arity(f!, deg, multiindex, multiplicity_external)

Throw an ArgumentError unless f! can be called with deg + 1 arguments.

hasmethod is the fast accept path, but it returns false for methods with concrete argument annotations, so a scan of the method table decides rejection. A varargs method that can absorb the arguments is accepted — that is what admits closures built by MORFESymbolicsExt and callable structs. A callable exposing no methods at all is trusted rather than rejected.

_counts_from_derivatives(derivatives) -> NTuple

Convert a per-slot list of 0-based derivative orders into a multiindex of counts: (0, 0, 1) → (2, 1). The list must be non-decreasing, since it is the argument order f! will be called with.

_definition_site(f!) -> String

Return " @ file.jl:12" for the first method of f!, or "" when f! exposes no methods. Used by show and by FullOrderModel._term_label to point a diagnostic at the definition of a term. Deliberately uses first rather than only: f! is allowed to carry several methods, and a display helper must never throw.

_infer_degree(f!) -> Int

Number of factors f! takes, i.e. its argument count less res. Requires a single fixed arity; varargs or conflicting arities are ambiguous and raise an ArgumentError telling the caller to state the degree explicitly.

_info_assumed(f!, multiindex, multiplicity_external, mi_source, assumed_order,
              assumed_me, forcing)

Emit an @info naming every value the constructor had to default, or nothing at all when the caller pinned them. mi_source is :stated, :arity or :degree, and names where an assumed multiindex came from.

fully_asymmetric is deliberately not reported here: it already has a dedicated diagnostic, FullOrderModel._info_implicit_symmetry, which fires at NthOrderModel construction exactly when the flag changes the result.

_method_arities(f!) -> (fixed, va)

Split the methods of f! by argument count, including res but excluding the callable itself. fixed lists the exact counts of the fixed-arity methods; va lists, for each varargs method, the minimum number of arguments it accepts.

_msg_assumed_compact(f!, multiindex, multiplicity_external, assumed) -> String

One-line report of the values the constructor had to default. assumed is the list of field names that were not stated; only those are named, so the closing sentence ("state them explicitly to silence this message") is always true.

_msg_assumed_forcing(f!, multiindex, multiplicity_external) -> String

Full report for the one inference that can silently reinterpret a linear term as a forcing term. Every call it presents as silent really is silent under the resolution rules — the testset "suggested silent forms are silent" pins that.

_pad_to_order(base, order) -> NTuple

Zero-pad base to length order. Padding is symmetry-neutral: trailing zeros change neither all(<=(1), mi) nor count(>(0), mi), so symmetry_type classifies (2,) and (2, 0, 0) identically. Truncation is refused — it would silently drop a slot.

_symmetry_label(multiindex, fully_asymmetric) -> String

Name the symmetry class a term will be given by the parametrisation solver.

Mirrors symmetry_type in src/ParametrisationMethod/RightHandSide/MultilinearTerms/Symmetry.jl, which cannot be called from here: MultilinearMaps is loaded twelve includes earlier. The testset "show agrees with symmetry_type" holds the two in step.

accumulate_qp!(Fe, ∇W_args, mult, element, q, dΩ, t::FEMMultilinearMap) -> nothing

Accumulate the integrand contribution at quadrature point q (weight dΩ) into the element residual vector Fe. ∇W_args is an NTuple of pre-scattered W columns; mult is the combinatorial multiplicity of the term.

assemble_element!(accum, Fe, element, t::FEMMultilinearMap) -> nothing

Scatter the element residual Fe into the global accumulator accum using the DOF map of element.

evaluate_term!(res, term, xs, r)

Evaluate a single MultilinearMap and accumulate (adds) the result into res.

Arguments

• res: output vector (modified in-place)

• term: multilinear term

• xs: tuple (x, x^(1), …, x^(ORD-1)) of state derivatives

• r: external state vector (or nothing if not used). If r is nothing but the term expects external arguments, an error is thrown.

evaluate_term!(res, t::FEMMultilinearMap{ORD}, xs, r)

Direct (uncached) evaluation of a FEM-backed multilinear term at state xs and external state r, accumulating the result into res. To be used in InvarianceError.jl.

Internal argument slots (determined by t.multiindex) are scattered to quadrature-point field quantities via scatter_qp! and assembled element-wise.

For me = 0: ∇W_args is a homogeneous NTuple{N_INT, QP_TYPE} — type-stable.

For me > 0: the external arg slots in ∇W_args receive r directly (not scattered). r is a small NEXT-dimensional external-state vector, not a FOM displacement field. The concrete `accumulateqp!evaluates the full multilinear mapF(∇u₁,…,∇uₙ, r₁,…,rₘₑ)` at the actual inputs.

fem_elements(t::FEMMultilinearMap) -> iterator

Return an iterator over the mesh elements (cells) for the FEM term t. Implement this for every concrete FEMMultilinearMap subtype.

fem_getdetJdV(element, q, t::FEMMultilinearMap) -> Real

Return the integration weight det(J) · w_q at quadrature point index q of element.

fem_n_qp(t::FEMMultilinearMap) -> Int

Return the number of quadrature points per element for the FEM term t.

fem_ndofs_per_cell(t::FEMMultilinearMap) -> Int

Return the number of degrees of freedom per element (cell) for the FEM term t.

fem_qp_buffer(t::FEMMultilinearMap)

Return a pre-allocated scratch buffer sized for one quadrature point of the FEM term t. Reused across calls to avoid allocation in the inner element loop.

fem_reinit!(element, t::FEMMultilinearMap) -> nothing

Reinitialise the FEM quadrature cache (e.g. CellValues) for element. Called once per element before any scatter_qp! or accumulate_qp! calls.

scatter_qp!(∇W_col, W_global, element, t::FEMMultilinearMap) -> nothing

Fill ∇W_col in-place with the field values (e.g. gradients) at all quadrature points of element for one column of the parametrisation W_global.

_evtype(::Type{T}) -> Type

Return the eigenvalue storage type for scalar type T: Complex{T} when T <: Real, or T itself when T <: Complex.

_rebase(poly, Q, Qinv) -> DensePolynomial

Re-express ṙ = f(r) in the coordinates r′ with r = Q r′, giving ṙ′ = Q⁻¹ f(Q r′).

The substitution itself is compose_linear (in Polynomials); the left multiplication by Q⁻¹ is a single matrix product on the coefficient array, since the polynomial is vector-valued with one row per component.

The promotion on the first line is load-bearing: compose_linear types its accumulators from the input polynomial, so a real polynomial composed with a complex Q cannot store the products back. promote_type rather than a blanket complex keeps a real A with an all-real spectrum in real arithmetic.

_subdiagonal_offenders(linear_matrix) -> Vector{Tuple{Int, Int}}

The (i, j) positions of the non-zero entries strictly below the diagonal.

Monomials are solved in GrLex order, so the |β| = 1 lower-order coupling can only read Λ[i, j] for i < j (LowerOrderCouplings._sum_degree_one_terms!); a strictly-lower entry of the external linear matrix would be silently discarded by the solver rather than producing a wrong-but-visible answer. Such entries are what triggers a re-basing, and this list is what the diagnostic reports. See the ExternalSystems module docstring.

_triangularising_basis(A) -> Union{Nothing, NamedTuple}

Choose a basis Q in which Q⁻¹ A Q is upper triangular, or nothing when A already is.

Returns (; Q, U, Qinv, route) with route ∈ (:eigen, :schur).

Two requirements, and why the route matters

Triangularity is the GrLex-causality constraint (module docstring); both routes give it.

Conjugate structure is the second, and only one route gives it. Realification.realify applies a single conj_map across all NVAR variables, external ones included, and fill_conjugate_monomial! implements W_{P·γ} = conj(W_γ) for a real FOM. Preserving that needs an involution σ with Q[:, σ(k)] = conj(Q[:, k]) and λ_{σ(k)} = conj(λ_k); a basis that breaks it silently invalidates any conjugate_permutation touching external indices — a wrong answer, not an error.

A Schur basis breaks it in general: in A = Q U Qᴴ the first Schur vector Q[:, 1] is an eigenvector for λ₁, but Q[:, 2] is a Schur vector, not an eigenvector for λ₂, so Q[:, 2] ≠ conj(Q[:, 1]) unless A is normal.

The eigenvector basis of a real A gives both. LAPACK's dgeev returns complex eigenvalues in adjacent conjugate pairs, each pair's eigenvectors stored as the real and imaginary parts of one column pair, so V[:, k+1] == conj(V[:, k]) holds bit-exactly; real eigenvalues get real, self-conjugate eigenvectors. And U = Diagonal(λ) is diagonal, which is strictly better than triangular: it lets the solver take the uncoupled fast path (coupled_external = !isdiag(...) in CohomologicalDriver).

So: real and diagonalisable → eigen; anything else (complex A, or real but near-defective, where the eigenvector basis would be catastrophically ill-conditioned) → complex schur. Note schur on a real matrix returns the real Schur form, which is only quasi-triangular — it must be given a complex matrix.

_unit(::Val{N}, j) -> SVector{N, Int}

The j-th unit multiindex / basis vector in N variables.

_zero_subdiagonal_linear_block!(poly; tol)

Zero the strictly-lower triangle of poly's linear block, after checking it is only round-off.

compose_linear produces sub-diagonal entries of size ~1e-17·‖A‖ where exact arithmetic would give zero, and istriu tests exact zeros — so without this the "triangularised" system is still not triangular by the package's own predicate, and linear_matrix would disagree with the polynomial it is supposed to be the Jacobian of.

linear_matrix_of_polynomial sets A[:, j] = coefficients[:, i_j] where i_j is the mset position of the unit multiindex eⱼ, so A[k, j] is coefficients[k, i_j] and the strict lower triangle is k > j within each unit-multiindex column.

A residual above tol means the basis is wrong; that must error rather than be truncated away, because truncating it would discard genuine coupling.

external_argument_vectors(sys, N_EXT) -> Vector{<:SVector{N_EXT}}

The external argument the solver passes for each external variable j.

During the cohomological solve the external factor of a monomial is a basis vector, not a state: the multilinear term is evaluated once per external direction and the coefficient is read off. In the model's own coordinates that vector is the unit vector eⱼ; after a re-basing it is the physical direction Q eⱼ = Q[:, j]. Multilinearity makes the substitution exact — f(…, Qeⱼ, Qe_k) = Σ_{a,b} Q[a,j] Q[b,k] f(…, e_a, e_b).

With basis === nothing this returns the integer unit vectors, i.e. exactly what the solver used before re-basing existed.

external_basis(sys) -> Union{Nothing, SMatrix}

The change of external coordinates Q with r = Q r′, or nothing when the system was left in the coordinates it was given in. nothing is the common case and is what every consumer branches on; it means "reduced external coordinates are the physical ones".

to_physical_external(sys, r) -> r_physical

Map a reduced external coordinate vector r′ to the physical r = Q r′.

The state-valued counterpart of external_argument_vectors, used where an external argument is a genuine state rather than a basis direction — notably the invariance-error residual, which samples reduced coordinates and must feed the model physical ones. Identity when the system was never re-based.

_check_external_terms(nonlinear_terms)

Reject terms that read the external state in a model built without one.

Without this the mismatch survives construction and only surfaces mid-solve, as "Term expects external arguments but no external state provided" from evaluate_term!.

evaluate_nonlinear_terms!(res, model, order, state_vectors, r = nothing)

Evaluate all nonlinear terms of a given polynomial degree for an NthOrderModel.

Arguments

• res: output vector (modified in-place)

• model: the NthOrderModel

• order: degree of the nonlinear terms to evaluate

• state_vectors: tuple (x, x^(1), …, x^(ORD-1)) of state derivatives

• r: external state vector (default nothing). Must be provided if any term uses external variables.

linear_first_order_matrices(model::NthOrderModel)

Construct the matrices A and B of the equivalent linear first-order system:

B Ẋ = A X

obtained from the ORD-th order model

B_ORD x^(ORD) + ... + B_1 x^(1) + B_0 x = F(...)

by introducing the augmented state vector

X = [x, x^(1), ..., x^(ORD-1)].

and the (ORDn_fom x ORDn_fom)-block matrices

B = [ I   0   0   ⋯   0
	  0   I   0   ⋯   0
	  ⋮       ⋱
	  0   0   0   ⋯  B_ORD ]

and

A = [ 0   I   0   ⋯   0
	  0   0   I   ⋯   0
	  ⋮       ⋱
	 -B₀ -B₁ -B₂ ⋯ -B_{ORD-1} ]

where I is the n_fom × n_fom identity matrix.

Base.show(io::IO, ::MIME"text/plain", m::NthOrderModel)

Print a summary of the model.

Without this, showing a model dumps every field — which for a FEM backend means the entire DofHandler behind each nonlinear term. Everything printed here is already on the type; nothing is computed.

The one-line method exists for the same reason: a model nested inside another container would otherwise expand its whole field tree there too.

_mode_numbers(σ, n_eigs) -> Vector{Int}

Number the physical modes of a spectrum: entry i gets out[i], the index of the conjugate pair it belongs to.

The pairs are the orbits of σ, numbered by first appearance — walk the spectrum in order and, on reaching an unvisited entry, assign the next number to it and to σ(i). A self-paired entry (σ(i) = i, a real eigenvalue) is its own mode.

This deliberately does not compute ⌈i/2⌉. That formula assumes the eigensolver emits conjugate partners adjacently, which is a convention of some solvers rather than a fact about spectra — a shift-invert or filtered solve can return {1, 5} as a pair. Deriving the numbering from σ agrees with the adjacent case and stays correct otherwise.

With no conjugate structure (σ === nothing), every entry is its own mode.

_stack_blocks(physical, companions, ROM, side) -> Array{ComplexF64, 3}

Assemble one side's order-blocks, applying the mirrored convention here and nowhere else: the right physical slice is block 1 with its derivatives in 2:ORD, the left physical slice is block ORD with its orthogonality blocks in 1:(ORD-1).

With companions === nothing the caller already holds whole blocks and they are used as given.

_structural_left_eigenmode_orders(λ, Y, mass, damping) -> Array{ComplexF64, 3}

Analytic companion left-eigenvector order-blocks for a proportionally damped second-order structure M ẍ + C ẋ + K x = 0 with real symmetric M, K and C = αM + βK. For the sesquilinear left eigenvector φ (solving φᴴ (λB − A) = 0 on the companion pencil) with real position mode ϕ = Y[:, 1, k]:

φ_2 = ϕ                    (physical slice)
φ_1 = (conj(λ) M + C) ϕ

These are exactly the blocks eigensolve_left would return — algebraically equal to -(1/conj(λ)) Kᵀϕ via the quadratic eigenrelation, but built from the moderate-norm M, C instead of K (which amplifies eigensolver noise by (ω_max/ω₁)²). No adjoint eigensolve is needed because the proportional damping makes the blocks analytic in ϕ.

check_biorthogonality(sd::SpectralData, model) -> Matrix{ComplexF64}

Return the master-block biorthogonality matrix G[r, s] = φᵣᴴ B ψₛ, which should be the identity for biorthogonally normalised eigenvectors.

This is the numerical guard on the mirrored right/left index convention: swapping the physical slice for a derivative block, or the two sides for each other, destroys G ≈ I loudly, whatever accessor was misused.

Diagnostic only — deliberately not called from parametrise or the solve, so it adds no cost to a normal run. Call it in tests, or when a result looks wrong.

detect_conjugate_permutation(lambda; atol = 1e-8) -> Union{Vector{Int}, Nothing}

Attempt to construct a conjugate-permutation vector from the eigenvalue vector lambda (length NVAR). Returns a Vector{Int} perm such that

lambda[perm[i]] ≈ conj(lambda[i])   for all i

and perm[perm[i]] == i (involution), or nothing if no such perfect pairing exists (e.g. an eigenvalue has no conjugate partner within atol).

Warning — necessary but not sufficient. Two eigenvalues forming a conjugate pair does not guarantee that the corresponding eigenvectors satisfy

master_modes[:, perm[r]] ≈ conj(master_modes[:, r]).

This condition can fail when:

• the eigenvalue is degenerate (eigenspace has dimension > 1),

• the solver returned a non-conjugate basis for a repeated eigenvalue,

• eigenvectors were post-processed with different phases or normalisation.

Always verify eigenvector conjugacy (e.g. check norm(master_modes[:, perm[r]] - conj(master_modes[:, r]))) before passing the returned vector to solve_parametrisation as conjugate_permutation. Passing an incorrect permutation silently corrupts W and R.

Arguments

• lambda — eigenvalue vector of length NVAR (master + external eigenvalues).

• atol — absolute tolerance for the conjugate-match test

		  `|lambda[j] - conj(lambda[i])| < atol`.

Returns

Vector{Int} (involution, 1-based) if a perfect pairing is found; nothing otherwise.

eigensolve(model::NthOrderModel, solver::DefaultEigensolver)

Solves right eigenproblem using eigen from LinearAlgebra. Let A and B be the first order matrices of model. Then it returns the eigenpairs

$ (A-\lambda_k B)y_k = 0$

eigensolve(model::NthOrderModel, solver::MorfeEigensolver)

Solves right eigenproblem using generalised_eigenpairs from MORFE.Eigensolvers. Let A and B be the first order matrices of model. Then it returns the eigenpairs

$ (A-\lambda_k B)y_k = 0$

eigensolve(mass::AbstractMatrix, stiffness::AbstractMatrix, solver::StructureModalDampingEigensolver)

Calculates right eigenvectors in the secod order form. Uses the relations:

$ (K-\omega_k^2*M)\y_{k,U}=0 $ $ y_{k,V}=\lambda_k y_{k,U} $

and assumes $C = \alpha*M + \beta*K$.

eigensolve(model::NthOrderModel, solver::StructureModalDampingEigensolver)

Thin wrapper: extracts M = model.linear_terms[end] and K = model.linear_terms[1] and delegates to eigensolve(M, K, solver).

eigensolve_left(model::NthOrderModel, solver::DefaultEigensolver)

Solves left eigenproblem using eigen from LinearAlgebra. let A and B be the first order matrices of model. Then it returns the eigenpairs

$ x_k^H(A-\lambda_k B) = 0$

eigensolve_left(model::NthOrderModel, solver::MorfeEigensolver)

Solves left eigenproblem using σ-shift to recover the correct eigenvalues.

external_conjugate_permutation(sys; atol = 1e-8) -> Union{Vector{Int}, Nothing}

The conjugate involution σ on the N_EXT external variables, or nothing when the external system has no conjugate structure to offer.

σ pairs external variable k with the one carrying conj(λ_k), and — when the system was re-based — additionally guarantees Q[:, σ(k)] == conj(Q[:, k]). Both conditions are needed for the reduction's conjugate symmetry: Realification.realify applies one conj_map across all NVAR variables, external ones included, and fill_conjugate_monomial! implements W_{P·γ} = conj(W_γ). With a real forcing F the external columns satisfy Φ[:, k] = L(λ_k)⁻¹ F Q[:, k], so conj(Φ[:, k]) = Φ[:, σ(k)] follows exactly from the two conditions together.

Eigenvalue pairing alone is not sufficient — see detect_conjugate_permutation's own warning — which is why the basis columns are verified rather than assumed. A system re-based onto a Schur basis generally fails that check and correctly returns nothing: its external variables simply are not conjugate pairs.

Use full_conjugate_permutation to assemble the full NVAR permutation.

full_conjugate_permutation(master_perm, sys) -> Vector{Int}

Assemble the full NVAR-length conjugate_permutation from the master block and the external system, appending ROM .+ σ where σ is the external involution.

Callers otherwise hand-write the whole vector ([2, 1, 3], [2, 1, 4, 3], …), which bakes in a pairing that a change of external coordinates can invalidate, and which has to special -case an odd number of external variables. Deriving the external block instead keeps it correct in both situations.

Throws when the external system has no conjugate structure — see external_conjugate_permutation.

generalised_eigenpairs(A, B; nev, shift=nothing, which=:LM, tol=0.0,
					   maxiter=3000, ncv=nothing, v0=nothing,
					   ritzvec=true, sort_largest_real=false)

Solve the generalised eigenproblem A x = lambda B x using Arpack.

Requires Arpack.jl and LinearMaps.jl. Load them to activate the MORFE extension.

indices(b::ModeBundle) -> Vector{Int}

The positions this bundle's modes occupy in the source spectrum.

Selecting master modes discards the original numbering, so it is recorded here. Two things need it: restricting the spectrum-wide conjugate involution to this bundle, and reporting a mode by the entry a user would index in their own spectrum. Conjugate partners need not be adjacent, so these are not recoverable by arithmetic.

left_eigenmode_orders_from_slice(linear_terms, left_slice, eigenvalues)
-> Array{ComplexF64, 3}

Reconstruct the full left eigenvector order-blocks from the physical-space slice, for callers that compute only the latter. For the sesquilinear left eigenvector φ = [φ_1; …; φ_ORD] of reported eigenvalue λ (physical slice ℓ = φ_ORD, satisfying L(λ)ᴴ ℓ = 0 with L(s) = Σ_k B_k s^k), the companion block equations give

φ_{ORD-1} = conj(λ) · (B_ORDᴴ ℓ) + B_{ORD-1}ᴴ ℓ
φ_j       = conj(λ) · φ_{j+1}   + B_jᴴ ℓ          j = ORD-2, …, 1

The eigenvalue is used only to define the eigenvector from its slice — the per-monomial cohomological solve reads the blocks and never touches it. Prefer eigensolvers that return the full blocks directly (eigensolve_left does); use this only when a slice is all you have.

The slice may carry an arbitrary per-mode scale: the blocks scale with it, and the orthogonality equations are invariant under per-mode row scaling.

apply — which operator each block is hit with

apply(B_j) is what multiplies ℓ; it defaults to adjoint, giving the B_jᴴ form derived above, which is correct for a general pencil.

Pass apply = identity when the pencil is self-adjoint — B_j real symmetric for every j, so L(s)ᵀ = L(s). Two things then hold simultaneously: B_jᴴ = B_j, and the left eigenvector is the conjugate of the right one, so a right position mode may legitimately be handed in as left_slice. That is exactly the structural case, and _structural_left_eigenmode_orders is the thin wrapper for it.

The keyword exists rather than hard-coding adjoint because exact symmetry is expected of assembled M, C but not guaranteed. When B is exactly symmetric the two routes agree bitwise (measured, dense and sparse alike: symmetry makes the scatter and gather traversals visit the same indices in the same order). But a single ulp of asymmetry — Ke[i, j] and Ke[j, i] are separate floating-point expressions in an element routine — makes B * x and B' * x differ at round-off. apply = identity therefore keeps the structural path bit-for-bit unconditionally, without resting on an assumption about the assembler, and states the self-adjointness claim at the call site.

Sharing the recurrence is the point: the φ_ORD = ℓ fill-downward index convention is the easy thing to get backwards, and it is now written once.

left_mode_blocks(b::ModeBundle) -> view or nothing

Blocks 1:(ORD-1) of the left eigenvectors — the ones feeding the orthogonality row operators. nothing when ORD == 1.

left_modes(b::ModeBundle) -> Matrix{ComplexF64}

The physical-space left eigenvectors, FOM × n — the highest-order block, mirroring right_modes, which is the lowest. Cached at construction.

master_by_sorting(nev::Integer) -> Vector{Int}

The first nev spectrum entries, as master indices — the eigenpairs were already sorted by spectrum.

Pure: it returns indices for SpectralData(model, sp; master = …) rather than marking them on the spectrum. The mutating select_master_modes_by_sorting it replaces wrote a mask into the Spectrum, which made the same object mean different things at different points in a script and left the selection invisible at the call site that used it.

master_by_target_frequency(sp::Spectrum, target_frequencies, tol) -> Vector{Int}

The spectrum entries within tol of any of target_frequencies, as master indices.

Distance used: dist(a, b) = abs(real(a - b)) + abs(imag(a - b)). A target that matches nothing warns and contributes no index — it is a mis-specified target, not an empty selection, and silently returning fewer masters than asked for is the failure mode worth naming.

master_conjugate_permutation(sd) -> Union{Nothing, Vector{Int}}

σ restricted to the master modes and re-indexed to 1:ROM — the form the cohomological solve consumes, before it is extended over the external variables.

Computed once in the constructor; this is a field read. Every solve consults it, so it is not something to re-derive per call.

normalise_biorthogonal!(
	model::NthOrderModel,
	eigenmodes::Matrix{T},
	left_eigenmodes::Matrix{T})

Normalise eigenmodes to fulfill the equation xi^H * B * yj = δ_{ij}

Both sides are scaled symmetrically: with s = sqrt(x_i^H B y_i), the right eigenmode is divided by s and the left eigenmode by conj(s), so the sesquilinear pairing becomes exactly 1.

outer_conjugate_permutation(sd) -> Union{Nothing, Vector{Int}}

σ restricted to the outer modes and re-indexed to 1:n_outer. Used to group conjugates when reporting off-manifold near-resonances.

Computed once in the constructor; this is a field read.

physical_mode(sd::SpectralData, i::Integer) -> Int

The physical mode number of spectrum entry i — conjugate partners share a number.

right_mode_derivatives(b::ModeBundle) -> view or nothing

Blocks 2:ORD of the right eigenvectors — the time derivatives ψ_{k+1} = λ ψ_k. nothing when ORD == 1.

right_modes(b::ModeBundle) -> Matrix{ComplexF64}

The physical-space right eigenvectors, FOM × n. Cached at construction; no indexing.

sort_by_magnitude!(eigenvalues, eigenmodes)

Sorts eigenpairs to resemble the order: |λ[1]| ≤ |λ[2]| ≤ ... where λ=eigenvalues.

sort_left_eigenmodes(eigenvalues, left_eigenvalues, left_eigenmodes)

Sorts the left eigenpairs to match right eigenpairs, by using the distance function dist(a, b) = abs(real(a - b)) + abs(imag(a - b))

spectrum(model::NthOrderModel, solver::StructureModalDampingEigensolver; sorter!)

Specialised path for StructureModalDampingEigensolver: mass-normalisation is built into eigensolve, and the left eigenvector order-blocks are analytic in the position mode (see _structural_left_eigenmode_orders). No adjoint solve or biorthogonal normalisation is needed.

The optional sorter! kwarg has the same semantics as in the general eigensolve: pass (args...) -> nothing to preserve the solver's natural ordering.

spectrum(stiffness, mass, solver::StructureModalDampingEigensolver; sorter!)

Convenience overload: pass K and M directly without constructing an NthOrderModel. The Rayleigh damping matrix C = αM + βK is rebuilt from the solver parameters for the left eigenvector order-blocks.

spectrum(
	model::NthOrderModel;
	solver::AbstractEigensolver = DefaultEigensolver(),
	sorter!::Function = sort_by_magnitude!,
	normaliser!::Function = normalise_biorthogonal!)

Computes left and right eigenpairs of the problem described in model by using the defined solver. Additionally sorter! sorts the eigenpairs and normaliser! is used to normalise the eigenmodes.

spectrum_entries(sd::SpectralData, p::Integer) -> Vector{Int}

The spectrum entries making up physical mode p — one for a real mode, two for a conjugate pair. Not necessarily consecutive, which is why they are looked up rather than computed.

_realify_term(exp_vec::SVector{N,Int}, coeff::C, n::Int)
	-> Dict{SVector{N,Int}, C} where {C,N}

Transform a single term (exponent vector exp_vec and coefficient coeff) of a polynomial in the canonical form (z, z̄, w) into a sum of real monomials. Returns a dictionary mapping new exponent vectors (in the real variables) to their coefficients.

Here N = 2n + m, with n conjugate pairs and m real variables.

_reorder_canonical(poly::DensePolynomial{C,N}, conj_map::Vector{Int})
	-> (DensePolynomial{C,N}, n, m)

Reorder variables according to a conjugation map conj_map of length N (where N = number of variables).

• conj_map[i] = j means variable i is conjugate to variable j.

• If variable i is real, then conj_map[i] = i.

The reordering groups variables as (z₁, …, zₙ, conj(z₁), …, conj(zₙ), w₁, …, wₘ) where n is the number of conjugate pairs and m the number of real variables. Terms with the same exponent after reordering are merged.

Returns the canonical polynomial (same concrete type as poly), n, and m.

realify(poly::DensePolynomial, conj_map::Vector{Int}) -> DensePolynomial

Transform a complex‑valued polynomial (with variables that may be conjugate pairs) into a polynomial in real variables.

Arguments

• poly: a polynomial in variables z₁, …, z_N.

• conj_map: a vector of length N where conj_map[i] = j means variable i is the conjugate of variable j; if i is real, then conj_map[i] = i.

Returns

A new polynomial in real variables x₁, …, x_n, y₁, …, y_n, w₁, …, w_m with n conjugate pairs and m real variables. The transformation uses the formulas z = x + i y, z̄ = x - i y. The returned polynomial has the same concrete type as the input poly (including the same number of variables).

realify_via_linear(poly::DensePolynomial, conj_map::Vector{Int}) -> DensePolynomial

Transform a complex‑valued polynomial into a polynomial in real variables by composing with the linear map that expresses complex variables in terms of real and imaginary parts. This is an alternative implementation to realify that uses the compose_linear function. The returned polynomial has the same concrete type as the input poly (real coefficients).

See also: realify, compose_linear

Build the n_int × NMON inner resonance matrix.

strategy applies graph/normal-form unconditional flags; inner_cond (optional) applies an eigenvalue-proximity check on the master eigenvalues.

Build the n_out × NMON outer resonance matrix using outer_cond.

_first_close_pair(eigenvalues, tol, n) -> (i, j, gap)

The first pair of eigenvalues separated by no more than tol, or (0, 0, 0.0) if there is none. Allocation-free, and stops at the first hit.

Split out from the caller so the @info that reports it sits outside the loop: string interpolation inside a loop body captures the loop variables and boxes them on every iteration, even when the branch is not taken and nothing is logged.

_outer_warn_tolerance(config, outer_eigenvalues)
	-> Union{Nothing, Float64, Vector{Float64}}

The threshold the off-manifold scan compares |λ_outer[j] - s| against, resolved without going through resolve_tolerances.

That function is wrong for this caller three times over: it re-emits every @info guard a second time, tol_relative makes it build n_monomials identical tolerance rows where the scan wants one, and for an explicitly per-target tol it returns nothing for the outer family — a vector sized for the inner targets cannot be indexed by an outer target number. That last case used to reach _resolve_outer_tol and throw, so a per-target tolerance combined with the default warn_outer = true aborted the solve. It now returns nothing, and the scan is skipped with a notice rather than taking the run down with it.

_resolve_outer_tol(tol, outer_tol, n_out, n_int)

Pick the tolerance for the outer targets, and reject the one combination that used to be a silent bounds error.

A per-target tolerance is read tol[k][local_idx] with local_idx local to its own condition — 1:n_int for the inner block, 1:n_out for the outer one. The two blocks are built from separate condition objects, so they can perfectly well carry separate tolerances; only the public constructors' single tol argument tied them together, and a per-target vector sized for n_int then overran whenever n_out > n_int.

outer_tol === nothing means "reuse tol", which is exactly right for a scalar (it applies to any number of targets) and impossible for a per-target vector — hence the error, which tells the caller what to pass instead of failing deep inside is_resonant.

_warn_outer_resonances(mset, master_eigs, outer_eigs, external_eigs, config, spectral)

Warn when a monomial is near-resonant with a physical mode that is not on the manifold.

That direction is then solved through a near-singular operator, so the ROM loses accuracy there regardless of how the load is shaped: solve_single_monomial! builds its operator from s = ⟨λ, α⟩ alone, which makes the conditioning independent of the right-hand side. Rounding injects a component along the near-null direction and 1/(λ_s - s) amplifies it.

Runs for autonomous and forced models alike — a monomial built purely from master coordinates that lands on an off-manifold eigenvalue is exactly as near-singular as a forced one, and that its cause is the chosen master set makes it more worth surfacing.

Why this does not build a ResonanceSet

It used to, and that made an on-by-default diagnostic cost O(NMON × n_outer) with a large constant: a whole second set including an n_int × NMON inner block that was discarded, _superharmonics allocating a temporary per monomial twice, and _local_index's findfirst inside is_resonant turning the outer build into O(NMON × n_outer²).

The criterion needs none of that — it is one distance test per (monomial, outer target) — so it is written out directly. A run that flags nothing now costs a constant 64 bytes (measured, unchanged from |mset| = 20 to 54) against the 246 400 the probe cost at |mset| = 35 with 58 outer modes. The test itself is unchanged: |λ_outer[j] - s| < tol, the same EigenvalueCondition comparison the probe applied, for every style. Reusing the already-built outer block instead was rejected deliberately: under :real_normal_form that block ORs each target with its conjugate, which would silently change what gets reported.

apply_internal_resonances!(mat, strategy, mi, n_int, k)

Set inner-resonance flags in column k of the inner matrix mat for the monomial with exponent vector mi.

• GraphInternal: all n_int rows flagged for degree ≥ 2; for a linear monomial eᵣ only row r (if r ≤ n_int) is flagged.

• NormalFormInternal: no-op.

build_resonance_set(model, mset, spectral::SpectralData, config::ResonanceConfig)
	-> ResonanceSet

Build the resonance set from a SpectralData bundle and a ResonanceConfig — the single entry point for resonance construction.

Reads the master and outer eigenvalues off spectral (replacing the eigenproblem plus master-mask plumbing), resolves the config's tolerances into correctly-sized inner and outer objects via resolve_tolerances, and warns about off-manifold near-resonances when config.warn_outer is set.

empty_resonance_set(multiindices, n_internal, n_outer=0) -> ResonanceSet

Construct a ResonanceSet with all resonance flags set to false. n_internal = number of master modes (ROM); n_outer = number of outer targets (0 = none).

is_resonant(rs, idx, target) -> Bool
is_resonant(rs, mi, target)  -> Bool

Return true when the monomial at position idx (or exponent vector mi) is resonant with target target. Targets 1:ROM query inner_resonances; targets > ROM query outer_resonances (returns false when outer is nothing).

n_internal(rs::ResonanceSet{ROM}) -> ROM

Return the number of internal master modes (compile-time constant from type parameter).

project_eigenvalues(λ, projection::Symbol) -> Vector

Apply a ResonanceConfig eigenvalue projection: :full is the identity, returning Vector{ComplexF64}; :imaginary_part_only returns the imaginary parts Im(λ) as Vector{Float64}.

A projected eigenvalue i·Im(λ) is carried by the real number Im(λ) — the two are isomorphic for this purpose, because for real a, b

|i·a − i·b| = |a − b|

and the whole of detection is abs(eig - s). So carrying the frequencies as Float64 reproduces the pure-imaginary complex computation exactly, with no i·Im reconstruction and one code path serving both modes.

Project every family the same way — master, external and outer. Mixing a projected family with an unprojected one promotes the superharmonic s = ⟨λ, α⟩ back to complex and reinstates the growth rates through the back door, which is precisely what the projection was asked to remove.

Always returns a plain Vector. The master eigenvalues arrive as a static vector, and a comprehension over one stays static — a SizedVector that the resonance_set_from_* signatures reject. The loops below are the same reason the caller used Vector{ComplexF64}(…) rather than collect.

resolve_tolerances(config, master_eigenvalues, outer_eigenvalues, n_monomials)
	-> (inner_tol, outer_tol)

Turn a ResonanceConfig into the two correctly-sized tolerance objects the resonance-set constructors want, and emit @info for settings that will have no effect.

With tol_relative, each family is sized for its own target count:

inner[k][r] = tol_relative * |λ_master[r]|      length n_int
outer[k][j] = tol_relative * |λ_outer[j]|       length n_out

which is why the two can now be combined at all.

resonance_set_from_complex_normal_form_style(
	multiindices, master_eigenvalues, tol;
	external_eigenvalues, outer_eigenvalues)

Build a ResonanceSet using the complex normal form style: inner resonances flagged by |λᵣ - s| < tol for each master mode r; outer resonances (if any) flagged by proximity to outer_eigenvalues.

Suitable for autonomous Invariant Manifolds (e.g. NNMs) with complex conjugate reduced variables; add external_eigenvalues for non-autonomous systems where the multiindex includes forcing directions.

resonance_set_from_condition_number_estimate(
	multiindices, master_eigenvalues, spectral_radius,
	target_condition_numbers, max_cond;
	external_eigenvalues, outer_eigenvalues,
	inner_target_indices, outer_target_indices, conjugacy_map)

Build a ResonanceSet using a condition-number criterion:

|λⱼ - s| * max_cond < spectral_radius * κ(λⱼ)

target_condition_numbers has n_int + n_out entries: the first n_int are for the master modes, the remainder for the outer modes.

inner_target_indices / outer_target_indices restrict which rows in each sub-matrix are populated (default: all).

resonance_set_from_graph_style(
	multiindices, master_eigenvalues, external_eigenvalues,
	outer_eigenvalues, tol)

Build a ResonanceSet using the graph style: every monomial of total degree ≥ 2 is marked resonant with all ROM master modes. Outer targets are flagged by eigenvalue proximity |λⱼ - s| < tol.

• master_eigenvalues: ROM eigenvalues; enter s and are inner targets.

• external_eigenvalues: enter s only (e.g. forcing frequencies in the multiindex, since external_eigenvalues cant be targets).

• outer_eigenvalues: outer targets (eigenvalues not included in master_eigenvalues, tested for near-resonance). Pass ComplexF64[] when there are no outer targets.

resonance_set_from_graph_style(multiindices, master_eigenvalues, external_eigenvalues, outer_condition)

Advanced overload that accepts a pre-built OuterResonanceCondition for the outer rows. outer_condition.target_indices must use local indices 1:n_out. n_out is inferred as maximum(outer_condition.target_indices) (or 0 if empty).

resonance_set_from_real_normal_form_style(
	multiindices, master_eigenvalues, conjugacy_map, tol;
	external_eigenvalues, outer_eigenvalues)

Build a ResonanceSet using the real normal form style: like CNF but conjugate pairs share the resonance flag — monomial k is resonant with target j when |λⱼ - s| < tol OR |λ_{conj(j)} - s| < tol.

conjugacy_map has length(master_eigenvalues) + length(outer_eigenvalues) entries; the first n_int cover inner targets, the remainder cover outer targets (re-indexed locally). conjugacy_map[i] is the local index of the conjugate of target i.

resonant_multiindices(rs, target) -> Vector{Int}

Return the positions of all monomials resonant with target. Targets 1:ROM query inner_resonances; targets > ROM query outer_resonances.

resonant_targets(rs, idx) -> AbstractVector{Bool}
resonant_targets(rs, mi)  -> Union{AbstractVector{Bool}, Nothing}

Return a boolean vector indicating which targets are resonant with the monomial at position idx (or exponent vector mi). Concatenates inner and outer rows. Returns nothing when mi is not in the multiindex set.

set_resonance!(rs, target, idx, value) -> rs
set_resonance!(rs, target, mi, value)  -> rs

Set the resonance flag for target target and monomial idx (or multiindex vector mi) to value. Targets 1:ROM address inner_resonances; targets > ROM address outer_resonances. Returns rs for chaining. Warns if mi is not found.

assemble_cohomological_matrix_and_rhs!(M, rhs, s, linear_terms, C_coeffs, E_coeffs,
										resonance, lower_order_couplings,
										external_dynamics, g_buffer) → nothing

In-place variant: writes the invariance-equation system matrix and RHS directly into the caller-supplied M and rhs buffers. No heap allocation occurs.

M must have size FOM × (FOM + ROM) — the border is constant width and the resonance mask selects its content rather than its size: column FOM + r holds C_r(s) when mode r is resonant and is zeroed otherwise. The zeroed columns are harmless because the matching orthogonality row pins R_{r,α} = 0 (see assemble_orthogonality_matrix_and_rhs!).

rhs must have length FOM. Both are overwritten on entry. g_buffer is the pre-allocated FOM-length scratch buffer for the external RHS.

build_sparse_L_and_rhs!(rhs, L_template, mappings, linear_terms, s, lower_order_couplings)
-> L_template

In-place Horner evaluation of the parametrisation operator using a pre-allocated union-pattern template and index mappings from precompute_sparse_L_template.

Avoids ALL per-monomial sparse arithmetic allocations regardless of whether the linear_terms share a common sparsity pattern. The template's nzval is overwritten on each call; its colptr/rowval are never modified.

evaluate_column!(c, s, r, C_coeffs) -> c

Evaluate the r-th reduced-dynamics operator column

C_r(s) = Σ_{L=1}^{ORD} C_coeffs[r][:, L] · s^{L-1}

in-place via Horner's method, overwriting the pre-allocated FOM-vector c.

c may be a plain Vector{T} or a column view view(M, :, col).

Horner recurrence

L runs from ORD-1 down to 1:

c  ←  C_coeffs[r][:, ORD]              (highest-degree coefficient)
for L = ORD-1, …, 1:
	c ← c · s + C_coeffs[r][:, L]

Column access C_coeffs[r][:, L] reads contiguous memory (Julia is column-major), so the loop touches sequential cache lines.

Arguments

• c :: AbstractVector{T} – output buffer (length FOM), overwritten.

• s :: T – evaluation frequency.

• r :: Int – 1-based master-mode index (1 ≤ r ≤ ROM).

• C_coeffs :: Vector{<:AbstractMatrix{T}} – pre-computed coefficients from precompute_column_polynomials; C_coeffs[r] is FOM × ORD.

Complexity

O(ORD · FOM)

evaluate_external_rhs!(rhs, s, external_dynamics, E_coeffs, g) -> rhs

Accumulate the external-forcing contribution to the cohomological right-hand side:

rhs += Σ_{e=1}^{N_EXT} external_dynamics[e] · E_e(s)

where E_e(s) = Σ_{L=1}^{ORD} E_coeffs[e][:, L] · s^{L-1} is the e-th external column polynomial (pre-computed by precompute_column_polynomials).

The sign flip was already absorbed into the pre‑computed coefficients E_coeffs, since they originate from the left‑hand side of the invariance equation and are moved to the right‑hand side, so the contribution is added to the RHS rather than subtracted.

Sparse exploitation

Only the non-zero entries of external_dynamics are processed. For periodic forcing of a few harmonics this is typically a small subset of N_EXT.

Combined Horner pass

Rather than evaluating each E_e(s) independently, the non-zero contributions are combined into a single degree-(ORD-1) vector polynomial

g(s) = Σ_{e active} external_dynamics[e] · E_e(s)

and evaluated in one Horner pass (ORD-1 scalar-vector multiplies plus ORD-1 saxpy operations over FOM), instead of N_EXT_active separate Horner passes.

Arguments

• rhs :: AbstractVector{T} – accumulator (length FOM), updated in-place.

• s :: T – evaluation frequency.

• external_dynamics :: AbstractVector{T} – known amplitudes of the N_EXT external forcing modes; typically sparse.

• E_coeffs :: Vector{<:AbstractMatrix{T}} – pre-computed external coefficients from precompute_column_polynomials; E_coeffs[e] is FOM × ORD.

• g :: AbstractVector{T} – pre-allocated FOM-length scratch buffer; zeroed on entry.

Complexity

• O(N_EXT_active · FOM · ORD) for combining coefficients.

• O(FOM · ORD) for the single Horner evaluation.

evaluate_system_matrix_and_lower_order_rhs!(parametrisation_operator, lower_order_rhs, s, lower_order_couplings, linear_terms)
-> parametrisation_operator

Evaluate the parametrisation operator L(s) and accumulate the lower-order right-hand-side contributions in a single Horner pass, reusing the transient intermediate matrices that are available only during the polynomial evaluation.

Mathematical context

At step j of the Horner recurrence (before the scalar multiply by s), the intermediate matrix

L[j](s) = Σ_{k=j+1}^{ORD+1} B[k] · s^{k-(j+1)}

is available. Multiplying by the pre-computed coupling vector ξ[j] = lower_order_couplings[j] gives the contribution of lower-order solution terms at derivative order j to the right-hand side:

contribution[j] = -L[j](s) · ξ[j]

The negative sign arises because these terms originate from the left-hand side of the cohomological equation and are transposed to the right-hand side.

Summed over j = 1, …, ORD, the full lower-order RHS is

lower_order_rhs = -Σ_{j=1}^{ORD} L[j](s) · ξ[j] = -Σ_{j=1}^{ORD} ( Σ_{k=j+1}^{ORD+1} B[k] · s^{k-(j+1)} ) · ξ[j]

This computation must share the Horner loop with L(s): the L[j] intermediates are transient, and recomputing them would double the O(ORD · FOM²) work.

The coupling vectors are obtained from MORFE.LowerOrderCouplings.compute_lower_order_couplings applied to the lower-order multi-indices associated with each Horner step.

Arguments

• parametrisation_operator :: AbstractMatrix{T} – output buffer (FOM × FOM), overwritten with L(s) = Σ_{k=1}^{ORD+1} B[k] · s^{k-1}.

• lower_order_rhs :: AbstractVector{T} – accumulator (length FOM), updated in-place. Must be initialised to zero (or the desired starting value) by the caller.

• s :: T – evaluation superharmonic.

• lower_order_couplings :: SVector{ORD, <:AbstractVector{T}} – coupling vectors ξ[j] for j = 1,…,ORD; each element is an AbstractVector{T} of length FOM.

• linear_terms :: NTuple{ORD+1, <:AbstractMatrix{T}} – linear_terms[k] = B[k].

Complexity

O(ORD · FOM²), shared with the L(s) evaluation.

precompute_column_polynomials(fom_matrices, generalised_right_eigenmodes,
							  reduced_dynamics_linear, ROM)
-> (C_coeffs, E_coeffs)

Pre-compute the polynomial coefficient arrays for the cohomological operator columns (master modes) and the external-forcing right-hand-side columns.

These arrays are computed once per polynomial order and reused for every multi-index at that order.

Return values

• C_coeffs :: Vector{Matrix{T}} of length ROM, where C_coeffs[r] is FOM × ORD. Column j of C_coeffs[r] is the degree-(j-1) coefficient of the r-th reduced-dynamics operator column:

  C_r(s) = Σ_{j=1}^{ORD} C_coeffs[r][:, j] · s^{j-1}

• E_coeffs :: Vector{Matrix{T}} of length N_EXT = NVAR - ROM, where E_coeffs[e] is FOM × ORD. Column j of E_coeffs[e] is the degree-(j-1) coefficient of the e-th external-forcing operator column:

  E_e(s) = Σ_{j=1}^{ORD} E_coeffs[e][:, j] · s^{j-1}

Arguments

• fom_matrices :: NTuple{ORD+1, <:AbstractMatrix{T}} – linear matrices of the full-order model; fom_matrices[k+1] corresponds to B[k] (0-indexed in the ODE).

• generalised_right_eigenmodes :: AbstractMatrix{T} of size FOM × NVAR – generalised eigenvectors; columns 1:ROM are the master modes, columns ROM+1:NVAR are the external forcing modes.

• reduced_dynamics_linear :: AbstractMatrix{T} of size NVAR × NVAR – Jordan-form matrix of the linear part of the reduced dynamics on the Invariant Manifold.

• ROM :: Int – number of master modes (dimension of the reduced-order model).

Recurrence

The output matrices are filled by a single downward Horner recurrence (j runs from ORD down to 1) using one FOM × NVAR working buffer D:

D ← B[ORD+1] * generalised_right_eigenmodes                               (j = ORD)
C_coeffs[r][:, j] ← D[:, r]       for r = 1…ROM
E_coeffs[e][:, j] ← D[:, ROM+e]   for e = 1…N_EXT

D ← D * reduced_dynamics_linear + B[j+1] * generalised_right_eigenmodes  (j = ORD-1, …, 1)
C_coeffs[r][:, j] ← D[:, r]       for r = 1…ROM
E_coeffs[e][:, j] ← D[:, ROM+e]   for e = 1…N_EXT

After step j, column j of every per-target matrix holds the exact degree-(j-1) coefficient

D[:, ·] = Σ_{k=j+1}^{ORD+1} B[k] * generalised_right_eigenmodes * reduced_dynamics_linear^{k-(j+1)}

Complexity

• Time: O(ORD · FOM · NVAR) (dominated by ORD matrix–matrix products)

• Storage: O(ORD · FOM · NVAR)

precompute_external_column_polynomials(fom_matrices, external_directions,
										reduced_dynamics_linear, D_master_steps)
-> E_coeffs

Φext-dependent half of [`precomputecolumnpolynomials](@ref). Computes the external-mode column polynomialsEe(s)reusing the pre-saved master Horner intermediatesDmastersteps(from [precomputemastercolumn_polynomials`](@ref)) instead of recomputing the master-column work.

Pass external_directions = zeros(FOM, N_EXT) to obtain the partial (Φext = 0) Ecoeffs needed for the initial external-forcing solve.

Arguments

• D_master_steps — from precompute_master_column_polynomials; length ORD, each FOM × ROM. D_master_steps[j] is the master Horner buffer at step j.

precompute_master_column_polynomials(fom_matrices, master_modes, Λ_master)
-> (C_coeffs, D_master_steps)

Φext-independent half of [`precomputecolumnpolynomials](@ref). Computes only the master-mode column polynomialsCr(s)and saves the intermediate FOM×ROM Horner buffer at every step for later reuse by [precomputeexternalcolumn_polynomials`](@ref).

Because Λ is upper-triangular, the master columns of the Horner buffer D form a closed subsystem under the recurrence D ← D·Λ + B[j+1]·Y: for column r ≤ ROM, (D·Λ)[:,r] = Σ_{k≤r} D[:,k]·Λ[k,r] depends only on master columns. Therefore C_coeffs is independent of the external directions Φ_ext.

Returns

• C_coeffs :: Vector{Matrix{T}} — same layout as from precompute_column_polynomials.

• D_master_steps :: Vector{Matrix{T}} — length ORD; D_master_steps[j] is the FOM×ROM master block of the Horner buffer at the step that wrote C_coeffs[:,j]. Used by precompute_external_column_polynomials to avoid recomputing master work.

precompute_sparse_L_template(linear_terms) -> (L_template, mappings)

Pre-allocate a SparseMatrixCSC{ComplexF64} with the union sparsity pattern of all linear_terms, and compute index mappings so that each entry of linear_terms[k] can be accumulated into the correct position of the template's nzval array.

Returns (L_template, mappings) where mappings[k][i] is the index into L_template.nzval that corresponds to position i of linear_terms[k].nzval.

Used once at context construction; the template is reused across all monomials by the in-place build_sparse_L_and_rhs! overload, eliminating per-monomial sparse arithmetic allocations even when the input matrices do not share a common pattern (e.g. when the damping matrix is C = α*M + β*K and Julia's sparse addition drops entries that cancel to exactly zero).

precompute_sparse_bordered_template(L_template, ROM) -> (M, border_row_base)

Allocate the constant-size (FOM+ROM) × (FOM+ROM) bordered cohomological matrix

	┌                              ┐
	│  L(s)      C(s) P            │   FOM rows  (invariance)
	│  P Ĵ(s)    P Ĉ(s) P + τ Q    │   ROM rows  (orthogonality / R_α = 0)
	└                              ┘

whose sparsity pattern depends only on the union pattern of L_template and on ROM — never on the resonance mask P = diag(ρ) of the monomial being solved. Non-resonant border entries are carried as numeric zeros in structural positions, which is exactly what allows one symbolic factorisation to be reused for every monomial. The system being represented is documented in the CohomologicalEquations module docstring; what follows is this function's own contract — where each block lands in the CSC arrays, which the assembly in _solve_monomial! writes to directly.

Block layout, in CSC order:

block	rows	cols	pattern
L	1:FOM	1:FOM	union pattern of L_template
C P	1:FOM	FOM+1:end	dense FOM × ROM
P Ĵ	FOM+1:end	1:FOM	dense ROM × FOM
P Ĉ P + τ Q	FOM+1:end	FOM+1:end	dense ROM × ROM

so nnz(M) = nnz(L_template) + 2·FOM·ROM + ROM². Appending the border rows to each of the first FOM columns preserves sorted rowval because every union row index is ≤ FOM.

Returns

• M :: SparseMatrixCSC — the template; only nzval is ever written afterwards.

• border_row_base :: Vector{Int} — length FOM; entry M[FOM+r, c] for c ≤ FOM lives at M.nzval[border_row_base[c] + r - 1].

No L → M index table is returned or needed. Because the L entries of column c are laid down as a contiguous prefix of that column, the map is affine within each column,

	L_template.nzval[p]  ↦  M.nzval[M.colptr[c] + (p - L_template.colptr[c])]

so scatter_L_into_bordered! is a per-column block copy rather than an indirect gather — which also saves an nnz(L)-length index vector.

Border column positions need no table either: column FOM+q starts at bq = M.colptr[FOM+q], so C_q(s) occupies the contiguous run M.nzval[bq : bq+FOM-1] and the corner entry M[FOM+m, FOM+q] sits at bq + FOM + m - 1.

scatter_L_into_bordered!(M, L_template) -> M

Copy the freshly evaluated L(s) from the standalone Horner workspace into the (1,1) block of the bordered template.

L's entries for column c sit contiguously in both matrices — at L_template.colptr[c] … and at M.colptr[c] … respectively (see precompute_sparse_bordered_template) — so each column is a single copyto! block move rather than an element-by-element gather through an index table. Keeping this a separate step is what leaves build_sparse_L_and_rhs! untouched: it needs its own square workspace for the transient Horner intermediates L[j](s) that accumulate the lower-order RHS.

assemble_orthogonality_matrix_and_rhs!(M, rhs, s, J_coeffs, C_coeffs, E_coeffs,
										resonance, lower_order_couplings,
										external_dynamics) → nothing

In-place variant: writes the orthogonality block and its RHS directly into the caller-supplied M (ROM × (FOM+ROM)) and rhs (length ROM) buffers. No heap allocation occurs.

The block has constant size — one row per master mode, resonant or not — and the resonance vector selects each row's content rather than the block's dimensions:

• resonant r: row r is the orthogonality condition Ĵ_r(s) W_α + Σ_m Ĉ_{rm}(s) R_{m,α} = g_{r,α}, with the corner entries masked to the resonant modes by evaluate_orthogonality_column_row!;

• non-resonant r: row r becomes the trivial equation τ R_{r,α} = 0 — everything zeroed except M[r, FOM+r] = τ = 1 — which encodes the style choice that non-resonant reduced-dynamics coefficients vanish.

Keeping the width constant is what lets the sparse path reuse one symbolic factorisation across every monomial. Note that τ is structural only: the solver never reads the trivial rows back out of the solution vector, it writes the hard zeros directly during unpacking.

evaluate_orthogonality_column_row!(c, s, r, C_coeffs, resonance) -> c

Evaluate the joint operator row C_r(s) in-place via Horner's method, overwriting the pre-allocated length-ROM vector c.

C_r(s) = Σ_{j=1}^{ORD-1} C_coeffs[r][j, :] · s^{j-1} is a 1 × ROM row polynomial. Entry c[j] holds C_r(s)[j] when master mode j is resonant and zero(T) otherwise — i.e. the layout is expanded and masked, indexed by the mode j itself rather than compacted into resonant rank order.

Masking rather than compacting is what keeps the bordered cohomological system at the constant size FOM + ROM for every monomial (one sparsity pattern, one cached symbolic factorisation). Dropping the non-resonant entries is lossless: the matching orthogonality row pins R_{j,α} = 0, so those coefficients multiply zero.

c may be a plain Vector{T} or a row view view(M, row, col_range).

Horner recurrence (column-wise, no allocation)

For each resonant column index j independently:

val  ←  C_coeffs[r][ORD-1, j]
for L = ORD-2, …, 1:
	val ← val · s + C_coeffs[r][L, j]
c[j] ← val

Column j of C_coeffs[r] is contiguous in memory (Julia is column-major), so each per-column Horner pass is cache-friendly.

Arguments

• c :: AbstractVector{T} – output buffer (length ROM), fully overwritten: resonant entries with C_r(s), the rest with zero.

• s :: T – evaluation frequency.

• r :: Int – 1-based master-mode index for the row equation (1 ≤ r ≤ ROM).

• C_coeffs :: Vector{<:AbstractMatrix{T}} – pre-computed coefficients from precompute_orthogonality_column_polynomials; C_coeffs[r] is (ORD-1) × ROM.

• resonance :: SVector{ROM, Bool} – resonance[j] is true iff master mode j is resonant at the current multi-index.

Complexity

O((ORD-1) · |R|), with no heap allocation.

evaluate_orthogonality_external_rhs(s, r, external_dynamics, E_coeffs) -> T

Compute the scalar external-forcing contribution to the right-hand side of the orthogonality equation for master mode r:

RHS_ext_r = -Σ_{e active} external_dynamics[e] · E_r_e(s)

where E_r_e(s) = Σ_{j=1}^{ORD-1} E_coeffs[r][j, e] · s^{j-1} is the scalar polynomial for forcing mode e in the row equation for mode r (pre-computed by precompute_orthogonality_column_polynomials).

The negative sign reflects that these terms are moved from the left-hand side of the cohomological equation to the right-hand side.

Sparse exploitation

Only the non-zero entries of external_dynamics are processed. For periodic forcing of a few harmonics this is typically a small subset of N_EXT.

Combined Horner pass

The non-zero contributions are combined into a single scalar polynomial

g(s) = Σ_{e active} external_dynamics[e] · E_r_e(s)

and evaluated in one Horner pass (ORD-2 scalar multiplies and ORD-2 · N_EXT_active scalar additions), instead of N_EXT_active separate Horner passes.

Arguments

• s :: T – evaluation frequency.

• r :: Int – 1-based master-mode index (1 ≤ r ≤ ROM).

• external_dynamics :: AbstractVector{T} – known amplitudes of the N_EXT external forcing modes; typically sparse.

• E_coeffs :: Vector{<:AbstractMatrix{T}} – pre-computed coefficients from precompute_orthogonality_column_polynomials; E_coeffs[r] is (ORD-1) × N_EXT.

Returns

The scalar RHS_ext_r = -g(s).

Complexity

O(N_EXT_active · (ORD-1)) for combining coefficients plus O(ORD-1) for the single Horner evaluation.

evaluate_orthogonality_row_and_lower_order_rhs!(row, s,
												lower_order_couplings,
												J_coeffs_r)
-> scalar_rhs :: T

Evaluate the orthogonality row operator Ĵ_r(s) and compute the scalar lower-order right-hand-side contribution for mode r in a single Horner pass, reusing the transient intermediate row vectors.

Mathematical context

At step j of the Horner recurrence (before the scalar multiply by s), the intermediate row vector

Ĵ_r[j](s) = Σ_{k=j+1}^{ORD} J_r[k, :] · s^{k-(j+1)}

is available. Contracting bilinearly with the pre-computed coupling vector ξ[j] gives the scalar contribution of lower-order solution terms at step j:

contribution[j] = -Ĵ_r[j](s) · ξ[j] = -Σᵢ Ĵ_r[j](s)ᵢ · ξ[j]ᵢ

The contraction must not conjugate: the sesquilinear conjugation of the orthogonality condition is already baked into the row coefficients J_r (see precompute_orthogonality_operator_coefficients), so row holds the Horner tail G_{r,j}(s) with the ᴴ applied. The negative sign arises because these terms originate from the left-hand side of the cohomological equation. Summed over j = 1, …, ORD-1:

RHS_lower_r = -Σ_{j=1}^{ORD-1} Ĵ_r[j](s) · ξ[j]

The sum runs to ORD-1 (one fewer than in InvarianceEquation) because the joint operator Q_r has one fewer degree. Sharing the loop with the Ĵ_r(s) evaluation avoids recomputing the Ĵ_r[j] intermediates.

Arguments

• row :: AbstractVector{T} – output buffer (length FOM), overwritten with Ĵ_r(s) = Σ_{j=1}^{ORD} J_r[j, :] · s^{j-1}.

• s :: T – evaluation superharmonic.

• lower_order_couplings :: SVector{ORD_M1, <:AbstractVector{T}} – coupling vectors ξ[j] for j = 1, …, ORD-1; each is a length-FOM vector.

• J_coeffs_r :: AbstractMatrix{T} – ORD × FOM matrix; row j is J_r[j, :], the degree-(j-1) coefficient of Ĵ_r. Obtained from precompute_orthogonality_operator_coefficients.

Returns

The scalar lower-order RHS accumulation RHS_lower_r = -Σ_{j=1}^{ORD-1} Ĵ_r[j](s) · ξ[j].

Complexity

O(ORD · FOM), shared with the Ĵ_r(s) evaluation.

precompute_orthogonality_column_polynomials(J_coeffs, right_master_blocks,
											external_directions,
											reduced_dynamics_linear)
-> (C_coeffs, E_coeffs)

Pre-compute the coefficient arrays of the operators C_r(s) (coupling to the unknown reduced dynamics) and E_r(s) (known external forcing) that appear in the orthogonality equation for master mode r:

J_r(s) W + C_r(s) f_m = − E_r(s) f_e − Σ_k G_{r,k}(s) ξ_k

Mathematical origin

With G_{r,k}(s) = Σ_{j=k+1}^{ORD} J_r[j, :] s^(j-1-k) the Horner tails of the row operator, the couplings are bilinear contractions against the right eigenmode order-blocks Y_k (Y_1 = physical mode, Y_{k+1} = next derivative block):

C_r(s) = Σ_{k=1}^{ORD-1} G_{r,k}(s) · Y_k^m        (master blocks, from the eigensolver)
E_r(s) = Σ_{k=1}^{ORD-1} G_{r,k}(s) · Y_k^e        (external blocks)

The master blocks are supplied directly (right_master_blocks). The external blocks are generalised — the reduced linear dynamics couples them back to the master modes — and are generated by the block recurrence

Y_1^e = Φ_ext,       Y_{k+1}^e = Y_k^e Λ_e + Y_k^m Λ_me

where Λ_me = Λ[1:ROM, ROM+1:NVAR] and Λ_e = Λ[ROM+1:NVAR, ROM+1:NVAR] are the master↔external and external blocks of reduced_dynamics_linear. This is the one place an eigenvalue-derived matrix remains, confined to the per-order precompute.

All contractions are bilinear (Σᵢ J_r[·,i] Y[i,·]): the sesquilinear conjugation of the orthogonality condition is already baked into J_coeffs.

Arguments

• J_coeffs :: Vector{<:AbstractMatrix{T}} – output of precompute_orthogonality_operator_coefficients; J_coeffs[r] is ORD × FOM.

• right_master_blocks :: AbstractArray{T,3} – right master-mode order-blocks, size FOM × ORD × ROM; right_master_blocks[:, k, m] = Y_k^m[:, m] (equal to the linear master monomials of the parametrisation W). Only blocks k ≤ ORD-1 are used.

• external_directions :: AbstractMatrix{T} – physical external directions Φ_ext, size FOM × N_EXT.

• reduced_dynamics_linear :: AbstractMatrix{T} – NVAR × NVAR linear reduced dynamics; only the Λ_me and Λ_e blocks are read.

Return values

• C_coeffs :: Vector{Matrix{T}} of length ROM; C_coeffs[r] is (ORD-1) × ROM, row p = degree-(p-1) coefficient of C_r(s).

• E_coeffs :: Vector{Matrix{T}} of length ROM; E_coeffs[r] is (ORD-1) × N_EXT, row p = degree-(p-1) coefficient of E_r(s).

When ORD == 1 both matrices have zero rows and the operators are identically zero (the corresponding blocks are absent from the assembled system).

Complexity

• Time: O(ROM² · ORD² · FOM) for the contractions plus O(ORD · FOM · NVAR · N_EXT) for the external block recurrence.

• Storage: O(ROM · ORD · NVAR)

precompute_orthogonality_operator_coefficients(fom_matrices, left_eigenmodes,
											   left_modes_derivatives = nothing)
-> Vector{Matrix{T}}

Pre-compute the polynomial coefficient arrays for the orthogonality row operators J_r(s) directly from the left-eigenvector order-blocks. No eigenvalue is used.

Mathematical origin

The orthogonality condition for master mode r at a monomial with superharmonic s is the sesquilinear B-orthogonality of the companion left eigenvector φ_r = [φ_{r,1}; …; φ_{r,ORD}] (defined by φ_rᴴ (λ_r B − A) = 0) against the first-order state 𝒲:

φ_rᴴ B 𝒲 = 0,     B = blockdiag(I, …, I, B_ORD)

Reducing block-by-block, the row acting on the physical unknown W is

J_r(s) = Σ_{j=1}^{ORD} J_r[j, :] · s^(j-1)

with coefficients read straight off the (conjugated) eigenvector blocks:

J_r[j, :]   = conj(φ_{r,j})          j = 1, …, ORD-1
J_r[ORD, :] = conj(B_ORDᴴ φ_{r,ORD}) = conj(B_ORDᴴ ℓ_r)

The conjugation is the sesquilinear ᴴ of the condition, stored so that the assembled matrix row acts bilinearly on W (row · W = Σᵢ rowᵢ Wᵢ).

Arguments

• fom_matrices :: NTuple{ORD+1, <:AbstractMatrix{T}} – linear matrices of the full-order model; fom_matrices[k+1] corresponds to B_k (0-indexed).

• left_eigenmodes :: AbstractMatrix{T} – physical-space (highest-order) left eigenvector slice; left_eigenmodes[:, r] is ℓ_r = φ_{r,ORD} (length FOM).

• left_modes_derivatives :: Union{Nothing, AbstractArray{T,3}} – lower-order left eigenvector blocks, size FOM × (ORD-1) × ROM; left_modes_derivatives[:, j, r] = φ_{r,j}. Required when ORD > 1 (the eigensolver returns them; see solve_left). May be nothing for ORD == 1.

Return value

A Vector{Matrix{T}} of length ROM; entry r is the ORD × FOM matrix J_coeffs[r] whose row j stores the degree-(j-1) coefficient of J_r(s).

Complexity

• Time: one B_ORDᴴ · ℓ_r product per mode — O(ROM · FOM²) dense (O(ROM · nnz) sparse); the remaining blocks are copies.

• Storage: O(ROM · ORD · FOM)

compute_higher_derivative_coefficients!(
	param_coeff, red_coeff, external_dynamics, superharmonic, global_index,
	generalised_eigenmodes, lower_order_couplings
) -> nothing

Compute the higher time‑derivative coefficients W^(j+1)[α] for j = 1 … ORD-1 using the superharmonic recurrence

W^(j+1)[α] = s · W^(j)[α]  +  Φ_master · R[α]  +  Φ_ext · e_dyn  +  ξ[j]

where:

• s = superharmonic is the frequency ⟨λ, α⟩,

• Φ = generalised_eigenmodes (FOM × NVAR) collects the right eigenmodes,

• R[α] = red_coeff[:, global_index] (ROM‑vector) contains the master‑mode reduced‑dynamics coefficients at the current monomial (already solved),

• e_dyn = external_dynamics (N_EXT‑vector) contains the known external dynamics at the current monomial,

• ξ[j] = lower_order_couplings[j] (FOM‑vector) contains the coupling from lower‑order monomials at derivative order j.

Modifies param_coeff in‑place. Does nothing when ORD = 1 (no higher derivatives exist for a first‑order ODE).

Arguments

• param_coeff :: AbstractArray{T, 3} — shape FOM × ORD × L; the coefficient tensor of the parametrisation polynomial.

• red_coeff :: AbstractMatrix{T} — shape ROM × L; master‑mode reduced‑dynamics coefficients.

• external_dynamics :: AbstractVector{T} — length N_EXT; known external dynamics at the current monomial.

• superharmonic :: T — scalar s = ⟨λ, α⟩.

• global_index :: Int — monomial index into the last axis of param_coeff and the last axis of red_coeff.

• generalised_eigenmodes :: AbstractMatrix{T} — shape FOM × NVAR; right generalised eigenvectors (master modes in columns 1:ROM, external modes in ROM+1:NVAR).

• lower_order_couplings :: AbstractVector{<:AbstractVector{T}} — length ORD; element j is a length‑FOM vector ξ[j] produced by LowerOrderCouplings.compute_lower_order_couplings.

create_parametrisation_method_objects(mset::MultiindexSet{NVAR}, ORD::Int, FOM::Int, ROM::Int, external_system_size::Int, ::Type{T}=Complex)

Create a consistent pair of polynomials:

• W: a Parametrisation{ORD, NVAR, T} with zero coefficients,

• R: a ReducedDynamics{ROM, NVAR, T} with zero coefficients.

Both polynomials share the same multiindex set mset and element type T. The total number of variables NVAR must satisfy NVAR == ROM + external_system_size. FOM is the full‑order dimension (size of the state vector). It is not stored but used to initialise the coefficient vectors correctly.

Arguments

• mset: multiindex set for NVAR variables.

• ORD: native order of the full ODE (1 or 2).

• FOM: dimension of the full‑order state in its native order.

• ROM: dimension of the reduced state.

• external_system_size: number of forcing variables (default 0).

• T: element type.

function restrict_Parametrisation_to_degree(W::Parametrisation, max_degree::Int) -> Parametrisation

Returns a new Parametrisation that contains all the monomials of poly that are of degree lower or equal than max_degree.

function restrict_ReducedDynamics_to_degree(R::ReducedDynamics, max_degree::Int) -> ReducedDynamics

Returns a new ReducedDynamics that contains all the monomials of poly that are of degree lower or equal than max_degree.

validate_multiindex_set(mset, nvar, rom; conjugate_permutation = nothing)

Check a custom multiindex set against everything the cohomological solve assumes, and throw an ArgumentError naming the offending exponent on the first violation.

The clauses, and why each one matters:

1. nvar variables, matching ROM + N_EXT.

2. Minimum total degree ≥ 1 — the expansion is centred on the fixed point, so the constant monomial has no coefficient to solve for.

3. Every unit multiindex eᵢ — the linear part of the parametrisation is initialised from the eigenvectors, one column per unit multiindex.

4. Downward closed — the graded solve reads W[α - β + eᵢ] while working on α and factorises α = β₁ + … + β_d over members of mset. A missing divisor is not an error at run time: it is read as zero, silently corrupting the right-hand side.

5. With a conjugate_permutation: that it is an involutive permutation of 1:nvar mapping 1:rom into itself, and that mset is closed under it. A member whose partner is absent is solved directly rather than filled by conjugation, so the result loses the conjugate structure the permutation asserts.

Returns nothing. Cost is O(nvar · |mset| · log|mset|) — negligible beside a solve, but parametrise and solve_parametrisation both take a validate_mset = false escape hatch for callers that have already checked.

See is_downward_closed and is_conjugate_closed for the two closure predicates on their own.

_accumulate_global_entries!(Fe_bufs, ∇W_qp, entries_by_deg, model, element, q, dΩ)

Recursive type-stable dispatch over the degree-grouped entry tuple. The base case (empty tuple) is a no-op. Each recursive step processes the head degree group — Julia specialises a method per DEG so Val(DEG) in ntuple and the accumulate_qp! dispatch are resolved at compile time — then recurses on the tail.

_accumulate_split!(accum, scratch, t, sym, W, set, rem, deg, candidate_indices, args_ext)

Accumulate contributions from one (monomial, external-split) pair into accum, dispatching on the SymmetryType tag sym.

• FullyAsymmetric: calls t.f! directly into accum for each factorisation.

• FullySymmetric / GroupwiseSymmetric: calls t.f! into scratch, then axpy!(multiplier, scratch, accum) to apply the symmetry count.

_build_fem_cached_split(::Val{DEG}, ::Val{ME}, cs) -> FEMCachedSplit{DEG, ME}

Convert a CachedSplit to a FEMCachedSplit{DEG, ME} for the FEM-batched replay path. Deduplicates (derivative_order, W_col_idx) pairs across all entries and remaps each entry's factor slots to local indices into the unique_cols table. DEG is the internal degree and ME the external multiplicity, carried as a type parameter so _replay_fem_split! can append the external arguments with a compile-time-known count. Called once at cache-build time; produces zero allocations in the hot path.

_build_global_fem_split(model, fem_splits_l, fem_term_indices) -> FEMGlobalSplit

Build the O4 combined-element-loop bookkeeping for one monomial, merging all me=0 FEM terms into a single FEMGlobalSplit.

• Deduplicates (order, col) pairs across all me=0 FEM splits into a global table.

• Remaps each FEMFactorisationEntry's local indices to global table indices.

• Groups entries by degree into a typed Tuple for type-stable compile-time dispatch.

• Returns an empty FEMGlobalSplit{Tuple{}} when no me=0 FEM terms are present.

_collect_entries(sym, t, mset, rem, deg, cands) -> Vector{FactorisationEntry}

Call the appropriate factorisation function for cache construction, routing FullyAsymmetric to factorisations_asymmetric, FullySymmetric to factorisations_fully_symmetric, and GroupwiseSymmetric to factorisations_groupwise_symmetric.

_derivative_orders(t) → NTuple

Map each factor slot to its 1-based derivative index. Example: multiindex = (2, 1) → (1, 1, 2).

_orders_for_cache(sym, t, deg) -> Vector{Int}

Return the per-slot derivative order vector for storage in a CachedSplit.

• FullySymmetric: all deg slots share one derivative index.

• Other: read from _derivative_orders(t).

_replay_all_fem_splits!(result, model, W, global_split, global_∇W_qp, global_Fe_buffers)

Execute the O4 combined element loop: traverse the mesh once, scatter all globally-unique W columns, then accumulate contributions from ALL me=0 FEM terms at each quadrature point.

fem_reinit! and scatter_qp! are called at most once per unique (element, W-column) pair; accumulate_qp! dispatches per degree group via _accumulate_global_entries!.

_replay_fem_split!(result, t, W, fem_split, Fe, temp, external_arguments)

Replay one FEMCachedSplit{DEG, ME} using the FEM-batched element loop.

For each element: calls fem_reinit! once, scatters each unique (order, col) W column to qp-level field values via scatter_qp!, accumulates all qp contributions via accumulate_qp!, and assembles the element residual into result via assemble_element!. The qp gradient buffer ∇W_qp is obtained from the term via fem_qp_buffer(t) and is owned by the term, not allocated here.

accumulate_qp! receives DEG internal arguments followed by ME external ones, so the tuple has length t.deg — matching what evaluate_term!'s direct FEM overload (_eval_fem_term_direct!) passes. Before this, the external slots were omitted here and ext_count was applied as a bare scalar multiplier, so an me > 0 FEM term could not tell one external direction from another and disagreed with the direct path. The external arguments are hoisted out of the element loop, exactly as in _replay_split!.

_replay_global_fem!(result, model, W, cache, exp_index)

Function barrier: retrieves cache.global_fem_splits[exp_index] (stored as Any) and delegates to _replay_all_fem_splits!. The barrier restores type stability because cache.global_∇W_qp::Matrix{QP} is concretely typed at the call site even though the global_split triggers runtime dispatch on ENTRIES_TUPLE.

_replay_split!(result, scratch, temp, t, W, split, deg, external_arguments)

Replay one CachedSplit into result using precomputed factorisation bookkeeping.

• me = 0 (args_ext_indices empty): accumulates directly into result.

• me > 0: accumulates into temp, then axpy!(ext_count, temp, result).

Dispatches asymmetric/symmetric accumulation via split.is_asymmetric.

_replay_term!(result, t, W, exp_index, t_idx, cache)

Replay all cached splits for one nonlinear term t and monomial exp_index into result. Dispatches on the concrete term type:

• MultilinearMap: replays via _replay_split! (closure path).

• FEMMultilinearMap: me=0 splits are handled by _replay_global_fem! (O4 combined loop); only me>0 fallback splits are processed here via _replay_fem_split!.

accumulate_multilinear_term!(result, scratch, temp, t, parametrisation,
							  exp, candidate_indices, external_exp, external_arguments)

Add the contribution of nonlinear term t for exponent exp to result.

When t has no external (forcing) slots (me = 0) there is exactly one split, so we skip bounded_index_tuples and write directly into result, saving one fill! and one axpy! per term. For me > 0 each split is accumulated into temp first, then scaled into result via axpy!.

build_multilinear_terms_cache(model, parametrisation[, skip_bits]) → MultilinearTermsCache

Precompute all factorisation data for every monomial and term. Valid as long as the multiindex set is unchanged.

When skip_bits[l] is true the cache entry for monomial l is left empty. This is safe for monomials that will never be replayed (linear monomials, conjugate secondaries): those entries are guarded by the same skip_bits check in the solve loop.

compute_multilinear_terms(model, exp, parametrisation) → Vector

Return the sum of all nonlinear-term contributions for exponent exp.

Scratch buffers and shared data (external_arguments, candidate_indices, external_exp) are allocated once and reused across all terms.

compute_multilinear_terms(model, exp_index, parametrisation, cache) → Vector

Cached variant of compute_multilinear_terms: replays precomputed factorisation data instead of calling the factorisation routines.

exp_index is the 1-based index into parametrisation's multiindex set.

compute_multilinear_terms!(result, model, exp_index, parametrisation, cache) → nothing

In-place variant: zeros result then accumulates all nonlinear contributions into it. Uses cache.scratch_buffer, cache.temp_buffer, and cache.external_arguments so that no heap allocation occurs during the inner solve loop.

symmetry_type(t) -> SymmetryType

Classify a AbstractMultilinearMap as FullyAsymmetric, FullySymmetric, or GroupwiseSymmetric based on its multiindex field.

Setting fully_asymmetric = true on any term overrides the multiindex-based classification and always returns FullyAsymmetric.

_is_zero_coeff(coeff) -> Bool

Return true when the FOM × ORD coefficient slice coeff is entirely zero. Used to skip trivial contributions without entering the accumulation loop.

_sum_degree_one_terms!(accumulator, upper_bound, multiindex_dict,
					   red_coefficients, param_coefficients, unit_vectors)

Accumulate the |β| = 1 part of the lower-order coupling sum into accumulator.

For each unit-vector eⱼ with eⱼ ≤ upper_bound, and for each i < j, adds R[i, eⱼ] * (upper_bound - eⱼ + eᵢ)[i] * W[:, :, upper_bound - eⱼ + eᵢ] to accumulator. Only reads precomputed multiindex_dict — no allocations.

The i < j restriction is causality, not an optimisation. upper_bound - eⱼ + eᵢ has the same total degree as upper_bound, and GrLex breaks ties by descending lexicographic order, so that exponent precedes upper_bound — and is therefore already solved — exactly when i < j. Consequently this branch reads only the strictly upper triangle of the reduced linear dynamics Λ, and a strictly-lower entry of Λ would be dropped without trace. That is why the external system's linear matrix is required to be upper triangular; see the ExternalSystems module docstring.

_sum_higher_degree_terms!(accumulator, upper_bound, mset, multiindex_dict,
						  red_coefficients, param_coefficients,
						  candidate_idxs, unit_vectors)

Accumulate the |β| ≥ 2 part of the lower-order coupling sum into accumulator.

Iterates over candidate_idxs (precomputed as exponents componentwise ≤ upper_bound with total degree in [2, sum(upper_bound)-1]) and adds the same weighted R[i, β] * W[:, :, upper_bound - β + eᵢ] contributions as _sum_degree_one_terms!.

compute_lower_order_couplings(upper_bound, parametrisation, reduced_dynamics,
							  multiindex_dict, accumulator, candidate_idxs,
							  unit_vectors) -> accumulator

Performance path. Evaluate the lower-order coupling vector ξ[α] for monomial α = upper_bound and accumulate the result into the pre-zeroed accumulator.

All auxiliary data (multiindex_dict, accumulator, candidate_idxs, unit_vectors) must be pre-allocated by the caller (typically via LowerOrderResources) to avoid heap allocation in the inner solve loop. Returns accumulator for convenience.

compute_lower_order_couplings(upper_bound, parametrisation, reduced_dynamics) -> Vector{Vector{T}}

Allocating convenience wrapper. Builds the multiindex dictionary, accumulator, and candidate index list internally, then delegates to the performance path.

Use the 7-argument form inside performance-critical loops; this overload is intended for interactive use or testing.

_backend_name(backend) -> Symbol

Return the stable diagnostic name of a sparse backend (:klu, :umfpack or :pardiso).

_backward_error(state, x, norm_r, norm_b, norm_A) -> Real

Return the normwise backward error norm_r / (norm_A * norm(x, Inf) + norm_b), guarded against a zero denominator in the state's real scalar type.

_bordered_solve!(state, solution, superharmonic, index, multiindex,
	resonance, resonance_set; reuse_factor = Val(false)) -> solution

Solve the bordered system in place through KLU, UMFPACK, or Pardiso. Exact grouped reuse skips numeric factorisation only when the caller proves that the assembled matrix is identical. Monomial data is carried solely to produce contextual failure diagnostics.

_cached_klu_factor(state, matrix) -> KLUFactorization

Return the cached KLU factorisation with value and index types derived from the active sparse matrix. Call only after the first successful factorisation has populated fact.

Return the selected backend's cached factor with its concrete type restored.

_cached_umfpack_factor(state, matrix) -> UmfpackLU

Return the cached UMFPACK factorisation with value and index types derived from the active sparse matrix. Call only after a successful factorisation has populated fact.

_configured_residual_tolerance(T, options) -> tolerance or nothing

Resolve the backward-error threshold in the real scalar type associated with T. options = nothing selects the standard checked default used by the low-level constructor.

_discard_klu_factor!(state, factorisation) -> nothing

Remove a failed KLU factorisation from state and release its native symbolic and numeric storage before another analysis begins. Merely assigning nothing to state.fact would leave that storage alive until garbage collection, which can retain damaged KLU state across the immediate recovery attempt on Windows.

Base.finalize is Julia's required external spelling; it runs the finaliser already owned by KLUFactorization and prevents it from running twice later.

Create and cache a checked fresh KLU factorisation of A.

Return whether a caught runtime failure must pass through without solver wrapping.

Return configured outer-resonance target indices for one monomial position.

Run Pardiso numeric factorisation followed by a solve.

Run Pardiso configuration and symbolic analysis for a bordered matrix.

Release extension-owned Pardiso factorisation storage.

Solve with the current Pardiso numeric factorisation.

Backward-compatible internal alias for _refactorise_klu!.

_refactorise_klu!(ss, A) -> KLU factorisation of `A`

Factorise the bordered matrix for the current monomial, reusing the symbolic analysis cached in ss.fact.

The first call analyses and factorises; every later call redoes only the numeric factorisation. Splitting them this way is what makes the constant-size formulation pay off: the sparsity pattern is identical for every monomial, so the ordering and symbolic phase — the expensive part — are computed once for the whole solve.

Why klu_factor! and not klu!

The numeric phase must re-pivot. The bordered matrix changes value on every monomial and its (1,1) block is near-singular at every resonance, where stability depends on pivoting being free to exchange rows across the border. klu_factor! re-pivots while reusing the cached symbolic analysis. KLU's exported klu! maps to klu_refactor, which replays the pivot sequence chosen at the first monomial and so loses accuracy exactly where it is needed.

Preconditions and aliasing

A's sparsity pattern must be identical on every call — SparseLinearSolverState guarantees this by construction, since only nzval is ever written.

klu takes A.nzval by reference rather than copying it, so the factorisation and the template share one value array: per-monomial assembly writes straight into what KLU reads, with no copy-in step. colptr/rowval stay KLU's own (0-based) copies, and are invariant in any case.

A factorisation is cached only if it succeeded. Both fresh and cached numeric factorisations use KLU's checked path because stored zeros in the constant border can make an unchecked status appear successful even when the numeric factor is singular. A failed cached factor has its native storage released immediately and is retried once from a fresh symbolic analysis.

Dispatch numeric refactorisation to the selected built-in sparse backend.

_refactorise_umfpack!(state, matrix) -> UMFPACK factorisation

Factorise the current bordered matrix through SuiteSparse UMFPACK. Successful symbolic analysis is reused by lu!. A failed cached numeric factorisation is discarded and retried once from a fresh analysis; only a successful factorisation is retained in state.fact.

_release_pardiso!(state) -> nothing

Finaliser: hand the Pardiso factorisation back. No-op on the KLU and UMFPACK paths, and never allowed to throw — a finaliser that raises would be reported out of context.

_sparse_backward_error!(state, factorisation, solution, norm_b) -> Real

Check the normwise backward error of a KLU or UMFPACK solution and, when necessary, perform at most state.max_refinement_steps correction solves with the supplied typed factorisation. Update the diagnostic maximum and refinement count and return the final backward error.

_sparse_inf_norm!(workspace, matrix) -> Real

Compute norm(matrix, Inf) for a CSC matrix using workspace for row sums and without allocating a dense intermediate.

_sparse_residual!(residual, matrix, solution) -> residual

Overwrite a vector containing b with matrix * solution - b using SparseArrays' five-argument mul! kernel.

Shared allocation-conscious solve for the built-in KLU and UMFPACK backends.

Construct and throw a contextual _BorderedSolveFailure.

_try_build_pardiso_solver() -> solver or nothing

Extension hook implemented by MORFEPardisoExt. Return a configured Pardiso solver when the extension is active, or nothing so automatic selection can fall back to KLU.

_assemble_bordered_system!(ctx, s, resonance, lower_order_couplings, external_dynamics)

Assemble the (FOM+ROM) × (FOM+ROM) bordered cohomological system into ctx.buffers.system_matrix and ctx.buffers.rhs.

• The first FOM rows come from the invariance equation (operator + nonlinear RHS).

• The last ROM rows come from the orthogonality conditions, with non-resonant modes contributing the trivial row R[r, α] = 0.

Called by the dense-path _solve_monomial!.

_assemble_nonlinear_rhs!(instrumentation, ctx, model, idx, W, ml_cache)

Instrumentation hook for nonlinear right-hand-side assembly. The production method returns nothing; benchmark instrumentation returns the corresponding @timed result.

_dense_backward_error(A, x, b) -> Real

Compute norm(A*x-b, Inf) / (norm(A, Inf)*norm(x, Inf) + norm(b, Inf)) without allocating a residual vector. The denominator is floored at the real scalar type's smallest positive normal value.

_finalise_monomial!(W, R, idx, ctx, s, resonance, external_dynamics,
	lower_order_couplings) -> nothing

Unpack the bordered solution into the primary coefficients of W and the resonant rows of R, write exact zeros to non-resonant master rows, and propagate the higher derivative coefficients using the generalised right eigenmodes.

_monomial_metrics(instrumentation, rhs_result, solve_result)

Convert instrumentation-specific phase results into the metrics delivered to solve observers. The production path returns nothing.

_resonance_vector(resonance_set, monomial_idx, ::Val{ROM}) -> SVector{ROM, Bool}

Return a compile-time-sized boolean vector indicating which of the ROM master modes are resonant with the monomial at position monomial_idx in the multiindex set. Using Val{ROM} enables the compiler to emit a fully unrolled ntuple loop.