First things first, I would check your code: there is no reason 400 unique integrals should take an hour, even with the most naive implementation. You are computing at a rate of an integral every ten seconds. While I agree with all the comments pointing to increasing parallelization, I think efficient scalar code is more important. Why scale up inefficient code? Considering GPUs is a whole 'nother beast, and from the problems, my coworkers have expressed trying to program them I'd recommend avoiding them for the time being.
The best thing you can do is to exploit molecular symmetry. Heads up though: this won't fix everything, because many interesting molecules lack symmetry. But for methane, as you tried, this should help a ton. Of course, this isn't as simple as just looking up some algorithm, because depending on your program you need to rotate the nuclear orientations to be symmetric, or develop routines to detect point group symmetry, etc. That being said, there are --- as far as I know --- two primary ways of exploiting point group symmetry for integrals. The first is to create symmetry-adapted basis functions, so your basis set is symmetric, which means your wave function/integrals/Fock matrix/whatever will be symmetric as well. Ernie Davidson did this in the seventies. A slightly more modern reference by Helgaker and Taylor is given here. Daniel Crawford's electronic structure tutorial explores this as well (it's really well written --- a great introduction).
The other way you can do it, and I believe this is more popular, is to create symmetry unique integrals and then build a partial Fock matrix only to take care of symmetrization later. If you see something like "generating petite list" in the output of your favorite electronic structure program, this is probably what it's doing. This method allows for integral screening more efficiently, I believe. Dupuis and King did this in the late seventies, based on earlier work by Dacre and Elder. These are the standard references you'll find. There may be some more recent work that I'm unaware of. The problem of point group symmetry had been effectively solved by 1980.
Of course, many interesting molecules don't have symmetry to take advantage of, so symmetry won't fix all your problems. The Schwarz equality is well addressed in Molecular Electronic Structure Theory (actually, most electronic structure theory is well addressed in that lovely pink book). If you haven't checked that out, please do. I think it will answer many of your questions about how modern codes work.
I don't want to repeat the other answers (Geoff gave you the papers I'd have suggested as well), so I'll finish with two lesser-considered references. First, there is an awesome (but tiny!) book called Parallel Computing in Quantum Chemistry that lays out the basics for efficient code including integral evaluation. It's not too hard to read, and I found the work on local MP2 to be pretty interesting. I'm often surprised more people haven't heard of it. The other place I would check out is Ed Valeev's LibInt package. He's a pretty talented programmer, and it's used in many modern packages like Psi and MPQC and ORCA, not to mention free for anyone to download and use. Some of my coworkers have played around with it for their own "homemade" electronic structure code, and speak highly of it. The user's manual is actually quite helpful in understanding integral evaluation. He has also written several exercises and tutorials, which you may find useful.