The matrix is square, so it's not sparse QR. You can ask MATLAB what it is using via:
spparms('spumoni',2)
Try this:
>> A=eye(5); A(1,5)=1; A(5,1)=1;
>> A
A =
1 0 0 0 1
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
1 0 0 0 1
>> inv(A)
Warning: Matrix is singular to working precision.
ans =
Inf Inf Inf Inf Inf
Inf Inf Inf Inf Inf
Inf Inf Inf Inf Inf
Inf Inf Inf Inf Inf
Inf Inf Inf Inf Inf
>> inv(sparse(A))
ans =
(1,1) 2
(5,1) 2
(2,2) 1
(3,3) 1
(4,4) 1
(1,5) 2
(5,5) 3
>> spparms('spumoni',2)
>> inv(sparse(A))
sp\: bandwidth = 4+1+4.
sp\: is A diagonal? no.
sp\: is band density (0.28) > 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.
sp\: A is not positive definite.
CHOLMOD version 1.7.0, Sept 20, 2008: : status: warning, matrix not positive definite
Architecture: unknown
sizeof(int): 4
sizeof(UF_long): 8
sizeof(void *): 8
sizeof(double): 8
sizeof(Int): 8 (CHOLMOD's basic integer)
sizeof(BLAS_INT): 8 (integer used in the BLAS)
Results from most recent analysis:
Cholesky flop count: 8
Nonzeros in L: 6
memory blocks in use: 13
memory in use (MB): 0.0
peak memory usage (MB): 0.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
nnz(L): 6
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.000000e02
Zero Pivot Tolerance: 1.000000e20
sp\: is MA57's symbolic analysis successful? yes.
Forecast number of Real entries in factors: 6
Forecast number of Integer entries in factors: 16
Forecast maximum frontal size: 2
Number of nodes in Assembly tree: 4
Forecast length of FACT array (real)
without numerical pivot: 39
Forecast length of ifact array (integer)
without numerical pivot 44
Length of fact required for a successful
completion of the numerical phase allowing
data compression (without numerical pivoting) 39
Length of ifact required for a successful
completion of the numerical phase allowing
data compression (without numerical pivoting) 44
Number of data compresses 0
Forecast number of floating point additions
for the assembly processes 6.000000e+00
Forecast number of floating point operations
to perform the elimination operations
counting multiplyadd pair as two operations 8.000000e+00
*** Return code from MA57AD: 1
sp\: is MA57's numeric Factorization successful? yes.
Number of entries in factors: 6
Storage for real data in factors: 6
Storage for integer data in factors: 19
Minimum length of fact required for a successful
completion of the numerical phase: 39
Minimum length of ifact required for a successful
completion of the numerical phase: 44
Minimum length of fact required for a successful
completion of the numerical phase without
data compression: 39
Minimum length of ifact required for a successful
completion of the numerical phase without
data compression: 44
Order of the largest frontal matrix: 2
Number of 2x2 numerical pivots: 0
Total number of fullysummed variables passed to
the father node because of numerical pivot: 1
Number of negative eigenvalues: 0
Rank of factorization of the matrix: 4
Pivot step where static pivot commences and
is set to zero if no modification performed: 0
Number of data compresses on real
data structures: 0
Number of data compresses on integer
data structures: 0
Number of block pivots in factors 5
Number of zeros in triangle of factors: 0
Number of zeros in rectangle of factors: 0
Number of zero columns in rectangle of factors: 0
Number of static pivots chosen: 0
Number of floating point additions
for the assembly processes: 6.000000e+00
Number of floating point operations
to perform the elimination operations
counting multiplyadd pair as two operations: 7.000000e+00
Minimum value of the scaling factor (MC64): 1.000000e+00
Maximum value of the scaling factor (MC64): 1.000000e+00
Maximum modulus of matrix entry: 1.000000e+00
*** Return code from MA57BD: 4
sp\: is MA57's solve successful? yes.
ans =
(1,1) 2
(5,1) 2
(2,2) 1
(3,3) 1
(4,4) 1
(1,5) 2
(5,5) 3
>> full(ans)
ans =
2 0 0 0 2
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
2 0 0 0 3
>>
From the MA57 spec sheet, at http://www.hsl.rl.ac.uk/specs/ma57.pdf this looks like a warning that the matrix is rank deficient:
+4 Matrix is rank deficient on exit from MA57B/BD. In this case, a decomposition will still have been produced that will enable the subsequent solution of consistent equations. INFO(25) will be set to the rank of the factorized matrix.
And sure enough, it is rank deficient:
>> rank(A)
ans =
4
>> svd(A)
ans =
2
1
1
1
0
(the rank(A) of 4 and the return code of 4 are coincidental). So the bug is not MA57, but the interface to it (MATLAB) that returns inv(A) in spite of the fact that the matrix is rank deficient.
