&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

"Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by Jim Demmel
and Guillaume Revy. See forum post 1783.

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

Fix problem in xTGEX2. The threshold value is too stringent and some matrices
are failing. Relax the threshold by a factor 2. More below.

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

See: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=1783
dgges issue >> by OndraKamenik >> Tue Mar 16, 2010 6:01 am

All,

I have the following problem with dgges. For version 3.1.1 and sooner, I get a
reasonable result, for version 3.2 and 3.2.1 I get info=n+2.

I am separating eigenvalues in the unit circle from one outside the unique
circle.

The two D and E matrices have relatively well separated null spaces, the
minimum angle is acos(0.97). However, if i calculate condition numbers of
E-lambda*D of lambda=[-1:0.01:1], they are quite bad.

The matrices are attached with a small c++ program which calls dgges and sorts
eigenvalues for a better comparison. There is also a Makefile, which links with
different version of lapack.

I use the reference blas.

My question is if it is a bug in Lapack introduced between 3.2. and 3.1.1 or
the matrix is just very bad and in version 3.1.1 I was just lucky to get a
reasonable solution.

Many thanks for any help.

Ondra Kamenik

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

Much more conversation on mailing list and forum [ skipped ]

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

Date: Wed, 31 Mar 2010 19:53:11 -0600
From: James Demmel

This bug was introduced by changing the routine that computes Householder
transformations in the last release of LAPACK from dlarfg to the new dlarfp,
which makes the diagonal of the R factor in QR nonnegative. The bug was a
failure (INFO = N+3) in the dgges routine for computing selected deflating
subspaces of a matrix pencil A - lambda*B. I will describe the bug, a quick
fix, and implications for the floating point debugging project that Guillaume
and others of us are working on.

The bug occurred in dtgex2, in the swapping of eigenvalues (adjacent 1x1 and/or
2x2 blocks on the diagonal of the generalized Schur form) in order to compute
a selected subspace. The swapping involves QR decompositions of  small (n <= 4)
matrices.

The code in dtgex2 performs two of its own internal correctness tests (a "weak"
one and a "strong" one) to see if the swapping has been performed stably. The
two tests compute a residual in slightly different ways. The weak test passed,
but the strong test  (which can be commented out by setting the internal
parameter WANDS to be .false.)  failed, leading to returning INFO = N+3.
However, it failed by exceeding the threshold THRESH only by a factor like 1.2
or less. THRESH is set to 10*macheps*dnorm, so if we changed the (probably
somewhat arbitrary) factor from 10 to 20, it would work. Or we could set WANDS
= false.

Since Bo's name is on this routine, his comments are particularly welcome.

Guillaume and I spent a while tracking this down, and will continue to find out
why dlarfp led to a (slightly!) larger residual than the old dlarfg.

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

Date: Thu, 1 Apr 2010 20:01:21 -0600
From: James Demmel <demmel@cs.berkeley.edu>

Actually, on reexamining the data, this example would be fixed by
changing 10 to 11,
but let's go with 20 to be on the safe side :) . I think it would be
most efficient
if Julie or you appropriately change the one line of code in DTGEX2:

      THRESH = MAX( TEN*EPS*DNORM, SMLNUM )

All versions (S/D/C/Z) have an analogous line of code with the same
constant TEN that I would change.

Sorry, I've lost track of the Mathworks bug report on dlarfp. Can you
remind me?

Thanks,
Jim

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

diff --git a/SRC/ctgex2.f b/SRC/ctgex2.f
index dbf9179c0b256a71eaf60949b6aded14dd8145d7..1f71a8bb78fed4cbffb112fa72f45a99c9db9249 100644 (file)
--- a/SRC/ctgex2.f
+++ b/SRC/ctgex2.f
@@ -118,8 +118,8 @@
       COMPLEX            CZERO, CONE
       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
-      REAL               TEN
-      PARAMETER          ( TEN = 10.0E+0 )
+      REAL               TWENTY
+      PARAMETER          ( TWENTY = 2.0E+1 )
       INTEGER            LDST
       PARAMETER          ( LDST = 2 )
       LOGICAL            WANDS
@@ -173,7 +173,16 @@
       CALL CLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M )
       CALL CLASSQ( 2*M*M, WORK, 1, SCALE, SUM )
       SA = SCALE*SQRT( SUM )
-      THRESH = MAX( TEN*EPS*SA, SMLNUM )
+*
+*     THRES has been changed from 
+*        THRESH = MAX( TEN*EPS*SA, SMLNUM )
+*     to
+*        THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
+*     on 04/01/10.
+*     "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
+*     Jim Demmel and Guillaume Revy. See forum post 1783.
+*
+      THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
 *
 *     Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks
 *     using Givens rotations and perform the swap tentatively.

diff --git a/SRC/dtgex2.f b/SRC/dtgex2.f
index 831783aee919ce0971ca1b4ea3079cc7a8df9c57..cb3763e8116700cd0e6506ac3e126f49b4aa465c 100644 (file)
--- a/SRC/dtgex2.f
+++ b/SRC/dtgex2.f
@@ -134,8 +134,8 @@
 *     .. Parameters ..
       DOUBLE PRECISION   ZERO, ONE
       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
-      DOUBLE PRECISION   TEN
-      PARAMETER          ( TEN = 1.0D+01 )
+      DOUBLE PRECISION   TWENTY
+      PARAMETER          ( TWENTY = 2.0D+01 )
       INTEGER            LDST
       PARAMETER          ( LDST = 4 )
       LOGICAL            WANDS
@@ -205,7 +205,16 @@
       CALL DLACPY( 'Full', M, M, T, LDST, WORK, M )
       CALL DLASSQ( M*M, WORK, 1, DSCALE, DSUM )
       DNORM = DSCALE*SQRT( DSUM )
-      THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
+*
+*     THRES has been changed from 
+*        THRESH = MAX( TEN*EPS*SA, SMLNUM )
+*     to
+*        THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
+*     on 04/01/10.
+*     "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
+*     Jim Demmel and Guillaume Revy. See forum post 1783.
+*
+      THRESH = MAX( TWENTY*EPS*DNORM, SMLNUM )
 *
       IF( M.EQ.2 ) THEN
 *

diff --git a/SRC/stgex2.f b/SRC/stgex2.f
index 99599e6352e4ddd7516596d54e3ffced81a4248a..ab4b5e468c5cc98c929a47e30d1ced5a0b601fc3 100644 (file)
--- a/SRC/stgex2.f
+++ b/SRC/stgex2.f
@@ -134,8 +134,8 @@
 *     .. Parameters ..
       REAL               ZERO, ONE
       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
-      REAL               TEN
-      PARAMETER          ( TEN = 1.0E+01 )
+      REAL               TWENTY
+      PARAMETER          ( TWENTY = 2.0E+01 )
       INTEGER            LDST
       PARAMETER          ( LDST = 4 )
       LOGICAL            WANDS
@@ -205,7 +205,16 @@
       CALL SLACPY( 'Full', M, M, T, LDST, WORK, M )
       CALL SLASSQ( M*M, WORK, 1, DSCALE, DSUM )
       DNORM = DSCALE*SQRT( DSUM )
-      THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
+*
+*     THRES has been changed from 
+*        THRESH = MAX( TEN*EPS*SA, SMLNUM )
+*     to
+*        THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
+*     on 04/01/10.
+*     "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
+*     Jim Demmel and Guillaume Revy. See forum post 1783.
+*
+      THRESH = MAX( TWENTY*EPS*DNORM, SMLNUM )
 *
       IF( M.EQ.2 ) THEN
 *

diff --git a/SRC/ztgex2.f b/SRC/ztgex2.f
index 95e780ef07108053ee4171f99b6fd8d16e95a036..3bd85ee4e99d649597105bfca0842727983e3d99 100644 (file)
--- a/SRC/ztgex2.f
+++ b/SRC/ztgex2.f
@@ -118,8 +118,8 @@
       COMPLEX*16         CZERO, CONE
       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
-      DOUBLE PRECISION   TEN
-      PARAMETER          ( TEN = 10.0D+0 )
+      DOUBLE PRECISION   TWENTY
+      PARAMETER          ( TWENTY = 2.0D+1 )
       INTEGER            LDST
       PARAMETER          ( LDST = 2 )
       LOGICAL            WANDS
@@ -173,7 +173,16 @@
       CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M )
       CALL ZLASSQ( 2*M*M, WORK, 1, SCALE, SUM )
       SA = SCALE*SQRT( SUM )
-      THRESH = MAX( TEN*EPS*SA, SMLNUM )
+*
+*     THRES has been changed from 
+*        THRESH = MAX( TEN*EPS*SA, SMLNUM )
+*     to
+*        THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
+*     on 04/01/10.
+*     "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
+*     Jim Demmel and Guillaume Revy. See forum post 1783.
+*
+      THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
 *
 *     Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks
 *     using Givens rotations and perform the swap tentatively.