Preconditioned Krylov methods are widely used and offer many advantages for solving sparse linear systems that do not have highly convergent, geometric multigrid solvers or specialized fast solvers. Unfortunately, however, Krylov methods encounter well-known difficulties in scaling beyond 10,000 processor cores because each iteration requires at least one vector inner product, which in turn requires a global synchronization that scales poorly because of internode latency. To help overcome these difficulties, we have developed hierarchical and nested Krylov methods in the PETSc library that reduce the number of global inner products required across the entire system (where they are expensive), while freely allowing inner products across smaller subsets of the entire system (where they are inexpensive) or using inner iterations that do not invoke vector inner products at all. We demonstrate that these methods significantly reduce overall simulation time on the Cray XK6 and Blue Gene/P for the PFLOTRAN subsurface flow application when using 10,000 through 224,000 cores.