It would be nice to simply define a Matrix of FixedEffects, i.e. maybe using the https://github.com/Jutho/LinearMaps.jl framework.
Users can then directly use any external solver to solve for Ax = b. The sparse stuff would be for free with convert(SparseMatrixCSC, A::LinearMap) method.
One potential issue is that I may lose the preconditioning in LSMR. Maybe define mul! with preconditionner and without preconditionner.