Operator splitting step

For a diagonal $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ the solution of Eq. 6 is just a division by $\left[1-\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi (1-\epsilon )\right]$ which requires very little effort. For the non-local $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ operator, a matrix equation has to be solved at each step in order to compute ${ J_{\rm new}}$. For the nearest neighbor operator the matrix structure is that of a sparse band matrix where the bandwidth of the matrix is proportional to the square of the maximum number of points along a coordinate and there are a total of 27 non-zero entries in every column of the matrix. This system can be solved by any suitable method. For example, we have implemented a solver based on LAPACK (Anderson et al., 1992) routines. Although this solver works fine for small grids [$(2*16+1)^3$ voxels with a bandwidth of 1123], the memory requirements rise beyond available single-CPU limits with already only $(2*32+1)^3$ voxels with a bandwidth of 4291. For a more realistic grid of $(2*256+1)^3$ voxels the bandwidth is more than 250,000. Therefore, a standard band matrix solver is only useful for comparison and testing on very small grids. Other sparse matrix solvers (for example, Xiaoye & Demmel, 2003) have similar memory scaling properties.

As an alternative to direct solutions of Eq. 6 we have implemented iterative solvers. These solvers have the huge advantage that their memory requirements are minimal to modest. In addition, they can be tailored to the special format of the matrix in Eq. 6 which makes coding the methods very simple. We obtained extremely good results with the Jordan and the Gauss-Seidel methods, which converge rapidly to a relative accuracy of $10^{-10}$. For either method, Ng acceleration turned out to be very useful, reducing the number of iterations significantly. In addition, the Rapido method (Zurmühl & Falk, 1986) which was constructed to solve systems of the form

$$(1-M)x = a$$ (27)

can be used. However, Rapido has strong constraints on the spectral radius of $M$ and this cannot be used for problems that involve substantial scattering. In the test calculations discussed below we have used the Jordan or Gauss-Seidel solvers as they are the fastest of the solvers we have implemented. We verified that all solvers give the same results.