Direct Solver

Hello everyone, I am trying to solve a linear system of equations A*x = b, where the matrix A is complex, symmetric, but NOT hermitian, and is of the order of 1e8 by 1e8 (or larger). The coefficient matrix A is assembled by using finite difference frequency domain method (on Yee grid) along with absorbing boundary condition (perfectly matched layer, about 30 layers) in 3D. I have attempted to solve this monster by using MATLAB efficient direct solver "\" on my work station that has 12 cores processor with 64 bit MAC OSX 10.6 operating system and 12 GB of RAM.

However, I have not had any success yet. The directly solver immediately uses 12 GB of RAM and then the computer starts transferring data to hard drive, I believe, which makes it run really slow, and due to this, it runs into the risk of being un- responsive after some time. In order to avoid this, I have never waited long enough to let the direct solver finish its job.

So, I am thinking there may be a better way since the matrix is highly sparse and most of the elements are located along the diagonal. I have also tried ilu + gimres iterative solver, and while this method of preconditioning and iteratively solving the equation converges, the solution is nothing but all zero, so it is not converging to the right solution. Also I can only use very low tolerance of 1e-4. On increasing the tolerance, gimres complains about the condition number of the matrix, which obviously is high, and then stagnates. I had the exact same problem with this combination of solver(gimres+ilu) in 2D also, where obviously the matrix size was smaller(1e6 square) than what I have here. In 2D, however, I was able to get correct solution in a minute or two using MATLAB's direct solver ("\") . So I really cannot trust ilu + gimres combo, and would like to get the direct solver work. With 12 GB RAM, I guess I should be fine memory wise for a sparse matrix of size I mentioned above, or is this not true? I will truly appreciate if any one here could suggest me a way to improve performance of MATLAB direct solver. I would also be interested to know what other iterative solvers (along with a preconditioner)may work on a monster like this. I have tried qmr without preconditioning, but that did not work.

Thanks for any help in advance.

[Walter Roberson]

Perhaps you could reformat the question to include some paragraph breaks? It is not easy to read that much continuous technical text.

[Efim]

Thanks. I hope it is slightly better now.

Answer 1

[Christine Tobler]

The direct sparse solver (\) unfortunately doesn't have methods for complex symmetric or hermitian problems, only for real symmetric and general nonsymmetric.

Most likely, the problem is with fill-in in the factorization. The order of the rows and columns of a sparse matrix has a strong influence on how many nonzeros its factorization will have, which translates directly to memory usage. Backslash makes a heuristic choice of some reordering, but this may not be the optimal one.

You can explicitly apply a reordering: Here's the list of available reordering algorithms (section "Reordering algorithms"), and example here's an additional example. Try applying sum(symbfact(A(p, p))) to get a quick estimate of the number of nonzeros needed in backslash when using the permutation p to reorder A. If one seems promising, you can apply x(p) = A(p, p)\b(p) instead of your original call x = A\b.

Another thing you can do: Represent your n-by-n complex symmetric system as a 2n-by-2n real symmetric system, solving for a larger sparse matrix but one that is symmetric. It's hard to predict if the symmetric will make up for the larger system, just something to try.

Here's a small example of how to do this:

>> A = randn(10) + 1i*randn(10); A = A + A.';
>> issymmetric(A)
ans =
  logical
   1
>> b = randn(10, 1) + 1i*randn(10, 1);
>> x = A \ b;

>> M = [imag(A) real(A); real(A) -imag(A)];
>> c = [imag(b); real(b)]
>> issymmetric(M)
ans =
  logical
   1
>> y = M \ c;
>> x2 = y(1:10) + 1i*y(11:end);
>> norm(x - xx)
ans =
   8.1146e-16

Note at this point, you would want to apply the reordering algorithm to M, not to the original matrix A, since M would not have a nice structure just because A has.

Answer 2

[Andrew Newell]

You need to make your matrix explicitly sparse, using commands like sparse or spdiags. The latter would take care of your finite difference formulae and then you could add the boundary conditions. That would reduce the number of elements to something like 3e8 or 5e8, and you could solve the matrix equation rapidly with the backslash operator.

[Efim]

Thanks for your input. I have not stored my matrix as a full. I do have sparse and spdiags everywhere. On top of that, my MATLAB code is vectorized, and does not have any for loop.

I will appreciate any other input in this regard.

Thanks for your time.