Provide API for initial guess on the iterative solvers

[lvella]

Hi, I am trying to use ViennaCL to solve a temporal sequence of linear systems, where the solution of the is very similar to the solution of the first. In linear algebra packages, it seems a common practice to start an iterative solver from an initial guess. Please provide an API for using a viennacl::vector as the initial guess (and possibly as output as well, so no copy is needed and the data can be updated in-place).

[karlrupp]

The usual approach for this is to write x = x_guess + x_update and to rewrite Ax = b into Ax_update = b - Ax_guess, and then recover x. This only requires a simple modification of the right hand side, which we thought is reasonable to do for a user.

As you're not the first user requesting this feature, it would be good to add this option. However, I see constraints by the way the iterative solvers are called now: If the initial guess is provided as an optional argument, it clashes with the specification of a preconditioner. We have some significant refactoring of the iterative solvers scheduled for ViennaCL 2.0.0, yet I'd prefer not to break the existing API in the ViennaCL 1.x.y series. Is a better documentation of how to provide initial guesses (including an example) an acceptable compromise in the meanwhile?

[lvella]

Ok, I think I understood how to do it now. I'll use it for now, but considering it involves one extra matrix multiplication and one extra vector addition, besides having to keep the previous x_guess allocated, it seems a suboptimal solution.

About the clash with the specification of preconditioner, it doesn't seem to be an issue. Consider the two functions below:

template<typename MatrixType, typename VectorType>
VectorType solve(
           const MatrixType & matrix,
           VectorType const & rhs,
           gmres_tag const & tag,
           VectorType & x) {...}

and

template<typename MatrixType, typename VectorType, typename PreconditionerType>
VectorType solve(
           const MatrixType & matrix,
           VectorType const & rhs,
           gmres_tag const & tag,
           PreconditionerType const & precond) {...}

A call to the first would not clash with a call to the later, because C++ compiler will try to match the most specific template type first, so if the types of my parameters were (A, B, gmres_tag, B), it would match with the first. Otherwise, if they were (A, B, gmres_tag, C), it would match with the second.

[lvella]

I am trying to use the Ax_update = b - Ax_guess technique, and I indeed get a much smaller number of iterations per solve() call, but the answer seems worse, and after a few time steps my simulation diverges.

Solving Ax=b directly I have a much higher iteration count, it is slower, but the answer better (my simulation slowly converges).

Considering this, I believe both techniques are not equivalent concerning the residual tolerance. So, how must I adjust the tolerance for solving Ax_update = b - Ax_guess in order to have the same error I would get if I was solving Ax=b directly?

PS: have you reviewed what I've said about adding this interface won't clash with current API? Do you agree?

[karlrupp]

Have you verified that the obtained solution is sufficiently accurate? You might want to tighten the solver tolerance and/or the number of iterations first, as both methods should ultimately result in similar behavior.

The API extension you proposed is, unfortunately, a bit dangerous, because there are situations where the 'rhs' and the 'initial guess' do not have the same type (e.g. viennacl::vector<> vs. viennacl::vector_proxy<>). I'll see whether I can make this sufficiently robust for 1.7.0, but I can't promise anything at this point.

[lvella]

How about a simpler solution is to use another name than solve, for instance solve_from. Instead of normal solve, solve_from would return void and the solution could be update in-place.

For how I know the solution was worse, I was calculating the residual myself, from the solve output, with the following code:

x = s->b - viennacl::linalg::prod(s->A, x);
std::vector<real> res(x.size());
viennacl::fast_copy(x.begin(), x.end(), res.begin());
real mean_res = 0.0;
for(auto r : res) {
    mean_res += fabs(r);
}
mean_res /= res.size();
std::cout << "Calculated mean residual: " << mean_res << std::endl;

I am solving pressure correction for Navier-Stokes, and at each time step I use the previously calculated pressure field as the X initial guess. By doing so, I see that after some time steps, the residual starts to grow, exponentially, getting very quickly to inf. If I start fresh in each time step, and don't provide any initial guess, the mean residual remains steady during all simulation.