Integral equation methods have been used with great success in electromagnetic scattering calculations. However, their application is in many cases limited by the storage requirements of dense matrices and also by the rapidly increasing computational time. However, the use of iterative solvers and special methods for computing the matrix-vector products can greatly reduce both the CPU and memory requirements. The application of iterative solvers for both volume and surface integral equations of electromagnetic scattering will be described. The complex symmetric version of QMR is working efficiently for both problems. Iterative solvers converge quickly even without preconditioner for the volume integral equation. We show how the eigenvalues of the coefficient matrix and the spectrum of the corresponding integral operator are related for a spherical scatterer. For the surface integral equation we have applied the sparse approximate inverse preconditioner. For the volume integral equation, the matrix-vector products can be computed with a 3-dimensional FFT. For the surface integral equation, the fast multipole method can be employed. For both approaches, problems with hundreds of thousands or even millions of unknowns can be solved with reasonable computational requirements.