LAPACK Archives

[Lapack] [lapackers] bug in DGEHRD

Forwarding David Day's post on the forum.
To whom it may concern,
Matlab worked for me. (R2009b)
The matrix is, modulo absolute perturbations the size of the machine 
precision, a Kronecker product.
A = numpy54;
amax = max( max( abs(A))), % approximately 122
A11 = A(1:7,1:7);
I = eye(4,4);
norm( kron(I,A11) - A )/amax, % approximately eps
[V11,D11] = eig(A11);
V = kron(I,V11);
D = kron(I,D11); % [V,D] = eig(A);
Each eigenvalue has multiplicity 4. The condition number of V
is sensitive to the algorithm used. Balance reduces A11 to
Schur (lower triangular) form.

On Thu, 22 Apr 2010, James Demmel wrote:

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.


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
does not work,
   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');.


<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