Path: news.mathworks.com!not-for-mail
From: "Tim Davis" <davis@cise.ufl.edu>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Solving Ax=b for large sparse matrix(1e+6 X 1e+5) and cholinc failure
Date: Tue, 13 Nov 2007 12:10:17 +0000 (UTC)
Organization: University of Florida
Lines: 39
Message-ID: <fhc479$lj5$1@fred.mathworks.com>
References: <fguq3f$ig1$1@fred.mathworks.com> <fgv6cb$2c5$1@fred.mathworks.com> <fgv898$6fk$1@fred.mathworks.com> <fgv96r$oov$1@fred.mathworks.com>
Reply-To: "Tim Davis" <davis@cise.ufl.edu>
NNTP-Posting-Host: webapp-03-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1194955817 22117 172.30.248.38 (13 Nov 2007 12:10:17 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Tue, 13 Nov 2007 12:10:17 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 45902
Xref: news.mathworks.com comp.soft-sys.matlab:437278


"Sterren Latsky" <stez17@hotmail.com> wrote in message 

> Honestly, it's not difficult at all. Let's say you have a
> set of simultaneous equations
...
> 
> >> x = A \ b
> 

That's what I meant by "use sparse QR", since if A is sparse
and rectangular, A\b uses sparse QR.  The problem is
manifold: (1) A\b for this case uses the colmmd preordering,
which for this matrix gives an R with 261e6 nonzeros (R
would require about 3GB of memory), (2) sparse QR in MATLAB
is slow (chol is much faster), and (3) the matrix is
ill-conditioned so chol(A'*A) would give a bad answer, (4)
the matrix is huge.  x=A\b just doesn't cut it here, for
those reasons.

With METIS(A,'col'), R has 115 million nonzeros.
METIS(A'*A) gives R with 128 million
COLAMD(A) gives R with 181 million
AMD(A'*A) gives R with 191 million
COLMMD(A) gives R with 261 million

so it would seem that metis(A,'col') gives the best column
ordering.  You could try

p = metis(A,'col') ;
C = A(:,p) ;
x = (C'*C)\(C'*b) ;
x (p) = x ;

to use a direct method.  This will be faster and use half
the memory as x=A\b.  However, it might be inaccurate since
A is ill-conditioned (I'm trying it now; details to follow).

I will post the matrix in my sparse matrix collection,
shortly.  I would need b as well, to try an iterative method.