diff --git a/src/bicgstab.jl b/src/bicgstab.jl
index 588f64838..b2421f1f7 100644
--- a/src/bicgstab.jl
+++ b/src/bicgstab.jl
@@ -16,14 +16,19 @@
 export bicgstab, bicgstab!
 
 """
-    (x, stats) = bicgstab(A, b::AbstractVector{FC}; c::AbstractVector{FC}=b,
-                          M=I, N=I, atol::T=√eps(T), rtol::T=√eps(T),
-                          itmax::Int=0, verbose::Int=0, history::Bool=false,
+    (x, stats) = bicgstab(A, b::AbstractVector{FC};
+                          c::AbstractVector{FC}=b, M=I, N=I,
+                          atol::T=√eps(T), rtol::T=√eps(T), itmax::Int=0,
+                          verbose::Int=0, history::Bool=false,
                           ldiv::Bool=false, callback=solver->false)
 
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, stats) = bicgstab(A, b, x0::AbstractVector; kwargs...)
+
+BICGSTAB can be warm-started from an initial guess `x0` where `kwargs` are the same keyword arguments as above.
+
 Solve the square linear system Ax = b of size n using BICGSTAB.
 BICGSTAB requires two initial vectors `b` and `c`.
 The relation `bᴴc ≠ 0` must be satisfied and by default `c = b`.
@@ -42,12 +47,6 @@ Information will be displayed every `verbose` iterations.
 
 This implementation allows a left preconditioner `M` and a right preconditioner `N`.
 
-BICGSTAB can be warm-started from an initial guess `x0` with
-
-    (x, stats) = bicgstab(A, b, x0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -56,6 +55,10 @@ and `false` otherwise.
 * `A`: a linear operator that models a matrix of dimension n;
 * `b`: a vector of length n.
 
+#### Optional argument
+
+* `x0`: a vector of length n that represents an initial guess of the solution x.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
diff --git a/src/bilq.jl b/src/bilq.jl
index 002858943..163b6339d 100644
--- a/src/bilq.jl
+++ b/src/bilq.jl
@@ -13,16 +13,20 @@
 export bilq, bilq!
 
 """
-    (x, stats) = bilq(A, b::AbstractVector{FC}; c::AbstractVector{FC}=b,
-                      atol::T=√eps(T), rtol::T=√eps(T), transfer_to_bicg::Bool=true,
-                      itmax::Int=0, verbose::Int=0, history::Bool=false,
-                      callback=solver->false)
+    (x, stats) = bilq(A, b::AbstractVector{FC};
+                      c::AbstractVector{FC}=b, atol::T=√eps(T),
+                      rtol::T=√eps(T), transfer_to_bicg::Bool=true,
+                      itmax::Int=0, verbose::Int=0,
+                      history::Bool=false, callback=solver->false)
 
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
-Solve the square linear system Ax = b of size n using BiLQ.
+    (x, stats) = bilq(A, b, x0::AbstractVector; kwargs...)
+
+BiLQ can be warm-started from an initial guess `x0` where `kwargs` are the same keyword arguments as above.
 
+Solve the square linear system Ax = b of size n using BiLQ.
 BiLQ is based on the Lanczos biorthogonalization process and requires two initial vectors `b` and `c`.
 The relation `bᴴc ≠ 0` must be satisfied and by default `c = b`.
 When `A` is symmetric and `b = c`, BiLQ is equivalent to SYMMLQ.
@@ -30,12 +34,6 @@ When `A` is symmetric and `b = c`, BiLQ is equivalent to SYMMLQ.
 An option gives the possibility of transferring to the BiCG point,
 when it exists. The transfer is based on the residual norm.
 
-BiLQ can be warm-started from an initial guess `x0` with
-
-    (x, stats) = bilq(A, b, x0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -44,6 +42,10 @@ and `false` otherwise.
 * `A`: a linear operator that models a matrix of dimension n;
 * `b`: a vector of length n.
 
+#### Optional argument
+
+* `x0`: a vector of length n that represents an initial guess of the solution x.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
diff --git a/src/bilqr.jl b/src/bilqr.jl
index 4677e9e9d..5e45cb372 100644
--- a/src/bilqr.jl
+++ b/src/bilqr.jl
@@ -21,6 +21,10 @@ export bilqr, bilqr!
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, y, stats) = bilqr(A, b, c, x0::AbstractVector, y0::AbstractVector; kwargs...)
+
+BiLQR can be warm-started from initial guesses `x0` and `y0` where `kwargs` are the same keyword arguments as above.
+
 Combine BiLQ and QMR to solve adjoint systems.
 
     [0  A] [y] = [b]
@@ -33,12 +37,6 @@ QMR is used for solving dual system `Aᴴy = c` of size n.
 An option gives the possibility of transferring from the BiLQ point to the
 BiCG point, when it exists. The transfer is based on the residual norm.
 
-BiLQR can be warm-started from initial guesses `x0` and `y0` with
-
-    (x, y, stats) = bilqr(A, b, c, x0, y0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -48,11 +46,16 @@ and `false` otherwise.
 * `b`: a vector of length n;
 * `c`: a vector of length n.
 
+#### Optional arguments
+
+* `x0`: a vector of length n that represents an initial guess of the solution x;
+* `y0`: a vector of length n that represents an initial guess of the solution y.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
 * `y`: a dense vector of length n;
-* `stats`: statistics collected on the run in a [`AdjointStats`](@ref) structure.
+* `stats`: statistics collected on the run in an [`AdjointStats`](@ref) structure.
 
 #### Reference
 
diff --git a/src/cg.jl b/src/cg.jl
index 09b91e64c..1c4b39cb5 100644
--- a/src/cg.jl
+++ b/src/cg.jl
@@ -26,6 +26,10 @@ export cg, cg!
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, stats) = cg(A, b, x0::AbstractVector; kwargs...)
+
+CG can be warm-started from an initial guess `x0` where `kwargs` are the same keyword arguments as above.
+
 The conjugate gradient method to solve the Hermitian linear system Ax = b of size n.
 
 The method does _not_ abort if A is not definite.
@@ -36,12 +40,6 @@ M also indicates the weighted norm in which residuals are measured.
 
 If `itmax=0`, the default number of iterations is set to `2 * n`.
 
-CG can be warm-started from an initial guess `x0` with
-
-    (x, stats) = cg(A, b, x0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -50,6 +48,10 @@ and `false` otherwise.
 * `A`: a linear operator that models a Hermitian positive definite matrix of dimension n;
 * `b`: a vector of length n.
 
+#### Optional argument
+
+* `x0`: a vector of length n that represents an initial guess of the solution x.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
diff --git a/src/cg_lanczos.jl b/src/cg_lanczos.jl
index b49cdd726..3efdcd90e 100644
--- a/src/cg_lanczos.jl
+++ b/src/cg_lanczos.jl
@@ -22,6 +22,10 @@ export cg_lanczos, cg_lanczos!
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, stats) = cg_lanczos(A, b, x0::AbstractVector; kwargs...)
+
+CG-LANCZOS can be warm-started from an initial guess `x0` where `kwargs` are the same keyword arguments as above.
+
 The Lanczos version of the conjugate gradient method to solve the
 Hermitian linear system Ax = b of size n.
 
@@ -30,12 +34,6 @@ The method does _not_ abort if A is not definite.
 A preconditioner M may be provided in the form of a linear operator and is
 assumed to be Hermitian and positive definite.
 
-CG-LANCZOS can be warm-started from an initial guess `x0` with
-
-    (x, stats) = cg_lanczos(A, b, x0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -44,6 +42,10 @@ and `false` otherwise.
 * `A`: a linear operator that models a Hermitian matrix of dimension n;
 * `b`: a vector of length n.
 
+#### Optional argument
+
+* `x0`: a vector of length n that represents an initial guess of the solution x.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
diff --git a/src/cgs.jl b/src/cgs.jl
index 9b4eb695f..4f770c5f3 100644
--- a/src/cgs.jl
+++ b/src/cgs.jl
@@ -11,14 +11,18 @@
 export cgs, cgs!
 
 """
-    (x, stats) = cgs(A, b::AbstractVector{FC}; c::AbstractVector{FC}=b,
-                     M=I, N=I, atol::T=√eps(T), rtol::T=√eps(T),
-                     itmax::Int=0, verbose::Int=0, history::Bool=false,
-                     ldiv::Bool=false, callback=solver->false)
+    (x, stats) = cgs(A, b::AbstractVector{FC};
+                     c::AbstractVector{FC}=b, M=I, N=I, atol::T=√eps(T),
+                     rtol::T=√eps(T), itmax::Int=0, verbose::Int=0,
+                     history::Bool=false, ldiv::Bool=false, callback=solver->false)
 
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, stats) = cgs(A, b, x0::AbstractVector; kwargs...)
+
+CGS can be warm-started from an initial guess `x0` where `kwargs` are the same keyword arguments as above.
+
 Solve the consistent linear system Ax = b of size n using CGS.
 CGS requires two initial vectors `b` and `c`.
 The relation `bᴴc ≠ 0` must be satisfied and by default `c = b`.
@@ -40,12 +44,6 @@ TFQMR and BICGSTAB were developed to remedy this difficulty.»
 
 This implementation allows a left preconditioner M and a right preconditioner N.
 
-CGS can be warm-started from an initial guess `x0` with
-
-    (x, stats) = cgs(A, b, x0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -54,6 +52,10 @@ and `false` otherwise.
 * `A`: a linear operator that models a matrix of dimension n;
 * `b`: a vector of length n.
 
+#### Optional argument
+
+* `x0`: a vector of length n that represents an initial guess of the solution x.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
diff --git a/src/cr.jl b/src/cr.jl
index 79f2cd289..4bb104534 100644
--- a/src/cr.jl
+++ b/src/cr.jl
@@ -16,13 +16,18 @@ export cr, cr!
 
 """
     (x, stats) = cr(A, b::AbstractVector{FC};
-                    M=I, atol::T=√eps(T), rtol::T=√eps(T), γ::T=√eps(T), itmax::Int=0,
-                    radius::T=zero(T), verbose::Int=0, linesearch::Bool=false, history::Bool=false,
+                    M=I, atol::T=√eps(T), rtol::T=√eps(T), γ::T=√eps(T),
+                    itmax::Int=0, radius::T=zero(T), verbose::Int=0,
+                    linesearch::Bool=false, history::Bool=false,
                     ldiv::Bool=false, callback=solver->false)
 
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, stats) = cr(A, b, x0::AbstractVector; kwargs...)
+
+CR can be warm-started from an initial guess `x0` where `kwargs` are the same keyword arguments as above.
+
 A truncated version of Stiefel’s Conjugate Residual method to solve the Hermitian linear system Ax = b
 of size n or the least-squares problem min ‖b - Ax‖ if A is singular.
 The matrix A must be Hermitian semi-definite.
@@ -34,12 +39,6 @@ In a linesearch context, 'linesearch' must be set to 'true'.
 
 If `itmax=0`, the default number of iterations is set to `2 * n`.
 
-CR can be warm-started from an initial guess `x0` with
-
-    (x, stats) = cr(A, b, x0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -48,6 +47,10 @@ and `false` otherwise.
 * `A`: a linear operator that models a Hermitian positive definite matrix of dimension n;
 * `b`: a vector of length n.
 
+#### Optional argument
+
+* `x0`: a vector of length n that represents an initial guess of the solution x.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
diff --git a/src/diom.jl b/src/diom.jl
index 58986ab47..23b6d48ca 100644
--- a/src/diom.jl
+++ b/src/diom.jl
@@ -11,15 +11,19 @@
 export diom, diom!
 
 """
-    (x, stats) = diom(A, b::AbstractVector{FC}; memory::Int=20,
-                      M=I, N=I, atol::T=√eps(T), rtol::T=√eps(T),
-                      reorthogonalization::Bool=false, itmax::Int=0,
-                      verbose::Int=0, history::Bool=false,
+    (x, stats) = diom(A, b::AbstractVector{FC};
+                      memory::Int=20, M=I, N=I, atol::T=√eps(T),
+                      rtol::T=√eps(T), reorthogonalization::Bool=false,
+                      itmax::Int=0, verbose::Int=0, history::Bool=false,
                       ldiv::Bool=false, callback=solver->false)
 
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, stats) = diom(A, b, x0::AbstractVector; kwargs...)
+
+DIOM can be warm-started from an initial guess `x0` where `kwargs` are the same keyword arguments as above.
+
 Solve the consistent linear system Ax = b of size n using DIOM.
 
 DIOM only orthogonalizes the new vectors of the Krylov basis against the `memory` most recent vectors.
@@ -34,12 +38,6 @@ and indefinite systems of linear equations can be handled by this single algorit
 
 This implementation allows a left preconditioner M and a right preconditioner N.
 
-DIOM can be warm-started from an initial guess `x0` with
-
-    (x, stats) = diom(A, b, x0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -48,6 +46,10 @@ and `false` otherwise.
 * `A`: a linear operator that models a matrix of dimension n;
 * `b`: a vector of length n.
 
+#### Optional argument
+
+* `x0`: a vector of length n that represents an initial guess of the solution x.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
diff --git a/src/dqgmres.jl b/src/dqgmres.jl
index 2dd5b0843..2d428b87c 100644
--- a/src/dqgmres.jl
+++ b/src/dqgmres.jl
@@ -11,15 +11,19 @@
 export dqgmres, dqgmres!
 
 """
-    (x, stats) = dqgmres(A, b::AbstractVector{FC}; memory::Int=20,
-                         M=I, N=I, atol::T=√eps(T), rtol::T=√eps(T),
-                         reorthogonalization::Bool=false, itmax::Int=0,
-                         verbose::Int=0, history::Bool=false,
+    (x, stats) = dqgmres(A, b::AbstractVector{FC};
+                         memory::Int=20, M=I, N=I, atol::T=√eps(T),
+                         rtol::T=√eps(T), reorthogonalization::Bool=false,
+                         itmax::Int=0, verbose::Int=0, history::Bool=false,
                          ldiv::Bool=false, callback=solver->false)
 
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, stats) = dqgmres(A, b, x0::AbstractVector; kwargs...)
+
+DQGMRES can be warm-started from an initial guess `x0` where `kwargs` are the same keyword arguments as above.
+
 Solve the consistent linear system Ax = b of size n using DQGMRES.
 
 DQGMRES algorithm is based on the incomplete Arnoldi orthogonalization process
@@ -34,12 +38,6 @@ Partial reorthogonalization is available with the `reorthogonalization` option.
 
 This implementation allows a left preconditioner M and a right preconditioner N.
 
-DQGMRES can be warm-started from an initial guess `x0` with
-
-    (x, stats) = dqgmres(A, b, x0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -48,6 +46,10 @@ and `false` otherwise.
 * `A`: a linear operator that models a matrix of dimension n;
 * `b`: a vector of length n.
 
+#### Optional argument
+
+* `x0`: a vector of length n that represents an initial guess of the solution x.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
diff --git a/src/fgmres.jl b/src/fgmres.jl
index b79c197ba..9fc53408a 100644
--- a/src/fgmres.jl
+++ b/src/fgmres.jl
@@ -11,8 +11,8 @@
 export fgmres, fgmres!
 
 """
-    (x, stats) = fgmres(A, b::AbstractVector{FC}; memory::Int=20,
-                        M=I, N=I, atol::T=√eps(T), rtol::T=√eps(T),
+    (x, stats) = fgmres(A, b::AbstractVector{FC};
+                        memory::Int=20, M=I, N=I, atol::T=√eps(T), rtol::T=√eps(T),
                         reorthogonalization::Bool=false, itmax::Int=0,
                         restart::Bool=false, verbose::Int=0, history::Bool=false,
                         ldiv::Bool=false, callback=solver->false)
@@ -20,6 +20,10 @@ export fgmres, fgmres!
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, stats) = fgmres(A, b, x0::AbstractVector; kwargs...)
+
+FGMRES can be warm-started from an initial guess `x0` where `kwargs` are the same keyword arguments as above.
+
 Solve the linear system Ax = b of size n using FGMRES.
 
 FGMRES computes a sequence of approximate solutions with minimum residual.
@@ -37,12 +41,6 @@ If `restart = true`, the restarted version FGMRES(k) is used with `k = memory`.
 If `restart = false`, the parameter `memory` should be used as a hint of the number of iterations to limit dynamic memory allocations.
 More storage will be allocated only if the number of iterations exceeds `memory`.
 
-FGMRES can be warm-started from an initial guess `x0` with
-
-    (x, stats) = fgmres(A, b, x0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -51,6 +49,10 @@ and `false` otherwise.
 * `A`: a linear operator that models a matrix of dimension n;
 * `b`: a vector of length n.
 
+#### Optional argument
+
+* `x0`: a vector of length n that represents an initial guess of the solution x.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
diff --git a/src/fom.jl b/src/fom.jl
index 4ea3fce92..fdd99708b 100644
--- a/src/fom.jl
+++ b/src/fom.jl
@@ -11,8 +11,8 @@
 export fom, fom!
 
 """
-    (x, stats) = fom(A, b::AbstractVector{FC}; memory::Int=20,
-                     M=I, N=I, atol::T=√eps(T), rtol::T=√eps(T),
+    (x, stats) = fom(A, b::AbstractVector{FC};
+                     memory::Int=20, M=I, N=I, atol::T=√eps(T), rtol::T=√eps(T),
                      reorthogonalization::Bool=false, itmax::Int=0,
                      restart::Bool=false, verbose::Int=0, history::Bool=false,
                      ldiv::Bool=false, callback=solver->false)
@@ -20,6 +20,10 @@ export fom, fom!
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, stats) = fom(A, b, x0::AbstractVector; kwargs...)
+
+FOM can be warm-started from an initial guess `x0` where `kwargs` are the same keyword arguments as above.
+
 Solve the linear system Ax = b of size n using FOM.
 
 FOM algorithm is based on the Arnoldi process and a Galerkin condition.
@@ -31,12 +35,6 @@ If `restart = true`, the restarted version FOM(k) is used with `k = memory`.
 If `restart = false`, the parameter `memory` should be used as a hint of the number of iterations to limit dynamic memory allocations.
 More storage will be allocated only if the number of iterations exceeds `memory`.
 
-FOM can be warm-started from an initial guess `x0` with
-
-    (x, stats) = fom(A, b, x0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -45,6 +43,10 @@ and `false` otherwise.
 * `A`: a linear operator that models a matrix of dimension n;
 * `b`: a vector of length n.
 
+#### Optional argument
+
+* `x0`: a vector of length n that represents an initial guess of the solution x.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
diff --git a/src/gmres.jl b/src/gmres.jl
index ac425054f..1af93328b 100644
--- a/src/gmres.jl
+++ b/src/gmres.jl
@@ -11,8 +11,8 @@
 export gmres, gmres!
 
 """
-    (x, stats) = gmres(A, b::AbstractVector{FC}; memory::Int=20,
-                       M=I, N=I, atol::T=√eps(T), rtol::T=√eps(T),
+    (x, stats) = gmres(A, b::AbstractVector{FC};
+                       memory::Int=20, M=I, N=I, atol::T=√eps(T), rtol::T=√eps(T),
                        reorthogonalization::Bool=false, itmax::Int=0,
                        restart::Bool=false, verbose::Int=0, history::Bool=false,
                        ldiv::Bool=false, callback=solver->false)
@@ -20,6 +20,10 @@ export gmres, gmres!
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, stats) = gmres(A, b, x0::AbstractVector; kwargs...)
+
+GMRES can be warm-started from an initial guess `x0` where `kwargs` are the same keyword arguments as above.
+
 Solve the linear system Ax = b of size n using GMRES.
 
 GMRES algorithm is based on the Arnoldi process and computes a sequence of approximate solutions with the minimum residual.
@@ -31,12 +35,6 @@ If `restart = true`, the restarted version GMRES(k) is used with `k = memory`.
 If `restart = false`, the parameter `memory` should be used as a hint of the number of iterations to limit dynamic memory allocations.
 More storage will be allocated only if the number of iterations exceeds `memory`.
 
-GMRES can be warm-started from an initial guess `x0` with
-
-    (x, stats) = gmres(A, b, x0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -45,6 +43,10 @@ and `false` otherwise.
 * `A`: a linear operator that models a matrix of dimension n;
 * `b`: a vector of length n.
 
+#### Optional argument
+
+* `x0`: a vector of length n that represents an initial guess of the solution x.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
diff --git a/src/gpmr.jl b/src/gpmr.jl
index f94f6e4ac..fac2270ba 100644
--- a/src/gpmr.jl
+++ b/src/gpmr.jl
@@ -12,16 +12,20 @@
 export gpmr, gpmr!
 
 """
-    (x, y, stats) = gpmr(A, B, b::AbstractVector{FC}, c::AbstractVector{FC}; memory::Int=20,
-                         C=I, D=I, E=I, F=I, atol::T=√eps(T), rtol::T=√eps(T),
-                         gsp::Bool=false, reorthogonalization::Bool=false,
-                         itmax::Int=0, λ::FC=one(FC), μ::FC=one(FC),
-                         verbose::Int=0, history::Bool=false,
-                         ldiv::Bool=false, callback=solver->false)
+    (x, y, stats) = gpmr(A, B, b::AbstractVector{FC}, c::AbstractVector{FC};
+                         memory::Int=20, C=I, D=I, E=I, F=I,
+                         atol::T=√eps(T), rtol::T=√eps(T), gsp::Bool=false,
+                         reorthogonalization::Bool=false, itmax::Int=0,
+                         λ::FC=one(FC), μ::FC=one(FC), verbose::Int=0,
+                         history::Bool=false, ldiv::Bool=false, callback=solver->false)
 
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, y, stats) = gpmr(A, B, b, c, x0::AbstractVector, y0::AbstractVector; kwargs...)
+
+GPMR can be warm-started from initial guesses `x0` and `y0` where `kwargs` are the same keyword arguments as above.
+
 Given matrices `A` of dimension m × n and `B` of dimension n × m,
 GPMR solves the unsymmetric partitioned linear system
 
@@ -59,12 +63,6 @@ Full reorthogonalization is available with the `reorthogonalization` option.
 Additional details can be displayed if verbose mode is enabled (verbose > 0).
 Information will be displayed every `verbose` iterations.
 
-GPMR can be warm-started from initial guesses `x0` and `y0` with
-
-    (x, y, stats) = gpmr(A, B, b, c, x0, y0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -75,6 +73,11 @@ and `false` otherwise.
 * `b`: a vector of length m;
 * `c`: a vector of length n.
 
+#### Optional arguments
+
+* `x0`: a vector of length m that represents an initial guess of the solution x;
+* `y0`: a vector of length n that represents an initial guess of the solution y.
+
 #### Output arguments
 
 * `x`: a dense vector of length m;
diff --git a/src/minres.jl b/src/minres.jl
index f2814403c..c63312f77 100644
--- a/src/minres.jl
+++ b/src/minres.jl
@@ -35,6 +35,10 @@ export minres, minres!
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, stats) = minres(A, b, x0::AbstractVector; kwargs...)
+
+MINRES can be warm-started from an initial guess `x0` where `kwargs` are the same keyword arguments as above.
+
 Solve the shifted linear least-squares problem
 
     minimize ‖b - (A + λI)x‖₂²
@@ -55,12 +59,6 @@ MINRES produces monotonic residuals ‖r‖₂ and optimality residuals ‖Aᴴr
 A preconditioner M may be provided in the form of a linear operator and is
 assumed to be Hermitian and positive definite.
 
-MINRES can be warm-started from an initial guess `x0` with
-
-    (x, stats) = minres(A, b, x0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -69,6 +67,10 @@ and `false` otherwise.
 * `A`: a linear operator that models a Hermitian matrix of dimension n;
 * `b`: a vector of length n.
 
+#### Optional argument
+
+* `x0`: a vector of length n that represents an initial guess of the solution x.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
diff --git a/src/minres_qlp.jl b/src/minres_qlp.jl
index a1c9fc01b..ba27efb44 100644
--- a/src/minres_qlp.jl
+++ b/src/minres_qlp.jl
@@ -26,6 +26,10 @@ export minres_qlp, minres_qlp!
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, stats) = minres_qlp(A, b, x0::AbstractVector; kwargs...)
+
+MINRES-QLP can be warm-started from an initial guess `x0` where `kwargs` are the same keyword arguments as above.
+
 MINRES-QLP is the only method based on the Lanczos process that returns the minimum-norm
 solution on singular inconsistent systems (A + λI)x = b of size n, where λ is a shift parameter.
 It is significantly more complex but can be more reliable than MINRES when A is ill-conditioned.
@@ -34,12 +38,6 @@ A preconditioner M may be provided in the form of a linear operator and is
 assumed to be Hermitian and positive definite.
 M also indicates the weighted norm in which residuals are measured.
 
-MINRES-QLP can be warm-started from an initial guess `x0` with
-
-    (x, stats) = minres_qlp(A, b, x0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -48,6 +46,10 @@ and `false` otherwise.
 * `A`: a linear operator that models a Hermitian matrix of dimension n;
 * `b`: a vector of length n.
 
+#### Optional argument
+
+* `x0`: a vector of length n that represents an initial guess of the solution x.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
diff --git a/src/qmr.jl b/src/qmr.jl
index cf3328649..a8db7e978 100644
--- a/src/qmr.jl
+++ b/src/qmr.jl
@@ -21,26 +21,24 @@
 export qmr, qmr!
 
 """
-    (x, stats) = qmr(A, b::AbstractVector{FC}; c::AbstractVector{FC}=b,
-                     atol::T=√eps(T), rtol::T=√eps(T),
-                     itmax::Int=0, verbose::Int=0, history::Bool=false,
-                     callback=solver->false)
+    (x, stats) = qmr(A, b::AbstractVector{FC};
+                     c::AbstractVector{FC}=b, atol::T=√eps(T),
+                     rtol::T=√eps(T), itmax::Int=0, verbose::Int=0,
+                     history::Bool=false, callback=solver->false)
 
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, stats) = qmr(A, b, x0::AbstractVector; kwargs...)
+
+QMR can be warm-started from an initial guess `x0` where `kwargs` are the same keyword arguments as above.
+
 Solve the square linear system Ax = b of size n using QMR.
 
 QMR is based on the Lanczos biorthogonalization process and requires two initial vectors `b` and `c`.
 The relation `bᴴc ≠ 0` must be satisfied and by default `c = b`.
 When `A` is symmetric and `b = c`, QMR is equivalent to MINRES.
 
-QMR can be warm-started from an initial guess `x0` with
-
-    (x, stats) = qmr(A, b, x0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -49,6 +47,10 @@ and `false` otherwise.
 * `A`: a linear operator that models a matrix of dimension n;
 * `b`: a vector of length n.
 
+#### Optional argument
+
+* `x0`: a vector of length n that represents an initial guess of the solution x.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
diff --git a/src/symmlq.jl b/src/symmlq.jl
index 10d903049..9d28dc07b 100644
--- a/src/symmlq.jl
+++ b/src/symmlq.jl
@@ -13,8 +13,8 @@ export symmlq, symmlq!
 
 
 """
-    (x, stats) = symmlq(A, b::AbstractVector{FC}; window::Int=0,
-                        M=I, λ::T=zero(T), transfer_to_cg::Bool=true,
+    (x, stats) = symmlq(A, b::AbstractVector{FC};
+                        window::Int=0, M=I, λ::T=zero(T), transfer_to_cg::Bool=true,
                         λest::T=zero(T), atol::T=√eps(T), rtol::T=√eps(T),
                         etol::T=√eps(T), itmax::Int=0, conlim::T=1/√eps(T),
                         verbose::Int=0, history::Bool=false,
@@ -23,6 +23,10 @@ export symmlq, symmlq!
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, stats) = symmlq(A, b, x0::AbstractVector; kwargs...)
+
+SYMMLQ can be warm-started from an initial guess `x0` where `kwargs` are the same keyword arguments as above
+
 Solve the shifted linear system
 
     (A + λI) x = b
@@ -35,12 +39,6 @@ SYMMLQ produces monotonic errors ‖x* - x‖₂.
 A preconditioner M may be provided in the form of a linear operator and is
 assumed to be Hermitian and positive definite.
 
-SYMMLQ can be warm-started from an initial guess `x0` with
-
-    (x, stats) = symmlq(A, b, x0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -49,6 +47,10 @@ and `false` otherwise.
 * `A`: a linear operator that models a Hermitian matrix of dimension n;
 * `b`: a vector of length n.
 
+#### Optional argument
+
+* `x0`: a vector of length n that represents an initial guess of the solution x.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
diff --git a/src/tricg.jl b/src/tricg.jl
index f4a3c35f2..1860d8492 100644
--- a/src/tricg.jl
+++ b/src/tricg.jl
@@ -22,6 +22,10 @@ export tricg, tricg!
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, y, stats) = tricg(A, b, c, x0::AbstractVector, y0::AbstractVector; kwargs...)
+
+TriCG can be warm-started from initial guesses `x0` and `y0` where `kwargs` are the same keyword arguments as above.
+
 Given a matrix `A` of dimension m × n, TriCG solves the symmetric linear system
 
     [ τE    A ] [ x ] = [ b ]
@@ -53,12 +57,6 @@ TriCG stops when `itmax` iterations are reached or when `‖rₖ‖ ≤ atol + 
 Additional details can be displayed if verbose mode is enabled (verbose > 0).
 Information will be displayed every `verbose` iterations.
 
-TriCG can be warm-started from initial guesses `x0` and `y0` with
-
-    (x, y, stats) = tricg(A, b, c, x0, y0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -68,6 +66,11 @@ and `false` otherwise.
 * `b`: a vector of length m;
 * `c`: a vector of length n.
 
+#### Optional arguments
+
+* `x0`: a vector of length m that represents an initial guess of the solution x;
+* `y0`: a vector of length n that represents an initial guess of the solution y.
+
 #### Output arguments
 
 * `x`: a dense vector of length m;
diff --git a/src/trilqr.jl b/src/trilqr.jl
index 7aa90d1c9..81e3f3bf3 100644
--- a/src/trilqr.jl
+++ b/src/trilqr.jl
@@ -21,6 +21,10 @@ export trilqr, trilqr!
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, y, stats) = trilqr(A, b, c, x0::AbstractVector, y0::AbstractVector; kwargs...)
+
+TriLQR can be warm-started from initial guesses `x0` and `y0` where `kwargs` are the same keyword arguments as above.
+
 Combine USYMLQ and USYMQR to solve adjoint systems.
 
     [0  A] [y] = [b]
@@ -32,12 +36,6 @@ USYMQR is used for solving dual system `Aᴴy = c` of size n × m.
 An option gives the possibility of transferring from the USYMLQ point to the
 USYMCG point, when it exists. The transfer is based on the residual norm.
 
-TriLQR can be warm-started from initial guesses `x0` and `y0` with
-
-    (x, y, stats) = trilqr(A, b, c, x0, y0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -47,11 +45,16 @@ and `false` otherwise.
 * `b`: a vector of length m;
 * `c`: a vector of length n.
 
+#### Optional arguments
+
+* `x0`: a vector of length n that represents an initial guess of the solution x;
+* `y0`: a vector of length m that represents an initial guess of the solution y.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
 * `y`: a dense vector of length m;
-* `stats`: statistics collected on the run in a [`AdjointStats`](@ref) structure.
+* `stats`: statistics collected on the run in an [`AdjointStats`](@ref) structure.
 
 #### Reference
 
diff --git a/src/trimr.jl b/src/trimr.jl
index 0905f351e..52bc1b8e9 100644
--- a/src/trimr.jl
+++ b/src/trimr.jl
@@ -22,6 +22,10 @@ export trimr, trimr!
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, y, stats) = trimr(A, b, c, x0::AbstractVector, y0::AbstractVector; kwargs...)
+
+TriMR can be warm-started from initial guesses `x0` and `y0` where `kwargs` are the same keyword arguments as above.
+
 Given a matrix `A` of dimension m × n, TriMR solves the symmetric linear system
 
     [ τE    A ] [ x ] = [ b ]
@@ -53,12 +57,6 @@ TriMR stops when `itmax` iterations are reached or when `‖rₖ‖ ≤ atol + 
 Additional details can be displayed if verbose mode is enabled (verbose > 0).
 Information will be displayed every `verbose` iterations.
 
-TriMR can be warm-started from initial guesses `x0` and `y0` with
-
-    (x, y, stats) = trimr(A, b, c, x0, y0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -68,6 +66,11 @@ and `false` otherwise.
 * `b`: a vector of length m;
 * `c`: a vector of length n.
 
+#### Optional arguments
+
+* `x0`: a vector of length m that represents an initial guess of the solution x;
+* `y0`: a vector of length n that represents an initial guess of the solution y.
+
 #### Output arguments
 
 * `x`: a dense vector of length m;
diff --git a/src/usymlq.jl b/src/usymlq.jl
index d61fc10b0..1d6d3e1d8 100644
--- a/src/usymlq.jl
+++ b/src/usymlq.jl
@@ -21,13 +21,17 @@ export usymlq, usymlq!
 
 """
     (x, stats) = usymlq(A, b::AbstractVector{FC}, c::AbstractVector{FC};
-                        atol::T=√eps(T), rtol::T=√eps(T), transfer_to_usymcg::Bool=true,
-                        itmax::Int=0, verbose::Int=0, history::Bool=false,
-                        callback=solver->false)
+                        atol::T=√eps(T), rtol::T=√eps(T),
+                        transfer_to_usymcg::Bool=true, itmax::Int=0,
+                        verbose::Int=0, history::Bool=false, callback=solver->false)
 
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, stats) = usymlq(A, b, c, x0::AbstractVector; kwargs...)
+
+USYMLQ can be warm-started from an initial guess `x0` where `kwargs` are the same keyword arguments as above.
+
 Solve the linear system Ax = b of size m × n using the USYMLQ method.
 
 USYMLQ is based on the orthogonal tridiagonalization process and requires two initial nonzero vectors `b` and `c`.
@@ -41,12 +45,6 @@ In all cases, problems must be consistent.
 An option gives the possibility of transferring to the USYMCG point,
 when it exists. The transfer is based on the residual norm.
 
-USYMLQ can be warm-started from an initial guess `x0` with
-
-    (x, stats) = usymlq(A, b, c, x0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -56,6 +54,10 @@ and `false` otherwise.
 * `b`: a vector of length m;
 * `c`: a vector of length n.
 
+#### Optional argument
+
+* `x0`: a vector of length n that represents an initial guess of the solution x.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;
diff --git a/src/usymqr.jl b/src/usymqr.jl
index 2cef2db0d..003a46dc5 100644
--- a/src/usymqr.jl
+++ b/src/usymqr.jl
@@ -21,13 +21,16 @@ export usymqr, usymqr!
 
 """
     (x, stats) = usymqr(A, b::AbstractVector{FC}, c::AbstractVector{FC};
-                        atol::T=√eps(T), rtol::T=√eps(T),
-                        itmax::Int=0, verbose::Int=0, history::Bool=false,
-                        callback=solver->false)
+                        atol::T=√eps(T), rtol::T=√eps(T), itmax::Int=0,
+                        verbose::Int=0, history::Bool=false, callback=solver->false)
 
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
 `FC` is `T` or `Complex{T}`.
 
+    (x, stats) = usymqr(A, b, c, x0::AbstractVector; kwargs...)
+
+USYMQR can be warm-started from an initial guess `x0` where `kwargs` are the same keyword arguments as above.
+
 Solve the linear system Ax = b of size m × n using USYMQR.
 
 USYMQR is based on the orthogonal tridiagonalization process and requires two initial nonzero vectors `b` and `c`.
@@ -38,12 +41,6 @@ It's considered as a generalization of MINRES.
 It can also be applied to under-determined and over-determined problems.
 USYMQR finds the minimum-norm solution if problems are inconsistent.
 
-USYMQR can be warm-started from an initial guess `x0` with
-
-    (x, stats) = usymqr(A, b, c, x0; kwargs...)
-
-where `kwargs` are the same keyword arguments as above.
-
 The callback is called as `callback(solver)` and should return `true` if the main loop should terminate,
 and `false` otherwise.
 
@@ -53,6 +50,10 @@ and `false` otherwise.
 * `b`: a vector of length m;
 * `c`: a vector of length n.
 
+#### Optional argument
+
+* `x0`: a vector of length n that represents an initial guess of the solution x.
+
 #### Output arguments
 
 * `x`: a dense vector of length n;