[Feature Request] Add Option for Initial Solution Vector

[RoyiAvital]

I want to use the package cg().
Yet it seems it is explicitly set the initial solution vector to zero:

Krylov.jl/src/cg.jl

Line 72 in b3ab247

x .= zero(T)

I think it would be great to have an option (In all solvers) to let the user set the initial condition.
I think you can let the user do it explicitly in the call for cg().
It will also play nicely with cg!(). Where you should allow set it in the constructor call: CgSolver(n, n, T, xInit).

Moreover, I think you should set the cg!() variant to return nothing.
Then the call for cg() will be like:

if isnothing(xinit)
  xinit = zeros(T, n);
end
solver = CgSolver(A, b, xinit = xinit);

return solver.x, solver.stats

Then in cg!() it will be shown nothing is returned.

The above is probably valid for all functions.

[dpo]

If you have a good initial vector, we recommend doing the following:

1. Compute r := b - A * xinit
2. solve A * ∆x = r with CG (which will initialize ∆x to zero)
3. compute x = xinit + ∆x.

The reason is that if ||xinit|| is "large", the successive steps generated by CG could be "small" in comparison and, when added to xinit, be absorbed because of numerical effects. That is less likely to happen if you perform the steps above as the final CG solution ∆x will be larger than individual steps.

That said, it's true that we could provide this kind of convenience interface. I seem to remember that this has come up before.

[RoyiAvital]

I see the logic. So indeed it will be great to have a wrapper to do just this.
As at the moment, I use cg() in a loop and I can convert from IterativeSolvers to Krylov because of that.
Also, some logic should be done when I call cg!() over and over. The logic would say more iterations will be added to the previous state. At the moment, the current implementation doesn't do that.

Also, I think it will be good to take the suggestion I made about the in place versions.
Id be even more extreme, let the user the option to supply buffers.

[dpo]

Also, some logic should be done when I call cg!() over and over. The logic would say more iterations will be added to the previous state. At the moment, the current implementation doesn't do that.

Not sure I understand. What's wrong with changing the max number of iterations?

Also, I think it will be good to take the suggestion I made about the in place versions.

It's already in place. The solution is stored in the solver object and can be reused.

[RoyiAvital]

OK, Let me try explain the current problems I have.
Assume we solve A x = b with cg().

I created a solver:

sCgSolver = CgSolver(100, 100, Float64);

Now, I run it for 100 iterations:

x, stats = cg!(sCgSolver, A, b, itmax = 100);

I check the solution and it is not good enough.
So I want to run 200 iterations more.
The intuition would say, do:

x, stats = cg!(sCgSolver, A, b, itmax = 200);

Yet it won't work since the solver won't use the last known state.
It will start from x .= zeros. Namely I "lost" my previous 100 iterations.
Instead of being at the iteration 100 + 200 = 300 I will be at 200.

I can do manually what you suggested above, but you understand it is not intuitive for the average user, right?

So my 3 suggestions are:

1. The function cg!() should use the current state of the solver as its starting point. So we have a constructor to get the vanilla solver but we can use it over and over with new iterations.
2. The constrictor itself and cg!() should allow defining init vector by the user.
3. The cg!() variant should not have explicit return (This is cosmetic, but it make sense).