What do people suggest for the linear system of this type? I know about trilinos, petsc, and sundials, but don't know the other alternatives or have exposure to them.
Generally speaking PETSc, Trilinos, and KINSOL (from SUNDIALS) are the best-of-breed when it comes to scalability. From an ecosystem standpoint, PETSc seems most flexible, since it does not try to take over
main or attempt to be a framework; many packages (SLEPc, libmesh, deal.II, PyClaw, MOOSE, TAO (before it was folded into PETSc), FEniCS) build on top of PETSc successfully.
Trilinos has more of a framework/ecosystem feel, as Trilinos consists of many packages; packages are still built on top of Trilinos, like FEniCS and deal.II. From mailing list traffic, the PETSc interface to FEniCS tends to get more use than the Trilinos interface.
KINSOL is the least flexible. It's basically a Newton solver or a fixed-point iteration solver, with a small number of interfaces to iterative linear solvers. It's probably the easiest to learn completely, but doesn't have the runtime flexibility that something like PETSc does. Choosing between KINSOL and PETSc, I'd pick PETSc; the code maintenance and quality alone are better, and PETSc is way more powerful.
I realize that the large solvers are overkill due to the relative lack of sparsity in the example, but can it hurt?
Probably not. There's some overhead due to PETSc initializing data structures, for instance, but I've only run into two people who have actually complained about it with good reason (and in one case, PETSc developers did a good job of reducing the PETSc memory footprint). JedBrown would know better; I'm just a PETSc user.