LAPACK Archives

[Lapack] dggev fails on non-singular system

Hello,

  I have done some testing and I have found that the following
patch fixes all the 3x3 and 4x4 matrix systems that I found failures
on:

Replace (in dhgeqz.f):

    629             IF( ( DBLE( MAXIT )*SAFMIN )*ABS( H( ILAST-1, ILAST
) ).LT.
    630      $          ABS( T( ILAST-1, ILAST-1 ) ) ) THEN
    631                ESHIFT = ESHIFT + H( ILAST-1, ILAST ) /
    632      $                  T( ILAST-1, ILAST-1 )
    633             ELSE
    634                ESHIFT = ESHIFT + ONE / ( SAFMIN*DBLE( MAXIT ) )
    635             END IF

with:

    629             IF( ( DBLE( MAXIT )*SAFMIN )*ABS( H( ILAST, ILAST-1
) ).LT.
    630      $          ABS( T( ILAST-1, ILAST-1 ) ) ) THEN
    631                ESHIFT = ESHIFT + 0.736 * H( ILAST, ILAST-1 ) /
    632      $                  T( ILAST-1, ILAST-1 )
    633             ELSE
    634                ESHIFT = ESHIFT + ONE / ( SAFMIN*DBLE( MAXIT ) )
    635             END IF

The idea here is the following: in the previous code, its certainly
possible for H(ILAST-1,ILAST) to equal 0, in which case the test
will pass but ESHIFT will not be modified at all. I have just
replaced it with H(ILAST,ILAST-1) which is guaranteed not to be
zero, since its a subdiagonal element and the previous tests will
search for zeros on the subdiagonal. The 0.736 coefficient is
arbitrary - without it, I still found a small number of convergence
failures on 4x4 integer systems.

Patrick Alken

On Tue, May 01, 2007 at 04:05:11PM +0100, Sven Hammarling wrote:
Dear Patrick,

Thank you for your problem report.  I can confirm that I get the same 
failure.

It is hoped to have improved, more efficient, versions of the 
generalized nonsymmetric solvers in a future release, but in the 
meantime we aim to look at this specific problem.  It does seem possible 
that the exceptional shift strategy needs some work.

Thanks you again.

Sven Hammarling.


<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