Choice of Iterative Algorithm

As discussed elsewhere (e.g., Draine 1988), the problem of electromagnetic scattering of an incident wave by an array of N point dipoles can be cast in the form
 
where is a 3N-dimensional (complex) vector of the incident electric field at the N lattice sites, is a 3N-dimensional (complex) vector of the (unknown) dipole polarizations, and is a complex matrix.

Because 3N is a large number, direct methods for solving this system of equations for the unknown vector are impractical, but iterative methods are useful: we begin with a guess (typically, ) for the unknown polarization vector, and then iteratively improve the estimate for until equation (10) is solved to some error criterion. The error tolerance may be specified as
 
where is the Hermitian conjugate of [], and h is the error tolerance. We typically use in order to satisfy eq.(10) to high accuracy. The error tolerance h can be specified by the user (see Appendix A).

A major change in going from DDSCAT.4b to 5a is the implementation of several different algorithms for iterative solution of the system of complex linear equations. DDSCAT.5a has been modified to permit solution algorithms to be treated in a fairly ``modular'' fashion, facilitating the testing of different algorithms. Many algorithms were compared by Flatau (1997); two of them performed well and are made available to the user in DDSCAT.5a. The choice of algorithm is made by specifying one of the options:

• PBCGST - Preconditioned BiConjugate Gradient with STabilitization method from the Parallel Iterative Methods (PIM) package created by R. Dias da Cunha and T. Hopkins.
• PETRKP - the complex conjugate gradient algorithm of Petravic & Kuo-Petravic (1979), as coded in the Complex Conjugate Gradient package (CCGPACK) created by P.J. Flatau. This is the algorithm discussed by Draine (1988) and used in previous versions of DDSCAT.

Both methods work well. Our experience suggests that PBCGST is often faster than PETRKP, by perhaps a factor of two. We therefore recommend it, but include the PETRKP method as an alternative.

The Parallel Iterative Methods (PIM) by Rudnei Dias da Cunha (rdd@ukc.ac.uk) and Tim Hopkins (trh@ukc.ac.uk) is a collection of Fortran 77 routines designed to solve systems of linear equations on parallel and scalar computers using a variety of iterative methods (available at http://www.mat.ufrgs.br/pim-e.html). PIM offers a number of iterative methods, including

• Conjugate-Gradients, CG (Hestenes 1952),
• Bi-Conjugate-Gradients, BICG (Fletcher 1976),
• Conjugate-Gradients squared, CGS (Sonneveld 1989),
• the stabilised version of Bi-Conjugate-Gradients, BICGSTAB (van der Vorst 1991),
• the restarted version of BICGSTAB, RBICGSTAB (Sleijpen & Fokkema 1992)
• the restarted generalized minimal residual, RGMRES (Saad 1986),
• the restarted generalized conjugate residual, RGCR (Eisenstat 1983),
• the normal equation solvers, CGNR (Hestenes 1952 and CGNE (Craig 1955),
• the quasi-minimal residual, QMR (highly parallel version due to Bucker & Sauren 1996),
• transpose-free quasi-minimal residual, TFQMR (Freund 1992),
• the Chebyshev acceleration, CHEBYSHEV (Young 1981).

The source code for these methods is distributed with DDSCAT but only PBCGST and PETRKP can be called directly via ddscat.par. It is possible (and was done) to add other options by changing the code in getfml.f . A helpful introduction to conjugate gradient methods is provided by the report ``Conjugate Gradient Method Without Agonizing Pain" by Jonathan R. Shewchuk, available as a postscript file: ftp://REPORTS.ADM.CS.CMU.EDU/usr0/anon/1994/CMU-CS-94-125.ps.