issue with LinearSolve[] when using SparseArray when size close
- To: mathgroup at smc.vnet.net
- Subject: [mg113622] issue with LinearSolve[] when using SparseArray when size close
- From: "Nasser M. Abbasi" <nma at 12000.org>
- Date: Fri, 5 Nov 2010 05:11:00 -0500 (EST)
- Reply-to: nma at 12000.org
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