Library

AbstractGlaOpr{T<:AbstractFloat}

Abstract base type for all operators in the Gila package. All concrete operator types must subtype this type. The operators are matrix free and are not AbstractArrays: they carry size, mul! and \ rather than an element at a time. T is the real storage precision of the operator, so the operator acts on Complex{T} data.

AbstractGlaVacOpr{T}

Abstract base type for the vacuum Green operators, the ones built from geometry alone with no material in them. Scattering operators are not vacuum operators; they subtype AbstractGlaOpr directly.

Traits shared by every vacuum operator are defined once on this type, so a subtype only overrides the ones where it differs.

Interface

A subtype has to implement:

• Base.size(opr): The number of rows and columns of the operator as a matrix
• mulAct!(opr, act::AbstractVector{Complex{T}}): The matrix-vector product on a flat vector (mutating)
• adjoint!(opr): The adjoint, with the value semantics of its docstring
• KernelAbstractions.get_backend(opr): The backend the operator computes on
• isadjoint(opr): Whether the operator is the adjoint of the Green operator
• isselfoperator(opr): Whether the source and target geometry are the same
• isexternaloperator(opr): Whether the source and target geometry are disjoint
• glaSze(opr): The target and source sizes in tensor form

isoverlappingoperator defaults to false here, eltype to Complex{T} on AbstractGlaOpr, and slv to a GCRODRSolver. Override any of the three when it does not hold.

BiCGStabSolver(maxItr = nothing, absTol = nothing, relTol = nothing; preCon = nothing)

Create a BiCGStabSolver. Unset values are resolved against the right hand side at solve time.

Returns

• BiCGStabSolver: A new solver instance

BiCGStabSolver

An iterative solver for linear systems of equations that uses the BiConjugate Gradient Stabilized Method (BiCGStab). This method is effective for non-symmetric linear systems and typically requires less memory than GMRES.

Fields

• maxItr::Union{Nothing, Int}: Maximum number of iterations (default: length(vec))
• absTol::Union{Nothing, Real}: Absolute tolerance for convergence (default: 0)
• relTol::Union{Nothing, Real}: Relative tolerance for convergence (default: √ε)
• preCon: Right preconditioner, anything with an ldiv!(out, preCon, inp) method (default: none)

The preconditioner is applied on the right: the Krylov space is that of opr * preCon⁻¹, and the residual driving the stopping test is inp - opr * out itself, not a preconditioned version of it.

GCRODRSolver(rstItr = nothing, maxItr = nothing, absTol = nothing, relTol = nothing; preCon = nothing, side = :left, rcyDim = 10, lrgFrc = 0.0)

Create a GCRODRSolver. Unset values are resolved against the right hand side at solve time.

Returns

• GCRODRSolver: A new solver instance

GCRODRSolver

An iterative solver for linear systems of equations that uses GCRO with deflated restarting (GCRO-DR). GCRO-DR keeps rcyDim directions of the Krylov space built in GMRES and reuses them in the following cycle and, through the rcy argument of solve, in a subsequent solve.

Fields

• rstItr::Union{Nothing, Int}: Iterations in a cycle, rcyDim included (default: min(20, length(vec)))
• maxItr::Union{Nothing, Int}: Maximum number of iterations (default: max(5000, length(vec)))
• absTol::Union{Nothing, Real}: Absolute tolerance for convergence (default: 0)
• relTol::Union{Nothing, Real}: Relative tolerance for convergence (default: √ε)
• preCon: Preconditioner, anything with an ldiv!(out, preCon, inp) method (default: none)
• side::Symbol: :left (default) or :right, the side preCon is applied on
• rcyDim::Int: Directions carried between cycles, 0 gives plain GMRES
• lrgFrc::Float64: Fraction of rcyDim taken from the large end of the spectrum

The recycled pair (U, C) satisfies A * U = C with C' * C = I, A being the preconditioned operator. The best correction inside range(U) therefore costs no operator applications, and the cycle's least squares problem decouples into the plain GMRES one.

lrgFrc * rcyDim of the directions approximate the eigenvectors of largest modulus (through ordinary Ritz values) and the rest those of smallest modulus (through harmonic Ritz values). Deflating the small end (making lrgFrc closer to 0) changes the rate of convergence.

There is no flexible field: flexible GMRES stores the preconditioned basis, which does not satisfy A * U = C.

GMRESSolver(rstItr = nothing, maxItr = nothing, absTol = nothing, relTol = nothing; preCon = nothing, flexible = false, side = flexible ? :right : :left)

Create a GMRESSolver. Unset values are resolved against the right hand side at solve time.

Returns

• GMRESSolver: A new solver instance

GMRESSolver

An iterative solver for linear systems of equations that uses the Generalized Minimal Residual Method (GMRES). This method is particularly effective for non-symmetric linear systems.

Fields

• rstItr::Union{Nothing, Int}: Number of iterations until restart (default: min(20, length(vec)))
• maxItr::Union{Nothing, Int}: Maximum number of iterations (default: max(5000, length(vec)))
• absTol::Union{Nothing, Real}: Absolute tolerance for convergence (default: 0)
• relTol::Union{Nothing, Real}: Relative tolerance for convergence (default: √ε)
• preCon: Preconditioner, anything with an ldiv!(out, preCon, inp) method (default: none)
• side::Symbol: :left (default) or :right, the side preCon is applied on
• flexible::Bool: Whether to store the preconditioned basis, admitting a varying preCon (default: false)

On the left the Krylov space is that of preCon⁻¹ * opr and the quantity driving the stopping test is preCon⁻¹ * (inp - opr * out), so two preconditioners reporting the same relative residual have not converged to the same thing: the discrepancy is the size of cond(preCon). On the right the space is that of opr * preCon⁻¹, the residual is inp - opr * out itself, and the price is one extra ldiv! and one extra pass over the solution per restart cycle.

Flexible GMRES keeps Z = [z₁ … z_k] with zₖ = preConₖ⁻¹ * vₖ and updates out += Z * y, so the applications never have to agree with each other. It implies side = :right, and it costs rstItr extra full length vectors.

GlaCmpVol

A set of disjoint GlaVol regions (potentially with different scales, and potentially separated) made to "feel" like a single volume.

Fields

• regions::Vector{GlaVol}: The regions of the tiling, in the order that fixes the flat degree of freedom layout

Flat degree of freedom layout

Fields over a composite volume are stored in one flat vector. Region r contributes a block of 3 * prod(r.cel) entries and the blocks are concatenated in regions order. Inside a block the ordering is vec of an array of size (r.cel..., 3), so the first Cartesian index runs fastest and the vector component index runs slowest. This matches the reshape that GlaOprVac's * applies to a flat vector, which means a region block can be reshaped into the 4-tensor an operator expects without any permutation.

coordinates walks the cells in this order but yields one entry per cell rather than per degree of freedom, so it is three times shorter than the flat vector. cellvolumes repeats each cell volume across the three vector components, so it lines up with the degrees of freedom entry by entry.

Invariants

The regions must be nonempty, have even cell counts, have pairwise commensurate scales, sit on a common grid, have disjoint interiors (face, edge, and corner contact is fine and expected), and satisfy the per-partition parity condition of chkParExtInf for every pair. Gaps between regions are allowed: a tiling may hold several separated bodies.

GlaCmpVol(vol::GlaVol)

Construct the one-region composite volume around a single GlaVol.

Arguments

• vol::GlaVol: The volume to wrap

Returns

• GlaCmpVol: A composite volume with vol as its only region

GlaFld{T}

A field (in the physics sense, i.e., a vector at each point in space) over a volume, stored as a flat buffer.

The volume is always a GlaCmpVol. A plain GlaVol is wrapped into a tiling of one region, so a field over an ordinary volume is the same type as a field over a refined one.

The buffer holds the density coefficient of each degree of freedom, the cell average of the current. Only dot and norm carry the cell volume ΔV of the region a degree of freedom belongs to: dot is the ΔV-weighted pairing and norm the ΔV-weighted L^p norm, with norm(fld, Inf) the plain maximum.

The layout inside the buffer is the flat degree of freedom layout of GlaCmpVol: region blocks in regions order, and inside a block the vec of an array of size (reg.cel..., 3).

Fields

• dat::AbstractVector{Complex{T}}: The flat buffer, on any KernelAbstractions backend, whose eltype fixes the storage precision T
• cvol::GlaCmpVol: The composite volume the field lives on
• off::Vector{Int}: Buffer offset of each region block, with the total length as a last entry

GlaOpr{T, S}

The full Green operator G₀(I - XG₀)⁻¹.

GlaOpr(opr::Union{GlaOprVac, GlaCmpOprVac}, sus; slv::GlaSlv=GCRODRSolver())

The full Green operator over an existing vacuum operator, which it shares. sus as for SctOpr over the same operator.

GlaOpr{T}(opr::GlaOpr)

Convert the storage precision of the operator, and of its susceptibility, to T. The solver is kept.

GlaOpr{T}(vol::Union{GlaVol, GlaCmpVol}, sus; dev, slv=GCRODRSolver(), kwargs...)

The full Green operator G₀(I - XG₀)⁻¹. Arguments as for SctOpr.

GlaOprVac{T}

The vacuum Green operator G₀, acting on Complex{T} data. srcMsk and trgMsk pick the two volumes out of their union when they overlap, and are nothing otherwise.

GlaOprVac(mem::GlaVacOprMem)

Wrap prepared operator memory, unmasked.

GlaOprVac{T}(cvol::GlaCmpVol; kwargs...)
GlaOprVac{T}(trgCvl::GlaCmpVol, srcCvl::GlaCmpVol; kwargs...)

Same as GlaCmpOprVac: a composite volume needs a block matrix rather than one circulant, so the result is a GlaCmpOprVac.

GlaOprVac{T}(opr::GlaOprVac)

Convert the storage precision of a vacuum Green function operator to T, without regenerating its Fourier coefficients.

GlaOprVac{T}(trgVol::GlaVol, srcVol::GlaVol; kwargs...)
GlaOprVac{T}(vol::GlaVol; kwargs...)

The vacuum Green operator from srcVol to trgVol, or of vol on itself, stored in Complex{T} (T defaults to dflPrc). Volumes that share interior are built on their union and masked, so they must sit on a common grid.

Keywords

• dev::Backend=CPU(): backend to compute on
• frqPhz=1.0+0.0im: complex frequency phase factor
• qssApx::Bool=false: use the quasistatic kernel
• shpCch::Bool=false: cache the far-field geometry table of each cell shape, see GlaVacOprMem. The in-memory caches of the far-field and contact fills are never evicted: a process keeps 20-26 MB per cell shape and one far-field set per shape and frequency of a sweep. The unexported GilaElectromagnetics.clrCch!() empties them
• prxWrn::Bool=true: warn when two separated volumes are closer than a third of a wavelength
• hrmPrt=nothing: :asy or :sym keeps only the anti-Hermitian or Hermitian part of a self operator, taken from the kernel before its transform; see asym and sym
• blc=nothing: Bloch lattice (ax = (1, 2), P = (32, 32), θ = (0.7, 1//2)), self operator only: vol is the unit cell, periodic along the axes ax with periods P cells and Bloch phases θ (radians, or exact twists θ/2π as Rationals), G(R + Pᵢ eᵢ) = e^{-iθᵢ} G(R); see setBlc!

GlaSlv

Abstract base type for all solvers in the Gila package. All concrete solver types must subtype this type.

GlaVol(cel::Array{<:Integer,1}, celScl::NTuple{3,Rational}, org::NTuple{3,Rational}, grdScl::NTuple{3,Rational}=celScl)

Constructor for Gila Volumes.

Arguments

• cel::Array{<:Integer,1}: Array of the number of cells in each dimension of the volume
• celScl::NTuple{3,Rational}: The size of each cell in each dimensions of the volume (in units of wavelength)
• org::NTuple{3,Rational}=(0//1, 0//1, 0//1): The origin of the volume in each dimension (in units of wavelength)
• grdScl::NTuple{3,Rational}=celScl: Spatial location of the center of each cell contained in the volume

Returns

• GlaVol: A new volume with the specified parameters

Throws

• ErrorException: If the cell scale is larger than the grid scale

GlaVol

Basic spatial memory structure for a rectangular prism volume.

Fields

• cel::NTuple{3,Int}: Number of cells in each dimension of the volume
• scl::NTuple{3,Rational{Int}}: Relative side length of a cuboid voxel (cell) compared to the wavelength
• org::NTuple{3,Rational{Int}}: Center position of the domain
• grd::NTuple{3,StepRange{Rational{Int},Rational{Int}}}: Spatial location of the center of each cell contained in the volume

InvSctOpr{T, O}

The inverse scattering operator (I - XG₀), with X over the source tiling of G₀.

InvSctOpr{T}(opr::InvSctOpr)

Convert the storage precision of the operator, and of its susceptibility, to T.

InvSctOpr{T}(vol::Union{GlaVol, GlaCmpVol}, sus; dev, kwargs...)

The inverse scattering operator (I - XG₀) of vol on itself, with sus in any form SusOpr takes and dev defaulting as there. Other keywords go to GlaOprVac, or to GlaCmpOprVac for a composite volume.

MixPrcRfn{Tlo}

Mixed precision iterative refinement (GMRES-IR): the residual and the solution update are formed in the precision of the operator, while the correction to the solution is solved for on a Tlo copy of that operator. One high precision matrix-vector product per outer step buys a high precision backward error at the matrix-vector cost of Tlo.

The operator passed to solve must be of higher precision than Tlo, and must match the precision of the right hand side.

Typical usage: GlaOpr{Float64}(vol, sus; slv=MixPrcRfn(Float32))

Fields

• innSlv::GlaSlv: The solver used for the low precision corrections
• cch::Base.RefValue{Any}: The (opr, oprLo) pair of the last solve, rebuilt when the operator changes
• maxItr::Union{Nothing, Int}: Maximum number of outer refinement steps (default: 20)
• absTol::Union{Nothing, Real}: Absolute tolerance for convergence (default: 0)
• relTol::Union{Nothing, Real}: Relative tolerance for convergence (default: √ε)

MixPrcRfn(::Type{Tlo}; innSlv::GlaSlv=GMRESSolver(), maxItr=nothing, absTol=nothing, relTol=nothing)

Create a MixPrcRfn refining in the precision Tlo. Unset tolerances are resolved against the right hand side at solve time.

Arguments

• Tlo::Type{<:AbstractFloat}: The precision of the correction solves
• innSlv::GlaSlv=GMRESSolver(): The solver for the correction solves
• maxItr, absTol, relTol: Outer loop settings, resolved at solve time when left unset

Returns

• MixPrcRfn{Tlo}: A new solver instance holding no operator copy yet

RcySpc()

The subspace a GCRODRSolver recycles, filled in place when passed as solve(opr, inp, slv; rcy = RcySpc()), and read by the next solve it is passed to.

Fields

• bas: U, the basis of the subspace, or nothing before a solve fills it
• img: C, equal to A * U with orthonormal columns, A being the preconditioned operator

SctOpr{T, I, S}

The scattering operator (I - XG₀)⁻¹, applied by solving with slv.

SctOpr{T}(opr::SctOpr)

Convert the storage precision of the operator, and of its susceptibility, to T. The solver is kept.

SctOpr(opr::Union{GlaOprVac, GlaCmpOprVac}, sus; slv::GlaSlv=GCRODRSolver())

The scattering operator over an existing vacuum operator, which it shares. A composite operator has to be a self operator.

SctOpr{T}(vol::Union{GlaVol, GlaCmpVol}, sus; dev, slv=GCRODRSolver(), kwargs...)

The scattering operator (I - XG₀)⁻¹, applied by solving with slv. Other arguments as for InvSctOpr.

SlvLog(; truSmp = :atrestart, hssSmp = :all)

A record of one solve, filled in place when passed as solve(opr, inp, slv; log = SlvLog()).

Arguments

• truSmp::Symbol: :atrestart (default) samples resTru at every GMRES restart, :never only at exit
• hssSmp::Symbol: :all (default) keeps every Hessenberg cycle, :first only the first

Fields

• status::Symbol: :converged, :maxiter, :breakdown, or :none before a solve fills it
• numItr::Int: Iterations taken
• resRec::Vector{Float64}: The relative recursive residual the algorithm tracks, before the first step and after each one
• resTru::Vector{Float64}: The true relative residual ‖inp - opr * out‖ / ‖inp‖, exit value last
• hss::Vector{Matrix{ComplexF64}}: The GMRES Hessenberg matrix per restart cycle, as built, augmented by the recycled directions for a GCRODRSolver
• oprApp::Int: Operator applications made by the solver
• pcnApp::Int: Preconditioner applications made by the solver
• oprAppLog::Int: Operator applications made to fill this record, kept out of oprApp
• oprAppPcn::Int: Operator applications made inside ldiv!, which only the preconditioner can count: left at zero here
• ρ, α, ω, dnm::Vector{ComplexF64}: BiCGStab's scalars and its dot(resShd, v) denominator, per iteration
• rcyErr::Vector{Float64}: ‖CᴴC - I‖ for a GCRODRSolver's recycled subspace, per cycle
• prm::NamedTuple: The settings the solve ran under, resolved and rescaled, and the iterations restarts fell on
• inner::Vector{SlvLog}: One record per inner solve, for a solver that wraps another

SusOpr{T, A}

The susceptibility operator X, taking a field to the polarization current it drives, over the tiling cvol.

sus is either a vector in the flat degree of freedom layout, applied entrywise, or a (celTot, 3, 3) array of cell tensors in coordinates(cvol) order, applied as a matrix in each cell (gyrotropy, rotated crystals). χ is dimensionless and diagonal in position, so it commutes with the ΔV weight of the GlaFld inner product. X \ f divides pointwise and throws where χ vanishes.

SusOpr{T}(opr::SusOpr)

Convert the storage precision of a susceptibility operator to T.

SusOpr{T}(cvol::GlaCmpVol, sus; dev::Backend)
SusOpr{T}(vol::GlaVol, sus; dev::Backend)
SusOpr(vol, sus; dev::Backend)

The susceptibility operator of a volume, plain or composite. dev defaults to the backend of sus when it is an array, else CPU().

sus is a number, a 3×3 matrix, an array of one value per cell, an array of one 3×3 tensor per cell ((cel..., 3, 3)), a function of the cell center returning either, a vector holding any of those per region, a vector of one value per cell in layout order, or a vector already in the degree of freedom layout. The shapes that name a cell grid only fit a tiling of one region.

asym(opr::Union{GlaOprVac, GlaCmpOprVac})
sym(opr::Union{GlaOprVac, GlaCmpOprVac})

The anti-Hermitian part (G₀ - G₀†) / 2i of a self vacuum operator, the radiative part that enters the fluctuation-dissipation theorem, or its Hermitian part (G₀ + G₀†) / 2. The part is rebuilt from the volume and options of opr with hrmPrt, which takes it from the kernel before the transform, so it costs a build rather than reusing the coefficients of opr, and leaves opr untouched.

cellvolumes(cvol::GlaCmpVol)

Cell volumes of a composite volume, one entry per degree of freedom.

The result has length sum(3 * prod(r.cel) for r in regions(cvol)) and lines up entry by entry with a flat field vector: each cell volume is repeated across the three vector components of that cell. A region has a single cell size, so the whole block of a region carries one value.

Arguments

• cvol::GlaCmpVol: The composite volume

Returns

• Vector{Float64}: The cell volumes in cubic wavelengths

coordinates(cvol::GlaCmpVol)

Iterate over the cells of a composite volume in flat layout order.

The iterator yields one (pos, vol, idx) triple per cell: pos is the cell center in wavelengths, vol is the cell volume in cubic wavelengths, and idx is the index of the region the cell belongs to. Cells come out region by region and, inside a region, in column major order over the cell grid, which is the same order the flat field layout uses. There is no repetition over the three vector components, so the iterator is a third as long as a field vector.

Arguments

• cvol::GlaCmpVol: The composite volume

Returns

• A generator of Tuple{NTuple{3,Float64},Float64,Int}

discretize!(fld::GlaFld, f)

Fill a field from a function of position.

f takes a cell center as an NTuple{3,Float64} in wavelengths and returns something indexable of length three, an SVector, a tuple, or a vector, holding the three components of the current density there. Sampling is setup code, so the values are always built on the CPU and copied to the buffer at the end.

Arguments

• fld::GlaFld: The field to fill
• f: The current density as a function of position

Returns

• GlaFld: The field, filled

eachregion(fld::GlaFld)

Iterate over the region views of a field, in regions order.

Arguments

• fld::GlaFld: The field

Returns

• An iterator of the regionview of every region

finest(cvol::GlaCmpVol)

The smallest cell scale present in a composite volume, taken elementwise.

Arguments

• cvol::GlaCmpVol: The composite volume

Returns

• NTuple{3,Rational}: The elementwise minimum of the region cell scales

isadjoint(opr::AbstractGlaOpr)

Whether opr is in adjoint mode, see adjoint!.

isexternaloperator(opr::AbstractGlaOpr)

Whether opr maps between two volumes that share no interior.

isoverlappingoperator(opr::AbstractGlaOpr)

Whether opr maps between two distinct volumes that share interior, built on their union and masked.

isquasistatic(opr::AbstractGlaOpr)

Whether opr was built from the quasistatic Green function rather than the full one.

isselfoperator(opr::AbstractGlaOpr)

Whether opr maps a volume to itself.

isvarying(preCon)

Whether preCon applies a different linear map on different applications.

The question is linearity in the vector handed over.

Defaults to false, so a preconditioner that varies has to say so.

nregions(cvol::GlaCmpVol)

The number of regions in a composite volume.

Arguments

• cvol::GlaCmpVol: The composite volume

Returns

• The region count

refine(cvol::GlaCmpVol, box; factor=2)

Refine a composite volume inside a box by carving.

Every region whose interior meets the box is replaced by a refined core plus the complement slabs that fill the rest of the region. The core takes factor times as many cells at factor times the resolution over the same physical extent. Regions the box misses are left alone and keep their position in the region list, so the pieces of region i appear where region i used to be, core first, then the nonempty slabs in the order xlo, xhi, ylo, yhi, zlo, zhi.

Snapping happens per region. The box is first intersected with the region, then grown outward to that region's cell boundaries, then grown outward again to even cell indices so that the core and every slab keep an even cell count. A box that crosses regions of different scale therefore snaps differently in each of them. The core always covers at least the requested box.

You may refine a volume multiple times, even in nested a fashion.

Arguments

• cvol::GlaCmpVol: The composite volume to refine
• box: The refinement box, either a GlaVol whose outer bounds are used, or a tuple (org, dims) of rational center and side lengths
• factor=2: The refinement factor, an Integer or an NTuple{3,Integer} with every entry at least 1. An entry of 1 leaves that direction at the parent resolution while still carving the region.

Returns

• GlaCmpVol: The refined composite volume, validated by the constructor

Throws

• ArgumentError: If any refinement factor is less than 1, or if the resulting tiling violates an invariant of GlaCmpVol

refine(vol::GlaVol, box; factor=2)

Refine a single volume inside a box, returning a composite volume.

Arguments

• vol::GlaVol: The volume to refine
• box: The refinement box, either a GlaVol or a tuple (org, dims)
• factor=2: The refinement factor, an Integer or an NTuple{3,Integer}

Returns

• GlaCmpVol: The refined composite volume

regions(cvol::GlaCmpVol)

The regions of a composite volume, in flat layout order.

Arguments

• cvol::GlaCmpVol: The composite volume

Returns

• Vector{GlaVol}: The regions

regionview(fld::GlaFld, idx::Integer)

View the block of one region as a 4-tensor.

The result has size (reg.cel..., 3) and shares memory with the field, so writing through it mutates the field. The entries are the density coefficients.

The shape and the index order are the ones GlaOprVac expects of a source or target array, so the view can be handed straight to an operator.

Arguments

• fld::GlaFld: The field
• idx::Integer: The region index, in the order of regions(fld.cvol)

Returns

• A 4-tensor view of size (reg.cel..., 3) into the field buffer

regrid(fld::GlaFld, scl::NTuple{3,Rational}=finest(fld.cvol))

Read a field off as current densities on a uniform grid.

The grid has cells of size scl and covers the bounding box of the composite volume, which the regions tile exactly. Every uniform cell takes the density of the source cell holding its center, so the result is piecewise constant and the regions coarser than scl come out as blocks of repeated values. Since scl divides every region scale, cell boundaries never cut across a source cell.

Arguments

• fld::GlaFld: The field to read
• scl::NTuple{3,Rational}=finest(fld.cvol): The uniform cell size

Returns

• Array{Complex{T},4} of size (cel..., 3) with cel the uniform cell counts, T the precision of fld

Throws

• ArgumentError: If scl is coarser than or incommensurate with the cell size of any region, or if it does not evenly divide the bounding box

setBlc!(mem::GlaVacOprMem, θ::Tuple, blk = nothing)

Refill the Fourier data of a Bloch operator in place at the Bloch phases θ, reusing its plans and its BlcSet, or the block blk = (H, A) blcCmp gave a composite pair.

setBlc!(opr::AbstractGlaOpr, θ::Tuple)

Set the Bloch phases θ of a Bloch self operator, or of the one inside a scattering operator, in place. On a self operator a phase costs one fill and one transform and reuses the plans and the per-frequency BlcSet; a composite operator keeps one per-frequency κ cache for all its fills, and each fine-mesh block pays one more fill.

setSus!(opr::Union{InvSctOpr, SctOpr, GlaOpr}, sus)

Replace the susceptibility of opr, given as a SusOpr or in any form SusOpr takes, converted to the precision of opr. Returns opr.

slv(opr::AbstractGlaOpr)

The solver opr uses for \ and for applying (I - XG₀)⁻¹; a default GCRODRSolver() for vacuum operators.

solve(opr::AbstractGlaOpr{Thi}, inp::AbstractVector{Complex{Thi}}, slv::MixPrcRfn{Tlo}; log = nothing) where {Thi, Tlo}

Solve the linear system opr * out = inp by iterative refinement, with corrections computed in the precision Tlo.

Arguments

• opr::AbstractGlaOpr{Thi}: The operator in the linear system, of precision Thi
• inp::AbstractVector{Complex{Thi}}: The right hand side vector
• slv::MixPrcRfn{Tlo}: The solver parameters
• log: A SlvLog to fill in place, or nothing (default) for no record. One inner record per refinement step lands in log.inner

Returns

• out::AbstractVector{Complex{Thi}}: The solution vector

Notes

• The Tlo copy of the operator is cached and reused across solves with the same operator, and rebuilt as soon as a different operator is passed
• Refinement in the precision of the operator is pointless, and throws

solve(opr, inp::AbstractArray{T}, slv::GCRODRSolver; x0 = nothing, log = nothing, rcy = nothing) where T

Solve the linear system opr * out = inp using the GCRO-DR method.

Arguments

• opr: The operator in the linear system, anything supporting mul! and size
• inp::AbstractArray{T}: The right-hand side vector
• slv::GCRODRSolver: The solver parameters
• x0: An initial guess, or nothing (default) to start from zero
• log: A SlvLog to fill in place, or nothing (default) for no record
• rcy: A RcySpc to start from and fill in place, or nothing (default) to recycle within this solve only

Returns

• out::AbstractArray{T}: The solution vector

Notes

• A supplied subspace is checked against the operator of this solve before it is trusted, and rebuilt when it does not match
• log.hss holds the augmented matrix [I B; 0 H] of the cycle, which is Hessenberg only when rcyDim is zero

solve(opr, inp::AbstractArray{T}, slv::GMRESSolver; x0 = nothing, log = nothing) where T

Solve the linear system opr * out = inp using the GMRES method.

Arguments

• opr: The operator in the linear system, anything supporting mul! and size
• inp::AbstractArray{T}: The right-hand side vector
• slv::GMRESSolver: The solver parameters
• x0: An initial guess, or nothing (default) to start from zero
• log: A SlvLog to fill in place, or nothing (default) for no record

Returns

• out::AbstractArray{T}: The solution vector

Notes

• The solver will iterate until either convergence is achieved or the maximum number of iterations is reached
• Convergence is determined by the absolute and relative tolerances specified in the solver parameters
• The solver uses restarts to manage memory usage, with the restart frequency specified in the solver parameters

solve(opr, inp::AbstractVector{T}, slv::BiCGStabSolver; x0 = nothing, log = nothing) where T

Solve the linear system opr * out = inp using the BiCGStab method.

Arguments

• opr: The operator in the linear system, anything supporting mul! and size
• inp::AbstractVector{T}: The right-hand side vector
• slv::BiCGStabSolver: The solver parameters
• x0: An initial guess, or nothing (default) to start from zero
• log: A SlvLog to fill in place, or nothing (default) for no record

Returns

• out::AbstractVector{T}: The solution vector

Notes

• The solver will iterate until either convergence is achieved or the maximum number of iterations is reached
• Convergence is determined by the absolute and relative tolerances specified in the solver parameters

sus(opr::AbstractGlaOpr)

The susceptibility operator X of a scattering operator, shared rather than copied, so an adjoint! on either is seen by both. setSus! replaces it, and an earlier handle keeps the old one.

vac(opr::AbstractGlaOpr)

The vacuum operator of opr, shared rather than copied. A vacuum operator is its own.

get_backend(opr::AbstractGlaOpr)

The KernelAbstractions backend the operator computes on.

adjoint!(opr::AbstractGlaOpr)

The adjoint of opr, reusing its memory. Use only the return value afterwards: the argument may have been mutated or replaced. adjoint!(adjoint!(opr)) restores the original; adjoint works on a copy instead. This is the adjoint in the ΔV inner product of GlaFld, which is the matrix adjoint whenever all cells are the same size.

GlaCmpOprVac{T}

The vacuum Green operator between two composite volumes, as the block matrix of pairwise vacuum operators over their regions: the Galerkin discretization on the non-uniform mesh, with no stitching at region boundaries.

Blocks act on density coefficients, as a GlaFld holds them, with no cell volume weight.

Cross-scale pairs in contact, and every cross-scale pair under a Bloch lattice, are a GlaSndOprVac; every other block is a GlaOprVac.