Skip to content

Commit

Permalink
Fix bound checks for gamma in xTGSJA
Browse files Browse the repository at this point in the history
  • Loading branch information
@weslleyspereira
weslleyspereira committed Apr 6, 2021
1 parent 6433162 commit 4f466dd
Show file tree
Hide file tree
Showing 4 changed files with 115 additions and 60 deletions.
36 changes: 23 additions & 13 deletions SRC/ctgsja.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -409,7 +409,7 @@ 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
$ RWK, SSMIN, SFMIN, HUGE
COMPLEX A2, B2, SNQ, SNU, SNV
* ..
* .. External Functions ..
Expand Down Expand Up @@ -469,6 +469,7 @@ SUBROUTINE CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
* Safe minimum
*
SFMIN = SLAMCH( 'Safe minimum' )
HUGE = SLAMCH( 'O' )
*
* Initialize U, V and Q, if necessary
*
Expand Down Expand Up @@ -616,21 +617,30 @@ SUBROUTINE CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
IF( ABS(A1).GE.SFMIN ) THEN
GAMMA = B1 / A1
*
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.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
*
CALL SLARTG( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ),
$ RWK )
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
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
*
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 )
ALPHA( K+I ) = ZERO
BETA( K+I ) = ONE
CALL CCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ),
$ LDA )
END IF
Expand Down
47 changes: 31 additions & 16 deletions SRC/dtgsja.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -405,15 +405,16 @@ 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
$ GAMMA, RWK, SNQ, SNU, SNV, SSMIN, SFMIN, HUGE
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
DOUBLE PRECISION DLAMCH
EXTERNAL DCOPY, DLAGS2, DLAPLL, DLARTG, DLASET, DROT,
$ DSCAL, XERBLA
$ DSCAL, XERBLA, DLAMCH
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN
Expand Down Expand Up @@ -460,6 +461,11 @@ 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 )
Expand Down Expand Up @@ -594,26 +600,35 @@ 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 )
*
IF( A1.NE.ZERO ) THEN
IF( ABS(A1).GE.SFMIN ) THEN
GAMMA = B1 / A1
*
* change sign if necessary
IF( GAMMA.LE.HUGE ) THEN
*
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
*
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
*
CALL DLARTG( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ),
$ RWK )
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
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
*
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 )
ALPHA( K+I ) = ZERO
BETA( K+I ) = ONE
CALL DCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ),
$ LDA )
END IF
Expand Down
47 changes: 31 additions & 16 deletions SRC/stgsja.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -405,15 +405,16 @@ 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
$ GAMMA, RWK, SNQ, SNU, SNV, SSMIN, SFMIN, HUGE
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
REAL SLAMCH
EXTERNAL SCOPY, SLAGS2, SLAPLL, SLARTG, SLASET, SROT,
$ SSCAL, XERBLA
$ SSCAL, XERBLA, SLAMCH
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN
Expand Down Expand Up @@ -460,6 +461,11 @@ 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 )
Expand Down Expand Up @@ -594,26 +600,35 @@ 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 )
*
IF( A1.NE.ZERO ) THEN
IF( ABS(A1).GE.SFMIN ) THEN
GAMMA = B1 / A1
*
* change sign if necessary
IF( GAMMA.LE.HUGE ) THEN
*
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
* change sign if necessary
*
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 ),
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
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
*
ELSE
ALPHA( K+I ) = ZERO
BETA( K+I ) = ONE
CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ),
$ LDA )
END IF
Expand Down
45 changes: 30 additions & 15 deletions SRC/ztgsja.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -409,16 +409,17 @@ 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
$ RWK, SSMIN, SFMIN, HUGE
COMPLEX*16 A2, B2, SNQ, SNU, SNV
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
DOUBLE PRECISION DLAMCH
EXTERNAL DLARTG, XERBLA, ZCOPY, ZDSCAL, ZLAGS2, ZLAPLL,
$ ZLASET, ZROT
$ ZLASET, ZROT, DLAMCH
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, DCONJG, MAX, MIN
Expand Down Expand Up @@ -465,6 +466,11 @@ 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 )
Expand Down Expand Up @@ -608,24 +614,33 @@ 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 ) )
*
IF( A1.NE.ZERO ) THEN
IF( ABS(A1).GE.SFMIN ) THEN
GAMMA = B1 / A1
*
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.HUGE ) THEN
*
CALL DLARTG( 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
*
CALL ZLARTG( 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
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
*
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 )
ALPHA( K+I ) = ZERO
BETA( K+I ) = ONE
CALL ZCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ),
$ LDA )
END IF
Expand Down

0 comments on commit 4f466dd

Please sign in to comment.