API Reference

rscg(A, b, shifts; kwargs...) -> (x, stats)
rscg(A, b, shifts, V; kwargs...) -> (Ξ, stats)

Out-of-place Reduced-Shifted Conjugate Gradient (RSCG) method.

Solves the family of shifted linear systems:

\[(σⱼ I + A) x(σⱼ) = b, \quad j = 1, 2, \ldots, N_{\text{shifts}}\]

This is a convenience wrapper that allocates a workspace internally. For repeated solves with the same dimensions, use rscg! with a pre-allocated workspace.

Arguments

• A: Hermitian positive semi-definite linear operator
• b::AbstractVector: Right-hand side vector
• shifts::AbstractVector: Shift values $σⱼ$
• V::AbstractMatrix: (optional) Reduction matrix for reduced mode

Keyword Arguments

See rscg! for the full list.

Returns

• x::Vector{Vector{FC}}: Solution vectors, x[j] solves $(σⱼ I + A) x = b$
• stats::ReducedShiftStats: Solver statistics including iteration count and convergence status

For reduced mode (with V):

• Ξ::Vector{Vector{FC}}: Reduced solutions, $Ξ[j] = V x(σⱼ) ∈ \mathbb{C}^m$
• stats::ReducedShiftStats: Solver statistics

Example

using LinearAlgebra
using ReducedShiftedKrylov

n = 100
A = Hermitian(rand(n, n) + rand(n, n)' + 2n * I)
b = rand(ComplexF64, n)
shifts = [0.1 + 0.01im, 0.2 + 0.01im, 0.3 + 0.01im]

# Solve all shifted systems
x, stats = rscg(A, b, shifts)

println("Converged: $(stats.solved)")
println("Iterations: $(stats.niter)")

Green's Function Application

For computing Green's function $G(z) = (zI - H)^{-1}$:

# Set A = -H, shifts = z_values
# Then (σI + A)x = b becomes (zI - H)x = b
H = Hermitian(...)  # Hamiltonian
z_values = [ω + im*η for ω in ω_range]  # complex frequencies

x, stats = rscg(-H, b, z_values)
# x[j] = G(z_values[j]) * b

See also: rscg!, RscgWorkspace, ReducedShiftStats

rscg!(workspace, A, b, shifts; kwargs...)
rscg!(workspace, A, b, shifts, V; kwargs...)

In-place Reduced-Shifted Conjugate Gradient (RSCG) method.

Solves the family of shifted linear systems simultaneously:

\[(σⱼ I + A) x(σⱼ) = b, \quad j = 1, 2, \ldots, N_{\text{shifts}}\]

The algorithm exploits the shift-invariance of Krylov subspaces: $K_k(σI + A, b) = K_k(A, b)$, requiring only one matrix-vector product per iteration regardless of the number of shifts.

Arguments

• workspace::RscgWorkspace{T,FC,S}: Pre-allocated workspace (created by RscgWorkspace)
• A: Hermitian positive semi-definite linear operator (matrix or LinearMap)
• b::AbstractVector{FC}: Right-hand side vector
• shifts::AbstractVector: Shift values $σⱼ$ (can be complex)
• V::AbstractMatrix: (optional) Reduction matrix $V ∈ \mathbb{C}^{m × n}$ for reduced mode. When provided, computes $Ξ(σⱼ) = V x(σⱼ)$ instead of full solutions.

Keyword Arguments

• atol::Real=√eps(T): Absolute tolerance for convergence $\|r_k(σⱼ)\| ≤ \text{atol}$
• rtol::Real=√eps(T): Relative tolerance $\|r_k(σⱼ)\| ≤ \text{rtol} \cdot \|b\|$
• itmax::Int=0: Maximum iterations (0 means 2n)
• timemax::Float64=Inf: Time limit in seconds
• verbose::Int=0: Verbosity level (0=silent, k>0 prints every k iterations)
• history::Bool=false: If true, store residual norm history in stats.residuals
• callback=workspace->false: Early termination callback. Return true to stop.

Returns

• workspace::RscgWorkspace: Updated workspace. Use results(workspace) to extract (x, stats).

Convergence

The residual for each shift is estimated as:

\[\|r_k(σⱼ)\| = |ρ_k(σⱼ)| \cdot \|r_k\|\]

where $r_k$ is the seed system residual and $ρ_k(σⱼ)$ is the shift coefficient.

Example

using LinearAlgebra
using ReducedShiftedKrylov

# Create Hermitian positive definite matrix
n = 100
H = rand(ComplexF64, n, n)
A = Hermitian(H + H' + 2n * I)
b = rand(ComplexF64, n)
shifts = [0.1im, 0.2im, 0.3im]

# Create workspace and solve
workspace = RscgWorkspace(A, b, length(shifts))
rscg!(workspace, A, b, shifts; verbose=1)
x, stats = results(workspace)

# Verify: x[j] solves (A + shifts[j]*I) x = b
for j in eachindex(shifts)
    r = (A + shifts[j] * I) * x[j] - b
    println("Shift $j: residual = $(norm(r))")
end

Reduced Mode Example

# Compute only V * x(σⱼ) for memory efficiency
m = 5
V = zeros(ComplexF64, m, n)
for i in 1:m
    V[i, i] = 1.0
end

workspace = RscgWorkspace(A, b, length(shifts), V)
rscg!(workspace, A, b, shifts, V)
Ξ, stats = results(workspace)
# Ξ[j] ∈ ℂᵐ instead of x[j] ∈ ℂⁿ

References

• Y. Nagai et al., "Reduced-Shifted Conjugate-Gradient Method for a Green's Function", J. Phys. Soc. Jpn. 86, 014708 (2017). DOI:10.7566/JPSJ.86.014708

See also: rscg, RscgWorkspace, results

krylov_solve(Val(:rscg), A, b, shifts; kwargs...) -> (x, stats)
krylov_solve(Val(:rscg), A, b, shifts, V; kwargs...) -> (Ξ, stats)

Generic interface for reduced-shifted Krylov solvers.

This provides a unified API for selecting different Krylov methods at runtime using a Val type selector. Currently supports:

• Val(:rscg): Reduced-Shifted Conjugate Gradient

Arguments

• Val(:rscg): Solver selector (use Val(:rscg) for RSCG)
• A: Hermitian positive semi-definite linear operator
• b::AbstractVector: Right-hand side vector
• shifts::AbstractVector: Shift values $σⱼ$
• V::AbstractMatrix: (optional) Reduction matrix for reduced mode

Keyword Arguments

See rscg! for the full list of keyword arguments.

Returns

Same as rscg:

• x::Vector{Vector{FC}}: Solution vectors (non-reduced mode)
• Ξ::Vector{Vector{FC}}: Reduced solutions (reduced mode)
• stats::ReducedShiftStats: Solver statistics

Example

using ReducedShiftedKrylov

# Select solver at runtime
solver = :rscg
x, stats = krylov_solve(Val(solver), A, b, shifts)

Extensibility

To add a new solver (e.g., :rsbicg), define:

krylov_solve(::Val{:rsbicg}, A, b, shifts; kwargs...) = rsbicg(A, b, shifts; kwargs...)

See also: rscg, krylov_workspace

krylov_workspace(Val(:rscg), A, b, nshifts) -> RscgWorkspace
krylov_workspace(Val(:rscg), A, b, nshifts, V) -> RscgWorkspace

Create a pre-allocated workspace for reduced-shifted Krylov solvers.

This provides a unified API for creating workspaces, enabling efficient memory reuse when solving multiple systems with the same dimensions.

Arguments

• Val(:rscg): Solver selector
• A: Linear operator (used to infer dimensions)
• b::AbstractVector: Right-hand side vector (used to infer element type)
• nshifts::Integer: Number of shift points
• V::AbstractMatrix: (optional) Reduction matrix for reduced mode

Returns

• workspace::RscgWorkspace: Pre-allocated workspace

Example

using ReducedShiftedKrylov
using LinearAlgebra

n = 100
A = Hermitian(rand(ComplexF64, n, n) + rand(ComplexF64, n, n)' + 2n * I)
b = rand(ComplexF64, n)
shifts = [0.1im, 0.2im, 0.3im]

# Create workspace once
workspace = krylov_workspace(Val(:rscg), A, b, length(shifts))

# Reuse for multiple solves
for i in 1:10
    b_i = rand(ComplexF64, n)
    rscg!(workspace, A, b_i, shifts)
    x_i, stats_i = results(workspace)
    # Process results...
end

Reduced Mode Example

# Create workspace for reduced mode
m = 5
V = zeros(ComplexF64, m, n)
for i in 1:m
    V[i, i] = 1.0
end

workspace = krylov_workspace(Val(:rscg), A, b, length(shifts), V)
rscg!(workspace, A, b, shifts, V)
Ξ, stats = results(workspace)

See also: RscgWorkspace, rscg!, krylov_solve

RscgWorkspace{T,FC,S} <: RSKrylovWorkspace{T,FC,S}

Pre-allocated workspace for the Reduced-Shifted Conjugate Gradient (RSCG) method.

Using a workspace allows efficient memory reuse when solving multiple systems with the same dimensions but different right-hand sides or shifts.

Type Parameters

• T <: AbstractFloat: Real floating-point type (e.g., Float64)
• FC: Float or Complex type (e.g., ComplexF64)
• S <: AbstractVector{FC}: Vector storage type (e.g., Vector{ComplexF64})

Constructors

RscgWorkspace(m, n, nshifts, S::Type; reduced=false, mreduced=0)
RscgWorkspace(A, b, nshifts)           # Non-reduced mode
RscgWorkspace(A, b, nshifts, V)        # Reduced mode

Fields

Problem Dimensions

• m::Int, n::Int: Matrix dimensions (must be square, m == n)
• nshifts::Int: Number of shift points
• reduced::Bool: Whether workspace is configured for reduced mode
• mreduced::Int: Dimension of reduced vectors (number of rows in V)

Seed System Vectors (dimension n)

These are used for the "seed" linear system (σ=0):

• x::S: Solution vector
• r::S: Residual vector $r_k = b - A x_k$
• p::S: Search direction
• Ap::S: Matrix-vector product $A p_k$

Non-Reduced Mode Arrays (dimension n × nshifts)

Used when reduced=false:

• x_shifts::Vector{S}: Solution vectors $x(σⱼ)$ for each shift
• p_shifts::Vector{S}: Search directions $p_k(σⱼ)$ for each shift

Reduced Mode Arrays (dimension mreduced × nshifts)

Used when reduced=true:

• Σ::S: Current reduced residual $Σ_k = V r_k$
• Ξ::Vector{S}: Reduced solutions $Ξ(σⱼ) = V x(σⱼ)$ for each shift
• Π::Vector{S}: Reduced search directions $Π_k(σⱼ) = V p_k(σⱼ)$

Scalar Coefficients

• ρ::Vector{FC}, ρ_prev::Vector{FC}: Shift coefficients $ρ_k(σⱼ)$, $ρ_{k-1}(σⱼ)$
• α_shifts::Vector{FC}, β_shifts::Vector{FC}: Per-shift CG coefficients
• α::Ref{FC}, α_prev::Ref{FC}: Seed system $α_k$, $α_{k-1}$
• β::Ref{FC}, β_prev::Ref{FC}: Seed system $β_k$, $β_{k-1}$

Convergence Tracking

• rNorms::Vector{T}: Current residual norms $\|r_k(σⱼ)\|$ for each shift
• converged::BitVector: Convergence flags for each shift
• not_cv::BitVector: Not-converged flags (inverse of converged)
• stagnated::BitVector: Flags indicating ρ coefficient instability (loss of significance)

Statistics

• stats::ReducedShiftStats{T}: Solver statistics (iterations, timing, status)

Example: Workspace Reuse

using LinearAlgebra
using ReducedShiftedKrylov

n = 100
A = Hermitian(rand(ComplexF64, n, n) + rand(ComplexF64, n, n)' + 2n * I)
shifts = [0.1im, 0.2im, 0.3im]

# Create workspace once
b = rand(ComplexF64, n)
workspace = RscgWorkspace(A, b, length(shifts))

# Solve multiple systems with different right-hand sides
for i in 1:10
    b_i = rand(ComplexF64, n)
    rscg!(workspace, A, b_i, shifts)
    x_i, stats_i = results(workspace)
    # Process results...
end

Memory Comparison

For Green's function calculations where m << n (e.g., m=10, n=10000), reduced mode provides significant memory savings.

See also: rscg!, rscg, results, ReducedShiftStats

ReducedShiftStats{T} <: RSKrylovStats{T}

Statistics returned by RSCG and related reduced-shifted Krylov methods.

This structure captures convergence information, timing, and optionally the residual history for each shift point.

Type Parameter

• T <: AbstractFloat: Real floating-point type (e.g., Float64)

Fields

• niter::Int: Total number of iterations performed
• solved::Bool: true if all shifts converged within tolerance
• residuals::Vector{Vector{T}}: Residual norm history for each shift (only populated if history=true was passed to the solver)
• indefinite::BitVector: Flags indicating indefinite matrix detection per shift (reserved for future use)
• timer::Float64: Elapsed wall-clock time in seconds
• status::String: Human-readable description of the termination reason

Status Messages

Possible values for status:

• "solution good enough given atol and rtol": All shifts converged
• "maximum number of iterations exceeded": Reached itmax without convergence
• "user-requested exit": Callback returned true
• "time limit exceeded": Exceeded timemax
• "x = 0 is exact solution": Zero right-hand side vector
• "unknown": Initial state before solving

Example

x, stats = rscg(A, b, shifts; history=true)

println("Converged: $(stats.solved)")
println("Iterations: $(stats.niter)")
println("Time: $(stats.timer) s")
println("Status: $(stats.status)")

# Access residual history (if history=true)
if !isempty(stats.residuals[1])
    using Plots
    for j in eachindex(shifts)
        plot!(stats.residuals[j], label="Shift $j", yscale=:log10)
    end
end

See also: rscg, rscg!, RscgWorkspace, reset!

results(workspace::RscgWorkspace) -> (solutions, stats)

Extract solutions and statistics from a solved workspace.

Arguments

• workspace::RscgWorkspace: Workspace after calling rscg!

Returns

For non-reduced mode (workspace.reduced == false):

• x::Vector{Vector{FC}}: Solution vectors, where x[j] solves $(σⱼ I + A) x = b$
• stats::ReducedShiftStats{T}: Solver statistics

For reduced mode (workspace.reduced == true):

• Ξ::Vector{Vector{FC}}: Reduced solutions, where $Ξ[j] = V x(σⱼ)$
• stats::ReducedShiftStats{T}: Solver statistics

Example

workspace = RscgWorkspace(A, b, length(shifts))
rscg!(workspace, A, b, shifts)
x, stats = results(workspace)

println("Solved: $(stats.solved)")
println("Iterations: $(stats.niter)")
println("Time: $(stats.timer) seconds")

Note

The returned vectors are references to the workspace's internal storage. If you need to keep the results while reusing the workspace, make a copy:

x_copy = deepcopy(x)

See also: rscg!, RscgWorkspace, ReducedShiftStats

reset!(stats::ReducedShiftStats) -> ReducedShiftStats

Reset statistics to initial state for reuse.

Clears all fields:

• niter → 0
• solved → false
• residuals → empty vectors
• indefinite → all false
• timer → 0.0
• status → "unknown"

Example

# After solving
x, stats = rscg(A, b, shifts)
println(stats.niter)  # e.g., 42

# Reset for reuse
reset!(stats)
println(stats.niter)  # 0

This is called automatically at the start of rscg!.