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. `````` ```
 Current Thread [Lapack] bug in DGEHRD, Julien Langou [Lapack] [lapackers] bug in DGEHRD, James Demmel [Lapack] [lapackers] bug in DGEHRD, James Demmel [Lapack] [lapackers] bug in DGEHRD, Julien Langou [Lapack] [lapackers] bug in DGEHRD, James Demmel <= [Lapack] [lapackers] bug in DGEHRD, Julien Langou [Lapack] [lapackers] bug in DGEHRD, Julien Langou

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