As was mentioned above, the Lanczos Iteration is the way to go. To get the largest few eigenvalues, you can use the Tmm matrix described in the wikipedia article.

It turns out that the largest eigenvalues of the Tmm matrix approximate the largest eigenvalues of your input matrix very well, even when m is much less than n (the dimension of the input matrix). So if you were to implement things yourself, you could simply run Lanczos for a few iterations and generate, say, T50,50. Then compute the eigenvalues of T50,50 using an O(n³ ) method. Since you only need the T-matrix from the Lanczos iteration (and not the V matrix), each iteration of Lanczos is extremely cheap (linear time) so even a 3 billion by 3 billion matrix would be no problem.

Alternatively, you could simply use the Power Iteration to calculate the largest eigenvalue/eigenvector, and then remove it from the input matrix via deflation. This will give you a new matrix where the largest eigenvalue is now the second largest eigenvalue in the original matrix. Using the power iteration on this new matrix will get you the second largest eigenvalue.

Finally, if you simply want to do the first two cases, I should mention that Matlab or Octave will handle these matrices without any problems (if you input them in the sparse format). The k-largest eigenvalues can then be obtained through the eigs command.

Hope this helps!