# Implementation details

## Orthogonalization

To denote a basis of vectors, e.g. to represent a given Krylov subspace, there is an abstract type `Basis{T}`

`KrylovKit.Basis`

— Type`abstract type Basis{T} end`

An abstract type to collect specific types for representing a basis of vectors of type `T`

.

Implementations of `Basis{T}`

behave in many ways like `Vector{T}`

and should have a `length`

, can be indexed (`getindex`

and `setindex!`

), iterated over (`iterate`

), and support resizing (`resize!`

, `pop!`

, `push!`

, `empty!`

, `sizehint!`

).

The type `T`

denotes the type of the elements stored in an `Basis{T}`

and can be any custom type that has vector like behavior (as defined in the docs of KrylovKit).

See `OrthonormalBasis`

for a specific implementation.

Many Krylov based algorithms use an orthogonal basis to parameterize the Krylov subspace. In that case, the specific implementation `OrthonormalBasis{T}`

can be used:

`KrylovKit.OrthonormalBasis`

— Type`OrthonormalBasis{T} <: Basis{T}`

A list of vector like objects of type `T`

that are mutually orthogonal and normalized to one, representing an orthonormal basis for some subspace (typically a Krylov subspace). See also `Basis`

Orthonormality of the vectors contained in an instance `b`

of `OrthonormalBasis`

(i.e. `all(dot(b[i],b[j]) == I[i,j] for i=1:length(b), j=1:length(b))`

) is not checked when elements are added; it is up to the algorithm that constructs `b`

to guarantee orthonormality.

One can easily orthogonalize or orthonormalize a given vector `v`

with respect to a `b::OrthonormalBasis`

using the functions `w, = orthogonalize(v,b,...)`

or `w, = orthonormalize(v,b,...)`

. The resulting vector `w`

of the latter can then be added to `b`

using `push!(b, w)`

. Note that in place versions `orthogonalize!(v, b, ...)`

or `orthonormalize!(v, b, ...)`

are also available.

Finally, a linear combination of the vectors in `b::OrthonormalBasis`

can be obtained by multiplying `b`

with a `Vector{<:Number}`

using `*`

or `mul!`

(if the output vector is already allocated).

We can orthogonalize or orthonormalize a given vector to another vector (assumed normalized) or to a given `KrylovKit.OrthonormalBasis`

.

`KrylovKit.orthogonalize`

— Function```
orthogonalize(v, b::OrthonormalBasis, [x::AbstractVector,] alg::Orthogonalizer]) -> w, x
orthogonalize!(v, b::OrthonormalBasis, [x::AbstractVector,] alg::Orthogonalizer]) -> w, x
orthogonalize(v, q, algorithm::Orthogonalizer]) -> w, s
orthogonalize!(v, q, algorithm::Orthogonalizer]) -> w, s
```

Orthogonalize vector `v`

against all the vectors in the orthonormal basis `b`

using the orthogonalization algorithm `alg`

of type `Orthogonalizer`

, and return the resulting vector `w`

and the overlap coefficients `x`

of `v`

with the basis vectors in `b`

.

In case of `orthogonalize!`

, the vector `v`

is mutated in place. In both functions, storage for the overlap coefficients `x`

can be provided as optional argument `x::AbstractVector`

with `length(x) >= length(b)`

.

One can also orthogonalize `v`

against a given vector `q`

(assumed to be normalized), in which case the orthogonal vector `w`

and the inner product `s`

between `v`

and `q`

are returned.

Note that `w`

is not normalized, see also `orthonormalize`

.

For more information on possible orthogonalization algorithms, see `Orthogonalizer`

and its concrete subtypes `ClassicalGramSchmidt`

, `ModifiedGramSchmidt`

, `ClassicalGramSchmidt2`

, `ModifiedGramSchmidt2`

, `ClassicalGramSchmidtIR`

and `ModifiedGramSchmidtIR`

.

`KrylovKit.orthonormalize`

— Function```
orthonormalize(v, b::OrthonormalBasis, [x::AbstractVector,] alg::Orthogonalizer]) -> w, β, x
orthonormalize!(v, b::OrthonormalBasis, [x::AbstractVector,] alg::Orthogonalizer]) -> w, β, x
orthonormalize(v, q, algorithm::Orthogonalizer]) -> w, β, s
orthonormalize!(v, q, algorithm::Orthogonalizer]) -> w, β, s
```

Orthonormalize vector `v`

against all the vectors in the orthonormal basis `b`

using the orthogonalization algorithm `alg`

of type `Orthogonalizer`

, and return the resulting vector `w`

(of norm 1), its norm `β`

after orthogonalizing and the overlap coefficients `x`

of `v`

with the basis vectors in `b`

, such that `v = β * w + b * x`

.

In case of `orthogonalize!`

, the vector `v`

is mutated in place. In both functions, storage for the overlap coefficients `x`

can be provided as optional argument `x::AbstractVector`

with `length(x) >= length(b)`

.

One can also orthonormalize `v`

against a given vector `q`

(assumed to be normalized), in which case the orthonormal vector `w`

, its norm `β`

before normalizing and the inner product `s`

between `v`

and `q`

are returned.

See `orthogonalize`

if `w`

does not need to be normalized.

For more information on possible orthogonalization algorithms, see `Orthogonalizer`

and its concrete subtypes `ClassicalGramSchmidt`

, `ModifiedGramSchmidt`

, `ClassicalGramSchmidt2`

, `ModifiedGramSchmidt2`

, `ClassicalGramSchmidtIR`

and `ModifiedGramSchmidtIR`

.

The expansion coefficients of a general vector in terms of a given orthonormal basis can be obtained as

`KrylovKit.project!`

— Function```
project!(y::AbstractVector, b::OrthonormalBasis, x,
[α::Number = 1, β::Number = 0, r = Base.OneTo(length(b))])
```

For a given orthonormal basis `b`

, compute the expansion coefficients `y`

resulting from projecting the vector `x`

onto the subspace spanned by `b`

; more specifically this computes

` y[j] = β*y[j] + α * dot(b[r[j]], x)`

for all $j ∈ r$.

whereas the inverse calculation is obtained as

`KrylovKit.unproject!`

— Function```
unproject!(y, b::OrthonormalBasis, x::AbstractVector,
[α::Number = 1, β::Number = 0, r = Base.OneTo(length(b))])
```

For a given orthonormal basis `b`

, reconstruct the vector-like object `y`

that is defined by expansion coefficients with respect to the basis vectors in `b`

in `x`

; more specifically this computes

` y = β*y + α * sum(b[r[i]]*x[i] for i = 1:length(r))`

An orthonormal basis can be transformed using a rank-1 update using

`KrylovKit.rank1update!`

— Function```
rank1update!(b::OrthonormalBasis, y, x::AbstractVector,
[α::Number = 1, β::Number = 1, r = Base.OneTo(length(b))])
```

Perform a rank 1 update of a basis `b`

, i.e. update the basis vectors as

` b[r[i]] = β*b[r[i]] + α * y * conj(x[i])`

It is the user's responsibility to make sure that the result is still an orthonormal basis.

Note that this changes the subspace. A mere rotation of the basis, which does not change the subspace spanned by it, can be computed using

`KrylovKit.basistransform!`

— Function`basistransform!(b::OrthonormalBasis, U::AbstractMatrix)`

Transform the orthonormal basis `b`

by the matrix `U`

. For `b`

an orthonormal basis, the matrix `U`

should be real orthogonal or complex unitary; it is up to the user to ensure this condition is satisfied. The new basis vectors are given by

` b[j] ← b[i] * U[i,j]`

and are stored in `b`

, so the old basis vectors are thrown away. Note that, by definition, the subspace spanned by these basis vectors is exactly the same.

## Dense linear algebra

KrylovKit relies on Julia's `LinearAlgebra`

module from the standard library for most of its dense linear algebra dependencies.

## Factorization types

The central ingredient in a Krylov based algorithm is a Krylov factorization or decomposition of a linear map. Such partial factorizations are represented as a `KrylovFactorization`

, of which `LanczosFactorization`

and `ArnoldiFactorization`

are two concrete implementations:

`KrylovKit.KrylovFactorization`

— Type`abstract type KrylovFactorization{T,S<:Number}`

Abstract type to store a Krylov factorization of a linear map `A`

of the form

`A * V = V * B + r * b'`

For a given Krylov factorization `fact`

of length `k = length(fact)`

, the basis `V`

is obtained via `basis(fact)`

and is an instance of some subtype of `Basis{T}`

, with also `length(V) == k`

and where `T`

denotes the type of vector like objects used in the problem. The Rayleigh quotient `B`

is obtained as `rayleighquotient(fact)`

and `typeof(B)`

is some subtype of `AbstractMatrix{S}`

with `size(B) == (k,k)`

, typically a structured matrix. The residual `r`

is obtained as `residual(fact)`

and is of type `T`

. One can also query `normres(fact)`

to obtain `norm(r)`

, the norm of the residual. The vector `b`

has no dedicated name and often takes a default form (see below). It should be a subtype of `AbstractVector`

of length `k`

and can be obtained as `rayleighextension(fact)`

(by lack of a better dedicated name).

A Krylov factorization `fact`

can be destructured as `V, B, r, nr, b = fact`

with `nr = norm(r)`

.

See also `LanczosFactorization`

and `ArnoldiFactorization`

for concrete implementations, and `KrylovIterator`

(with in particular `LanczosIterator`

and `ArnoldiIterator`

) for iterators that construct progressively expanding Krylov factorizations of a given linear map and a starting vector.

`KrylovKit.LanczosFactorization`

— Type`mutable struct LanczosFactorization{T,S<:Real} <: KrylovFactorization{T,S}`

Structure to store a Lanczos factorization of a real symmetric or complex hermitian linear map `A`

of the form

`A * V = V * B + r * b'`

For a given Lanczos factorization `fact`

of length `k = length(fact)`

, the basis `V`

is obtained via `basis(fact)`

and is an instance of `OrthonormalBasis{T}`

, with also `length(V) == k`

and where `T`

denotes the type of vector like objects used in the problem. The Rayleigh quotient `B`

is obtained as `rayleighquotient(fact)`

and is of type `SymTridiagonal{S<:Real}`

with `size(B) == (k,k)`

. The residual `r`

is obtained as `residual(fact)`

and is of type `T`

. One can also query `normres(fact)`

to obtain `norm(r)`

, the norm of the residual. The vector `b`

has no dedicated name but can be obtained via `rayleighextension(fact)`

. It takes the default value $e_k$, i.e. the unit vector of all zeros and a one in the last entry, which is represented using `SimpleBasisVector`

.

A Lanczos factorization `fact`

can be destructured as `V, B, r, nr, b = fact`

with `nr = norm(r)`

.

`LanczosFactorization`

is mutable because it can `expand!`

or `shrink!`

. See also `LanczosIterator`

for an iterator that constructs a progressively expanding Lanczos factorizations of a given linear map and a starting vector. See `ArnoldiFactorization`

and `ArnoldiIterator`

for a Krylov factorization that works for general (non-symmetric) linear maps.

`KrylovKit.ArnoldiFactorization`

— Type`mutable struct ArnoldiFactorization{T,S} <: KrylovFactorization{T,S}`

Structure to store an Arnoldi factorization of a linear map `A`

of the form

`A * V = V * B + r * b'`

For a given Arnoldi factorization `fact`

of length `k = length(fact)`

, the basis `V`

is obtained via `basis(fact)`

and is an instance of `OrthonormalBasis{T}`

, with also `length(V) == k`

and where `T`

denotes the type of vector like objects used in the problem. The Rayleigh quotient `B`

is obtained as `rayleighquotient(fact)`

and is of type `B::PackedHessenberg{S<:Number}`

with `size(B) == (k,k)`

. The residual `r`

is obtained as `residual(fact)`

and is of type `T`

. One can also query `normres(fact)`

to obtain `norm(r)`

, the norm of the residual. The vector `b`

has no dedicated name but can be obtained via `rayleighextension(fact)`

. It takes the default value $e_k$, i.e. the unit vector of all zeros and a one in the last entry, which is represented using `SimpleBasisVector`

.

An Arnoldi factorization `fact`

can be destructured as `V, B, r, nr, b = fact`

with `nr = norm(r)`

.

`ArnoldiFactorization`

is mutable because it can `expand!`

or `shrink!`

. See also `ArnoldiIterator`

for an iterator that constructs a progressively expanding Arnoldi factorizations of a given linear map and a starting vector. See `LanczosFactorization`

and `LanczosIterator`

for a Krylov factorization that is optimized for real symmetric or complex hermitian linear maps.

`KrylovKit.GKLFactorization`

— Type`mutable struct GKLFactorization{TU,TV,S<:Real}`

Structure to store a Golub-Kahan-Lanczos (GKL) bidiagonal factorization of a linear map `A`

of the form

```
A * V = U * B + r * b'
A' * U = V * B'
```

For a given GKL factorization `fact`

of length `k = length(fact)`

, the two bases `U`

and `V`

are obtained via `basis(fact, :U)`

and `basis(fact, :V)`

. Here, `U`

and `V`

are instances of `OrthonormalBasis{T}`

, with also `length(U) == length(V) == k`

and where `T`

denotes the type of vector like objects used in the problem. The Rayleigh quotient `B`

is obtained as `rayleighquotient(fact)`

and is of type `Bidiagonal{S<:Number}`

with `size(B) == (k,k)`

. The residual `r`

is obtained as `residual(fact)`

and is of type `T`

. One can also query `normres(fact)`

to obtain `norm(r)`

, the norm of the residual. The vector `b`

has no dedicated name but can be obtained via `rayleighextension(fact)`

. It takes the default value $e_k$, i.e. the unit vector of all zeros and a one in the last entry, which is represented using `SimpleBasisVector`

.

A GKL factorization `fact`

can be destructured as `U, V, B, r, nr, b = fact`

with `nr = norm(r)`

.

`GKLFactorization`

is mutable because it can `expand!`

or `shrink!`

. See also `GKLIterator`

for an iterator that constructs a progressively expanding GKL factorizations of a given linear map and a starting vector `u₀`

.

A `KrylovFactorization`

or `GKLFactorization`

can be destructured into its defining components using iteration, but these can also be accessed using the following functions

`KrylovKit.basis`

— Function` basis(fact::KrylovFactorization)`

Return the list of basis vectors of a `KrylovFactorization`

, which span the Krylov subspace. The return type is a subtype of `Basis{T}`

, where `T`

represents the type of the vectors used by the problem.

`KrylovKit.rayleighquotient`

— Function`rayleighquotient(fact::KrylovFactorization)`

Return the Rayleigh quotient of a `KrylovFactorization`

, i.e. the reduced matrix within the basis of the Krylov subspace. The return type is a subtype of `AbstractMatrix{<:Number}`

, typically some structured matrix type.

`KrylovKit.residual`

— Function`residual(fact::KrylovFactorization)`

Return the residual of a `KrylovFactorization`

. The return type is some vector of the same type as used in the problem. See also `normres(F)`

for its norm, which typically has been computed already.

`KrylovKit.normres`

— Function`normres(fact::KrylovFactorization)`

Return the norm of the residual of a `KrylovFactorization`

. As this has typically already been computed, it is cheaper than (but otherwise equivalent to) `norm(residual(F))`

.

`KrylovKit.rayleighextension`

— Function`rayleighextension(fact::KrylovFactorization)`

Return the vector `b`

appearing in the definition of a `KrylovFactorization`

; often it is simply the last coordinate unit vector, which can be represented using `SimpleBasisVector`

.

As the `rayleighextension`

is typically a simple basis vector, we have created a dedicated type to represent this without having to allocate an actual vector, i.e.

`KrylovKit.SimpleBasisVector`

— Type`SimpleBasisVector(m, k)`

Construct a simple struct `SimpleBasisVector <: AbstractVector{Bool}`

representing a coordinate basis vector of length `m`

in the direction of `k`

, i.e. for `e_k = SimpleBasisVector(m, k)`

we have `length(e_k) = m`

and `e_k[i] = (i == k)`

.

Furthermore, to store the Rayleigh quotient of the Arnoldi factorization in a manner that can easily be expanded, we have constructed a custom matrix type to store the Hessenberg matrix in a packed format (without zeros):

`KrylovKit.PackedHessenberg`

— Type```
struct PackedHessenberg{T,V<:AbstractVector{T}} <: AbstractMatrix{T}
data::V
n::Int
end
```

A custom struct to store a Hessenberg matrix in a packed format (without zeros). Hereto, the non-zero entries are stored sequentially in vector `data`

of length `n(n+1)/2`

.

## Factorization iterators

Given a linear map $A$ and a starting vector $x₀$, a Krylov factorization is obtained by sequentially building a Krylov subspace ${x₀, A x₀, A² x₀, ...}$. Rather then using this set of vectors as a basis, an orthonormal basis is generated by a process known as Lanczos or Arnoldi iteration (for symmetric/hermitian and for general matrices, respectively). These processes are represented as iterators in Julia:

`KrylovKit.KrylovIterator`

— Type`abstract type KrylovIterator{F,T}`

Abstract type for iterators that take a linear map of type `F`

and an initial vector of type `T`

and generate an expanding `KrylovFactorization`

thereof.

When iterating over an instance of `KrylovIterator`

, the values being generated are subtypes of `KrylovFactorization`

, which can be immediately destructured into a `basis`

, `rayleighquotient`

, `residual`

, `normres`

and `rayleighextension`

.

See `LanczosIterator`

and `ArnoldiIterator`

for concrete implementations and more information.

`KrylovKit.LanczosIterator`

— Type```
struct LanczosIterator{F,T,O<:Orthogonalizer} <: KrylovIterator{F,T}
LanczosIterator(f, v₀, [orth::Orthogonalizer = KrylovDefaults.orth, keepvecs::Bool = true])
```

Iterator that takes a linear map `f::F`

(supposed to be real symmetric or complex hermitian) and an initial vector `v₀::T`

and generates an expanding `LanczosFactorization`

thereof. In particular, `LanczosIterator`

uses the Lanczos iteration scheme to build a successively expanding Lanczos factorization. While `f`

cannot be tested to be symmetric or hermitian directly when the linear map is encoded as a general callable object or function, it is tested whether the imaginary part of `dot(v, f(v))`

is sufficiently small to be neglected.

The argument `f`

can be a matrix, or a function accepting a single argument `v`

, so that `f(v)`

implements the action of the linear map on the vector `v`

.

The optional argument `orth`

specifies which `Orthogonalizer`

to be used. The default value in `KrylovDefaults`

is to use `ModifiedGramSchmidtIR`

, which possibly uses reorthogonalization steps. One can use to discard the old vectors that span the Krylov subspace by setting the final argument `keepvecs`

to `false`

. This, however, is only possible if an `orth`

algorithm is used that does not rely on reorthogonalization, such as `ClassicalGramSchmidt()`

or `ModifiedGramSchmidt()`

. In that case, the iterator strictly uses the Lanczos three-term recurrence relation.

When iterating over an instance of `LanczosIterator`

, the values being generated are instances of `LanczosFactorization`

, which can be immediately destructured into a `basis`

, `rayleighquotient`

, `residual`

, `normres`

and `rayleighextension`

, for example as

```
for (V, B, r, nr, b) in LanczosIterator(f, v₀)
# do something
nr < tol && break # a typical stopping criterion
end
```

Note, however, that if `keepvecs=false`

in `LanczosIterator`

, the basis `V`

cannot be extracted.

Since the iterator does not know the dimension of the underlying vector space of objects of type `T`

, it keeps expanding the Krylov subspace until the residual norm `nr`

falls below machine precision `eps(typeof(nr))`

.

The internal state of `LanczosIterator`

is the same as the return value, i.e. the corresponding `LanczosFactorization`

. However, as Julia's Base iteration interface (using `Base.iterate`

) requires that the state is not mutated, a `deepcopy`

is produced upon every next iteration step.

Instead, you can also mutate the `KrylovFactorization`

in place, using the following interface, e.g. for the same example above

```
iterator = LanczosIterator(f, v₀)
factorization = initialize(iterator)
while normres(factorization) > tol
expand!(iterator, factorization)
V, B, r, nr, b = factorization
# do something
end
```

Here, `initialize(::KrylovIterator)`

produces the first Krylov factorization of length 1, and `expand!(::KrylovIterator, ::KrylovFactorization)`

(@ref) expands the factorization in place. See also `initialize!(::KrylovIterator, ::KrylovFactorization)`

to initialize in an already existing factorization (most information will be discarded) and `shrink!(::KrylovFactorization, k)`

to shrink an existing factorization down to length `k`

.

`KrylovKit.ArnoldiIterator`

— Type```
struct ArnoldiIterator{F,T,O<:Orthogonalizer} <: KrylovIterator{F,T}
ArnoldiIterator(f, v₀, [orth::Orthogonalizer = KrylovDefaults.orth])
```

Iterator that takes a general linear map `f::F`

and an initial vector `v₀::T`

and generates an expanding `ArnoldiFactorization`

thereof. In particular, `ArnoldiIterator`

iterates over progressively expanding Arnoldi factorizations using the Arnoldi iteration.

The argument `f`

can be a matrix, or a function accepting a single argument `v`

, so that `f(v)`

implements the action of the linear map on the vector `v`

.

The optional argument `orth`

specifies which `Orthogonalizer`

to be used. The default value in `KrylovDefaults`

is to use `ModifiedGramSchmidtIR`

, which possibly uses reorthogonalization steps.

When iterating over an instance of `ArnoldiIterator`

, the values being generated are instances of `ArnoldiFactorization`

, which can be immediately destructured into a `basis`

, `rayleighquotient`

, `residual`

, `normres`

and `rayleighextension`

, for example as

```
for (V, B, r, nr, b) in ArnoldiIterator(f, v₀)
# do something
nr < tol && break # a typical stopping criterion
end
```

Since the iterator does not know the dimension of the underlying vector space of objects of type `T`

, it keeps expanding the Krylov subspace until the residual norm `nr`

falls below machine precision `eps(typeof(nr))`

.

The internal state of `ArnoldiIterator`

is the same as the return value, i.e. the corresponding `ArnoldiFactorization`

. However, as Julia's Base iteration interface (using `Base.iterate`

) requires that the state is not mutated, a `deepcopy`

is produced upon every next iteration step.

Instead, you can also mutate the `ArnoldiFactorization`

in place, using the following interface, e.g. for the same example above

```
iterator = ArnoldiIterator(f, v₀)
factorization = initialize(iterator)
while normres(factorization) > tol
expand!(iterator, factorization)
V, B, r, nr, b = factorization
# do something
end
```

Here, `initialize(::KrylovIterator)`

produces the first Krylov factorization of length 1, and `expand!(::KrylovIterator, ::KrylovFactorization)`

(@ref) expands the factorization in place. See also `initialize!(::KrylovIterator, ::KrylovFactorization)`

to initialize in an already existing factorization (most information will be discarded) and `shrink!(::KrylovFactorization, k)`

to shrink an existing factorization down to length `k`

.

Similarly, there is also an iterator for the Golub-Kahan-Lanczos bidiagonalization proces:

`KrylovKit.GKLIterator`

— Type```
struct GKLIterator{F,TU,O<:Orthogonalizer}
GKLIterator(f, u₀, [orth::Orthogonalizer = KrylovDefaults.orth, keepvecs::Bool = true])
```

Iterator that takes a general linear map `f::F`

and an initial vector `u₀::TU`

and generates an expanding `GKLFactorization`

thereof. In particular, `GKLIterator`

implements the Golub-Kahan-Lanczos bidiagonalization procedure. Note, however, that this implementation starts from a vector `u₀`

in the codomain of the linear map `f`

, which will end up (after normalisation) as the first column of `U`

.

The argument `f`

can be a matrix, a tuple of two functions where the first represents the normal action and the second the adjoint action, or a function accepting two arguments, where the first argument is the vector to which the linear map needs to be applied, and the second argument is either `Val(false)`

for the normal action and `Val(true)`

for the adjoint action. Note that the flag is thus a `Val`

type to allow for type stability in cases where the vectors in the domain and the codomain of the linear map have a different type.

The optional argument `orth`

specifies which `Orthogonalizer`

to be used. The default value in `KrylovDefaults`

is to use `ModifiedGramSchmidtIR`

, which possibly uses reorthogonalization steps.

When iterating over an instance of `GKLIterator`

, the values being generated are instances `fact`

of `GKLFactorization`

, which can be immediately destructured into a `basis(fact, :U)`

, `basis(fact, :V)`

, `rayleighquotient`

, `residual`

, `normres`

and `rayleighextension`

, for example as

```
for (U, V, B, r, nr, b) in GKLIterator(f, u₀)
# do something
nr < tol && break # a typical stopping criterion
end
```

Since the iterator does not know the dimension of the underlying vector space of objects of type `T`

, it keeps expanding the Krylov subspace until the residual norm `nr`

falls below machine precision `eps(typeof(nr))`

.

The internal state of `GKLIterator`

is the same as the return value, i.e. the corresponding `GKLFactorization`

. However, as Julia's Base iteration interface (using `Base.iterate`

) requires that the state is not mutated, a `deepcopy`

is produced upon every next iteration step.

Instead, you can also mutate the `GKLFactorization`

in place, using the following interface, e.g. for the same example above

```
iterator = GKLIterator(f, u₀)
factorization = initialize(iterator)
while normres(factorization) > tol
expand!(iterator, factorization)
U, V, B, r, nr, b = factorization
# do something
end
```

Here, `initialize(::GKLIterator)`

produces the first GKL factorization of length 1, and `expand!(::GKLIterator, ::GKLFactorization)`

(@ref) expands the factorization in place. See also `initialize!(::GKLIterator, ::GKLFactorization)`

to initialize in an already existing factorization (most information will be discarded) and `shrink!(::GKLIterator, k)`

to shrink an existing factorization down to length `k`

.

As an alternative to the standard iteration interface from Julia Base (using `iterate`

), these iterative processes and the factorizations they produce can also be manipulated using the following functions:

`KrylovKit.expand!`

— Function`expand!(iter::KrylovIterator, fact::KrylovFactorization)`

Expand the Krylov factorization `fact`

by one using the linear map and parameters in `iter`

.

`KrylovKit.shrink!`

— Function`shrink!(fact::KrylovFactorization, k)`

Shrink an existing Krylov factorization `fact`

down to have length `k`

. Does nothing if `length(fact)<=k`

.

`KrylovKit.initialize`

— Function`initialize(iter::KrylovIterator)`

Initialize a length 1 Krylov factorization corresponding to `iter`

.

`KrylovKit.initialize!`

— Function`initialize!(iter::KrylovIterator, fact::KrylovFactorization)`

Initialize a length 1 Krylov factorization corresponding to `iter`

in the already existing factorization `fact`

, thereby destroying all the information it currently holds.