try using iterative solvers and linear maps (particularly for large-scale problems)

[chriscoey]

currently we only use direct linsys solvers (cholesky factorization, since we have a symmetric PSD matrix).

can try indirect/iterative (eg JuliaLinearAlgebra/IterativeSolvers.jl#1), especially for large problems

• can initialize with previous direction (warm-start)
• during early iterations, can use loose tolerances, then decrease tolerance to eps as we seem to converge

libraries like https://github.com/JuliaSparse/Pardiso.jl have both direct and iterative solvers

for a matrix-free approach, see https://github.com/Jutho/LinearMaps.jl

[chriscoey]

closely related: #65

[chriscoey]

closed by #74