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.
|