With ``SparseArray`KrylovLinearSolve[A,b, "Preconditioner" -> f]`` you can use an arbitrary function `f` as preconditioner. Here an example: precdata = SparseArray`SparseMatrixILU[A, "Method" -> "ILU0"] f = x \[Function] SparseArray`SparseMatrixApplyILU[precdata, x] Other accepted values for the option `"Method"` are "ILUT" and "ILUTP". More or less, SparseArray`KrylovLinearSolve[A,b, "Method" -> "BiCGSTAB", "Preconditioner" -> f] should then be equivalent to LinearSolve[A,b, Method -> {"Krylov", "Method" -> "BiCGSTAB", "Preconditioner" -> "ILU0"}] Other supported Krylov methods are `"ConjugateGradient"` (only for symmetric positive-definite matrices) and `"GMRES"` See the documentation of `LinearSolve`, section Options, subsection Methods subsubsection "Krylov". Using ``SparseArray`KrylovLinearSolve`` allows you to build your own preconditioners. For example, f = x \[Function] Evaluate[1/Diagonal[A] x] is the well-known [Jacobi preconditioner](https://en.wikipedia.org/wiki/Preconditioner#Jacobi_(or_diagonal)_preconditioner) ( -- or rather its inverse. In constrast to the mathematical nomenclature, a preconditioner itself gets applied, not its inverse, in Mathematica. That's very plausible from an algorithmic point pf view. See [this](https://en.wikipedia.org/wiki/Preconditioner#Jacobi_(or_diagonal)_preconditioner) and the standard literature about numerical linear equations for more information about what a preconditioner is. Basically, a preconditioner is meant to speed up the convergence of iterative linear solvers. Expanding on that would be out of scope for this site.