bug fix and docs

docs/make.jl

@@ -10,6 +10,7 @@ makedocs(; modules=[KrylovKit],
                              "man/eig.md",
                              "man/svd.md",
                              "man/matfun.md",
+                             "man/reallinear.md",
                              "man/algorithms.md",
                              "man/implementation.md"]],
          format=Documenter.HTML(; prettyurls=get(ENV, "CI", nothing) == "true"))

docs/src/index.md

@@ -20,6 +20,11 @@ The high level interface of KrylovKit is provided by the following functions:
 *   [`expintegrator`](@ref): exponential integrator for a linear non-homogeneous ODE
     (generalization of `exponentiate`)
 
+Furthermore, for specialised use cases, there are functions that can deal with so-called
+"real linear maps", which arise e.g. in the context of differentiable programming:
+*   [`reallinsolve`](@ref) and [`realeigsolve`](@ref)
+
+
 ## Package features and alternatives
 This section could also be titled "Why did I create KrylovKit.jl"?
 

docs/src/man/eig.md

@@ -41,11 +41,7 @@ corresponding to the largest magnitude eigenvalue of `A`.
 
 More generally, if you want to compute several eigenvalues of a real linear map, and you know
 that all of them are real, so that also the associated eigenvectors will be real, then you
-can use the `realeigsolve` method, which is also restricted to the 'expert' method call and
-which will error if any of the requested eigenvalues turn out to be complex
-```@docs
-realeigsolve
-```
+can use the [`realeigsolve`](@ref) method.
 
 ## Automatic differentation
 

docs/src/man/linear.md

@@ -1,5 +1,7 @@
 # Linear problems
 
+## Linear systems
+
 Linear systems are of the form `A*x=b` where `A` should be a linear map that has the same
 type of output as input, i.e. the solution `x` should be of the same type as the right hand
 side `b`. They can be solved using the function `linsolve`:

docs/src/man/matfun.md

@@ -1,5 +1,5 @@
-# Functions of matrices and linear maps
-Applying a function of a matrix or linear map to a given vector can in some cases also be
+# Functions of matrices and linear operator
+Applying a function of a matrix or linear operator to a given vector can in some cases also be
 computed using Krylov methods. One example is the inverse function, which exactly
 corresponds to what `linsolve` computes: ``A^{-1} * b``. There are other functions ``f``
 for which ``f(A) * b`` can be computed using Krylov techniques, i.e. where ``f(A) * b`` can

docs/src/man/reallinear.md

@@ -0,0 +1,46 @@
+# Real linear maps
+
+A map $$f: V \to V$$ from some vector space $$V$$ to itself is said to be real linear map if
+it satisfies $$f(\alpha x + \beta y) = \alpha f(x) + \beta f(y)$$ for all $$x, y \in V$$ and
+all $$\alpha, \beta \in \mathbb{R}$$. When $$V$$ is itself a real vector space, this is just
+the natural concept of a linear map. However, this definition can be used even if $$x$$ and
+$$y$$ are naturally represented using complex numbers and arithmetic and also admit complex linear
+combinations, i.e. if $$V$$ is a complex vector space.
+
+Such real linear maps arise whenever `f(x)` involves calling `conj(x)`, and are for example
+obtained in the context of Jacobians (pullbacks) of complex valued functions that are not
+holomorphic.
+
+To deal with real linear maps, one should reinterpret $$V$$ as a real vector space, by
+restricting the possible linear combinations to those with real scalar coefficients, and by
+using the real part of the inner product. When the vectors are explictly represented as
+some `AbstractVector{Complex{T}}`, this could be obtained by explicitly splitting
+them in their real and imaginary parts and stacking those into `AbstractVector{T}` objects
+with twice the original length.
+
+However, KrylovKit.jl admits a different approach, where the original representation of
+vectors is kept, and the inner product is simply replaced by its real part. KrylovKit.jl
+offers specific methods for solving linear systems and eigenvalue systems in this way. For
+linear problems, this is implemented using `reallinsolve`:
+
+```@docs
+reallinsolve
+```
+
+In the case of eigenvalue systems, a similar method `realeigsolve` is available. In this
+context, only real eigenvalues are meaningful, as the corresponding eigenvectors should be
+built from real linear combinations of the vectors that span the (real) Krylov subspace.
+This approach can also be applied to linear maps on vectors that were naturally real to
+begin with, if it is guaranteed that the targetted eigenvalues are real. In that case, also
+the associated eigenvectors will be returned using only real arithmic. This is contrast
+with `eigsolve`, which will always turn to complex arithmetic if the linear map is real but
+not symmetric. An error will be thrown if complex eigenvalues are encountered within the
+targetted set.
+
+```@docs
+realeigsolve
+```
+
+Note that both `reallinsolve` and `realeigsolve` currently only exist with the "expert" mode
+interface, where the user has to manually specify the underlying Krylov algorithm and its
+parameters, i.e. `GMRES` or `BiCGStab` for `reallinsolve` and `Arnoldi` for `realeigsolve`.

docs/src/man/svd.md

@@ -1,4 +1,6 @@
 # Singular value problems
+
+## Singular values and singular vectors
 It is possible to iteratively compute a few singular values and corresponding left and
 right singular vectors using the function `svdsolve`:
 

test/testsetup.jl

@@ -133,6 +133,7 @@ end
 
 if VERSION < v"1.9"
     stack(f, itr) = mapreduce(f, hcat, itr)
+    stack(itr) = reduce(hcat, itr)
 end
 
 end