Fast Solvers

Implicit and semi-implicit time integrators (backward Euler, Crank-Nicolson, SICN, CNAB) all solve a linear system of the form

\[\left(\mathbf{I} - \alpha\mathbf{A}\right)\mathbf{u}(t_{k+1}) = \mathbf{r}_k\]

at every time step, where $\mathbf{A}$ is the linear operator of the discretized PDE and $\alpha$ depends on the scheme:

Scheme	$\alpha$
Backward Euler (`:BE`)	$\Delta t$
Crank-Nicolson / SICN / CNAB (`:CN`)	$\Delta t / 2$

Crank-Nicolson-type schemes also need the forward product $\left(\mathbf{I} + \alpha\mathbf{A}\right)\mathbf{u}(t_k)$. Because $\mathbf{A}$ and $\alpha$ are fixed during a simulation, the FastSolvers submodule pre-processes $\mathbf{A}$ once (eigendecomposition, FFT diagonalization, or LU factorization) and bakes $\alpha$ into a cached factor. Each time step then reduces to a few structured multiplies with no allocations.

Usage with the models

Every model provides a build_fast_solver method that picks the right solver for its boundary conditions, and an integrate_model_fast method that uses it. The scheme keyword sets $\alpha$ and must match the integrator you use.

using PolynomialModelReductionDataset
const pomoreda = PolynomialModelReductionDataset

heat1d = pomoreda.Heat1DModel(
    spatial_domain=(0.0, 1.0), time_domain=(0.0, 1.0),
    Δx=1/2^7, Δt=1e-3, diffusion_coeffs=0.1
)
μ = heat1d.diffusion_coeffs
A, B = heat1d.finite_diff_model(heat1d, μ; same_on_both_ends=true)
Ubc = ones(heat1d.time_dim)           # boundary input

solver = pomoreda.build_fast_solver(heat1d, μ; scheme=:CN, same_on_both_ends=true)   # one-time precompute
U = pomoreda.integrate_model_fast(heat1d, solver, heat1d.tspan, heat1d.IC, Ubc;
                                  control_matrix=B, scheme=:CN)

See each model's page for the keyword arguments its integrate_model_fast method accepts.

Using the solvers directly

The solvers can also be built from an operator matrix and used in your own time-stepping loop, for example to integrate a reduced-order model.

solver = FastSymTridiagSolver(A, Δt/2)   # α = Δt/2 for Crank-Nicolson

linsolve!(u_new, solver, rhs)            # solve (I - α A) u_new = rhs
mulIpA!(out, solver, u_old)              # out = (I + α A) u_old
update_timestep!(solver, Δt_new/2)       # change α without rebuilding

Solver	Structure of $\mathbf{A}$	Typical use
`FastSymTridiagSolver`	1D symmetric tridiagonal	1D Dirichlet, Neumann, or mixed boundaries
`FastCirculant1DSolver`	1D circulant	1D periodic boundaries
`FastKronSumSolver`	$(\mathbf{A}_y \otimes \mathbf{I}) + (\mathbf{I} \otimes \mathbf{A}_x)$, symmetric tridiagonal factors	2D Dirichlet boundaries
`FastFFT2DSolver`	Kronecker sum of two circulants	2D periodic boundaries
`FastDenseSolver`	dense, unstructured	reduced-order models
`FactorizedSolver`	any (sparse) matrix	fallback when no structure applies

API

PolynomialModelReductionDataset.FastSolvers — Module

Fast implicit linear solvers for PDE time integration.

This submodule provides reusable, allocation-free solvers for systems of the form (I - α A) u_new = rhs (and the companion (I + α A) x product) that arise repeatedly in implicit / semi-implicit time stepping. The factor α is baked into the solver at construction time:

Backward Euler : α = Δt
Crank-Nicolson : α = Δt / 2

A typical workflow is therefore

solver = FastSymTridiagSolver(A, Δt/2)        # CN-flavoured
linsolve!(u_new, solver, rhs)                  # (I - α A) u_new = rhs
mulIpA!(out,    solver, u_old)                 # out = (I + α A) u_old
update_timestep!(solver, Δt_new/2)             # rebuild cached factor in O(N)

Available concrete solver types and the structure they exploit:

Type	Structure of A	Per-step cost
`FastSymTridiagSolver`	1D symmetric tridiagonal (Dirichlet/Neumann)	O(N²) (BLAS-2)
`FastCirculant1DSolver`	1D circulant (periodic)	O(N log N) (FFT)
`FastKronSumSolver`	2D `(Ay ⊗ I) + (I ⊗ Ax)` with sym-tridiag	O(Nx Ny (Nx+Ny))
`FastFFT2DSolver`	2D circulant Kronecker sum (periodic, periodic)	O(N log N) (FFT)
`FastDenseSolver`	dense unstructured (e.g. reduced-order)	O(r²) (BLAS-2)
`FactorizedSolver`	any sparse A (LU factorization once)	O(nnz) per solve

source

PolynomialModelReductionDataset.FastSolvers.AbstractFastSolver — Type

abstract type AbstractFastSolver

Abstract supertype for all fast linear solvers that pre-process A once so that each time step reduces to a structured solve / multiply.

source

PolynomialModelReductionDataset.FastSolvers.FastSymTridiagSolver — Type

struct FastSymTridiagSolver{T<:Real} <: AbstractFastSolver

Fast implicit solver for a 1-D symmetric tridiagonal operator A. Eigen- decomposes A = V Λ Vᵀ once and caches both 1/(1 − α λ_i) and 1 + α λ_i so that each solve / multiply is a pair of matrix-vector products plus an elementwise scaling.

Fields

V::Matrix{T} where T<:Real: Eigenvectors of A (real, N×N)
Vt::Matrix{T} where T<:Real: Transpose of V, materialized for BLAS
λ::Vector{T} where T<:Real: Eigenvalues of A (length N)
inv_denom::Vector{T} where T<:Real: 1 / (1 - α λ_i) — cached inverse
mul_factor::Vector{T} where T<:Real: 1 + α λ_i — cached forward (CN) factor
N::Int64: Spatial dimension
tmp::Vector{T} where T<:Real: Workspace (length N)

source

PolynomialModelReductionDataset.FastSolvers.FastCirculant1DSolver — Type

struct FastCirculant1DSolver{P, IP} <: AbstractFastSolver

Fast implicit solver for a 1-D circulant operator A (periodic BCs). The eigenvalues are λ = fft(first column of A); the time loop reduces to a pair of in-place FFTs plus an elementwise scaling.

Fields

plan_f::Any: In-place forward FFT plan
plan_if::Any: In-place inverse FFT plan
λ::Vector{ComplexF64}: Eigenvalues of A (complex, length N)
inv_denom::Vector{ComplexF64}: 1 / (1 - α λ) — cached inverse
mul_factor::Vector{ComplexF64}: 1 + α λ — cached forward (CN) factor
N::Int64: Spatial dimension
buffer::Vector{ComplexF64}: Complex workspace of length N (FFT scratchpad)

source

PolynomialModelReductionDataset.FastSolvers.FastKronSumSolver — Type

struct FastKronSumSolver{T<:Real} <: AbstractFastSolver

Fast implicit solver for a 2-D operator with Kronecker-sum structure

A = (Ay ⊗ I_Nx) + (I_Ny ⊗ Ax)

where both 1-D factors are symmetric tridiagonal. Eigendecomposing Ax and Ay once turns each solve / multiply into four dense matrix-matrix products (BLAS-3, very cache friendly) plus an elementwise scaling. Total per-step cost is O(Nx Ny (Nx + Ny)).

Fields

Vx::Matrix{T} where T<:Real: Eigenvectors of Ax
Vy::Matrix{T} where T<:Real: Eigenvectors of Ay
Vxt::Matrix{T} where T<:Real: Vxᵀ materialized for BLAS
Vyt::Matrix{T} where T<:Real: Vyᵀ materialized for BLAS
λx::Vector{T} where T<:Real: Eigenvalues of Ax (length Nx)
λy::Vector{T} where T<:Real: Eigenvalues of Ay (length Ny)
inv_denom::Matrix{T} where T<:Real: 1 ./ (1 .- α .* (λx .+ λy')) — cached inverse
mul_factor::Matrix{T} where T<:Real: 1 .+ α .* (λx .+ λy') — cached forward (CN) factor
Nx::Int64: x-dimension
Ny::Int64: y-dimension
tmp1::Matrix{T} where T<:Real: Workspace 1 (Nx × Ny)
tmp2::Matrix{T} where T<:Real: Workspace 2 (Nx × Ny)

source

PolynomialModelReductionDataset.FastSolvers.FastFFT2DSolver — Type

struct FastFFT2DSolver{P, IP} <: AbstractFastSolver

Fast implicit solver for a 2-D doubly-periodic operator (Kronecker sum of two 1-D circulants). Each step is a single in-place 2-D FFT, an elementwise scaling, and an inverse FFT.

Fields

plan_f::Any: In-place 2-D forward FFT plan
plan_if::Any: In-place 2-D inverse FFT plan
λx::Vector{ComplexF64}: 1-D eigenvalues in x
λy::Vector{ComplexF64}: 1-D eigenvalues in y
inv_denom::Matrix{ComplexF64}: 1 ./ (1 .- α .* (λx .+ λy')) — cached inverse
mul_factor::Matrix{ComplexF64}: 1 .+ α .* (λx .+ λy') — cached forward (CN) factor
Nx::Int64: x-dimension
Ny::Int64: y-dimension
buffer::Matrix{ComplexF64}: Complex (Nx × Ny) workspace

source

PolynomialModelReductionDataset.FastSolvers.FastDenseSolver — Type

struct FastDenseSolver <: AbstractFastSolver

Fast implicit solver for a dense, unstructured matrix A (typically from a reduced-order model). Eigendecomposes A once and pre-computes the full inverse (I - α A)⁻¹ (and I + α A for CN-style steps) so that each step is a single dense matrix-vector multiply.

If A is nearly defective (cond(V) > cond_threshold), falls back to a direct LU-based inverse for robustness; a warning is issued.

Fields

M_inv::Matrix{Float64}: (I - α A)⁻¹, real r×r matrix applied via mul! each step
IpA::Matrix{Float64}: I + α A, real r×r matrix used by mulIpA!
V::Matrix{ComplexF64}: Eigenvectors of A (complex, stored for update_timestep!)
Vinv::Matrix{ComplexF64}: Inverse of V (complex)
λ::Vector{ComplexF64}: Eigenvalues of A (complex)
A::Matrix{Float64}: Stored copy of A (real, for update_timestep! rebuilding IpA)
r::Int64: Dimension of the system
eigen_based::Bool: Whether the solver was constructed via eigendecomposition (false = LU fallback)

source

PolynomialModelReductionDataset.FastSolvers.FactorizedSolver — Type

mutable struct FactorizedSolver{TF, TA<:(AbstractMatrix)} <: AbstractFastSolver

Fallback fast solver for arbitrary A. Forms M = I - α A once and stores its LU factorization; each linsolve! is a single triangular back-solve. Use this when none of the structured solvers apply (e.g. asymmetric tridiag, matrices with Δt baked into boundary rows).

Fields

fact::Any: LU factorization of (I - α A)
A::AbstractMatrix: Stored A (used for mulIpA! and update_timestep!)
α::Float64: Baked-in α
N::Int64: System dimension

source

PolynomialModelReductionDataset.FastSolvers.linsolve! — Function

linsolve!(out, solver::AbstractFastSolver, rhs)

Solve (I - α A) out = rhs in place using the factor cached in solver, where α is the value baked in at construction (or by update_timestep!). Returns out.

source

PolynomialModelReductionDataset.FastSolvers.mulIpA! — Function

mulIpA!(out, solver::AbstractFastSolver, x)

Compute the forward product out = (I + α A) x in place, as needed by the explicit half of Crank-Nicolson-type schemes. Returns out.

source

PolynomialModelReductionDataset.FastSolvers.update_timestep! — Function

update_timestep!(solver::AbstractFastSolver, α_new)

Replace the baked-in α with α_new (e.g. after a change of time step) and rebuild the cached factors. The structured solvers only rescale cached eigenvalue factors; FactorizedSolver refactorizes I - α_new A. Returns solver.

source

PolynomialModelReductionDataset.FastSolvers.backward_euler_solve! — Function

backward_euler_solve!(unew, F, rhs)

Backward compatibility alias for linsolve!. Identical behaviour: it applies the cached factor of (I - α A)⁻¹ baked into the solver. Prefer linsolve! in new code.

source