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