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: Thu, 8 Nov 2007 14:27:23 +0000 (UTC)
Organization: University of Florida
Lines: 236
Message-ID: <fgv6cb$2c5$1@fred.mathworks.com>
References: <fguq3f$ig1$1@fred.mathworks.com>
Reply-To: "Tim Davis" <davis@cise.ufl.edu>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1194532043 2437 172.30.248.35 (8 Nov 2007 14:27:23 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Thu, 8 Nov 2007 14:27:23 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 45902
Xref: news.mathworks.com comp.soft-sys.matlab:436578


"Alessio Rucci" <alessio.rucci@libero.it> wrote in message
<fguq3f$ig1$1@fred.mathworks.com>...
> I have to solve Ax=b where A is a  large sparse matrix (1e+6
> x 1e+5) with about
> 
> 7e+6 non-zero elements.
> 
>  
> 
> I tried to solve A'Ax=A'b using pcg, but the new matrix
> (A'A) is very bad conditioned. I tried to use cholinc to
> precondition the system, but A&#8217;A has about 1e+7 non-zero
> elements, and Matlab gives me out of memory error. 
> 
>  
> 
> I tried also to solve the system using:
> 
> x=(A&#8217;A)\(A&#8217;b)
> 
>  
> 
> but here it is what I get:
> 
>  
> 
> sp\: bandwidth = 5864+1+5864.
> 
> sp\: is A diagonal? no.
> 
> sp\: is band density (0.01) > bandden (0.50) to try banded
> solver? no.
> 
> sp\: is A triangular? no.
> 
> sp\: is A morally triangular? no.
> 
> sp\: is A a candidate for Cholesky (symmetric, real positive
> diagonal)? yes.
> 
> sp\: is CHOLMOD's symbolic Cholesky factorization (with
> automatic reordering) successful? yes.
> 
> sp\: is CHOLMOD's numeric Cholesky factorization
successful? no.
> 
>  
> 
> CHOLMOD version 1.1.0, May 5, 2006: : status: ERROR, out of
> memory
> 
>   Architecture: Linux
> 
>     sizeof(int):      4
> 
>     sizeof(UF_long):  4
> 
>     sizeof(void *):   4
> 
>     sizeof(double):   8
> 
>     sizeof(Int):      4 (CHOLMOD's basic integer)
> 
>     sizeof(BLAS_INT): 4 (integer used in the BLAS)
> 
>   Results from most recent analysis:
> 
>     Cholesky flop count: 8.0627e+11
> 
>     Nonzeros in L:       1.9098e+08
> 
>   memory blocks in use:          11
> 
>   memory in use (MB):           8.8
> 
>   peak memory usage (MB):    1961.0
> 
>   maxrank:    update/downdate rank:   8
> 
>   supernodal control: 1 40 (supernodal if flops/lnz >= 40)
> 
>   nmethods=0: default strategy:  Try user permutation if
> given.  Try AMD.
> 
>     Select best ordering tried.
> 
>     method 0: user permutation (if given)
> 
>     method 1: AMD (or COLAMD if factorizing AA')
> 
>         prune_dense: for pruning dense nodes:   10
> 
>         a dense node has degree >= max(16,(10)*sqrt(n))
> 
>         flop count: 8.0627e+11
> 
>         nnz(L):     1.9098e+08
> 
>   final_asis: TRUE, leave as is
> 
>   dbound:  LDL' diagonal threshold:  0
> 
>     Entries with abs. value less than dbound are replaced
> with +/- dbound.
> 
>   grow0: memory reallocation:  1.2
> 
>   grow1: memory reallocation:  1.2
> 
>   grow2: memory reallocation: 5
> 
>   nrelax, zrelax:  supernodal amalgamation rule:
> 
>     s = # columns in two adjacent supernodes
> 
>     z = % of zeros in new supernode if they are merged.
> 
>     Two supernodes are merged if (s <= 4) or (no new zero
> entries) or
> 
>     (s <= 16 and z < 80%) or (s <= 48 and z < 10%) or (z < 5%)
> 
>   OK
> 
> sp\: But A is real symmetric; try MA57.
> 
> MA57 Control Parameters
> 
>     Ordering used:    AMD ordering using MC47.
> 
>     Block Size for Level 3 BLAS:                    16
> 
>     Merge assembly tree nodes if number of eliminations
> 
>     is less than:                                   16
> 
>     Matrix is scaled using symmetrized version of HSL code
MC64.
> 
>     Maximum iterative refinement steps:             10
> 
>     Pivot Threshold:                               
1.000000e-02
> 
>     Zero Pivot Tolerance:                          
1.000000e-20
> 
> sp\: is MA57's symbolic analysis successful?  yes.
> 
>     Forecast number of Real entries in factors:     191653054
> 
>     Forecast number of Integer entries in factors:  1466442
> 
>     Forecast maximum frontal size:                  8071
> 
>     Number of nodes in Assembly tree:               5001
> 
>     Forecast length of FACT array (real)
> 
>     without numerical pivot:                        288897514
> 
>     Forecast length of ifact array (integer)
> 
>     without numerical pivot                         15335967
> 
>     Length of fact required for a successful
> 
>     completion of the numerical phase allowing
> 
>     data compression (without numerical pivoting)   284102031
> 
>     Length of ifact required for a successful
> 
>     completion of the numerical phase allowing
> 
>     data compression (without numerical pivoting)   15335967
> 
>     Number of data compresses                       0
> 
>     Forecast number of floating point additions
> 
>     for the assembly processes                     
6.036817e+08
> 
>     Forecast number of floating point operations
> 
>     to perform the elimination operations
> 
>     counting multiply-add pair as two operations   
7.745417e+11
> 
>     *** Return code from MA57AD:  1
> 
> ??? Error using ==> mldivide
> 
> Out of memory. Type HELP MEMORY for your options.
> 
>  
> 
> Can anyone suggest me a solution or a different strategy to
> solve the original problem Ax=b?
> 
> 
> Thank you
>  

My last post seems to have gotten lost.  Or not.
I clicked "submit"  but I got sent to the "update your
profile" page.  So this might be a duplicate.

Upload the matrix to http://www.cise.ufl.edu/~web-gfs
and let me take a crack at it.  I'd like to add it to
my collection (http://www.cise.ufl.edu/research/sparse).
Send A itself, not A'*A, as a mat file.

x=(A'*A)\(A'*b) is not a good idea, since it is very
ill conditioned.  The failure here is CHOLMOD running
out of memory, however.  It's puzzling that backslash
then tries MA57.  CHOLMOD and MA57 both use AMD, and
the latter tends to use a tad bit more memory.  So if
CHOLMOD+AMD runs out of memory, MA57+AMD almost certainly
will.

You should use sparse QR, although that doesn't use
the BLAS and is thus slower than sparse Cholesky.

Try using METIS.  There is an interface to the METIS
ordering method in SuiteSparse on the File Exchange.
METIS can give better orderings than AMD, although it
can segfault sometimes :-(