diff --git a/SRC/ctgsja.f b/SRC/ctgsja.f index 4976489b23..22725e86ce 100644 --- a/SRC/ctgsja.f +++ b/SRC/ctgsja.f @@ -398,7 +398,7 @@ SUBROUTINE CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, * .. Parameters .. INTEGER MAXIT PARAMETER ( MAXIT = 40 ) - REAL ZERO, ONE + REAL ZERO, ONE, HUGENUM PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) COMPLEX CZERO, CONE PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), @@ -409,20 +409,20 @@ SUBROUTINE CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LOGICAL INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV INTEGER I, J, KCYCLE REAL A1, A3, B1, B3, CSQ, CSU, CSV, ERROR, GAMMA, - $ RWK, SSMIN, SFMIN, HUGE + $ RWK, SSMIN COMPLEX A2, B2, SNQ, SNU, SNV * .. * .. External Functions .. LOGICAL LSAME - REAL SLAMCH - EXTERNAL LSAME, SLAMCH + EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL CCOPY, CLAGS2, CLAPLL, CLASET, CROT, CSSCAL, $ SLARTG, XERBLA * .. * .. Intrinsic Functions .. - INTRINSIC ABS, CONJG, MAX, MIN, REAL + INTRINSIC ABS, CONJG, MAX, MIN, REAL, HUGE + PARAMETER ( HUGENUM = HUGE(ZERO) ) * .. * .. Executable Statements .. * @@ -466,11 +466,6 @@ SUBROUTINE CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, RETURN END IF * -* Safe minimum -* - SFMIN = SLAMCH( 'Safe minimum' ) - HUGE = SLAMCH( 'O' ) -* * Initialize U, V and Q, if necessary * IF( INITU ) @@ -613,34 +608,25 @@ SUBROUTINE CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, * A1 = REAL( A( K+I, N-L+I ) ) B1 = REAL( B( I, N-L+I ) ) + GAMMA = B1 / A1 * - IF( ABS(A1).GE.SFMIN ) THEN - GAMMA = B1 / A1 -* - IF( GAMMA.LE.HUGE ) THEN -* - IF( GAMMA.LT.ZERO ) THEN - CALL CSSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB ) - IF( WANTV ) - $ CALL CSSCAL( P, -ONE, V( 1, I ), 1 ) - END IF + IF( (GAMMA.LE.HUGENUM).AND.(GAMMA.GE.-HUGENUM) ) THEN * - CALL SLARTG( ABS( GAMMA ), ONE, BETA( K+I ), - $ ALPHA( K+I ), RWK ) + IF( GAMMA.LT.ZERO ) THEN + CALL CSSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB ) + IF( WANTV ) + $ CALL CSSCAL( P, -ONE, V( 1, I ), 1 ) + END IF * - IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN - CALL CSSCAL( L-I+1, ONE / ALPHA( K+I ), - $ A( K+I, N-L+I ), LDA ) - ELSE - CALL CSSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ), - $ LDB ) - CALL CCOPY( L-I+1, B( I, N-L+I ), LDB, - $ A( K+I, N-L+I ), LDA ) - END IF + CALL SLARTG( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ), + $ RWK ) * + IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN + CALL CSSCAL( L-I+1, ONE / ALPHA( K+I ), A( K+I, N-L+I ), + $ LDA ) ELSE - ALPHA( K+I ) = ZERO - BETA( K+I ) = ONE + CALL CSSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ), + $ LDB ) CALL CCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ), $ LDA ) END IF diff --git a/SRC/dtgsja.f b/SRC/dtgsja.f index d6d38a1da1..37668c7391 100644 --- a/SRC/dtgsja.f +++ b/SRC/dtgsja.f @@ -397,7 +397,7 @@ SUBROUTINE DTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, * .. Parameters .. INTEGER MAXIT PARAMETER ( MAXIT = 40 ) - DOUBLE PRECISION ZERO, ONE + DOUBLE PRECISION ZERO, ONE, HUGENUM PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) * .. * .. Local Scalars .. @@ -405,19 +405,19 @@ SUBROUTINE DTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LOGICAL INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV INTEGER I, J, KCYCLE DOUBLE PRECISION A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, ERROR, - $ GAMMA, RWK, SNQ, SNU, SNV, SSMIN, SFMIN, HUGE + $ GAMMA, RWK, SNQ, SNU, SNV, SSMIN * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. - DOUBLE PRECISION DLAMCH EXTERNAL DCOPY, DLAGS2, DLAPLL, DLARTG, DLASET, DROT, - $ DSCAL, XERBLA, DLAMCH + $ DSCAL, XERBLA * .. * .. Intrinsic Functions .. - INTRINSIC ABS, MAX, MIN + INTRINSIC ABS, MAX, MIN, HUGE + PARAMETER ( HUGENUM = HUGE(ZERO) ) * .. * .. Executable Statements .. * @@ -461,11 +461,6 @@ SUBROUTINE DTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, RETURN END IF * -* Safe minimum -* - SFMIN = DLAMCH( 'Safe minimum' ) - HUGE = DLAMCH( 'O' ) -* * Initialize U, V and Q, if necessary * IF( INITU ) @@ -599,36 +594,27 @@ SUBROUTINE DTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, * A1 = A( K+I, N-L+I ) B1 = B( I, N-L+I ) + GAMMA = B1 / A1 * - IF( ABS(A1).GE.SFMIN ) THEN - GAMMA = B1 / A1 -* - IF( GAMMA.LE.HUGE ) THEN + IF( (GAMMA.LE.HUGENUM).AND.(GAMMA.GE.-HUGENUM) ) THEN * -* change sign if necessary -* - IF( GAMMA.LT.ZERO ) THEN - CALL DSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB ) - IF( WANTV ) - $ CALL DSCAL( P, -ONE, V( 1, I ), 1 ) - END IF +* change sign if necessary * - CALL DLARTG( ABS( GAMMA ), ONE, BETA( K+I ), - $ ALPHA( K+I ), RWK ) + IF( GAMMA.LT.ZERO ) THEN + CALL DSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB ) + IF( WANTV ) + $ CALL DSCAL( P, -ONE, V( 1, I ), 1 ) + END IF * - IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN - CALL DSCAL( L-I+1, ONE / ALPHA( K+I ), - $ A( K+I, N-L+I ), LDA ) - ELSE - CALL DSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ), - $ LDB ) - CALL DCOPY( L-I+1, B( I, N-L+I ), LDB, - $ A( K+I, N-L+I ), LDA ) - END IF + CALL DLARTG( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ), + $ RWK ) * + IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN + CALL DSCAL( L-I+1, ONE / ALPHA( K+I ), A( K+I, N-L+I ), + $ LDA ) ELSE - ALPHA( K+I ) = ZERO - BETA( K+I ) = ONE + CALL DSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ), + $ LDB ) CALL DCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ), $ LDA ) END IF diff --git a/SRC/stgsja.f b/SRC/stgsja.f index 28f5b3da3e..5a81b7e4fd 100644 --- a/SRC/stgsja.f +++ b/SRC/stgsja.f @@ -397,7 +397,7 @@ SUBROUTINE STGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, * .. Parameters .. INTEGER MAXIT PARAMETER ( MAXIT = 40 ) - REAL ZERO, ONE + REAL ZERO, ONE, HUGENUM PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) * .. * .. Local Scalars .. @@ -405,19 +405,19 @@ SUBROUTINE STGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LOGICAL INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV INTEGER I, J, KCYCLE REAL A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, ERROR, - $ GAMMA, RWK, SNQ, SNU, SNV, SSMIN, SFMIN, HUGE + $ GAMMA, RWK, SNQ, SNU, SNV, SSMIN * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. - REAL SLAMCH EXTERNAL SCOPY, SLAGS2, SLAPLL, SLARTG, SLASET, SROT, - $ SSCAL, XERBLA, SLAMCH + $ SSCAL, XERBLA * .. * .. Intrinsic Functions .. - INTRINSIC ABS, MAX, MIN + INTRINSIC ABS, MAX, MIN, HUGE + PARAMETER ( HUGENUM = HUGE(ZERO) ) * .. * .. Executable Statements .. * @@ -461,11 +461,6 @@ SUBROUTINE STGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, RETURN END IF * -* Safe minimum -* - SFMIN = SLAMCH( 'Safe minimum' ) - HUGE = SLAMCH( 'O' ) -* * Initialize U, V and Q, if necessary * IF( INITU ) @@ -599,36 +594,27 @@ SUBROUTINE STGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, * A1 = A( K+I, N-L+I ) B1 = B( I, N-L+I ) + GAMMA = B1 / A1 * - IF( ABS(A1).GE.SFMIN ) THEN - GAMMA = B1 / A1 -* - IF( GAMMA.LE.HUGE ) THEN + IF( (GAMMA.LE.HUGENUM).AND.(GAMMA.GE.-HUGENUM) ) THEN * -* change sign if necessary +* change sign if necessary * - IF( GAMMA.LT.ZERO ) THEN - CALL SSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB ) - IF( WANTV ) - $ CALL SSCAL( P, -ONE, V( 1, I ), 1 ) - END IF -* - CALL SLARTG( ABS( GAMMA ), ONE, BETA( K+I ), - $ ALPHA( K+I ),RWK ) + IF( GAMMA.LT.ZERO ) THEN + CALL SSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB ) + IF( WANTV ) + $ CALL SSCAL( P, -ONE, V( 1, I ), 1 ) + END IF * - IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN - CALL SSCAL( L-I+1, ONE / ALPHA( K+I ), - $ A( K+I, N-L+I ),LDA ) - ELSE - CALL SSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ), - $ LDB ) - CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, - $ A( K+I, N-L+I ),LDA ) - END IF + CALL SLARTG( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ), + $ RWK ) * + IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN + CALL SSCAL( L-I+1, ONE / ALPHA( K+I ), A( K+I, N-L+I ), + $ LDA ) ELSE - ALPHA( K+I ) = ZERO - BETA( K+I ) = ONE + CALL SSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ), + $ LDB ) CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ), $ LDA ) END IF diff --git a/SRC/ztgsja.f b/SRC/ztgsja.f index e248c4c617..c1c5abdb18 100644 --- a/SRC/ztgsja.f +++ b/SRC/ztgsja.f @@ -398,7 +398,7 @@ SUBROUTINE ZTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, * .. Parameters .. INTEGER MAXIT PARAMETER ( MAXIT = 40 ) - DOUBLE PRECISION ZERO, ONE + DOUBLE PRECISION ZERO, ONE, HUGENUM PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) COMPLEX*16 CZERO, CONE PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), @@ -409,7 +409,7 @@ SUBROUTINE ZTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LOGICAL INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV INTEGER I, J, KCYCLE DOUBLE PRECISION A1, A3, B1, B3, CSQ, CSU, CSV, ERROR, GAMMA, - $ RWK, SSMIN, SFMIN, HUGE + $ RWK, SSMIN COMPLEX*16 A2, B2, SNQ, SNU, SNV * .. * .. External Functions .. @@ -417,12 +417,12 @@ SUBROUTINE ZTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, EXTERNAL LSAME * .. * .. External Subroutines .. - DOUBLE PRECISION DLAMCH EXTERNAL DLARTG, XERBLA, ZCOPY, ZDSCAL, ZLAGS2, ZLAPLL, - $ ZLASET, ZROT, DLAMCH + $ ZLASET, ZROT * .. * .. Intrinsic Functions .. - INTRINSIC ABS, DBLE, DCONJG, MAX, MIN + INTRINSIC ABS, DBLE, DCONJG, MAX, MIN, HUGE + PARAMETER ( HUGENUM = HUGE(ZERO) ) * .. * .. Executable Statements .. * @@ -466,11 +466,6 @@ SUBROUTINE ZTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, RETURN END IF * -* Safe minimum -* - SFMIN = DLAMCH( 'Safe minimum' ) - HUGE = DLAMCH( 'O' ) -* * Initialize U, V and Q, if necessary * IF( INITU ) @@ -613,34 +608,25 @@ SUBROUTINE ZTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, * A1 = DBLE( A( K+I, N-L+I ) ) B1 = DBLE( B( I, N-L+I ) ) + GAMMA = B1 / A1 * - IF( ABS(A1).GE.SFMIN ) THEN - GAMMA = B1 / A1 -* - IF( GAMMA.LE.HUGE ) THEN -* - IF( GAMMA.LT.ZERO ) THEN - CALL ZDSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB ) - IF( WANTV ) - $ CALL ZDSCAL( P, -ONE, V( 1, I ), 1 ) - END IF + IF( (GAMMA.LE.HUGENUM).AND.(GAMMA.GE.-HUGENUM) ) THEN * - CALL ZLARTG( ABS( GAMMA ), ONE, BETA( K+I ), - $ ALPHA( K+I ), RWK ) + IF( GAMMA.LT.ZERO ) THEN + CALL ZDSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB ) + IF( WANTV ) + $ CALL ZDSCAL( P, -ONE, V( 1, I ), 1 ) + END IF * - IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN - CALL ZDSCAL( L-I+1, ONE / ALPHA( K+I ), - $ A( K+I, N-L+I ), LDA ) - ELSE - CALL ZDSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ), - $ LDB ) - CALL ZCOPY( L-I+1, B( I, N-L+I ), LDB, - $ A( K+I, N-L+I ), LDA ) - END IF + CALL DLARTG( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ), + $ RWK ) * + IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN + CALL ZDSCAL( L-I+1, ONE / ALPHA( K+I ), A( K+I, N-L+I ), + $ LDA ) ELSE - ALPHA( K+I ) = ZERO - BETA( K+I ) = ONE + CALL ZDSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ), + $ LDB ) CALL ZCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ), $ LDA ) END IF