From b4873c19f60844641cceacfb8281abe13e2bd4f0 Mon Sep 17 00:00:00 2001 From: thijs Date: Thu, 4 Feb 2021 15:57:39 +0100 Subject: [PATCH] i remember thinking to myself, i don't need to handle this case, it should never happen... i was wrong --- SRC/claqr5.f | 62 +++++++++++++++++++++++++++------------------------- SRC/dlaqr5.f | 62 ++++++++++++++++++++++++++++------------------------ SRC/slaqr5.f | 62 ++++++++++++++++++++++++++++------------------------ SRC/zlaqr5.f | 54 +++++++++++++++++++++++---------------------- 4 files changed, 126 insertions(+), 114 deletions(-) diff --git a/SRC/claqr5.f b/SRC/claqr5.f index 314382982..71f26d8c9 100644 --- a/SRC/claqr5.f +++ b/SRC/claqr5.f @@ -462,37 +462,39 @@ * . is zero (as done here) is traditional but probably * . unnecessary. ==== * - IF( H( K+1, K ).NE.ZERO ) THEN - TST1 = CABS1( H( K, K ) ) + CABS1( H( K+1, K+1 ) ) - IF( TST1.EQ.RZERO ) THEN - IF( K.GE.KTOP+1 ) - $ TST1 = TST1 + CABS1( H( K, K-1 ) ) - IF( K.GE.KTOP+2 ) - $ TST1 = TST1 + CABS1( H( K, K-2 ) ) - IF( K.GE.KTOP+3 ) - $ TST1 = TST1 + CABS1( H( K, K-3 ) ) - IF( K.LE.KBOT-2 ) - $ TST1 = TST1 + CABS1( H( K+2, K+1 ) ) - IF( K.LE.KBOT-3 ) - $ TST1 = TST1 + CABS1( H( K+3, K+1 ) ) - IF( K.LE.KBOT-4 ) - $ TST1 = TST1 + CABS1( H( K+4, K+1 ) ) - END IF - IF( CABS1( H( K+1, K ) ).LE.MAX( SMLNUM, ULP*TST1 ) ) - $ THEN - H12 = MAX( CABS1( H( K+1, K ) ), - $ CABS1( H( K, K+1 ) ) ) - H21 = MIN( CABS1( H( K+1, K ) ), - $ CABS1( H( K, K+1 ) ) ) - H11 = MAX( CABS1( H( K+1, K+1 ) ), - $ CABS1( H( K, K )-H( K+1, K+1 ) ) ) - H22 = MIN( CABS1( H( K+1, K+1 ) ), - $ CABS1( H( K, K )-H( K+1, K+1 ) ) ) - SCL = H11 + H12 - TST2 = H22*( H11 / SCL ) + IF( K.GE.KTOP) THEN + IF( H( K+1, K ).NE.ZERO ) THEN + TST1 = CABS1( H( K, K ) ) + CABS1( H( K+1, K+1 ) ) + IF( TST1.EQ.RZERO ) THEN + IF( K.GE.KTOP+1 ) + $ TST1 = TST1 + CABS1( H( K, K-1 ) ) + IF( K.GE.KTOP+2 ) + $ TST1 = TST1 + CABS1( H( K, K-2 ) ) + IF( K.GE.KTOP+3 ) + $ TST1 = TST1 + CABS1( H( K, K-3 ) ) + IF( K.LE.KBOT-2 ) + $ TST1 = TST1 + CABS1( H( K+2, K+1 ) ) + IF( K.LE.KBOT-3 ) + $ TST1 = TST1 + CABS1( H( K+3, K+1 ) ) + IF( K.LE.KBOT-4 ) + $ TST1 = TST1 + CABS1( H( K+4, K+1 ) ) + END IF + IF( CABS1( H( K+1, K ) ) + $ .LE.MAX( SMLNUM, ULP*TST1 ) ) THEN + H12 = MAX( CABS1( H( K+1, K ) ), + $ CABS1( H( K, K+1 ) ) ) + H21 = MIN( CABS1( H( K+1, K ) ), + $ CABS1( H( K, K+1 ) ) ) + H11 = MAX( CABS1( H( K+1, K+1 ) ), + $ CABS1( H( K, K )-H( K+1, K+1 ) ) ) + H22 = MIN( CABS1( H( K+1, K+1 ) ), + $ CABS1( H( K, K )-H( K+1, K+1 ) ) ) + SCL = H11 + H12 + TST2 = H22*( H11 / SCL ) * - IF( TST2.EQ.RZERO .OR. H21*( H12 / SCL ).LE. - $ MAX( SMLNUM, ULP*TST2 ) )H( K+1, K ) = ZERO + IF( TST2.EQ.RZERO .OR. H21*( H12 / SCL ).LE. + $ MAX( SMLNUM, ULP*TST2 ) )H( K+1, K ) = ZERO + END IF END IF END IF * diff --git a/SRC/dlaqr5.f b/SRC/dlaqr5.f index db1fbf87d..12e7db637 100644 --- a/SRC/dlaqr5.f +++ b/SRC/dlaqr5.f @@ -483,36 +483,40 @@ * . is zero (as done here) is traditional but probably * . unnecessary. ==== * - IF( K.GE.KTOP .AND. H( K+1, K ).NE.ZERO ) THEN - TST1 = ABS( H( K, K ) ) + ABS( H( K+1, K+1 ) ) - IF( TST1.EQ.ZERO ) THEN - IF( K.GE.KTOP+1 ) - $ TST1 = TST1 + ABS( H( K, K-1 ) ) - IF( K.GE.KTOP+2 ) - $ TST1 = TST1 + ABS( H( K, K-2 ) ) - IF( K.GE.KTOP+3 ) - $ TST1 = TST1 + ABS( H( K, K-3 ) ) - IF( K.LE.KBOT-2 ) - $ TST1 = TST1 + ABS( H( K+2, K+1 ) ) - IF( K.LE.KBOT-3 ) - $ TST1 = TST1 + ABS( H( K+3, K+1 ) ) - IF( K.LE.KBOT-4 ) - $ TST1 = TST1 + ABS( H( K+4, K+1 ) ) - END IF - IF( ABS( H( K+1, K ) ).LE.MAX( SMLNUM, ULP*TST1 ) ) - $ THEN - H12 = MAX( ABS( H( K+1, K ) ), ABS( H( K, K+1 ) ) ) - H21 = MIN( ABS( H( K+1, K ) ), ABS( H( K, K+1 ) ) ) - H11 = MAX( ABS( H( K+1, K+1 ) ), - $ ABS( H( K, K )-H( K+1, K+1 ) ) ) - H22 = MIN( ABS( H( K+1, K+1 ) ), - $ ABS( H( K, K )-H( K+1, K+1 ) ) ) - SCL = H11 + H12 - TST2 = H22*( H11 / SCL ) + IF( K.GE.KTOP ) THEN + IF( H( K+1, K ).NE.ZERO ) THEN + TST1 = ABS( H( K, K ) ) + ABS( H( K+1, K+1 ) ) + IF( TST1.EQ.ZERO ) THEN + IF( K.GE.KTOP+1 ) + $ TST1 = TST1 + ABS( H( K, K-1 ) ) + IF( K.GE.KTOP+2 ) + $ TST1 = TST1 + ABS( H( K, K-2 ) ) + IF( K.GE.KTOP+3 ) + $ TST1 = TST1 + ABS( H( K, K-3 ) ) + IF( K.LE.KBOT-2 ) + $ TST1 = TST1 + ABS( H( K+2, K+1 ) ) + IF( K.LE.KBOT-3 ) + $ TST1 = TST1 + ABS( H( K+3, K+1 ) ) + IF( K.LE.KBOT-4 ) + $ TST1 = TST1 + ABS( H( K+4, K+1 ) ) + END IF + IF( ABS( H( K+1, K ) ) + $ .LE.MAX( SMLNUM, ULP*TST1 ) ) THEN + H12 = MAX( ABS( H( K+1, K ) ), + $ ABS( H( K, K+1 ) ) ) + H21 = MIN( ABS( H( K+1, K ) ), + $ ABS( H( K, K+1 ) ) ) + H11 = MAX( ABS( H( K+1, K+1 ) ), + $ ABS( H( K, K )-H( K+1, K+1 ) ) ) + H22 = MIN( ABS( H( K+1, K+1 ) ), + $ ABS( H( K, K )-H( K+1, K+1 ) ) ) + SCL = H11 + H12 + TST2 = H22*( H11 / SCL ) * - IF( TST2.EQ.ZERO .OR. H21*( H12 / SCL ).LE. - $ MAX( SMLNUM, ULP*TST2 ) ) THEN - H( K+1, K ) = ZERO + IF( TST2.EQ.ZERO .OR. H21*( H12 / SCL ).LE. + $ MAX( SMLNUM, ULP*TST2 ) ) THEN + H( K+1, K ) = ZERO + END IF END IF END IF END IF diff --git a/SRC/slaqr5.f b/SRC/slaqr5.f index 042df5431..d60a1d3c0 100644 --- a/SRC/slaqr5.f +++ b/SRC/slaqr5.f @@ -483,36 +483,40 @@ * . is zero (as done here) is traditional but probably * . unnecessary. ==== * - IF( K.GE.KTOP .AND. H( K+1, K ).NE.ZERO ) THEN - TST1 = ABS( H( K, K ) ) + ABS( H( K+1, K+1 ) ) - IF( TST1.EQ.ZERO ) THEN - IF( K.GE.KTOP+1 ) - $ TST1 = TST1 + ABS( H( K, K-1 ) ) - IF( K.GE.KTOP+2 ) - $ TST1 = TST1 + ABS( H( K, K-2 ) ) - IF( K.GE.KTOP+3 ) - $ TST1 = TST1 + ABS( H( K, K-3 ) ) - IF( K.LE.KBOT-2 ) - $ TST1 = TST1 + ABS( H( K+2, K+1 ) ) - IF( K.LE.KBOT-3 ) - $ TST1 = TST1 + ABS( H( K+3, K+1 ) ) - IF( K.LE.KBOT-4 ) - $ TST1 = TST1 + ABS( H( K+4, K+1 ) ) - END IF - IF( ABS( H( K+1, K ) ).LE.MAX( SMLNUM, ULP*TST1 ) ) - $ THEN - H12 = MAX( ABS( H( K+1, K ) ), ABS( H( K, K+1 ) ) ) - H21 = MIN( ABS( H( K+1, K ) ), ABS( H( K, K+1 ) ) ) - H11 = MAX( ABS( H( K+1, K+1 ) ), - $ ABS( H( K, K )-H( K+1, K+1 ) ) ) - H22 = MIN( ABS( H( K+1, K+1 ) ), - $ ABS( H( K, K )-H( K+1, K+1 ) ) ) - SCL = H11 + H12 - TST2 = H22*( H11 / SCL ) + IF( K.GE.KTOP ) THEN + IF( H( K+1, K ).NE.ZERO ) THEN + TST1 = ABS( H( K, K ) ) + ABS( H( K+1, K+1 ) ) + IF( TST1.EQ.ZERO ) THEN + IF( K.GE.KTOP+1 ) + $ TST1 = TST1 + ABS( H( K, K-1 ) ) + IF( K.GE.KTOP+2 ) + $ TST1 = TST1 + ABS( H( K, K-2 ) ) + IF( K.GE.KTOP+3 ) + $ TST1 = TST1 + ABS( H( K, K-3 ) ) + IF( K.LE.KBOT-2 ) + $ TST1 = TST1 + ABS( H( K+2, K+1 ) ) + IF( K.LE.KBOT-3 ) + $ TST1 = TST1 + ABS( H( K+3, K+1 ) ) + IF( K.LE.KBOT-4 ) + $ TST1 = TST1 + ABS( H( K+4, K+1 ) ) + END IF + IF( ABS( H( K+1, K ) ).LE.MAX( SMLNUM, ULP*TST1 ) ) + $ THEN + H12 = MAX( ABS( H( K+1, K ) ), + $ ABS( H( K, K+1 ) ) ) + H21 = MIN( ABS( H( K+1, K ) ), + $ ABS( H( K, K+1 ) ) ) + H11 = MAX( ABS( H( K+1, K+1 ) ), + $ ABS( H( K, K )-H( K+1, K+1 ) ) ) + H22 = MIN( ABS( H( K+1, K+1 ) ), + $ ABS( H( K, K )-H( K+1, K+1 ) ) ) + SCL = H11 + H12 + TST2 = H22*( H11 / SCL ) * - IF( TST2.EQ.ZERO .OR. H21*( H12 / SCL ).LE. - $ MAX( SMLNUM, ULP*TST2 ) ) THEN - H( K+1, K ) = ZERO + IF( TST2.EQ.ZERO .OR. H21*( H12 / SCL ).LE. + $ MAX( SMLNUM, ULP*TST2 ) ) THEN + H( K+1, K ) = ZERO + END IF END IF END IF END IF diff --git a/SRC/zlaqr5.f b/SRC/zlaqr5.f index 3a79a0455..c12f4b780 100644 --- a/SRC/zlaqr5.f +++ b/SRC/zlaqr5.f @@ -462,37 +462,39 @@ * . is zero (as done here) is traditional but probably * . unnecessary. ==== * - IF( H( K+1, K ).NE.ZERO ) THEN - TST1 = CABS1( H( K, K ) ) + CABS1( H( K+1, K+1 ) ) - IF( TST1.EQ.RZERO ) THEN - IF( K.GE.KTOP+1 ) - $ TST1 = TST1 + CABS1( H( K, K-1 ) ) - IF( K.GE.KTOP+2 ) - $ TST1 = TST1 + CABS1( H( K, K-2 ) ) - IF( K.GE.KTOP+3 ) - $ TST1 = TST1 + CABS1( H( K, K-3 ) ) - IF( K.LE.KBOT-2 ) - $ TST1 = TST1 + CABS1( H( K+2, K+1 ) ) - IF( K.LE.KBOT-3 ) - $ TST1 = TST1 + CABS1( H( K+3, K+1 ) ) - IF( K.LE.KBOT-4 ) - $ TST1 = TST1 + CABS1( H( K+4, K+1 ) ) - END IF - IF( CABS1( H( K+1, K ) ).LE.MAX( SMLNUM, ULP*TST1 ) ) - $ THEN - H12 = MAX( CABS1( H( K+1, K ) ), + IF( K.GE.KTOP ) THEN + IF( H( K+1, K ).NE.ZERO ) THEN + TST1 = CABS1( H( K, K ) ) + CABS1( H( K+1, K+1 ) ) + IF( TST1.EQ.RZERO ) THEN + IF( K.GE.KTOP+1 ) + $ TST1 = TST1 + CABS1( H( K, K-1 ) ) + IF( K.GE.KTOP+2 ) + $ TST1 = TST1 + CABS1( H( K, K-2 ) ) + IF( K.GE.KTOP+3 ) + $ TST1 = TST1 + CABS1( H( K, K-3 ) ) + IF( K.LE.KBOT-2 ) + $ TST1 = TST1 + CABS1( H( K+2, K+1 ) ) + IF( K.LE.KBOT-3 ) + $ TST1 = TST1 + CABS1( H( K+3, K+1 ) ) + IF( K.LE.KBOT-4 ) + $ TST1 = TST1 + CABS1( H( K+4, K+1 ) ) + END IF + IF( CABS1( H( K+1, K ) ) + $ .LE.MAX( SMLNUM, ULP*TST1 ) ) THEN + H12 = MAX( CABS1( H( K+1, K ) ), $ CABS1( H( K, K+1 ) ) ) - H21 = MIN( CABS1( H( K+1, K ) ), + H21 = MIN( CABS1( H( K+1, K ) ), $ CABS1( H( K, K+1 ) ) ) - H11 = MAX( CABS1( H( K+1, K+1 ) ), + H11 = MAX( CABS1( H( K+1, K+1 ) ), $ CABS1( H( K, K )-H( K+1, K+1 ) ) ) - H22 = MIN( CABS1( H( K+1, K+1 ) ), + H22 = MIN( CABS1( H( K+1, K+1 ) ), $ CABS1( H( K, K )-H( K+1, K+1 ) ) ) - SCL = H11 + H12 - TST2 = H22*( H11 / SCL ) + SCL = H11 + H12 + TST2 = H22*( H11 / SCL ) * - IF( TST2.EQ.RZERO .OR. H21*( H12 / SCL ).LE. - $ MAX( SMLNUM, ULP*TST2 ) )H( K+1, K ) = ZERO + IF( TST2.EQ.RZERO .OR. H21*( H12 / SCL ).LE. + $ MAX( SMLNUM, ULP*TST2 ) )H( K+1, K ) = ZERO + END IF END IF END IF *