I found the permutation that makes the matrix blocktriangular by eyeball,
based on the Kroneckerproductlike spyplot of the matrix:
prm1 =
Columns 1 through 16
1 8 15 22 2 9 16 23 3 10 17
24 4 11 18 25
Columns 17 through 28
5 12 19 26 6 13 20 27 7 14 21 28
This makes it block lower triangular where the first block is 24x24 and
the second 4x4.
Let me suggest a quick fix, perhaps Beresford can comment. Right now
the inner loop of the code basically does this:
for i = 1 to n
compute C = norm of column i *omitting the diagonal*
compute R = norm of row i *omitting the diagonal*
choose a scale factor F so that multiplying row i by F
and column i by 1/F makes their new norms more nearly equal
end for
If we changed the definition of C and R to include the diagonal entry
(very little change to the code), then we would stop balancing when
the off diagonal parts were both somewhat smaller than the diagonal,
which is enough for improved backward error.
For this matrix, where the diagonal is O(1), the problem would probably
go away (unless there is something more subtle about the stopping
criterion that I don't see). If the diagonal were really tiny or 0, then
this suggested trick would not help. More aggressively, we could
add the maximum diagonal entry to C and R. But this would "turn off"
balancing in the case of graded matrices, when it may be most useful.
The most general fix would be strongly connected components, a lot
more work.
Jim
Julien Langou wrote:
On Thu, 22 Apr 2010, James Demmel wrote:
There is a permutation that makes this matrix 2x2 block triangular,
Mmm, I can not find this permutation ... I tried symamd and colamd and
they can not either.
Note that while
[V,D]=eig(A);
does not work,
[Q,H]=hess(A);
[V,D]=eig(H);
V = Q*V;
works ...
Definetley balancing ... I've just remembered that actually you can
desactivate balancing in eig of matlab, all relative backward errors
are back to 1e16 with [V,D]=eig(A,'nobalance');.
Julien.
