Skip to content

Commit 5a7bde2

Browse files
weslleyspereiraweslleyspereira
committed
Simplifies the bug fix
1 parent 4f466dd commit 5a7bde2

4 files changed

Lines changed: 78 additions & 134 deletions

File tree

‎SRC/ctgsja.f‎

Lines changed: 19 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -398,7 +398,7 @@ SUBROUTINE CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
398398
* .. Parameters ..
399399
INTEGER MAXIT
400400
PARAMETER ( MAXIT = 40 )
401-
REAL ZERO, ONE
401+
REAL ZERO, ONE, HUGENUM
402402
PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
403403
COMPLEX CZERO, CONE
404404
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,
409409
LOGICAL INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV
410410
INTEGER I, J, KCYCLE
411411
REAL A1, A3, B1, B3, CSQ, CSU, CSV, ERROR, GAMMA,
412-
$ RWK, SSMIN, SFMIN, HUGE
412+
$ RWK, SSMIN
413413
COMPLEX A2, B2, SNQ, SNU, SNV
414414
* ..
415415
* .. External Functions ..
416416
LOGICAL LSAME
417-
REAL SLAMCH
418-
EXTERNAL LSAME, SLAMCH
417+
EXTERNAL LSAME
419418
* ..
420419
* .. External Subroutines ..
421420
EXTERNAL CCOPY, CLAGS2, CLAPLL, CLASET, CROT, CSSCAL,
422421
$ SLARTG, XERBLA
423422
* ..
424423
* .. Intrinsic Functions ..
425-
INTRINSIC ABS, CONJG, MAX, MIN, REAL
424+
INTRINSIC ABS, CONJG, MAX, MIN, REAL, HUGE
425+
PARAMETER ( HUGENUM = HUGE(ZERO) )
426426
* ..
427427
* .. Executable Statements ..
428428
*
@@ -466,11 +466,6 @@ SUBROUTINE CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
466466
RETURN
467467
END IF
468468
*
469-
* Safe minimum
470-
*
471-
SFMIN = SLAMCH( 'Safe minimum' )
472-
HUGE = SLAMCH( 'O' )
473-
*
474469
* Initialize U, V and Q, if necessary
475470
*
476471
IF( INITU )
@@ -613,34 +608,25 @@ SUBROUTINE CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
613608
*
614609
A1 = REAL( A( K+I, N-L+I ) )
615610
B1 = REAL( B( I, N-L+I ) )
611+
GAMMA = B1 / A1
616612
*
617-
IF( ABS(A1).GE.SFMIN ) THEN
618-
GAMMA = B1 / A1
619-
*
620-
IF( GAMMA.LE.HUGE ) THEN
621-
*
622-
IF( GAMMA.LT.ZERO ) THEN
623-
CALL CSSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB )
624-
IF( WANTV )
625-
$ CALL CSSCAL( P, -ONE, V( 1, I ), 1 )
626-
END IF
613+
IF( (GAMMA.LE.HUGENUM).AND.(GAMMA.GE.-HUGENUM) ) THEN
627614
*
628-
CALL SLARTG( ABS( GAMMA ), ONE, BETA( K+I ),
629-
$ ALPHA( K+I ), RWK )
615+
IF( GAMMA.LT.ZERO ) THEN
616+
CALL CSSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB )
617+
IF( WANTV )
618+
$ CALL CSSCAL( P, -ONE, V( 1, I ), 1 )
619+
END IF
630620
*
631-
IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN
632-
CALL CSSCAL( L-I+1, ONE / ALPHA( K+I ),
633-
$ A( K+I, N-L+I ), LDA )
634-
ELSE
635-
CALL CSSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ),
636-
$ LDB )
637-
CALL CCOPY( L-I+1, B( I, N-L+I ), LDB,
638-
$ A( K+I, N-L+I ), LDA )
639-
END IF
621+
CALL SLARTG( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ),
622+
$ RWK )
640623
*
624+
IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN
625+
CALL CSSCAL( L-I+1, ONE / ALPHA( K+I ), A( K+I, N-L+I ),
626+
$ LDA )
641627
ELSE
642-
ALPHA( K+I ) = ZERO
643-
BETA( K+I ) = ONE
628+
CALL CSSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ),
629+
$ LDB )
644630
CALL CCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ),
645631
$ LDA )
646632
END IF

‎SRC/dtgsja.f‎

Lines changed: 20 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -397,27 +397,27 @@ SUBROUTINE DTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
397397
* .. Parameters ..
398398
INTEGER MAXIT
399399
PARAMETER ( MAXIT = 40 )
400-
DOUBLE PRECISION ZERO, ONE
400+
DOUBLE PRECISION ZERO, ONE, HUGENUM
401401
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
402402
* ..
403403
* .. Local Scalars ..
404404
*
405405
LOGICAL INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV
406406
INTEGER I, J, KCYCLE
407407
DOUBLE PRECISION A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, ERROR,
408-
$ GAMMA, RWK, SNQ, SNU, SNV, SSMIN, SFMIN, HUGE
408+
$ GAMMA, RWK, SNQ, SNU, SNV, SSMIN
409409
* ..
410410
* .. External Functions ..
411411
LOGICAL LSAME
412412
EXTERNAL LSAME
413413
* ..
414414
* .. External Subroutines ..
415-
DOUBLE PRECISION DLAMCH
416415
EXTERNAL DCOPY, DLAGS2, DLAPLL, DLARTG, DLASET, DROT,
417-
$ DSCAL, XERBLA, DLAMCH
416+
$ DSCAL, XERBLA
418417
* ..
419418
* .. Intrinsic Functions ..
420-
INTRINSIC ABS, MAX, MIN
419+
INTRINSIC ABS, MAX, MIN, HUGE
420+
PARAMETER ( HUGENUM = HUGE(ZERO) )
421421
* ..
422422
* .. Executable Statements ..
423423
*
@@ -461,11 +461,6 @@ SUBROUTINE DTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
461461
RETURN
462462
END IF
463463
*
464-
* Safe minimum
465-
*
466-
SFMIN = DLAMCH( 'Safe minimum' )
467-
HUGE = DLAMCH( 'O' )
468-
*
469464
* Initialize U, V and Q, if necessary
470465
*
471466
IF( INITU )
@@ -599,36 +594,27 @@ SUBROUTINE DTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
599594
*
600595
A1 = A( K+I, N-L+I )
601596
B1 = B( I, N-L+I )
597+
GAMMA = B1 / A1
602598
*
603-
IF( ABS(A1).GE.SFMIN ) THEN
604-
GAMMA = B1 / A1
605-
*
606-
IF( GAMMA.LE.HUGE ) THEN
599+
IF( (GAMMA.LE.HUGENUM).AND.(GAMMA.GE.-HUGENUM) ) THEN
607600
*
608-
* change sign if necessary
609-
*
610-
IF( GAMMA.LT.ZERO ) THEN
611-
CALL DSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB )
612-
IF( WANTV )
613-
$ CALL DSCAL( P, -ONE, V( 1, I ), 1 )
614-
END IF
601+
* change sign if necessary
615602
*
616-
CALL DLARTG( ABS( GAMMA ), ONE, BETA( K+I ),
617-
$ ALPHA( K+I ), RWK )
603+
IF( GAMMA.LT.ZERO ) THEN
604+
CALL DSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB )
605+
IF( WANTV )
606+
$ CALL DSCAL( P, -ONE, V( 1, I ), 1 )
607+
END IF
618608
*
619-
IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN
620-
CALL DSCAL( L-I+1, ONE / ALPHA( K+I ),
621-
$ A( K+I, N-L+I ), LDA )
622-
ELSE
623-
CALL DSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ),
624-
$ LDB )
625-
CALL DCOPY( L-I+1, B( I, N-L+I ), LDB,
626-
$ A( K+I, N-L+I ), LDA )
627-
END IF
609+
CALL DLARTG( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ),
610+
$ RWK )
628611
*
612+
IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN
613+
CALL DSCAL( L-I+1, ONE / ALPHA( K+I ), A( K+I, N-L+I ),
614+
$ LDA )
629615
ELSE
630-
ALPHA( K+I ) = ZERO
631-
BETA( K+I ) = ONE
616+
CALL DSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ),
617+
$ LDB )
632618
CALL DCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ),
633619
$ LDA )
634620
END IF

‎SRC/stgsja.f‎

Lines changed: 20 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -397,27 +397,27 @@ SUBROUTINE STGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
397397
* .. Parameters ..
398398
INTEGER MAXIT
399399
PARAMETER ( MAXIT = 40 )
400-
REAL ZERO, ONE
400+
REAL ZERO, ONE, HUGENUM
401401
PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
402402
* ..
403403
* .. Local Scalars ..
404404
*
405405
LOGICAL INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV
406406
INTEGER I, J, KCYCLE
407407
REAL A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, ERROR,
408-
$ GAMMA, RWK, SNQ, SNU, SNV, SSMIN, SFMIN, HUGE
408+
$ GAMMA, RWK, SNQ, SNU, SNV, SSMIN
409409
* ..
410410
* .. External Functions ..
411411
LOGICAL LSAME
412412
EXTERNAL LSAME
413413
* ..
414414
* .. External Subroutines ..
415-
REAL SLAMCH
416415
EXTERNAL SCOPY, SLAGS2, SLAPLL, SLARTG, SLASET, SROT,
417-
$ SSCAL, XERBLA, SLAMCH
416+
$ SSCAL, XERBLA
418417
* ..
419418
* .. Intrinsic Functions ..
420-
INTRINSIC ABS, MAX, MIN
419+
INTRINSIC ABS, MAX, MIN, HUGE
420+
PARAMETER ( HUGENUM = HUGE(ZERO) )
421421
* ..
422422
* .. Executable Statements ..
423423
*
@@ -461,11 +461,6 @@ SUBROUTINE STGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
461461
RETURN
462462
END IF
463463
*
464-
* Safe minimum
465-
*
466-
SFMIN = SLAMCH( 'Safe minimum' )
467-
HUGE = SLAMCH( 'O' )
468-
*
469464
* Initialize U, V and Q, if necessary
470465
*
471466
IF( INITU )
@@ -599,36 +594,27 @@ SUBROUTINE STGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
599594
*
600595
A1 = A( K+I, N-L+I )
601596
B1 = B( I, N-L+I )
597+
GAMMA = B1 / A1
602598
*
603-
IF( ABS(A1).GE.SFMIN ) THEN
604-
GAMMA = B1 / A1
605-
*
606-
IF( GAMMA.LE.HUGE ) THEN
599+
IF( (GAMMA.LE.HUGENUM).AND.(GAMMA.GE.-HUGENUM) ) THEN
607600
*
608-
* change sign if necessary
601+
* change sign if necessary
609602
*
610-
IF( GAMMA.LT.ZERO ) THEN
611-
CALL SSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB )
612-
IF( WANTV )
613-
$ CALL SSCAL( P, -ONE, V( 1, I ), 1 )
614-
END IF
615-
*
616-
CALL SLARTG( ABS( GAMMA ), ONE, BETA( K+I ),
617-
$ ALPHA( K+I ),RWK )
603+
IF( GAMMA.LT.ZERO ) THEN
604+
CALL SSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB )
605+
IF( WANTV )
606+
$ CALL SSCAL( P, -ONE, V( 1, I ), 1 )
607+
END IF
618608
*
619-
IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN
620-
CALL SSCAL( L-I+1, ONE / ALPHA( K+I ),
621-
$ A( K+I, N-L+I ),LDA )
622-
ELSE
623-
CALL SSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ),
624-
$ LDB )
625-
CALL SCOPY( L-I+1, B( I, N-L+I ), LDB,
626-
$ A( K+I, N-L+I ),LDA )
627-
END IF
609+
CALL SLARTG( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ),
610+
$ RWK )
628611
*
612+
IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN
613+
CALL SSCAL( L-I+1, ONE / ALPHA( K+I ), A( K+I, N-L+I ),
614+
$ LDA )
629615
ELSE
630-
ALPHA( K+I ) = ZERO
631-
BETA( K+I ) = ONE
616+
CALL SSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ),
617+
$ LDB )
632618
CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ),
633619
$ LDA )
634620
END IF

‎SRC/ztgsja.f‎

Lines changed: 19 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -398,7 +398,7 @@ SUBROUTINE ZTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
398398
* .. Parameters ..
399399
INTEGER MAXIT
400400
PARAMETER ( MAXIT = 40 )
401-
DOUBLE PRECISION ZERO, ONE
401+
DOUBLE PRECISION ZERO, ONE, HUGENUM
402402
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
403403
COMPLEX*16 CZERO, CONE
404404
PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
@@ -409,20 +409,20 @@ SUBROUTINE ZTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
409409
LOGICAL INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV
410410
INTEGER I, J, KCYCLE
411411
DOUBLE PRECISION A1, A3, B1, B3, CSQ, CSU, CSV, ERROR, GAMMA,
412-
$ RWK, SSMIN, SFMIN, HUGE
412+
$ RWK, SSMIN
413413
COMPLEX*16 A2, B2, SNQ, SNU, SNV
414414
* ..
415415
* .. External Functions ..
416416
LOGICAL LSAME
417417
EXTERNAL LSAME
418418
* ..
419419
* .. External Subroutines ..
420-
DOUBLE PRECISION DLAMCH
421420
EXTERNAL DLARTG, XERBLA, ZCOPY, ZDSCAL, ZLAGS2, ZLAPLL,
422-
$ ZLASET, ZROT, DLAMCH
421+
$ ZLASET, ZROT
423422
* ..
424423
* .. Intrinsic Functions ..
425-
INTRINSIC ABS, DBLE, DCONJG, MAX, MIN
424+
INTRINSIC ABS, DBLE, DCONJG, MAX, MIN, HUGE
425+
PARAMETER ( HUGENUM = HUGE(ZERO) )
426426
* ..
427427
* .. Executable Statements ..
428428
*
@@ -466,11 +466,6 @@ SUBROUTINE ZTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
466466
RETURN
467467
END IF
468468
*
469-
* Safe minimum
470-
*
471-
SFMIN = DLAMCH( 'Safe minimum' )
472-
HUGE = DLAMCH( 'O' )
473-
*
474469
* Initialize U, V and Q, if necessary
475470
*
476471
IF( INITU )
@@ -613,34 +608,25 @@ SUBROUTINE ZTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
613608
*
614609
A1 = DBLE( A( K+I, N-L+I ) )
615610
B1 = DBLE( B( I, N-L+I ) )
611+
GAMMA = B1 / A1
616612
*
617-
IF( ABS(A1).GE.SFMIN ) THEN
618-
GAMMA = B1 / A1
619-
*
620-
IF( GAMMA.LE.HUGE ) THEN
621-
*
622-
IF( GAMMA.LT.ZERO ) THEN
623-
CALL ZDSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB )
624-
IF( WANTV )
625-
$ CALL ZDSCAL( P, -ONE, V( 1, I ), 1 )
626-
END IF
613+
IF( (GAMMA.LE.HUGENUM).AND.(GAMMA.GE.-HUGENUM) ) THEN
627614
*
628-
CALL ZLARTG( ABS( GAMMA ), ONE, BETA( K+I ),
629-
$ ALPHA( K+I ), RWK )
615+
IF( GAMMA.LT.ZERO ) THEN
616+
CALL ZDSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB )
617+
IF( WANTV )
618+
$ CALL ZDSCAL( P, -ONE, V( 1, I ), 1 )
619+
END IF
630620
*
631-
IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN
632-
CALL ZDSCAL( L-I+1, ONE / ALPHA( K+I ),
633-
$ A( K+I, N-L+I ), LDA )
634-
ELSE
635-
CALL ZDSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ),
636-
$ LDB )
637-
CALL ZCOPY( L-I+1, B( I, N-L+I ), LDB,
638-
$ A( K+I, N-L+I ), LDA )
639-
END IF
621+
CALL DLARTG( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ),
622+
$ RWK )
640623
*
624+
IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN
625+
CALL ZDSCAL( L-I+1, ONE / ALPHA( K+I ), A( K+I, N-L+I ),
626+
$ LDA )
641627
ELSE
642-
ALPHA( K+I ) = ZERO
643-
BETA( K+I ) = ONE
628+
CALL ZDSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ),
629+
$ LDB )
644630
CALL ZCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ),
645631
$ LDA )
646632
END IF

0 commit comments

Comments
 (0)