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