Solving large linear systems

Hello,

I want to solve a large system Ax=b, the size of A is 300,000, the matrix is sparse and diagonal dominant.

Backslash is not more viable, so I'm trying to use "ilu" and "luinc" as a preconditioner within an iterative solver (bicg, gmres, pcg) but I get "out of memory" in the ilu/luinc computation.

Is there any suggestion you could give me? Is there any other function?

Thank you!

[Wolfgang Schwanghart]

What makes you think that backslash is not viable?

[Hugo]

Because I got "out of memory" and the error seemed to be in: "mldivide"

[Andrew Newell]

By "size", do you mean each dimension?

Answer 1

[Andrew Newell]

With the information you have provided, it is difficult to diagnose your problem. Do you get an "out of memory" error if you run the following commands?

n = 300000;
e = rand(n,1);
A = spdiags([e -2*e e], -1:1, n, n);
b = rand(n,1);
y = A\b;

EDIT: Here is code based on the new details you have provided:

N = 700;
B = 1e6*ones(N*N,3);
d = [-N 0 N]; % This might be [-N-1 0 N+1], depending on what "separate" means
Z = spdiags(B,d,N*N,N*N);
H = rand(N*N,1);
Y=Z\H;

This took my computer 2.3 seconds and I had no memory problems.

[Hugo]

Thanks for your answer, I tried your example and I don't have any problem solving it, but I don't know why with my system I get "out of memory". Let me tell you what I have:

N = 700

Z = spalloc(N*N,N*N,nnz) % Z almost [500,000 x 500,000]

H = full vector [500,000 x 1]

for x = 1:N

for y = 1:N

Z = Computation of Z;

end

end

% The matrix Z is diagonal, symmetric, sparse,

% all of its values are "10e+6", its form is three diagonal lines in the diagonal, one in the right side and one in the left side, both last ones separate from the center diagonal 700 points (hope this is clear)

Y = Z\H : Out of memory

[L,U] = luinc(Z,1e-7); : Out of memory

[H,flag,relres,iter,resvec] = pcg(Z,H,1e-6,100,L,U);

Is there any other preconditioner I could use?

Thank you!!!

[John D'Errico]

By the way, if you are generating a tridiagonal banded matrix using a spalloc and a loop, then you need to learn how to use spdiags.

[Hugo]

Hi Andrew,

I tried your example and you're right, I had no problem at all, but I made a mistake and I forgot to say that my matrix is complex, so using your example I ran this one and I got "out of memory".

N = 700;

B = 1e6*ones(N*N,5) + 1i*1e6*ones(N*N,5);

d = [-N -1 0 1 N];

Z = spdiags(B,d,N*N,N*N);

figure(1), spy(Z)

H = 1e6*ones(N*N,1);

Y=Z\H;

The matrix is around 50MB. Do you have any other idea of what can I do? I would really appreciate it!! I tried iterative but building the preconditioner "luinc" gave me again "out of memory"

Thanks!!

[Hugo]

BTW, the matrix is tridiagonal with two band off-diagonals separate N (N=700)