Bugfix: numerical_setup! for PETScLinearSolver

[JordiManyer]

We previously threw away PETSc's matrix ns.B, and re-set the KSP object (commented code).
This is not only unnecessary, but can also cause some issues with MUMPS.

I am not completely sure why, but here are some notes on the issue:

• It is matrix-dependent, and only happens in parallel (nprocs > 1).
• It has to do with the re-use of the symmetric permutation created by MUMPS to
  find the pivots.
  I think it probably re-orders the matrix internally, and does not re-order it again when we swap it
  using KSPSetOperators. So when accessing the new (non-permuted) matrix using the old permutation,
  it throws a stack overflow error.

In fact, when updating the SNES setups in the nonlinear solvers, we do not re-set the matrix,
but use copy! instead.

[amartinhuertas]

We previously threw away PETSc's matrix ns.B, and re-set the KSP object (commented code).
This is not only unnecessary

I agree there is no need to do so (threw away ns.B), as the (implicit) assumption is that A has the same sparsity pattern as ns.B. I would not check that this is the case (it is costly), but let the code run even if not fulfilled. I guess that the copy! method will fail because of that?

[amartinhuertas]

are we sure that we do no longer need to call KSPSetOperators?

[JordiManyer]

This is what the manual says:

To solve successive linear systems that have different preconditioner matrices 
(i.e., the matrix elements and/or the matrix data structure change), the user 
must call KSPSetOperators and KSPSolve for each solve.

which would indicate we do. However, my tests point to this not being the case. This works, for instance:

solver = PETScLinearSolver(mykspsetup)
println("First call")
ns = numerical_setup(symbolic_setup(solver,A),A)
x = allocate_in_domain(A)
fill!(x,0.0)
solve!(x,ns,b)

println("Second call")
map(Ai -> rmul!(Ai,2.0), partition(A))
b *= 2.0
numerical_setup!(ns,A)
x2 = allocate_in_domain(A)
fill!(x2,0.0)
solve!(x2,ns,b)

e = norm(x-x2)
println("Error =", e)

[amartinhuertas]

In this copy! call, which are the possible types for A?

[JordiManyer]

In theory, all supported types