LAPACK Archives

[Lapack] [lapackers] bug in DGEHRD

I found the permutation that makes the matrix block-triangular by eyeball,
based on the Kronecker-product-like 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 1e-16 with [V,D]=eig(A,'nobalance');.

Julien.



<Prev in Thread] Current Thread [Next in Thread>


For additional information you may use the LAPACK/ScaLAPACK Forum.
Or one of the mailing lists, or