linear_solve_design.jl

#=
Replace calcJ, calcW, linsolve, do_newJW to [u, p, t, newJ, newW, tol, error_estimate]
foo!(z, integrator, alg, nlsolver) -> z .= W\z
foo(z, integrator, alg, nlsolver) -> W\z

foo!(z, integrators, alg, nlsolver)

- User control inverse of W(_t)
- user defined factorization (object with ldiv! defined) (Info: u, p, t, newW,
  tol, error_estimate)
- User control Jacobian (structured/sparse) (newJ)
- partial and real Jacobian
- multiple Jacobians (J = sum(Js))
=#

#Replace calcJ, calcW, linsolve, do_newJW to [u, p, t, newJ, newW, tol, error_estimate]
#foo!(z, integrator, alg, nlsolver) -> z .= W\z
#foo(z, integrator, alg, nlsolver) -> W\z

mutable struct LinSolveFactorize...
  factorizeW
  F... # factorized
  Js...
  Ws... # mutated in factor
end
#Rodas5(linsolve=LinSolve(qr!))

function (linearsolver::LinSolveFactorize)(z, integrators, alg, nlsolver) # so does Ros
  newJ, newW = do_newJW(...)

  if newJ
    calc_J!(linearsolver, ...)
  end
  if newW
    # prec
    calc_W!(linearsolver, ...)
    linearsolver.factorizeW(F, W) # lu -> qr from alg
  end
  ldiv!(F, z)
end

# GMRES
#Rodas5(linsolve=LinSolveKrylov(solver=gmres!, left_pre = ...))
mutable struct LinSolveKrylov...
  factorizeW
  Js...
  krylovcache...
end
function (linearsolver::LinSolveKrylov)(z, integrators, alg, nlsolver) # so does Ros
  newJ, newW = do_newJW(...)
  if newJ
    calc_J!(linearsolver, ...)
  end
  if newW
    calc_W!(linearsolver, ...)
    linearsolver.factorizeW(F, W) # lu -> qr from alg
  end
  #integrator.opts.reltol default preconditioner
  gmres!(action(f), z, cache=linearsolver.krylovcache, precond = F)
end


# partical factorization
mutable struct LinSolveApproximateFactorization...
  factorizeWs
  F... # factorized
  Js...
  Ws... # mutated in factor
end

function (linearsolver::LinSolveApproximateFactorization)(z, integrators, alg, nlsolver) # so does Ros
  newJ, newW = do_newJW(...)

  if newJ
    calc_J!(linearsolver.Js[1], integrator.f[1], ...)
    calc_J!(linearsolver.Js[2], integrator.f[2], ...)
  end
  if newW
    # prec
    calc_W!(linearsolver.Ws[1], linearsolver.Js[1], ...)
    calc_W!(linearsolver.Ws[2], linearsolver.Js[2], ...)
    linearsolver.factorizeW(F, Ws) # lu -> qr from alg
  end
  ldiv!(F, z)
end