Solving multiple linear systems where A (the matrix to be inverted) is always the same

I have the following code excerpt:

%suppose A is a 40000 by 40000 matrix
phi = ones(40000,1);  
for t = 1:360
      phivec = -phi;
      phi = full(A\sparse(phivec));
end

For each t, I solve the system A*x = phivec. For t=1, phivec is given. For t>1, phivec comes from the phivec at t-1. Is there a way to speed up the process, since A is the same for each t. I did try to invert A before the loop starts and do a matrix multiplication, but since A is huge (40000*40000), it's taking a lot of memory. Is there a way to invert A for just once to save time?

Answer 1

[John D'Errico]

This is incredibly silly. Sorry, but it is.

On the first iteration through, you solve the linear system

phi = A\zeros(40000,1);

Then you essentially solve the problem

phi = A\(-phi);

Sorry, but simply ridiculous. What is the solution to the first iteration? ALL ZEROS. What is -1 times a vector of zeros? ZERO.

It matters not how many times you solve the same homogeneous linear system. The solution will be all zeros.

Now, if you decide to start with a non-zero vector phi, it starts to look like the inverse power method for eigenvalues. Not quite. But close.

[jh2011]

Sorry, I meant to say a vector of all ones.

[John D'Errico]

So factorize A, OUTSIDE of the loop. Use a pivoted LU factorization. Or, I don't care, invert A. Yes, it will require a massive amount of memory in either case, because any such factorization of a sparse system tends to be pretty full.

[jh2011]

The excessive use of memory makes my computer really slow. What can I do to get around it? I only want to invert the matrix once since it's the same for the entire loop.

[John D'Errico]

Yep. And so? Big problems require big computers. Get more memory. Find a bigger, faster computer. Or solve smaller problems. Do you really expect more than that? Seriously. There is little magic here.

You might try a sparse solver, something like lsqr, gmres, etc. Because these iterative solvers never factorize the system at all, each iteration requires not much more than some matrix multiplies. And since you can provide an initial guess, if each iteration is converging on a solution, you can gain because the solver has a relatively small distance to go on each iteration.

A problems is you will probably need to learn about preconditioning tools like ilu to make this work well. When I was doing this type of work, I found that choosing a good preconditioner was quite valuable.