MathGroup Archive 2010

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: issue with LinearSolve[] when using SparseArray when

  • To: mathgroup at
  • Subject: [mg113739] Re: issue with LinearSolve[] when using SparseArray when
  • From: Oliver Ruebenkoenig <ruebenko at>
  • Date: Wed, 10 Nov 2010 06:27:28 -0500 (EST)

Dear Nasser,

On Mon, 8 Nov 2010, Nasser M. Abbasi wrote:

> On 11/8/2010 5:13 AM, Oliver Ruebenkoenig wrote:
>> 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
> Thanks you so much Oliver.

You are most welcome. I am glad I could help out.

> You do not know how much this helped. I have given up on using Mathematica 
> for my school project due to the sparse matrix performance, and had to 
> implement everything in that other system (starts with M and ends with B) to 
> be able to complete the assignment.
> Using your method to make A was the breakthrough needed. The A matrix from 
> your method and what I was doing are the same in shape and content, but your 
> method is very fast and efficient, not sure why, and all the RAM problem went 
> away.

The RAM problem actually came up during the construction of A. That's what 
I noted and then wrote a different version to construct A which did not 
have this memory problem. Check your method with MaxMemoryUsed[] before 
and after the construction of A.

I quite often define


to keep track of where and when memory is consumed.

> I am still not sure why the way I was making A is slow. May be it is the use 
> of ConstantArray to make the band? I need to find out, so I can avoid this 
> again. If you have more insight on this, it will be good to know for all of 
> us.

I have not really looked at it, since the memory consumption was too large 
in the first place.

> To make it easy to show the problem again, in the hope someone will learn 
> from it, I posted the complete notebook, see below for link.
> The speed improvement was amazing.

I am glad.

> I compared the speed of the Mathematica LinearSolver with the other system 
> "\" solver for sparse A, for different size Matrices, and Mathematica easily 
> now is 2 times as fast: (both systems use the UMFPACK multifrontal direct 
> solver methods when A is sparse in the direct solvers from my understanding)

Spread the word....

> note: n below means A is of size n^ by n^2 and not n by n, but these are 
> mostly zero elements, hence sparse matrix is used) (this is a problem to 
> solve Poisson pde on 2D on unit square)
> The other system, used cputime() to time the direct solver cpu time used, and 
> in Mathematics I used Timing[]

I do not know what cputime() does. Also, my favorite is AbsoluteTiming 
since this is the time I loose in my life while waiting for something. In 
my opinion the only thing relevant.

> --------------CPU speed -------------
> n=31,   CPU(other)=0.015   CPU(Mathematica)=0.016
> n=63,   CPU(other)=0.13    CPU(Mathematica)=0.015
> n=127,  CPU(other)=0.25    CPU(Mathematica)=0.063
> n=255,  CPU(other)=1.54    CPU(Mathematica)=0.36
> n=511,  CPU(other)=5.54    CPU(Mathematica)=2
> n=1023, CPU(other)=27.1    CPU(Mathematica)=14.62
> -----------------------
> Here is the output of the run for Mathematica
> -------------------------------------------------
> size of A={961,961}
> CPU time for making A=0.015999999999998238
> before solver, MemoryInUse[].634727478027344
> direct solver, n=31
> cpu time used=0.016000000000005343 sec
> -------------------------------------------------
> size of A={3969,3969}
> CPU time for making A=0.046999999999997044
> before solver, MemoryInUse[].893272399902344
> direct solver, n=63
> cpu time used=0.014999999999993463 sec
> -------------------------------------------------
> size of A={16129,16129}
> CPU time for making A=0.2810000000000059
> before solver, MemoryInUse[]=21.933853149414062
> direct solver, n=127
> cpu time used=0.06299999999999528 sec
> -------------------------------------------------
> size of A={65025,65025}
> CPU time for making A=1.3100000000000165
> before solver, MemoryInUse[]=26.124313354492188
> direct solver, n=255
> cpu time used=0.3580000000000041 sec
> -------------------------------------------------
> size of A={261121,261121}
> CPU time for making A=5.492000000000004
> before solver, MemoryInUse[]=42.942710876464844
> direct solver, n=511
> cpu time used=2.0120000000000005 sec
> -------------------------------------------------
> size of A={1046529,1046529}
> CPU time for making A=21.278999999999996
> before solver, MemoryInUse[]=110.32942962646484
> direct solver, n=1023
> cpu time used=14.61800000000001 sec
> -------------------------------------------
> 	And finally, here are the 2 functions to make sparse A, one is very 
> fast and one was slow.
> ------------- FAST sparse maker, due to Oliver Ruebenkoenig -------
> getSparseALaplace2DfastMethod[n_?(IntegerQ[#]&&Positive[#]&)]:=Module[{i,tmp,A},
> A=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}];
> A[[tmp,tmp+1]]=0.;
> A[[tmp+1,tmp]]=0.;

This is very fast. I first create the general structure of the 
matrix and then assign to the matrix with Part. Part[] in general is a 
"good" function to use, since it deals well with packed arrays and is 
compilable, accepts vectorized input.

Speaking about packed arrays, you might want to read this


Also On["Packing"] and Off["Packing"] are useful friends to track down 
performance / memory issues.

> A
> ]
> ------------- SLOW sparse maker, but why slow, thanks to me --------
> getSparseALaplace2DslowMethod[n_?(IntegerQ[#]&&Positive[#]&)]:=Module[{r,off,block},
> (* Make the block along the A matrix diagonal*)
> 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
> }
> ]
> ]
> -------------------------------------------------------
> If someone like to have the full notebook itself, you can download it from 
> here:
> Thanks again Oliver
> --Nasser

If you do PDE numerics you might want to read this talk I gave a while 
back. It has some other fast code that is useful in the PDE world.

Have fun solving PDEs,

  • Prev by Date: Re: Balance point of a solid
  • Next by Date: Repost (Boundary Condition NDSolve)
  • Previous by thread: Re: issue with LinearSolve[] when using SparseArray when
  • Next by thread: 2-D Fourier Transform?