Re: issue with LinearSolve[] when using SparseArray when
- To: mathgroup at smc.vnet.net
- Subject: [mg113726] Re: issue with LinearSolve[] when using SparseArray when
- From: "Nasser M. Abbasi" <nma at 12000.org>
- Date: Tue, 9 Nov 2010 03:53:17 -0500 (EST)
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 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.
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.
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 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)
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[]
--------------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.;
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