Rouzaud Denis

Hello again, thanks again.
 
There is an indexing used for the displaymode 'graph' that is used to display progressively the diagonal of the inverse while computation. There was an error in the use of the unique function, which is now corrected. It was leading to missing values in the computed inverse!
 
It seems that the computation works for non symetric matrices. Here are the result I get now.
 
Let me know what you think.
 
Thanks for your help.
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
load bcsstk28
N=Problem.A;
sm = smartinv(N,1,'tri');
sq = inv(N);
dm = sq-sm;
dmax = abs(dm(1,1));
P = [ 1 1];
for i = 2:size(dm)
if abs(dm(i,i))>dmax,dmax=abs(dm(i,i));p=[i i];end
if abs(dm(i-1,i))>dmax,dmax=abs(dm(i-1,i));p=[i-1 i];end
if abs(dm(i,i-1))>dmax,dmax=abs(dm(i,i-1));p=[i i-1];end
end
dm(p(1),p(2))
sq(p(1),p(2))
 
ans =
    2.664812814856532e-013
 
ans =
   0.033834948807814
 
clear
load west0479
N = west0479;
sm = smartinv(N,1,'penta');
sq = inv(N);
dm = sq-sm;
dmax = abs(dm(1,1));
P = [ 1 1];
if abs(dm(1,2))>dmax,dmax=abs(dm(1,2));p=[1 2];end
if abs(dm(2,1))>dmax,dmax=abs(dm(2,1));p=[2 1];end
for i = 3:size(dm)
if abs(dm(i,i))>dmax,dmax=abs(dm(i,i));p=[i i];end
if abs(dm(i-1,i))>dmax,dmax=abs(dm(i-1,i));p=[i-1 i];end
if abs(dm(i,i-1))>dmax,dmax=abs(dm(i,i-1));p=[i i-1];end
if abs(dm(i-2,i))>dmax,dmax=abs(dm(i-2,i));p=[i-2 i];end
if abs(dm(i,i-2))>dmax,dmax=abs(dm(i,i-2));p=[i i-2];end
end
dm(p(1),p(2))
sq(p(1),p(2))
 
ans =
    6.366462912410498e-011
 
ans =
    5.218763111339072e+003

There's still something wrong with this submission. I computed the exact inverse of this matrix:

http://www.cise.ufl.edu/research/sparse/matrices/HB/bcsstk28.html

and then extracted the tridiagonal entries of the inversion (entries on the diagonal, subdiagonal, and +1 diagonal).

I then computed S = smartinv(A,1,'tri'), excepting to get the same thing.

But I don't I get a matrix with zeros on the +1 and -1 diagonals, every 15 entries or so. So this revised submission is still not computing what it says it computes.

This function attempts to do the right thing when computing selected entries of the inverse, by solving a sequence of linear equations after doing an LU factorization. However, the implementation is flawed.

One suggestion: the statement [L,U,P]=lu(A) does not allow LU to permute A to improve sparsity. It would be better to use the 4-output LU, [L,U,P,Q]=lu(A), except in your case you'd need to use another name for Q since you have a Q variable already.

You would then need to use the column permutation Q in the forward/backsolve steps.

Including a "waitbar" in a function like this is not a good idea.

There seems to be a bug. The documentation says that s=smartinv(A) returns the diagonal of the inverse, but in fact I get the whole inverse back.

If I do s=smartinv(A,1), I get a diagonal matrix, but not the diagonal entries of the inverse, which is what the documentatioon said I would get. (I tried with this matrix, which is included in MATLAB):

load west0479
A=west0479;

With smartinv(A,1) or smartinv(A,1,'mono'), I should get the diagonal of inv(A) as the result, but I don't. SOmething else is returned.

Finally, the method is very slow. Computing the entire inverse takes 0.13 seconds using S=inv(A). Using S=smartinv(A) I get the same result, but it takes 15 seconds.

In the profiler, all the time is spent in this one line of code:

 Q(q0+1:q0+w,q0+1:q0+w) = insblock;

that's a very bad idea. You should never use sparse submatrix assignment, if at all possible. A better approach is to build the matrix all at once, from the blocks, using the "sparse" command.

So there appears to be at least three serious bugs in the code:
(1) does not perform according to the documentation (smartinv(A) returns inv(A), not diag(inv(A))).
(2) very slow as compared with just "inv"
(3) cannot compute the diagonal of the inverse, in spite of the fact that the documentation says it can do so. Instead it returns something else altogether.

And 2 serious drawbacks:
(1) Not really a bug, but a serious performance issue: doesn't allow for pivoting to maintain sparsity, and is thus unsuitable for large matrices.
(2) Remove the waitbar. If your function is used within other codes, the end user isn't interested in the "Inverting matrix..." message.