### DQGMRES Solver with Various Preconditioning Strategies

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/examples/dqgmres.md

This example demonstrates solving a linear system Ax = b using DQGMRES. It showcases solving without preconditioning, and with split, left, and right preconditioning using an ILU(0) preconditioner. The setup involves loading a sparse matrix from MatrixMarket and constructing LinearOperator objects for the preconditioner.

```julia
using Krylov, LinearOperators, ILUZero, MatrixMarket
using LinearAlgebra, Printf, SuiteSparseMatrixCollection

ssmc = ssmc_db(verbose=false)
matrix = ssmc_matrices(ssmc, "Simon", "raefsky1")
path = fetch_ssmc(matrix, format="MM")

n = matrix.nrows[1]
A = MatrixMarket.mmread(joinpath(path[1], "$(matrix.name[1]).mtx"))
b = A * ones(n)

F = ilu0(A)

@printf("nnz(ILU) / nnz(A): %7.1e\n", nnz(F) / nnz(A))

# Solve Ax = b with DQGMRES and an ILU(0) preconditioner
# Remark: DIOM, FOM and GMRES can be used in the same way
opM = LinearOperator(Float64, n, n, false, false, (y, v) -> forward_substitution!(y, F, v))
opN = LinearOperator(Float64, n, n, false, false, (y, v) -> backward_substitution!(y, F, v))
opP = LinearOperator(Float64, n, n, false, false, (y, v) -> ldiv!(y, F, v))

# Without preconditioning
x, stats = dqgmres(A, b, memory=50, history=true)
r = b - A * x
@printf("[Without preconditioning] Residual norm: %8.1e\n", norm(r))
@printf("[Without preconditioning] Number of iterations: %3d\n", length(stats.residuals) - 1)

# Split preconditioning
x, stats = dqgmres(A, b, memory=50, history=true, M=opM, N=opN)
r = b - A * x
@printf("[Split preconditioning] Residual norm: %8.1e\n", norm(r))
@printf("[Split preconditioning] Number of iterations: %3d\n", length(stats.residuals) - 1)

# Left preconditioning
x, stats = dqgmres(A, b, memory=50, history=true, M=opP)
r = b - A * x
@printf("[Left preconditioning] Residual norm: %8.1e\n", norm(r))
@printf("[Left preconditioning] Number of iterations: %3d\n", length(stats.residuals) - 1)

# Right preconditioning
x, stats = dqgmres(A, b, memory=50, history=true, N=opP)
r = b - A * x
@printf("[Right preconditioning] Residual norm: %8.1e\n", norm(r))
@printf("[Right preconditioning] Number of iterations: %3d\n", length(stats.residuals) - 1)
```

--------------------------------

### Example MPI Krylov Solver Script

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/custom_workspaces.md

An example script demonstrating the usage of MPIVector and MPIOperator with Krylov.jl. This script can be executed with multiple MPI processes.

```julia
MPI.Init()

comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
n = 10  # global_len in MPIVector
T = Float64
method = :cg

# Version without MPI
A_global = Diagonal(1:n) .|> T
b_global = collect(n:-1:1) .|> T
x_global, _ = krylov_solve(Val(method), A_global, b_global)

# Version with MPI
A_mpi = MPIOperator{T}(n, n)
b_mpi = MPIVector(b_global)

kc = KrylovConstructor(b_mpi)
workspace = krylov_workspace(Val(method), kc)
krylov_solve!(workspace, A_mpi, b_mpi)
x_mpi = Krylov.solution(workspace)

```

--------------------------------

### Solve Linear Systems with Shifted Matrices using cg_lanczos_shift

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/examples/cg_lanczos_shift.md

This example demonstrates how to solve a linear system (A + αI)x = b for multiple shift values α using the cg_lanczos_shift function. It includes setup for matrix loading, defining shifts, solving the system, and calculating/printing relative residuals.

```julia
using Krylov, MatrixMarket, SuiteSparseMatrixCollection
using LinearAlgebra, Printf

function residuals(A, b, shifts, x)
  nshifts = length(shifts)
  r = [ (b - A * x[i] - shifts[i] * x[i]) for i = 1 : nshifts ]
  return r
end
ssmc = ssmc_db(verbose=false)
matrix = ssmc_matrices(ssmc, "HB", "1138_bus")
path = fetch_ssmc(matrix, format="MM")

A = MatrixMarket.mmread(joinpath(path[1], "$(matrix.name[1]).mtx"))
n, m = size(A)
b = ones(n)

# Solve (A + αI)x = b.
shifts = [1.0, 2.0, 3.0, 4.0]
(x, stats) = cg_lanczos_shift(A, b, shifts)
show(stats)
r = residuals(A, b, shifts, x)
resids = map(norm, r) / norm(b)
@printf("Relative residuals with shifts:\n")
for resid in resids
  @printf(" %8.1e", resid)
end
@printf("\n")

```

--------------------------------

### Explicit Warm-start with cg_lanczos!

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/warm-start.md

Manually perform a warm start for `cg_lanczos!` by calculating the residual and solving for the update.

```julia
workspace = CgLanczosWorkspace(A, b)

cg_lanczos!(workspace, A, b)
x₀ = workspace.x  # Ax₀ ≈ b
r = b - A * x₀    # r = b - Ax₀

cg_lanczos!(workspace, A, r)
Δx = workspace.x  # AΔx = r
x = x₀ + Δx       # Ax = b
```

--------------------------------

### Install and Test Krylov.jl Package

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/README.md

Install the Krylov package using the Julia package manager and run its tests to ensure proper functionality.

```julia
julia> ]
pkg> add Krylov
pkg> test Krylov
```

--------------------------------

### MINRES-QLP Example for Symmetric Systems

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/examples/minres_qlp.md

Demonstrates how to use the MINRES-QLP solver to find the minimum-norm solution for a symmetric, singular linear system. Requires Krylov.jl and LinearAlgebra packages.

```julia
using Krylov
using LinearAlgebra, Printf

A = diagm([1.0; 2.0; 3.0; 0.0])
n = size(A, 1)
b = [1.0; 2.0; 3.0; 4.0]
b_norm = norm(b)

# MINRES-QLP returns the minimum-norm solution of symmetric, singular and inconsistent systems
(x, stats) = minres_qlp(A, b);
r = b - A * x;

@printf("Residual r: %s\n", Krylov.vec2str(r))
@printf("Relative residual norm ‖r‖: %8.1e\n", norm(r) / b_norm)
@printf("Solution x: %s\n", Krylov.vec2str(x))
@printf("Minimum-norm solution? %s\n", x ≈ [1.0; 1.0; 1.0; 0.0])
```

--------------------------------

### 2D Poisson Equation Solver Setup

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/custom_workspaces.md

Sets up parameters for solving a 2D Poisson equation on a square domain using Krylov methods. This includes defining domain size and grid spacing.

```julia
using Krylov, OffsetArrays
import Krylov: solution, statistics

# Parameters
L = 1.0            # Length of the square domain
Nx = 200           # Number of interior grid points in x
Ny = 200           # Number of interior grid points in y
Δx = L / (Nx + 1)  # Grid spacing in x
Δy = L / (Ny + 1)  # Grid spacing in y
```

--------------------------------

### Implementing GMRES(k) with Warm Starts

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/warm-start.md

Use warm starts to implement the restarted GMRES(k) method, allowing for a fixed number of iterations per restart.

```julia
k = 50
workspace = GmresWorkspace(A, b, k)  # FomWorkspace(A, b, k)
workspace.x .= 0                     # workspace.x .= x₀
nrestart = 0
while !issolved(workspace) || nrestart ≤ 10
  solve!(workspace, A, b, workspace.x, itmax=k)
  nrestart += 1
end
```

--------------------------------

### BICGSTAB Solver with Various Preconditioning Strategies

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/examples/bicgstab.md

This example demonstrates solving a linear system Ax = b using the BICGSTAB algorithm. It showcases solving without preconditioning, and with split, left, and right preconditioning using an incomplete LU factorization.

```julia
using Krylov, LinearOperators, KrylovPreconditioners, HarwellRutherfordBoeing
using LinearAlgebra, Printf, SuiteSparseMatrixCollection, SparseArrays

ssmc = ssmc_db(verbose=false)
matrix = ssmc_matrices(ssmc, "HB", "sherman5")
path = fetch_ssmc(matrix, format="RB")

n = matrix.nrows[1]
A = RutherfordBoeingData(joinpath(path[1], "$(matrix.name[1]).rb")).data
b = A * ones(n)

F = ilu(A, τ = 0.05)

@printf("nnz(ILU) / nnz(A): %7.1e\n", nnz(F) / nnz(A))

# Solve Ax = b with BICGSTAB and an incomplete LU factorization
# Remark: CGS can be used in the same way
opM = LinearOperator(Float64, n, n, false, false, (y, v) -> forward_substitution!(y, F, v))
opN = LinearOperator(Float64, n, n, false, false, (y, v) -> backward_substitution!(y, F, v))
opP = LinearOperator(Float64, n, n, false, false, (y, v) -> ldiv!(y, F, v))

# Without preconditioning
x, stats = bicgstab(A, b, history=true)
r = b - A * x
@printf("[Without preconditioning] Residual norm: %8.1e\n", norm(r))
@printf("[Without preconditioning] Number of iterations: %3d\n", length(stats.residuals) - 1)

# Split preconditioning
x, stats = bicgstab(A, b, history=true, M=opM, N=opN)
r = b - A * x
@printf("[Split preconditioning] Residual norm: %8.1e\n", norm(r))
@printf("[Split preconditioning] Number of iterations: %3d\n", length(stats.residuals) - 1)

# Left preconditioning
x, stats = bicgstab(A, b, history=true, M=opP)
r = b - A * x
@printf("[Left preconditioning] Residual norm: %8.1e\n", norm(r))
@printf("[Left preconditioning] Number of iterations: %3d\n", length(stats.residuals) - 1)

# Right preconditioning
x, stats = bicgstab(A, b, history=true, N=opP)
r = b - A * x
@printf("[Right preconditioning] Residual norm: %8.1e\n", norm(r))
@printf("[Right preconditioning] Number of iterations: %3d\n", length(stats.residuals) - 1)

```

--------------------------------

### Warm-start for Block Method (TriMR)

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/warm-start.md

Implement a warm start for the TriMR block method by calculating initial residuals and solving for updates.

```julia
# [E  A] [x] = [b]
# [Aᴴ F] [y]   [c]
M = inv(E)
N = inv(F)
x₀, y₀, stats = trimr(A, b, c, M=M, N=N)

# E and F are not available inside TriMR
b₀ = b -  Ex₀ - Ay
c₀ = c - Aᴴx₀ - Fy

Δx, Δy, stats = trimr(A, b₀, c₀, M=M, N=N)
x = x₀ + Δx
y = y₀ + Δy
```

--------------------------------

### MPI Krylov Solver Script

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/custom_workspaces.md

This script demonstrates a basic distributed Krylov solver setup using MPI. It checks the solution across different MPI ranks and finalizes the MPI environment.

```julia
for i = 0:MPI.Comm_size(comm)-1
    if i == rank
        println("Local data: ", x_mpi.data)
        println("Global data: ", x_global)
    end
    MPI.Barrier(comm)
end

MPI.Finalize()
```

--------------------------------

### Newton's Method using In-place CG Solver

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/inplace.md

An example of implementing Newton's method for convex optimization using an in-place Conjugate Gradient (CG) solver. It shows how to create a CG workspace and use the `cg!` function to solve the linear system at each iteration.

```julia
using Krylov
import Krylov: solution

function newton(∇f, ∇²f, x₀; itmax = 200, tol = 1e-8)

    n = length(x₀)
    x = copy(x₀)
    gx = ∇f(x)
    
    iter = 0
    S = typeof(x)
    workspace = CgWorkspace(n, n, S)

    solved = false
    tired = false

    while !(solved || tired) 

        Hx = ∇²f(x)              # Compute ∇²f(xₖ)
        cg!(workspace, Hx, -gx)  # Solve ∇²f(xₖ)Δx = -∇f(xₖ)
        Δx = solution(workspace) # Recover Δx from the workspace
        x = x + Δx               # Update xₖ₊₁ = xₖ + Δx
        gx = ∇f(x)               # ∇f(xₖ₊₁)
        
        iter += 1
        solved = norm(gx) ≤ tol
        tired = iter ≥ itmax
    end
    return x
end
```

--------------------------------

### CGNE Solver Example

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/examples/cgne.md

Solves a linear system using the CGNE method. It fetches a matrix from the SuiteSparseMatrixCollection, constructs a system, solves it, and prints the relative residual and the error in the solution.

```julia
using Krylov, HarwellRutherfordBoeing, SuiteSparseMatrixCollection
using LinearAlgebra, Printf

ssmc = ssmc_db(verbose=false)
matrix = ssmc_matrices(ssmc, "HB", "wm2")
path = fetch_ssmc(matrix, format="RB")

A = RutherfordBoeingData(joinpath(path[1], "$(matrix.name[1]).rb")).data
(m, n) = size(A)
@printf("System size: %d rows and %d columns\n", m, n)

x_exact = A' * ones(m)
x_exact_norm = norm(x_exact)
x_exact /= x_exact_norm
b = A * x_exact
(x, stats) = cgne(A, b)
show(stats)
resid = norm(A * x - b) / norm(b)
@printf("CGNE: Relative residual: %7.1e\n", resid)
@printf("CGNE: ‖x - x*‖₂: %7.1e\n", norm(x - x_exact))
```

--------------------------------

### ILU Preconditioner with GMRES (Left Preconditioning)

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/preconditioners.md

Example of using an Incomplete LU factorization preconditioner from KrylovPreconditioners.jl with GMRES, specifying left preconditioning.

```julia
using KrylovPreconditioners, Krylov
Pℓ = ilu(A)
x, stats = gmres(A, b, M=Pℓ, ldiv=true)  # left preconditioning
```

--------------------------------

### Sparse LU Preconditioner for Gpmr

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/preconditioners.md

Example of using a sparse LU factorization from SuiteSparse.jl as a preconditioner (C) for the Gpmr solver on a block system.

```julia
using SuiteSparse, Krylov

C = lu(M)

# [M  A] [x] = [b]
# [B  0] [y]   [c]
x, y, stats = gpmr(A, B, b, c, C=C, gsp=true, ldiv=true)
```

--------------------------------

### Block-Jacobi Preconditioner with GMRES

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/preconditioners.md

Example of using a Block-Jacobi preconditioner from KrylovPreconditioners.jl with the GMRES solver.

```julia
using KrylovPreconditioners, Krylov

P⁻¹ = BlockJacobiPreconditioner(A)  # Block-Jacobi preconditioner
x, stats = gmres(A, b, M=P⁻¹)
```

--------------------------------

### Solve Sparse System with ILU(0) Preconditioner on NVIDIA GPU

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/gpu.md

Demonstrates solving a sparse linear system using bicgstab with an ILU(0) preconditioner on an NVIDIA GPU. Requires CUDA.jl and cuSPARSE/cuSOLVER.

```julia
using SparseArrays, Krylov, LinearOperators
using CUDA, CUDA.CUSPARSE, CUDA.CUSOLVER

if CUDA.functional()
  # Optional -- Compute a permutation vector p such that A[:,p] has no zero diagonal
  p = zfd(A_cpu)
  p .+= 1
  A_cpu = A_cpu[:,p]

  # Transfer the linear system from the CPU to the GPU
  A_gpu = CuSparseMatrixCSR(A_cpu)  # A_gpu = CuSparseMatrixCSC(A_cpu)
  b_gpu = CuVector(b_cpu)

  # ILU(0) decomposition LU ≈ A for CuSparseMatrixCSC or CuSparseMatrixCSR matrices
  P = ilu02(A_gpu)

  # Additional vector required for solving triangular systems
  n = length(b_gpu)
  T = eltype(b_gpu)
  z = CUDA.zeros(T, n)

  # Solve Py = x
  function ldiv_ilu0!(P::CuSparseMatrixCSR, x, y, z)
    ldiv!(z, UnitLowerTriangular(P), x)  # Forward substitution with L
    ldiv!(y, UpperTriangular(P), z)      # Backward substitution with U
    return y
  end

  function ldiv_ilu0!(P::CuSparseMatrixCSC, x, y, z)
    ldiv!(z, LowerTriangular(P), x)      # Forward substitution with L
    ldiv!(y, UnitUpperTriangular(P), z)  # Backward substitution with U
    return y
  end

  # Operator that model P⁻¹
  symmetric = hermitian = false
  opM = LinearOperator(T, n, n, symmetric, hermitian, (y, x) -> ldiv_ilu0!(P, x, y, z))

  # Solve a non-Hermitian system with an ILU(0) preconditioner on GPU
  x̄, stats = bicgstab(A_gpu, b_gpu, M=opM)

  # Recover the solution of Ax = b with the solution of A[:,p]x̄ = b
  invp = invperm(p)
  x = x̄[invp]
end
```

--------------------------------

### Get Maximum Available CPU Threads

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/tips.md

Retrieve the maximum number of threads available on the system.

```julia
NMAX = Sys.CPU_THREADS
```

--------------------------------

### Initialize Krylov.jl and LinearAlgebra for Helmholtz in Julia

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/matrix_free.md

Initializes the necessary packages, Krylov.jl and LinearAlgebra, for solving the 3D Helmholtz equation. This is a prerequisite for setting up and applying the matrix-free Helmholtz operator.

```julia
using Krylov, LinearAlgebra
```

--------------------------------

### GPU Computation with Krylov.jl and Metal.jl

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/gpu.md

Demonstrates converting CPU arrays to GPU arrays and solving a dense least-norm problem on an Apple M1 GPU using `craig` solver. Requires `Krylov` and `Metal` packages.

```julia
using Krylov, Metal

T = Float32  # Metal.jl also works with ComplexF32
n = 10
m = 20

# CPU Arrays
A_cpu = rand(T, n, m)
b_cpu = rand(T, n)

# GPU Arrays
A_gpu = MtlMatrix(A_cpu)
b_gpu = MtlVector(b_cpu)

# Solve a dense least-norm problem on an Apple M1 GPU
x, stats = craig(A_gpu, b_gpu)
```

--------------------------------

### Get Total Bytes of Workspace

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/storage.md

Calculates the total number of bytes used by a Krylov workspace. This value can then be formatted for human readability.

```julia
nbytes = sizeof(workspace)
```

--------------------------------

### Enable Julia Multithreading at Startup

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/tips.md

Configure the number of Julia threads to use when launching the Julia executable. This can be done via environment variables or command-line arguments.

```shell
julia -t auto  # alternative: --threads auto
julia -t N     # alternative: --threads N

JULIA_NUM_THREADS=N julia
```

--------------------------------

### Solve SPD System with IC(0) Preconditioner on GPU

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/gpu.md

Demonstrates solving a Hermitian positive-definite system on an NVIDIA GPU using an IC(0) preconditioner. This involves creating a CuSparseMatrixCSR, computing the IC(0) factorization, and defining a LinearOperator for the preconditioner.

```julia
using SparseArrays, Krylov, LinearOperators
using CUDA, CUDA.CUSPARSE

if CUDA.functional()
  # Transfer the linear system from the CPU to the GPU
  A_gpu = CuSparseMatrixCSR(A_cpu)  # A_gpu = CuSparseMatrixCSC(A_cpu)
  b_gpu = CuVector(b_cpu)

  # IC(0) decomposition LLᴴ ≈ A for CuSparseMatrixCSC or CuSparseMatrixCSR matrices
  P = ic02(A_gpu)

  # Additional vector required for solving triangular systems
  n = length(b_gpu)
  T = eltype(b_gpu)
  z = CUDA.zeros(T, n)

  # Solve Py = x
  function ldiv_ic0!(P::CuSparseMatrixCSR, x, y, z)
    ldiv!(z, LowerTriangular(P), x)   # Forward substitution with L
    ldiv!(y, LowerTriangular(P)', z)  # Backward substitution with Lᴴ
    return y
  end

  function ldiv_ic0!(P::CuSparseMatrixCSC, x, y, z)
    ldiv!(z, UpperTriangular(P)', x)  # Forward substitution with L
    ldiv!(y, UpperTriangular(P), z)   # Backward substitution with Lᴴ
    return y
  end

  # Operator that model P⁻¹
  symmetric = hermitian = true
  opM = LinearOperator(T, n, n, symmetric, hermitian, (y, x) -> ldiv_ic0!(P, x, y, z))

  # Solve a Hermitian positive definite system with an IC(0) preconditioner on GPU
  x, stats = cg(A_gpu, b_gpu, M=opM)
end
```

--------------------------------

### Set and Get BLAS Threads

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/tips.md

Programmatically set and retrieve the number of BLAS threads. The number of threads N must be between 1 and NMAX.

```julia
BLAS.set_num_threads(N)  # 1 ≤ N ≤ NMAX
BLAS.get_num_threads()
```

--------------------------------

### LDL Factorization for Tridiagonal Systems with TriCG

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/preconditioners.md

Example of using LDL factorizations from LDLFactorizations.jl as preconditioners (M and N) for the TriCG solver on a block tridiagonal system.

```julia
using LDLFactorizations, Krylov

M = ldl(E)
N = ldl(F)

# [E   A] [x] = [b]
# [Aᴴ -F] [y]   [c]
x, y, stats = tricg(A, b, c, M=M, N=N, ldiv=true)
```

--------------------------------

### In-Place Krylov Solver Interface with Workspace

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/generic_interface.md

Illustrates the in-place interface using `krylov_workspace` and `krylov_solve!`. This is useful for iterative refinement, performance-critical applications, or when detailed statistics and control over the solving process are needed. It provides access to convergence status, iteration counts, timers, and operator-vector product counts.

```julia
using Krylov, SparseArrays, LinearAlgebra

# Define a square nonsymmetric matrix A and a right-hand side vector b
n = 100
A = sprand(n, n, 0.05) + I
b = rand(n)

# In-place interface
for method in (:bicgstab, :gmres, :qmr)
    # Create a workspace for the Krylov method
    workspace = krylov_workspace(Val(method), A, b)

    # Solve the system in-place
    krylov_solve!(workspace, A, b)

    # Get the statistics
    stats = Krylov.statistics(workspace)

    # Retrieve the solution
    x = Krylov.solution(workspace)

    # Check if the solver converged
    solved = Krylov.issolved(workspace)
    println("Convergence of $method: ", solved)

    # Display the number of iterations
    niter = Krylov.iteration_count(workspace)
    println("Number of iterations for $method: ", niter)

    # Display the allocation timer
    allocation_timer = Krylov.elapsed_allocation_time(workspace)
    println("Elapsed allocation time for $method: ", allocation_timer, " seconds")

    # Display the computation timer
    computation_timer = Krylov.elapsed_time(workspace)
    println("Elapsed time for $method: ", computation_timer, " seconds")

    # Display the number of operator-vector products with A and A'
    nAprod = Krylov.Aprod_count(workspace)
    nAtprod = Krylov.Atprod_count(workspace)
    println("Number of operator-vector products with A and A' for $method: ", (nAprod, nAtprod))

    println()
end
```

--------------------------------

### Define MPIOperator for Distributed Linear Operators

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/custom_workspaces.md

Defines an MPIOperator struct to represent a distributed diagonal linear operator. This is a basic example for implementing distributed operators.

```julia
using Krylov

struct MPIOperator{T}
    m::Int
    n::Int
end

Base.size(A::MPIOperator) = (A.m, A.n)
Base.eltype(A::MPIOperator{T}) where T = T

function Krylov.kmul!(y::MPIVector{Float64}, A::MPIOperator{Float64}, u::MPIVector{Float64})
    y.data .= u.data_range .* u.data
    return y
end
```

--------------------------------

### Solving Linear System with Craig Method

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/examples/craig.md

This snippet demonstrates how to use the `craig` function to solve a linear system with a regularization parameter λ. It also shows how to display the statistics of the solution process and verify the minimum-norm property of the solution.

```julia
using Krylov
using LinearAlgebra, Printf

m = 5
n = 8
λ = 1.0e-3
A = rand(m, n)
b = A * ones(n)
xy_exact = [A  λ*I] \ b # In Julia, this is the min-norm solution!

(x, y, stats) = craig(A, b, λ=λ, atol=0.0, rtol=1.0e-20, verbose=1)
show(stats)

# Check that we have a minimum-norm solution.
# When λ > 0 we solve min ‖(x,s)‖  s.t. Ax + λs = b, and we get s = λy.
@printf("Primal feasibility: %7.1e\n", norm(b - A * x - λ^2 * y) / norm(b))
@printf("Dual   feasibility: %7.1e\n", norm(x - A' * y) / norm(x))
@printf("Error in x: %7.1e\n", norm(x - xy_exact[1:n]) / norm(xy_exact[1:n]))
if λ > 0.0
  @printf("Error in y: %7.1e\n", norm(λ * y - xy_exact[n+1:n+m]) / norm(xy_exact[n+1:n+m]))
end
```

--------------------------------

### Warm-start CG with CgWorkspace

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/warm-start.md

Use the approximate solution from a previous CG run as the initial guess for a new run.

```julia
workspace = CgWorkspace(A, b)
cg!(workspace, A, b, itmax=100)
if !issolved(workspace)
  # Use the approximate solution `workspace.x` as starting point
  cg!(workspace, A, b, workspace.x, itmax=100)
end
```

--------------------------------

### Krylov Solver Workspace Constructors

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/inplace.md

Demonstrates the standard constructors for Krylov workspaces, which store intermediate data for repeated solves with consistent dimensions and precision. The workspace is always the first argument for in-place methods.

```julia
minres_workspace = MinresWorkspace(m, n, Vector{Float64})
minres!(minres_workspace, A1, b1)

bicgstab_workspace = BicgstabWorkspace(m, n, Vector{ComplexF64})
bicgstab!(bicgstab_workspace, A2, b2)

gmres_workspace = GmresWorkspace(m, n, Vector{BigFloat})
gmres!(gmres_workspace, A3, b3)

lsqr_workspace = LsqrWorkspace(m, n, CuVector{Float32})
lsqr!(lsqr_workspace, A4, b4)
```

--------------------------------

### Direct Warm-start CG

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/warm-start.md

Provide an initial guess `x0` directly to the `cg` function.

```julia
cg(A, b, x0)
```

--------------------------------

### Storing GMRES Iterates with a Callback

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/callbacks.md

An example of using a callback function with the GMRES method to store all intermediate iterates ($x_k$). The callback accesses workspace variables like `z`, `inner_iter`, `V`, and `R` to reconstruct and store each iterate.

```julia
S = Krylov.ktypeof(b)
global X = S[]  # Storage for GMRES iterates

function gmres_callback(workspace)
  z = workspace.z
  k = workspace.inner_iter
  nr = sum(1:k)
  V = workspace.V
  R = workspace.R
  y = copy(z)

  # Solve Rk * yk = zk
  for i = k : -1 : 1
    pos = nr + i - k
    for j = k : -1 : i+1
      y[i] = y[i] - R[pos] * y[j]
      pos = pos - j + 1
    end
    y[i] = y[i] / R[pos]
  end

  # xk = Vk * yk
  xk = sum(V[i] * y[i] for i = 1:k)
  push!(X, xk)

  return false  # We don't want to add new stopping conditions
end

(x, stats) = gmres(A, b, callback = gmres_callback)
```

--------------------------------

### ILU0 Preconditioner with BiCGStab (Right Preconditioning)

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/preconditioners.md

Applies an ILU0 preconditioner from ILUZero.jl to the BiCGStab solver, demonstrating right preconditioning.

```julia
using ILUZero, Krylov
Pᵣ = ilu0(A)
x, stats = bicgstab(A, b, N=Pᵣ, ldiv=true)  # right preconditioning
```

--------------------------------

### Solve Sparse Systems on AMD GPU with Different Formats

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/gpu.md

Shows how to solve sparse linear systems on an AMD GPU using ROCSparseMatrixCSC, ROCSparseMatrixCSR, and ROCSparseMatrixCOO formats. Requires AMDGPU.jl and rocSPARSE.

```julia
using AMDGPU, Krylov
using AMDGPU.rocSPARSE, SparseArrays

if AMDGPU.functional()
  # CPU Arrays
  A_cpu = sprand(100, 200, 0.3)
  b_cpu = rand(100)

  # GPU Arrays
  A_csc_gpu = ROCSparseMatrixCSC(A_cpu)
  A_csr_gpu = ROCSparseMatrixCSR(A_cpu)
  A_coo_gpu = ROCSparseMatrixCOO(A_cpu)
  b_gpu = ROCVector(b_cpu)

  # Solve a rectangular and sparse system on an AMD GPU
  x_csc, y_csc, stats_csc = lnlq(A_csc_gpu, b_gpu)
  x_csr, y_csr, stats_csr = craig(A_csr_gpu, b_gpu)
  x_coo, y_coo, stats_coo = craigmr(A_coo_gpu, b_gpu)
end
```

--------------------------------

### Index Documentation

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/reference.md

Generates an index of the package's public API.

```julia
```@index
```
```

--------------------------------

### Solve Symmetric, Singular, Consistent System with Symmlq

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/examples/symmlq.md

Use symmlq to find the minimum-norm solution for a symmetric, singular, and consistent system. Ensure LinearAlgebra and Printf are imported for norm calculations and formatted output.

```julia
using Krylov
using LinearAlgebra, Printf

A = diagm([1.0; 2.0; 3.0; 0.0])
n = size(A, 1)
b = [1.0; 2.0; 3.0; 0.0]
b_norm = norm(b)

# SYMMLQ returns the minimum-norm solution of symmetric, singular and consistent systems
(x, stats) = symmlq(A, b, transfer_to_cg=false);
r = b - A * x;

@printf("Residual r: %s\n", Krylov.vec2str(r))
@printf("Relative residual norm ‖r‖: %8.1e\n", norm(r) / b_norm)
@printf("Solution x: %s\n", Krylov.vec2str(x))
@printf("Minimum-norm solution? %s\n", x ≈ [1.0; 1.0; 1.0; 0.0])

```

--------------------------------

### Set OpenBLAS/MKL Threads at Runtime

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/tips.md

Configure the number of threads for OpenBLAS or MKL at runtime. Ensure the number of threads is between 1 and the maximum available threads.

```shell
OPENBLAS_NUM_THREADS=N julia  # Set the number of OpenBLAS threads to N
MKL_NUM_THREADS=N julia       # Set the number of MKL threads to N
```

--------------------------------

### Jacobi Preconditioner with Symmlq

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/preconditioners.md

Demonstrates creating and using a Jacobi preconditioner (diagonal matrix of inverse absolute diagonal elements) with the Symmlq solver.

```julia
using Krylov
n, m = size(A)
d = [A[i,i] ≠ 0 ? 1 / abs(A[i,i]) : 1 for i=1:n]  # Jacobi preconditioner
P⁻¹ = diagm(d)
x, stats = symmlq(A, b, M=P⁻¹)
```

--------------------------------

### Solving Saddle-Point System with Default Parameters

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/examples/usymlqr.md

Solves a saddle-point system with the default parameters for `usymlqr`. It then verifies the solution by calculating the residual.

```julia
using LinearAlgebra, Printf, SparseArrays
using Krylov

# Identity matrix.
eye(n::Int) = sparse(1.0 * I, n, n)

# Saddle-point systems
n = m = 5
A = [2^(i/j)*j + (-1)^(i-j) * n*(i-1) for i = 1:n, j = 1:n]
b = ones(n)
D = diagm(0 => [2.0 * i for i = 1:n])
m, n = size(A)
c = -3*b

# [I   A] [x] = [b]
# [Aᴴ  0] [y]   [c]
(x, y, stats) = usymlqr(A, b, c)
K = [I A; A' zeros(n,n)]
d = [b; c]
r = d - K * [x; y]
resid = norm(r)
@printf("USYMLQR: Relative residual: %8.1e\n", resid)
```

--------------------------------

### Solving 1D Poisson Equation with FFT and Krylov

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/matrix_free.md

Demonstrates solving the 1D Poisson equation on a periodic domain using FFTW for matrix-free Laplacian application within a Krylov solver. This approach is efficient for structured grids and avoids explicit matrix storage. Requires FFTW, Krylov, and LinearAlgebra.

```julia
using FFTW, Krylov, LinearAlgebra

# Define the problem size and domain
n = 32768                         # Number of grid points (2^15)
L = 4π                            # Length of the domain
x = LinRange(0, L, n+1)[1:end-1]  # Periodic grid (excluding the last point)

# Define the source term f(x)
f = sin.(x)
```

--------------------------------

### Solve Linear System with craigmr

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/examples/craigmr.md

This snippet demonstrates solving a linear system Ax = b using the craigmr solver. It fetches a matrix from the SuiteSparseMatrixCollection, constructs a known solution, and then solves the system, reporting the relative residual and solution accuracy.

```julia
using Krylov, HarwellRutherfordBoeing, SuiteSparseMatrixCollection
using LinearAlgebra, Printf

ssmc = ssmc_db(verbose=false)
matrix = ssmc_matrices(ssmc, "HB", "wm1")
path = fetch_ssmc(matrix, format="RB")

A = RutherfordBoeingData(joinpath(path[1], "$(matrix.name[1]).rb")).data
(m, n) = size(A)
@printf("System size: %d rows and %d columns\n", m, n)

x_exact = A' * ones(m)
x_exact_norm = norm(x_exact)
x_exact /= x_exact_norm
b = A * x_exact
(x, y, stats) = craigmr(A, b)
show(stats)
resid = norm(A * x - b) / norm(b)
@printf("CRAIGMR: Relative residual: %7.1e\n", resid)
@printf("CRAIGMR: ‖x - x*‖₂: %7.1e\n", norm(x - x_exact))
@printf("CRAIGMR: %d iterations\n", length(stats.residuals))
```

--------------------------------

### BasicLU Maxvolbasis for Least-Norm with Craigmr

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/preconditioners.md

Employs BasicLU.jl's maxvolbasis to construct a preconditioner for solving minimum norm problems with the Craigmr solver.

```julia
using BasicLU, LinearOperators, Krylov

# Least-norm problem
m, n = size(A)
basis, B = maxvolbasis(A)
opA = LinearOperator(A)
B⁻¹ = LinearOperator(Float64, m, m, false, false, (y, v) -> (y .= v ; BasicLU.solve!(B, y, 'N')),
                                                  (y, v) -> (y .= v ; BasicLU.solve!(B, y, 'T')),
                                                  (y, v) -> (y .= v ; BasicLU.solve!(B, y, 'T')))

x, y, stats = craigmr(B⁻¹ * opA, B⁻¹ * b)  # min ‖x‖₂  s.t.  B⁻¹Ax = B⁻¹b
```

--------------------------------

### BiLQ Solver

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/solvers/unsymmetric.md

Provides the BiLQ (Biconjugate Low-rank) solver for non-Hermitian linear systems. Includes the function `bilq` and its in-place variant `bilq!`.

```APIDOC
## BiLQ Solver

### Description
Provides the BiLQ (Biconjugate Low-rank) solver for non-Hermitian linear systems. Includes the function `bilq` and its in-place variant `bilq!`.

### Functions
- `bilq(A, b, ...)`: Solves the linear system Ax = b using the BiLQ method.
- `bilq!(A, x, b, ...)`: Solves the linear system Ax = b in-place using the BiLQ method, updating x.

### Workspace
- `BilqWorkspace`: Structure to hold the state for the BiLQ solver.
```

--------------------------------

### Solve Dense Hermitian System on AMD GPU

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/gpu.md

Demonstrates solving a dense Hermitian linear system on an AMD GPU using minres. Requires AMDGPU.jl.

```julia
using Krylov, AMDGPU

if AMDGPU.functional()
  # CPU Arrays
  A_cpu = rand(ComplexF64, 20, 20)
  A_cpu = A_cpu + A_cpu'
  b_cpu = rand(ComplexF64, 20)

  A_gpu = ROCMatrix(A_cpu)
  b_gpu = ROCVector(b_cpu)

  # Solve a dense Hermitian system on an AMD GPU
  x, stats = minres(A_gpu, b_gpu)
end
```

--------------------------------

### Diagonal Preconditioner with Minres

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/preconditioners.md

Shows how to construct a diagonal preconditioner using the inverse of the column norms and apply it with the Minres solver.

```julia
using Krylov
n, m = size(A)
d = [1 / norm(A[:,i]) for i=1:m]  # diagonal preconditioner
P⁻¹ = diagm(d)
x, stats = minres(A, b, M=P⁻¹)
```

--------------------------------

### Out-of-Place Krylov Solver Interface

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/generic_interface.md

Demonstrates the out-of-place interface for Krylov methods like CG, CR, and CAR. This is suitable when you don't need to modify an existing workspace or want a simpler API for direct solving.

```julia
using Krylov, SparseArrays, LinearAlgebra

# Define a symmetric positive definite matrix A and a right-hand side vector b
n = 1000
A = sprandn(n, n, 0.005)
A = A * A' + I
b = randn(n)

# Out-of-place interface
for method in (:cg, :cr, :car)
    x, stats = krylov_solve(Val(method), A, b)
    r = b - A * x
    println("Residual norm for $(method): ", norm(r))
end
```

--------------------------------

### Solve Dense Least-Squares Problem on Intel GPU

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/gpu.md

Demonstrates solving a dense least-squares problem on an Intel GPU using lsqr. Requires oneAPI.jl.

```julia
using Krylov, oneAPI

if oneAPI.functional()
  T = Float32  # oneAPI.jl also works with ComplexF32
  m = 20
  n = 10

  # CPU Arrays
  A_cpu = rand(T, m, n)
  b_cpu = rand(T, m)

  # GPU Arrays
  A_gpu = oneMatrix(A_cpu)
  b_gpu = oneVector(b_cpu)

  # Solve a dense least-squares problem on an Intel GPU
  x, stats = lsqr(A_gpu, b_gpu)
end
```

--------------------------------

### Solve Dense System on NVIDIA GPU

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/gpu.md

Demonstrates solving a dense linear system using Krylov.jl on an NVIDIA GPU. Ensure CUDA is functional and convert CPU arrays to CuMatrix and CuVector.

```julia
using CUDA, Krylov

if CUDA.functional()
  # CPU Arrays
  A_cpu = rand(20, 20)
  b_cpu = rand(20)

  # GPU Arrays
  A_gpu = CuMatrix(A_cpu)
  b_gpu = CuVector(b_cpu)

  # Solve a square and dense system on an Nivida GPU
  x, stats = bilq(A_gpu, b_gpu)
end
```

--------------------------------

### Solving a Regularized Least Squares Problem with LSQR

Source: https://github.com/juliasmoothoptimizers/krylov.jl/blob/main/docs/src/examples/lsqr.md

This snippet demonstrates solving a linear least squares problem Ax = b with Tikhonov regularization (λ) using the `lsqr` function. It fetches a matrix from the SuiteSparse Matrix Collection, sets up the problem, and then solves it. The `atol` and `btol` parameters are set to 0.0 for demonstration purposes, meaning convergence is determined by the iteration limit or other internal criteria. The results, including the relative residual and the norm of the solution vector, are printed.

```julia
using MatrixMarket, SuiteSparseMatrixCollection
using Krylov, LinearOperators
using LinearAlgebra, Printf

ssmc = ssmc_db(verbose=false)
matrix = ssmc_matrices(ssmc, "HB", "illc1033")
path = fetch_ssmc(matrix, format="MM")

A = MatrixMarket.mmread(joinpath(path[1], "$(matrix.name[1]).mtx"))
b = MatrixMarket.mmread(joinpath(path[1], "$(matrix.name[1])_b.mtx"))[:]
(m, n) = size(A)
@printf("System size: %d rows and %d columns\n", m, n)

# Define a regularization parameter.
λ = 1.0e-3

(x, stats) = lsqr(A, b, λ=λ, atol=0.0, btol=0.0)
show(stats)
resid = norm(A' * (A * x - b) + λ * x) / norm(b)
@printf("LSQR: Relative residual: %8.1e\n", resid)
@printf("LSQR: ‖x‖: %8.1e\n", norm(x))
```