## In-place methods

All solvers in Krylov.jl have an in-place variant implemented in a method whose name ends with `!`

. A workspace (`KrylovSolver`

), which contains the storage needed by a Krylov method, can be used to solve multiple linear systems with the same dimensions and the same floating-point precision. The section storage requirements specifies the memory needed for each Krylov method. Each `KrylovSolver`

has two constructors:

```
XyzSolver(A, b)
XyzSolver(m, n, S)
```

`Xyz`

is the name of the Krylov method with lowercase letters except its first one (`Cg`

, `Minres`

, `Lsmr`

, `Bicgstab`

, ...). Given an operator `A`

and a right-hand side `b`

, you can create a `KrylovSolver`

based on the size of `A`

and the type of `b`

or explicitly give the dimensions `(m, n)`

and the storage type `S`

.

For example, use `S = Vector{Float64}`

if you want to solve linear systems in double precision on the CPU and `S = CuVector{Float32}`

if you want to solve linear systems in single precision on an Nvidia GPU.

`DiomSolver`

, `FomSolver`

, `DqgmresSolver`

, `GmresSolver`

, `BlockGmresSolver`

, `FgmresSolver`

, `GpmrSolver`

, `CgLanczosShiftSolver`

and `CglsLanczosShiftSolver`

require an additional argument (`memory`

or `nshifts`

).

The workspace is always the first argument of the in-place methods:

```
minres_solver = MinresSolver(n, n, Vector{Float64})
minres!(minres_solver, A1, b1)
dqgmres_solver = DqgmresSolver(n, n, memory, Vector{BigFloat})
dqgmres!(dqgmres_solver, A2, b2)
lsqr_solver = LsqrSolver(m, n, CuVector{Float32})
lsqr!(lsqr_solver, A3, b3)
```

A generic function `solve!`

is also available and dispatches to the appropriate Krylov method.

`Krylov.solve!`

— Function`solve!(solver, args...; kwargs...)`

Generic function that dispatches to the appropriate in-place Krylov method based on the type of `solver`

.

In-place methods return an updated `solver`

workspace. Solutions and statistics can be recovered via `solver.x`

, `solver.y`

and `solver.stats`

. Functions `solution`

, `statistics`

and `results`

can be also used.

`Krylov.nsolution`

— Function`nsolution(solver)`

Return the number of outputs of `solution(solver)`

.

`Krylov.issolved`

— Function`issolved(solver)`

Return a boolean that determines whether the Krylov method associated to `solver`

succeeded.

`Krylov.solution`

— Function`solution(solver)`

Return the solution(s) stored in the `solver`

. Optionally you can specify which solution you want to recover, `solution(solver, 1)`

returns `x`

and `solution(solver, 2)`

returns `y`

.

`Krylov.statistics`

— Function`statistics(solver)`

Return the statistics stored in the `solver`

.

`Krylov.results`

— Function`results(solver)`

Return a tuple containing the solution(s) and the statistics associated with the `solver`

. Allows retrieving the output arguments of an out-of-place method from the in-place method.

For example, instead of `x, stats = cg(A, b)`

, you can use:

```
solver = CgSolver(A, b)
cg!(solver, A, b)
x, stats = results(solver)
```

## Examples

We illustrate the use of in-place Krylov solvers with two well-known optimization methods. The details of the optimization methods are described in the section about Factorization-free operators.

### Example 1: Newton's method for convex optimization without linesearch

```
using Krylov
function newton(∇f, ∇²f, x₀; itmax = 200, tol = 1e-8)
n = length(x₀)
x = copy(x₀)
gx = ∇f(x)
iter = 0
S = typeof(x)
solver = CgSolver(n, n, S)
Δx = solver.x
solved = false
tired = false
while !(solved || tired)
Hx = ∇²f(x) # Compute ∇²f(xₖ)
cg!(solver, Hx, -gx) # Solve ∇²f(xₖ)Δx = -∇f(xₖ)
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
```

### Example 2: The Gauss-Newton method for nonlinear least squares without linesearch

```
using Krylov
function gauss_newton(F, JF, x₀; itmax = 200, tol = 1e-8)
n = length(x₀)
x = copy(x₀)
Fx = F(x)
m = length(Fx)
iter = 0
S = typeof(x)
solver = LsmrSolver(m, n, S)
Δx = solver.x
solved = false
tired = false
while !(solved || tired)
Jx = JF(x) # Compute J(xₖ)
lsmr!(solver, Jx, -Fx) # Minimize ‖J(xₖ)Δx + F(xₖ)‖
x = x + Δx # Update xₖ₊₁ = xₖ + Δx
Fx_old = Fx # F(xₖ)
Fx = F(x) # F(xₖ₊₁)
iter += 1
solved = norm(Fx - Fx_old) / norm(Fx) ≤ tol
tired = iter ≥ itmax
end
return x
end
```