Re: issue with LinearSolve[] when using SparseArray when
- To: mathgroup at smc.vnet.net
- Subject: [mg113739] Re: issue with LinearSolve[] when using SparseArray when
- From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
- 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
Mem:={MaxMemoryUsed[]/1024.^2,MemoryInUse[]/1024.^2}
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
tutorial/LinearAlgebraPerformance
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:
>
> http://12000.org/tmp/sparse_matrix_issue/
>
> 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.
http://library.wolfram.com/infocenter/Conferences/7549/
Have fun solving PDEs,
Oliver