Current implementation of iterative solver in DDSCAT
DDSCAT numerical solution is based on iterative solution involving matrix times vector multiplication proposed by Goodman et al. (1990) [3]. The basic idea of this algorithm, in case of one dimensional problem, is succintly described in
http://www.cs.utk.edu/~dongarra/etemplates/node384.html
The method is based on calculating matrix times vector by extending problem to larger one.
Other proposals
Barrowes et al. (2001) (PDF) [1] proposed a new FFT-based method to expedite matrix-vector multiplies by reducing calculations to a one dimensional case. They describe their method as having a similar purpose to the FFT-based method of Goodman et al [3], but more general in implementation due to its ability to multiply arbitrary blocked Toeplitz matrices times a vector. This approach is similar to that of Lee (1986) [4].
Flatau (2004) looked at algorithm of 1d case(PDF) [2]. His proposal was based on writting down matrix times vector multiplication in terms of the symmetric and skew-symmetric matrices. Extension to more dimensions are possible, see [5] and http://www.liv.ac.uk/maths/ETC/mpta/
Test Fortran codes related to this paper are available.