How can I successfully solve a certain huge sparse matrix; i.e., 8e6 x 8e6 entries?

Hello everybody,

I have implemented a numerical model which consists--among others--of a Finite Difference scheme that is used to solve the quasi-stationary magnetic diffusion equation in three dimensions; i.e., a modified version of the Laplacian operator in Cartesian coordinates. For solving the generated system of linear equations (A)x = b I use the SOR algorithm. The model ran well so far, no matter what input I chose. I should mention that the system matrix A consists of complex-valued elements.

However, for some input parameters the SOR solver does not deliver useful results, or more specifically, it does not converge to the given level of accuracy at all. The system matrix has 8e6x8e6 entries which shouldn't be a problem because, for most input parameters, the SOR algorithm delivers physically correct results. It has to have something to do with the entries of the system matrix per se.

For solving the equation on the Cartesian grid I use the 7-point-stencil. If the unmodified version of the Laplacian operator is solved, the weighing factor at the actual node (i,j,k) is -6. In my modified version, it is -6 + i*x, where x denotes the imaginary part. I even tried SPQR from the SuiteSparse package provided by Tim Davis but since this is not an iterative method MATLAB throws an "out of memory" error (I can use up to 192 GiB (!)). It would really be great if someone encountered a similar problem and wants to share their solution. If interested, I could upload the system matrix as well as the boundary condition vector.

Thank you very much in anticipation!

Best regards,

Chris Volkmar

Answer 1

[Fabio Freschi]

Sorry for late reply, I am on holiday!

The method is described here

O. Axelsson, A. Kucherov, "Real valued iterative methods for solving complex symmetric linear systems", Numerical Linear Algebra with applications, 2000, 7:197--218

Check Algorithm 2.1 on page 202. I have the draft of a Matlab function if you want a black-box solver

Answer 2

[Fabio Freschi]

Chris, is the matrix complex symmetric (not-hermitian)? If so I think I have a possible solution to try. Fabio

[Chris]

Fabio, the matrix indeed is non-hermitian. I'm really interested in this solution of yours. -Chris

Answer 3

[Chris]

Thanks for your reply! I downloaded the paper and will have a look to the suggested algorithm. Thanks again for this useful advice. I think I'll be able to implement it. Have a nice holiday!

Chris