Selecting a linear solver for integration of stiff problems

• If the problem size is fairly small (say N < 100), then using the dense solver is probably best; it is the simplest to use and reasonably inexpensive for small N. For larger N, it is important to take advantage of sparsity (zero-nonzero) structure within the problem.

• If there is local (nearest-neighbor) coupling, or if the coupling is local after a suitable reordering of y, then use the banded linear solver. Local coupling means that the i-th component of the RHS or residual function depends only on components $y_j$ for which $|i-j|$ is small relative to N. (Note that the dense and band solvers are only applicable for the serial version of the solver.)

• For even larger problems, consider one of the Krylov iterative methods. These are hardest to use, because for best results they usually require preconditioning. However they offer the best opportunity to exploit the sparsity structure in the problem. The preconditioner is a matrix which, at least crudely, approximates the actual matrix in the linear system to be solved, and is typically built from an approximation of the relevant Jacobian matrix. Typically, that approximation uses only part of the true Jacobian, but as a result is much less expensive to solve. If the Jacobian can be approximated by a matrix that is banded (serial case) or block-diagonal with banded blocks (parallel case), SUNDIALS includes preconditioner modules for such cases. In each of the user guides, the section 'Linear solver specification functions' and the section on preconditioner modules contain more detailed comments on preconditioning.

• On the construction of preconditioners for problems arising from the spatial discretization of time-dependent partial differential equation systems, there is considerable discussion in the paper P. N. Brown and A. C. Hindmarsh, Reduced Storage Matrix Methods in Stiff ODE Systems, J. Appl. Math. & Comp., 31 (1989), pp.40-91.