Re: issue with LinearSolve[] when using SparseArray when
- To: mathgroup at smc.vnet.net
- Subject: [mg113737] Re: issue with LinearSolve[] when using SparseArray when
- From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
- Date: Tue, 9 Nov 2010 03:55:21 -0500 (EST)
On Fri, 5 Nov 2010, Nasser M. Abbasi wrote:
> There is some serious problem I am getting when trying to use
> Mathematica LinearSolve[] to solve Au=f when using SparseArray when the
> size of non-zero elements gets close to 1 GB.
>
> The non-zero elements in A that I am building is of the order n^3.
>
> Therefore, for n=255, memory needed should be around 130 Mega Bytes
> (255^2*8) (using 8 bytes for each value, I am using numerical everything).
>
> For n=511, memory used should be a little over 1 Giga Bytes.
>
> I have 8 GB Ram. Windows 7, new Intel CPU.
>
> When I run the solver for n=31 or 63 or 127 or 255, it all works, and
> Mathematica is fast solving Au=f, takes few seconds, and no problem.
> These are the CPU times reported by Mathemtica Timing command for the
> LinearSolve[] call
>
> n=63, cpu=0.03
> n=127, cpu=0.078
> n=255, cpu=0.437
> n=511, been running for many hrs, memory problem.
>
> When I changed to n=511, I see memory of the Mathematica Kernel going up
> to almost 100% of the PC memory, using almost 7 GB, and I have waited
> for 9 hrs, and Mathematica is still not done.
>
> There seems, on the face of it, something really wrong here. It seems
> SparseArray behavior or how LinearSolve[] uses it, does not scale well
> at all? or is this a windows OS issue? if it is a bug in my code, but it
> works so well for all the other n values?
>
> I am posting the code, it is really small code, to see if someone can
> please try it on their PC and see if they get the same behavior.
>
> ------------ code -----------------
> $MinPrecision=$MachinePrecision;$MaxPrecision=$MachinePrecision;
> Share[];
>
> (* make sparse A*)
> makeA[n_?(IntegerQ[#]&&Positive[#]&)]:=Module[{r,off,block},
> r=Table[-4,{i,n}];
> off=Table[1,{i,n-1}];
>
> block =DiagonalMatrix[r,0]+DiagonalMatrix[off,1]+DiagonalMatrix[off,-1];
>
> SparseArray[{Band[{1,1}]->ConstantArray[block,{n}],Band[{1,n+1}]->1,Band[{n+1,1}]->1}]
> ];
>
> (* f(x,y) *)
> force[i_?(IntegerQ[#]&),j_?(IntegerQ[#]&),h_?(NumericQ[#]&&Positive[#]
> &)]:=Module[{x=i*h,y=j*h},
> N@Exp[-(x-0.25)^2-(y-0.6)^2]
> ];
>
> (*n=127; h=2^-7;*) (*these values are OK *)
> (*n=255; h=2^-8;*) (*these values are OK *)
> n=511; h=2^-9; (*these cause problem *)
>
> A=N[makeA[n]];
>
> (* fill in f vector, in correct order for problem*)
> f=Table[0,{i,n^2}];
> For[j=1,j<=n,j++,
> For[i=1,i<=n,i++,
> f[[j+n*(i-1)]]=force[i,j,h]
> ]
> ];
>
> Print["before solver, MemoryInUse[]=",MemoryInUse[]];
> {cpu,sol}=Timing[LinearSolve[A,f]];
> Print["after solver, MemoryInUse[]=",MemoryInUse[]];
> Print["after solver, cpu=",cpu];
>
> ---------------- end code -----------------
>
> thanks
> --Nasser
>
>
here is my version:
$HistoryLength = 0;
(*f (x,y)*)
force[i_?(IntegerQ[#] &), j_?(IntegerQ[#] &),
h_?(NumericQ[#] && Positive[#] &)] :=
Module[{x = i*h, y = j*h}, N@Exp[-(x - 0.25)^2 - (y - 0.6)^2]];
n = 127; h =2^-7;(*these values are OK*)
(*n=255;h=2^-8;*)(*these values are OK*)
(*n=511;h=2^-9;*)(*these cause problem*)
AbsoluteTiming[
AA = SparseArray[{Band[{1, 1}] -> -4., Band[{2, 1}] -> 1.,
Band[{1, 2}] -> 1., Band[{1, n + 1}] -> 1.,
Band[{n + 1, 1}] -> 1.}, {n^2, n^2}, 0.]]
tmp = Table[i*n, {i, n - 1}];
AA[[tmp, tmp + 1]] = 0.;
AA[[tmp + 1, tmp]] = 0.;
MaxMemoryUsed[]/1024.^2
(*fill in f vector,in correct order for problem*)
f = Table[0, {i, n^2}];
For[j = 1, j <= n, j++,
For[i = 1, i <= n, i++, f[[j + n*(i - 1)]] = force[i, j, h]]];
Print["before solver, MemoryInUse[]=", MemoryInUse[]];
{cpu, sol} = AbsoluteTiming[LinearSolve[AA, f]];
Print["after solver, MemoryInUse[]=", MemoryInUse[]];
Print["after solver, cpu=", cpu];
MaxMemoryUsed[]/1024.^2
Nasser, please check that the Matrix AA is really the same as A. But this
gives you some ideas what to try. If you want to solve even bigger
problems, you could think about Method->"Krylov" or other for LinearSolve.
I hope this helps,
Oliver