This would get better answers at [scicomp.se], anyway:
- if your matrix happens to be positive semidefinite, then after a suitable shift the smallest eigenvalues become the largest ones in modulus and can be computed with the power method. Similar tricks apply in cases when the sought eigenvalues are "at the border of the spectrum".
- For all other cases, shifted Arnoldi (implemented in ARPACK which is suggested in the comments) is the go-to algorithm. It requires solving several linear systems with $A$, possibly shifted (i.e., $(A-tI)x=b$).
- If your matrix is singular, you can get by inverting $A-tI$ for a small $t$ rather than $A$, with similar results.
- If the eigenvalues you are looking for are somewhere in the middle of the spectrum, in general you are going to need to solve shifted linear systems with your matrix inside the algorithm. There aren't many ways around this. Luckily, for "medium-large" matrices, this is feasible with sparse direct methods, and there are off-the-shelf libraries that do it for you taking care of the low-level details.
- A recent family of algorithms that aims, among other things, to compute eigenvalues from inside the spectrum without the need to solve linear systems is the FEAST one. I have never used them personally, but I think they are worth exploring after Arnoldi if you determine that it does not meet your needs. Another algorithm aimed at larger-than-Arnoldi scale is called Jacobi-Davidson (JD).