Thread Subject: Large sparse matrix

Hello,

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

Backslash is not more viable for a large problem since I got "out of memory", 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 I could use?

I have 4G in RAM

Thank you!

"Hugo" wrote in message <isrkdc$9p0$1@newscl01ah.mathworks.com>...
> Hello,
>
> I want to solve a large system Ax=b, the size of A is 300,000, the matrix is sparse, square, symmetrical and diagonal dominant.
======================

The following iterative method sometimes works well for diagonally dominant systems

absA=abs(A);
C=absA*sum(absA,2);
x=initialValues;

for ii=1:numIter

 x=x-A.'*(A*x-b)./C;

end
  
 

"Matt J" wrote in message <iss37m$f27$1@newscl01ah.mathworks.com>...
> "Hugo" wrote in message <isrkdc$9p0$1@newscl01ah.mathworks.com>...
> > Hello,
> >
> > I want to solve a large system Ax=b, the size of A is 300,000, the matrix is sparse, square, symmetrical and diagonal dominant.
> ======================
>
> The following iterative method sometimes works well for diagonally dominant systems
>
> absA=abs(A);
> C=absA*sum(absA,2);
> x=initialValues;
>
> for ii=1:numIter
>
> x=x-A.'*(A*x-b)./C;
>
> end
>

Thanks for your answer, Let me explain more about my problem so you can see where I'm getting "out of memory":

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!!!

"Hugo" wrote in message <it9l40$8se$1@newscl01ah.mathworks.com>...
>
> % 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)
=================

No, it's not clear. First you said Z is diagonal. Then you seem to say it has off-diagonal elements. If Z is truly diagonal then the solution is simply H./diag(Z) and you're done.




>
> 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?
======================


I already suggested a method in my previous post.

"Matt J" wrote in message <itab50$4si$1@newscl01ah.mathworks.com>...
> "Hugo" wrote in message <it9l40$8se$1@newscl01ah.mathworks.com>...
> >
> > % 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)
> =================
>
> No, it's not clear. First you said Z is diagonal. Then you seem to say it has off-diagonal elements. If Z is truly diagonal then the solution is simply H./diag(Z) and you're done.
>
>
>
>
> >
> > 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?
> ======================
>
>
> I already suggested a method in my previous post.


Yes, you're right, is not completely diagonal since it has two elements off the diagonal, I think those elements are the problem when I try to get "luinc", I know the method you suggested works but not for this case..., do you have any other idea?

Thank you very much!!

"Hugo " <hugo@radtronics.com.au> wrote in message
news:it9l40$8se$1@newscl01ah.mathworks.com...
> "Matt J" wrote in message <iss37m$f27$1@newscl01ah.mathworks.com>...
>> "Hugo" wrote in message <isrkdc$9p0$1@newscl01ah.mathworks.com>...
>> > Hello,
>> >
>> > I want to solve a large system Ax=b, the size of A is 300,000, the
>> > matrix is sparse, square, symmetrical and diagonal dominant.
>> ======================
>>
>> The following iterative method sometimes works well for diagonally
>> dominant systems
>>
>> absA=abs(A);
>> C=absA*sum(absA,2);
>> x=initialValues;
>>
>> for ii=1:numIter
>>
>> x=x-A.'*(A*x-b)./C;
>>
>> end
>
> Thanks for your answer, Let me explain more about my problem so you can
> see where I'm getting "out of memory":
>
> N = 700
> Z = spalloc(N*N,N*N,nnz) % Z almost [500,000 x 500,000]

What is nnz?

> H = full vector [500,000 x 1]
> for x = 1:N
> for y = 1:N
> Z = Computation of Z;

I assume you're filling in elements in Z rather than simply overwriting?

Have you considered constructing vectors R, C, and D in the loop and
constructing Z from those vectors later on with sparse(R, C, D)?

> 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?

Depending on what nnz is, it might be easier (and less memory intensive) to
write a function that accepts a vector x and computes Z*x (and possibly
Z'*x) without explicitly creating Z and pass a function handle to that
function into the iterative solvers like GMRES or LSQR. Given the
tridiagonal form you described above, that would be my next step. Look at
example 2 in the reference page for LSQR for an example you can adapt.

--
Steve Lord
slord@mathworks.com
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

"Hugo" wrote in message <itabv8$7pk$1@newscl01ah.mathworks.com>...
>
>do you have any other idea?
===================

That would depend on what you mean by "I know the method you suggested works but not for this case". What exactly happens when you tried it?