Improvements in the discrete-dipole approximation method of computing scattering and absorption

August 15, 1997 / Vol. 22, No. 16 / OPTICS LETTERS 1205 Improvements in the discrete-dipole approximation method of computing scattering and absorption Piotr J. Flatau Scripps Institution of Oceanography, Univeristy of California, San Diego, La Jolla, California 92093-0221 Received April 14, 1997 Improvements in complex-conjugate gradient algorithms applied to the discrete-dipole approximation are reported. It is shown that computational time is reduced by use of the stabilized version of the biconjugate gratings algorithm, with diagonal left preconditioning.  1997 Optical Society of America The discrete-dipole approximation (DDA) is a f lexible technique for computing scattering and absorption by targets of arbitrary geometry. In Ref. 1 the DDA for scattering calculations is reviewed. Rather than direct methods for solving linear systems of equations that arise in the DDA problem, iterative methods for finding the solution have proved effective and efficient. Here we systematically study various nonstation- ary iterative (conjugate gradient) methods in a search for the most efficient one. We document implemen- tation of these methods in the public-domain code DDSCAT.5a. 1 Numerical aspects of the DDA continue to be of great interest. Yung 2 applied a conjugate gradient method in the DDA approach. Hoekstra 3 identified Yung's scheme as the conjugate gradient (CG) algo- rithm proposed by Hestenes and Stiefel. 4 Rahola5 discussed solution of dense systems of linear equa- tions in the DDA and the choice of the best iterative method for this application. Draine 6 implemented a CG method based on research of Petravic and Kuo- Petravic.7 This implementation is quite robust and has been used for many years.1 However, Lumme and Rahola8 applied the quasi-minimal residual (QMR) CG algorithm to the system of linear equations that arises in the DDA applications. They claim that the QMR method is approximately three times faster than the one employed in the DDSCAT code. 1 This study prompted me to perform the analysis re- ported here. PIM9 is a collection of FORTRAN 77 routines designed to solve systems of linear equations on parallel and sequential computers by a variety of iterative meth- ods. PIM contains implementations of various meth- ods: CG; conjugate gradients for normal equations with minimization of the residual norm (CGNR); conjugate gradients for normal equations with minimization of the error norm (CGNE); biconjugate gradients (Bi-CG); conjugate gradients squared (CGS); the stabilized ver- sion of Bi-CG (Bi-CGSTAB); the restarted, stabilized version of Bi-CG (RBi-CGSTAB); the restarted, gener- alized minimal residual (RGMRES); the restarted, gen- eralized conjugate residual (RGCR), the quasi-minimal residual with coupled two-term recurrences (QMR); the transpose-free QMR (TFQMR); and Chebyshev acceler- ation. The routines permit the use of preconditioners; the user may choose to use left, right, or symmetric preconditioning. The convergence rate of an iterative method depends on the coefficient matrix. Hence one can attempt to transform the linear system into one that is equiva- lent (in the sense that it has the same solution) but is easier to solve. A preconditioner is a matrix M that effects such a transformation. It is possible to intro- duce left and right preconditioners. 10 The simplest preconditioner consists of just the diagonal of the co- efficient matrix. This is known as the (point) Jacobi preconditioner. To compare these different algorithms, I used them to find solutions to the problem of scattering by a homogeneous sphere. The scattering problem is specified by the usual size parameter x ­ 2payl, where a is the radius. Tables 1 and 2 list the numbers of iterations and their CPU times for size parameters x ­ 0.1 and x ­ 1 and for several values of the refractive indices. The CG methods are defined as above. In the tables the label (L) indicates left Jacobi preconditioning. For ex- ample, CGNE(L) is the CG method for normal equa- tions with minimization of the error norm and left Jacobi preconditioning. Similarly, (R) indicates right Jacobi preconditioning. CPU time (on a sequential Silicon Graphics workstation) is normalized to the best method. An asterisk indicates that the method did not converge in the maximum allowed number of itera- tions or that the method failed to converge. Fractional error 10 25 was used as the stopping criterion. The DDSCAT.5a code1 with the newly implemented GPFA fast-Fourier-transform method was used. For the Bi- CGSTAB and CGNE algorithms, left and right Neu- mann polynomial preconditioners truncated after the first term were used. Thus Bi-CGSTAB(N) (L) indi- cates the stabilized version of the Bi-CG method with a left Neumann polynomial preconditioner. Table 1 lists results for size parameter x ­ 0.1 and real refractive indices n ­ 1.33, 2, 3, 5 as well as one case with a small complex part of refractive index n ­ s5, 0.0001d and size parameter x ­ 0.1. For Table 1 the CPU times are normalized to the CG(L) method, which was found to be the best method. For example, it is 4.0 times faster than the CGNE method for n ­ s1.33, 0d. For larger values of real refractive index the CGNE method is almost an order of magnitude slower than CG. This is so because more iterations are needed for the same convergence and because the cost of one CG iteration is less than cost of one CGNE 0146-9592/97/161205-03$10.00/0  1997 Optical Society of America

1206 OPTICS LETTERS / Vol. 22, No. 16 / August 15, 1997 Table 1. CPU Time (Normalized) and Number of Iterations for x5 0.1a Method n ­ s1.33, 0 d (2, 0) (3, 0) (5, 0) (5, 0.0001) CGNE 4.0(9) 4.9(24) 8.7(76) *(540) *(540) CGNE(L) 3.3(7) 4.0(19) 7.8(67) *(540) *(540) CGNE(R) 4.1(9) 4.9(24) 8.8(76) *(540) *(540) CGNE(N) (L) 3.1(4) 4.1(13) 19.3(113) *(540) *(540) CGNE(N) (R) 4.4(6) 7.7(25) *(140) *(540) *(540) CGNR 4.0(9) 4.7(23) 8.0(69) *(540) *(540) CGNR(L) 3.3(7) 4.0(19) 5.9(50) 4.6(329) 3.5(330) CGNR(R) 4.1(9) 4.7(23) 8.0(69) *(540) *(540) QMR 3.7(6) 3.3(11) 3.3(19) *(111) *(78) QMR(L) 2.6(4) 2.8(9) 2.7(15) *(75) *(92) QMR(R) 3.8(6) 3.4(11) 3.4(19) *(268) *(540) CG 1.4(6) 1.2(11) 1.2(20) 1.2(163) 1.1(213) CG(L) 1.0(4) 1.0(9) 1.0(16) 1.0(138) 1.0(182) CG(R) 1.4(6) 1.2(11) 1.2(20) 1.1(157) 1.2(213) BiCG 2.3(6) 2.1(11) *(140) *(540) *(540) BiCG(L) 1.6(4) 1.8(9) *(140) *(540) *(540) BiCG(R) 2.4(6) 2.2(11) *(140) *(540) *(540) Bi-CGSTAB 1.8(4) 1.5(7) 1.5(13) *(540) *(540) Bi-CGSTAB(L) 1.4(3) 1.3(6) 1.3(11) 4.0(281) 4.2(388) Bi-CGSTAB(R) 1.8(4) 1.5(7) 1.6(13) *(540) *(540) Bi-CGSTAB(N) (L) 1.9(2) 6.8(17) 14.9(65) *(540) *(540) Bi-CGSTAB(N) (R) 2.1(2) 13.1(33) 27.9(122) *(540) *(540) TFQMR 3.8(5) 3.4(9) 4.3(19) *(540) *(540) TFQMR(L) 3.1(4) 3.1(8) 4.1(18) *(540) *(540) TFQMR(R) 3.9(5) 3.5(9) 4.4(19) *(540) *(540) CGS 1.7(4) 1.5(7) 1.4(12) *(540) *(540) CGS(L) 1.4(3) 1.3(6) 1.2(10) *(540) *(540) CGS(R) 1.8(4) 1.5(70) 1.4(12) *(540) *(540) aThe labels (L), (R), and (N) and the asterisk are defined in text. iteration. The QMR algorithm is never competitive and actually fails to converge for large real refractive indices. For a small refractive index the Bi-CGSTAB algorithm is comparable with the CG method and re- quires fewer iterations. However, the cost per itera- tion is larger than for a CG iteration, which offsets the advantage of fewer iterations. The Petravic±Kuo- Petravic 7 algorithm that has been used for many years6 is similar to CGNR and CGNE. However, on occa- sion, slightly different convergence rates were observed because of stabilization of Petravic±Kuo-Petravic al- gorithm every 10th time step. 6 This is true for all other cases. The storage requirement for CG, CGNE, and CGNR is 6 3 N, for BiCG it is 8 3 N, for CGS, Bi-CGSTAB, and TFQMR it is 10 3 N; QMR re- quires 11 3 N. Thus, for a pure real refractive in- dex, the CG algorithm not only is the fastest method but also requires the least amount of temporary stor- age. It can be seen that left preconditioning by the inverse of the diagonal of the DDA matrix 1 reduces the number of iterations needed. The added time needed for division by diagonal elements is generally negli- gible compared with the time saved by the smaller number of iterations. It can be seen that for Bi- CGSTAB, Bi-CGSTAB(L), and Bi-CGSTAB(R) the left Jacobi preconditioning is the only method that con- verges for larger refractive index. Restarted methods (RBi-CBSTAB and RGCR) appear not to be compet- itive, but further study may be needed (I used the orthogonal base of 10 vectors for all restarted meth- ods). The CG method is also competitive in cases with small absorption (column 6 of Table 1). I also calcu- lated (not presented) results for the size parameter x ­ 0.1 and an increasingly complex part of refractive indices n ­ s1.33, 0d, (1.33, 0.01), (1.33, 0.1), (1.33, 1), (1.33, 2), (1.33, 3). The BiCGSTAB(L), method proved to be the most robust. However, CGS(L) is competi- tive and faster for n ­ s1.33, 3d. Both CGS(L) and BiCGSTAB(L) require the same amount of iteration for convergence, and their costs are similar. These meth- ods are 1.6±2.9 times faster than CGNR, the method used in the DDSCAT code. QMR and TFQMR, which Lumme and Rahola 8 claim is faster than CGNR and the DDSCAT implementations, do not converge on oc- casion, and when they work they are only slightly better in this case. As before, left Jacobi precondi- tioning is almost always beneficial. The CG(L) al- gorithm is faster than BiCGSTAB(L) for refractive indices n ­ s1.33, 0d, (1.33, 0.01), (1.33, 0.1). Table 2 is for size parameter x ­ 1. All the re- sults are normalized to Bi-CGSTAB(L). This method is clearly superior to the CGNR method, and it is 2±4.3 times faster. It can be seen that CGNR con- verges slowly and has not satisfied the stopping cri- terion in 140 iterations for n ­ s3, 0.0001d. For this larger value of size parameter the QMR algorithm does not converge well, but its smooth version TFQMR does. However, TFQMR is slower than Bi-CGSTAB(L) and comparable with CGNR. The CG(L) method for re- fractive index n ­ s1.33, 0d and n ­ 1.33, 0.01 is faster than the reference scheme Bi-CGSTAB(L). It can be seen that the Neumann polynomial preconditioning

August 15, 1997 / Vol. 22, No. 16 / OPTICS LETTERS 1207 Table 2. CPU Time (Normalized) and Number of Iterations for x5 1a Method n ­ s1.33, 0 d (1.33, 0.01) (1.33, 1) (2, 0) (3, 0.0001) CGNE 3.2(10) 3.2(10) 2.0(16) 4.5(33) *(140) CGNE (L) 2.7(8) 2.6(8) 1.7(13) 3.8(27) *(140) CGNE(R) 3.2(10) 3.2(10) 2.0(16) 4.6(33) *(140) CGNE(N) (L) 2.6(5) 2.6(5) *(140) *(140) *(140) CGNE(N) (R) 3.6(7) 3.6(7) *(140) *(140) *(140) CGNR 3.5(11) 3.5(11) 2.0(16) 4.3(32) *(140) CGNR(L) 2.7(8) 2.6(8) 1.7(13) 3.7(27) *(140) CGNR(R) 3.5(11) 3.5(11) 2.0(16) 4.4(32) *(140) QMR *(47) *(58) *(25) *(76) *(50) QMR(L) 5.3(12) *(59) *(21) *(71) *(39) QMR(R) *(52) *(63) *(22) *(70) *(37) CG 1.3(8) 1.3(8) *(140) *(140) *(140) CG(L) 0.9(5) 0.9(5) *(140) *(140) *(140) CG(R) 1.3(8) 1.3(8) *(140) *(140) *(140) BiCG *(140) *(140) *(140) *(140) *(140) BiCG(L) *(140) *(140) *(140) *(140) *(140) BiCG(R) *(140) *(140) *(140) *(140) *(140) Bi-CGSTAB 1.3(4) 1.3(4) 1.2(10) 1.2(9) 1.1(24) Bi-CGSTAB(L) 1.0(3) 1.0(3) 1.0(8) 1.0(7) 1.0(21) Bi-CGSTAB(R) 1.3(4) 1.3(4) 1.2(10) 1.3(9) 1.1(24) Bi-CGSTAB(N) (L) 1.4(2) 1.4(2) 1.5(6) 4.6(17) *(140) Bi-CGSTAB(N) (R) 1.5(2) 1.5(2) 1.6(6) *(140) *(140) TFQMR 3.3(6) 3.3(6) 3.3(14) 3.4(13) 3.8(42) TFQMR(L) 2.8(5) 2.8(5) 3.0(13) 3.0(11) 3.7(40) TFQMR(R) 3.4(6) 3.4(6) 3.3(14) 3.4(13) 3.9(42) CGS 1.3(4) 1.3(4) 1.3(11) 1.4(10) 1.6(34) CGS(L) 1.0(3) 1.0(3) 1.2(10) 1.3(9) 1.1(23) CGS(R) 1.3(4) 1.3(4) 1.3(11) 1.4(10) 1.6(33) aThe labels (L), (R), and (N) and the asterisk are defined in text. Bi-CGSTAB(N) (L) or Bi-CGSTAB(N) (R) does reduce the number of iterations needed for certain refractive indices. However, the cost associated with the addi- tional calculations always offsets this improved conver- gence rate. As before, the left Jacobi preconditioner is superior to right- or no-preconditioner cases. CG(L) works well for small refractive index but is comparable with Bi-CGSTAB(L). The QMR algorithm fails to con- verge, but the TFQMR algorithm converges well and is comparable with CGNR. The CG method is theoreti- cally valid for Hermitian positive definite matrices. The matrix that arises in the DDA is not Hermitian but symmetric. Therefore, strictly speaking, the CG method is not valid for use in the DDA. Users are ad- vised to test the CG method when extrapolating results presented here to different size parameters, particle sizes, and refractive-index values. I recommend use of the stabilized version of Bi- CGSTAB algorithm with left Jacobi preconditioning. This algorithm requires 67% greater storage capacity than the CGNR algorithm but is typically 2±3 times faster. The recent version of DDA code DDSCAT5a devel- oped by Draine and Flatau contains improvements documented in this Letter. The code is available by anonymous ftp from ftp.astro.princeton.edu or from the Light Scattering and Radiative Transfer Codes Li- brary, SCATTERLIB (http://atol.ucsd.edu/,pf latau). I have been supported in part by the U.S. Office of Naval Research Y oung Investigator Program and in part by DuPont Corporate Educational Assistance. I thank M. J. Wolff and A. E. Ilin, who helped with com- puter tests. Bruce Draine checked the manuscript. I thank R. J. Riegert of DuPont for his continuing inter- est in DDSCAT developments. References 1. B. Draine and P. Flatau, J. Opt. Soc. Am. A 11, 1491 (1994). 2. Y. L. Yung, Appl. Opt. 17, 3707 (1978). 3. A. G. Hoekstra, ``Computer simulations of elastic light scattering,'' (University of Amsterdam, Amsterdam, 1994). 4. M. R. Hestenes and E. L. Stiefel, J. Res. Natl. Bur. Stand. (U.S.) 49, 435 (1952). 5. J. Rahola, SIAM J. Sci. Comput. 17, 78 (1996). 6. B. Draine, Astrophys. J. 333, 848 (1988). 7. M. Petravic and G. Kuo-Petravic, J. Comput. Phys. 32, 263 (1979). 8. K. Lumme and J. Rahola, Astrophys. J. 425, 653 (1994). 9. R. D. da Cunha and T. R. Hopkins, Appl. Numer. Math. 19, 33 (1995). 10. R. Barrett, M. Berry, T. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, and H. van der Vorst,Templates for the Solution of Linear Sys- tems: Building Blocks for Iterative Methods (Society for Industrial and Applied Mathematics, Philadelphia, Pa., 1994).