Is there a reason that a direct method (e.g. cyclic reduction) can be less accurate than ADI in solving the linear algebraic equations from a separable elliptic differential equation (whose exact analytical solution is used in calculating the error).. e.g. the larger operation cound and matrix conditioning issues affecting the round-off error in cyclic reduction and offsetting the advantage of having no iterative convergence error...

Can this observation, my observations on compiler dependence and architecture dependence in previous messages in this forum stem from a bug in programming my code..

Incidentally I am using a fishpack subroutine for cyclic reduction (I don't understand CR to the degree of coding it) from :

http://www.scd.ucar.edu/css/software/fishpack/

As you speculate, direct methods can indeed have larger roundoff errors than iterative methods. Often, the diagonal dominance of the effective iteration matrix of an iterative method is better than that of the matrix being directly solved. Also, the iterative method involves fewer arithmetic operations. Therefore, the iterative method upon being sufficiently converged, can yield a solution with less roundoff error than the direct method.

With either method, you can improve the solution by the technique of iterative improvement, where the system is put in the so-called delta form, with the residual vector as the right-hand side. Of course, this residual vector must be computed using higher-precision arithmetic. As an aside, given the floating-point architecture of modern cpus, it costs you more clock cycles to do 32-bit floating point computations than to do 64-bit ones.

You should note that it is incorrect to use the exact analytical solution to the differential equation to calculate the numerical error in order to draw any conclusions about roundoff or iterative convergence error levels. You must use the exact solution to the discretized equations, which of course begs the question of how you should obtain it. In general, you can only use a better approximation to the exact discrete solution, perhaps from a 128-bit precision calculation. As an alternative you can create an exact solution yourself (perhaps the exact analytical solution you have right now), but you must then create a modified system of discretized equations by subjecting that exact solution to the discretized differential operator. This will yield a non-zero source term in your discretized system, representing the discretization error of the discrete operator.

Thanks for the excellent explanation.