Put options to have log in Iterative solvers

[rveltz]

When using Iterative solvers, it is very useful to have the log of convergence to attest the usefullness / efficiency of a preconditionner. The option is not present in the current file linear_nonlinear.jl.

I suggest to do the following unless you have a better way of triggering the verbose option:

function (f::LinSolveIterativeSolvers)(x,A,b,update_matrix=false;verbose = false, Pl=nothing, Pr=nothing, tol=nothing, kwargs...)
  if f.iterable === nothing
    Pl = ComposePreconditioner(f.Pl, Pl, true)
    Pr = ComposePreconditioner(f.Pr, Pr, false)
    @warn "Romain modified here"
    f.iterable = f.generate_iterator(x,A,b,f.args...;
                                     initially_zero=true,restart=5,
                                     maxiter=5,tol=1e-16,Pl=Pl,Pr=Pr,
                                     f.kwargs...,kwargs...)
    tol′ = get(f.kwargs, :tol, nothing)
    tol′ !== nothing && (tol = tol′)
    f.iterable.reltol = tol
  end
  x .= false
  iter = f.iterable
  purge_history!(iter, x, b)

  verbose && @printf("=== gmres ===\n%4s\t%4s\t%7s\n","rest","iter","resnorm")
  restart = get(f.kwargs, :restart, 5)

  for (iteration, residual) = enumerate(iter)
    verbose && @printf("%3d\t%3d\t%1.2e\n", 1 + div(iteration - 1, restart), 1 + mod(iteration - 1, restart), residual)
  end
  return nothing
end

and also:

function (p::DefaultLinSolve)(x,A,b,update_matrix=false;tol=nothing,verbose=false, kwargs...)
  if p.iterable isa Vector && eltype(p.iterable) <: LinearAlgebra.BlasInt # `iterable` here is the pivoting vector
    F = LU{eltype(A)}(A, p.iterable, zero(LinearAlgebra.BlasInt))
    ldiv!(x, F, b)
    return nothing
  end
  if update_matrix
    if typeof(A) <: Matrix
      blasvendor = BLAS.vendor()
      if (blasvendor === :openblas || blasvendor === :openblas64) && size(A,1) <= 500 && ArrayInterface.can_setindex(x) # if the user doesn't use OpenBLAS, we assume that is a much better BLAS implementation like MKL
        p.A = RecursiveFactorization.lu!(A)
      else
        p.A = lu!(A)
      end
    elseif typeof(A) <: Tridiagonal
      p.A = lu!(A)
    elseif typeof(A) <: Union{SymTridiagonal}
      p.A = ldlt!(A)
    elseif typeof(A) <: Union{Symmetric,Hermitian}
      p.A = bunchkaufman!(A)
    elseif typeof(A) <: SparseMatrixCSC
      p.A = lu(A)
    elseif ArrayInterface.isstructured(A)
      p.A = factorize(A)
    elseif !(typeof(A) <: AbstractDiffEqOperator)
      # Most likely QR is the one that is overloaded
      # Works on things like CuArrays
      p.A = qr(A)
    end
  end

  if typeof(A) <: Union{Matrix,SymTridiagonal,Tridiagonal,Symmetric,Hermitian} # No 2-arg form for SparseArrays!
    x .= b
    ldiv!(p.A,x)
  # Missing a little bit of efficiency in a rare case
  #elseif typeof(A) <: DiffEqArrayOperator
  #  ldiv!(x,p.A,b)
  elseif ArrayInterface.isstructured(A) || A isa SparseMatrixCSC
    ldiv!(x,p.A,b)
  elseif typeof(A) <: AbstractDiffEqOperator
    # No good starting guess, so guess zero
    if p.iterable === nothing
      p.iterable = IterativeSolvers.gmres_iterable!(x,A,b;initially_zero=true,restart=15,maxiter=15,tol=1e-16,kwargs...)
      p.iterable.reltol = tol
    end
    x .= false
    iter = p.iterable
    purge_history!(iter, x, b)
    verbose && @printf("=== gmres ===\n%4s\t%4s\t%7s\n","rest","iter","resnorm")

    for residual in iter
    end
  else
    ldiv!(x,p.A,b)
  end
  return nothing
end

[rveltz]

Any chance to have this improved? The issue is that generate_iterator does not accept the verbose keyword. It has to be either shaved from f.kwargs , or passed as above, or put the field verbose in LinSolveIterativeSolvers

Also, the defaults should be handled more carefully, like (maybe)

f.iterable = f.generate_iterator(x,A,b,f.args...;
                                     initially_zero=true,restart=get(f.kwargs, :restart, 5),
                                     maxiter=get(f.kwargs, :maxiter, 5),tol=get(f.kwargs, :tol, 1e-16),Pl=Pl,Pr=Pr,
                                     f.kwargs...,kwargs...)

[ChrisRackauckas]

Yeah that looks good to me.

[vpuri3]

would be great to allow for a user-defined function within the Krylov iteration.

[ChrisRackauckas-Claude]

Triage: LinSolveGMRES, DefaultLinSolve, and related iterative solver infrastructure have been removed from DiffEqBase. This functionality now lives in LinearSolve.jl. If the issue persists there, please file it at https://github.com/SciML/LinearSolve.jl. Recommend closing.