Initial infrastructure for Krylov-based linear solvers

[mtanneau]

cc @amontoison

This PR introduces some infrastructure for Krylov-based linear solvers.
There are no Krylov-specific algorithmic modifications to the IPM at this point.

TODO:

• Identify relevant data structures and interfaces for the task
• Decide on how to best implement preconditioners [left for future PR]
• Add unit tests
• Update docs

Currently, a Krylov solver can be selected as follows:

model.params.KKTOptions = Tulip.KKT.SolverOptions(
    Tulip.KKT.KrylovPDSolver,       # Solve the normal equations system with a Krylov method
    method=Tulip.KKT.Krylov.minres  # Select the Krylov method
)

the syntax also allows for passing a callback:

model.params.KKTOptions = Tulip.KKT.SolverOptions(
    Tulip.KKT.KrylovPDSolver,       # Solve the normal equations system with a Krylov method
    method=(A, b) -> Tulip.KKT.Krylov.cg(A, b)  # "custom" Krylov method
)

[codecov-commenter]

Codecov Report

Merging #56 into master will decrease coverage by 0.04%.
The diff coverage is 87.75%.

@@            Coverage Diff             @@
##           master      #56      +/-   ##
==========================================
- Coverage   89.72%   89.67%   -0.05%     
==========================================
  Files          30       31       +1     
  Lines        2073     2122      +49     
==========================================
+ Hits         1860     1903      +43     
- Misses        213      219       +6     

Impacted Files	Coverage Δ
src/KKTSolver/KKTSolver.jl	77.77% <ø> (ø)
src/KKTSolver/krylov.jl	87.75% <87.75%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 68a3b96...2003908. Read the comment docs.

[amontoison]

You should move v1 and v2 in KrylovPDSolver structure and allocate them only once.

[amontoison]

You could use v1 and v2 here to avoid 2 new allocations at each solve! pass:

mul!(v1, D, ξd)
mul!(v2, kkt.A, v1)
v2 .+= ξp

[amontoison]

You could reuse v1 and v2 :

mul!(v1, kkt.A', dy)
v1 .-= ξd
mul!(v2, D, v1)

[amontoison]

You don't need to form d and D explicitly. If you want to keep it because it's more user-friendly, I recommend to preallocate a vector d inside the KrylovPDSolver structure and reuse the memory.

[amontoison]

v3 = zeros(Tv, n)
opD = LinearOperator(Tv, n, n, true, true, v -> v3 .= v ./ (kkt.θ .+ kkt.regP))

[mtanneau]

I don't want to optimize the code too much for now.
I'll update the docs.

[mtanneau]

I'm happy with the current state.

@amontoison I will introduce Jacobi preconditioners in a forthcoming PR.

[amontoison]

Let's merge thIs PR. 👍

[mtanneau]

One last issue: unified syntax for the resolution of SQD systems.

For functions like minres or minres_qlp, things look like

K = [-(Θ⁻¹ + Rp) Aᵀ; A Rd]  # obviously this would be a linear operator
Δ, stats = minres(K, ξ)  # solve 2x2 system

# Recover search directions
dx .= Δ[1:n]
dy .= Δ[n+1:m]

while for trimr we have

_dx, _dy, stats = trimr(A, ξp, ξd)  # solve 2x2 system

# Recover search directions
dx .= _dx
dy .= _dy

We could have an internal flag for KrylovSQDSolvers that indicates whether the method exploits the 2x2 block structure or not, but I would prefer to use dispatch instead.

So, I suggest to have a KrylovSIDSolver type for methods that solve the Symmetric InDefinite augmented system without exploiting the SQD structure explicitly.

[amontoison]

I prefer the second solution with a KrylovSIDSolver type.