Hi Chris,

the SOR method can diverge if the matrix is not self-adjoint positive definite. However the convergence can be very slow even for self-adjoint matrices.

The first thing that you should try is one of the Krylov subspace method, i.e., the Conjugate Gradient method (for self-adjoint positive definite matrices) or the GMRES method (for general matrices). These methods only require the implementation of the matrix-vector product and the scalar product. Thus are easy to implement.

If these methods do not yield a desirable speed of convergence you can still go for more elaborated methods like preconditioned GMRES, Multigrid, AMG, AMLI, etc. But try the simple methods first. ;)

You can find the methods e.g. in Saads book (see link below).

Regards

http://www-users.cs.umn.edu/~saad/books.html

Try compressed row storage.

For the divergence of SOR, try more complicated iterative method, such as ILU or AMG (hypre, LLNL), or their combinations.

Dear Chris,

We dealt with analogous problem and suggest an approach based on using of Woodbury formula. I attached two papers. May be, it will be useful. It is not for the Laplace equation, but for the Stokes equation. But for the Laplace equation , it is also appropriate and even simpler.

Hi Chris,

the SOR method can diverge if the matrix is not self-adjoint positive definite. However the convergence can be very slow even for self-adjoint matrices.

The first thing that you should try is one of the Krylov subspace method, i.e., the Conjugate Gradient method (for self-adjoint positive definite matrices) or the GMRES method (for general matrices). These methods only require the implementation of the matrix-vector product and the scalar product. Thus are easy to implement.

If these methods do not yield a desirable speed of convergence you can still go for more elaborated methods like preconditioned GMRES, Multigrid, AMG, AMLI, etc. But try the simple methods first. ;)

You can find the methods e.g. in Saads book (see link below).

Regards

http://www-users.cs.umn.edu/~saad/books.html

Wow - I'd like to thank all of you for your nice answers! Now I have to start some reading I guess :)

As Rittich writes, I also strongly recommend iterative algorithms such as multigrid or a Krylov method.  For Laplacian-like discretisations you could try my "Simple and fast multigrid solution of Poisson’s equation using diagonally oriented grids". ANZIAM J., 43(E):E1–E36, July 2001. http://journal.austms.org.au/ojs/index.php/ANZIAMJ/ article/view/465.

In Matlab store your matrix as a sparse matrix and then implicitly invoke the Lapack sparse solvers that underlie Matlab's simple "A\b" operator for sparse matrices.

Why you do not try two term iterative scheme based on incomplete LU-factorization. Split your matrix to three matrices instead of two A=M-N-E, where M nonsingular which resemble A to maintain the sparsity of A since you have a structured matrix and to minimize the fill-ins. May you can use a parameter on N and E to accelerate the convergence of the scheme.

Hi Chris, 

I had the same problem one year ago. Just follow the following steps:

1- Decompose the original system Ax=b into real and imaginary parts, i.e., A=a+jb, then the new system matrix becomes; Bx=b

See O. Z. Mehdizadeh and M. Paraschivoiu, “Investigation of a two-dimensional spectral element method for Helmholtz’s equation,” Journal of Computational Physics, vol. 189, pp. 111-129, 2003.

2- in order to make sure that the system is positive definite, multiply by transpose of B.

3- use the iterative method you like.

4- compose the solution into real and imaginary.

5- Enjoy your solution.

I also recommend using preconditioned (!) iterative Krylov-subspace solvers (GMRES, BiCGstab, IDR,....) or sparse-direct approaches. Typically they are perfectly capable of handling complex data, but this might depend on your actual computing environment.

You can either form a matrix A and then compress it using standard techniques for example the CSR (Compressed Sparse Row) or you can directly solve it without matrix storage because you are using a finite difference technique your matrix-vector product 'Ax' is actually your discretised equation. In both forms I would recommend not to go for SOR but for better iterative techniques like BiCGSTAB or GMRES. SOR can be used as a preconditioner in both these techniques for accelerated convergence and SOR preconditioners are easier to implement.

http://netlib.org/linalg/html_templates/node91.html

http://www-users.cs.umn.edu/~saad/books.html

http://www.netlib.org/linalg/html_templates/Templates.html