Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Refactor xGEBAL #808

Merged
merged 2 commits into from
Mar 28, 2023
Merged

Refactor xGEBAL #808

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
21231d8
Refactored GEBAL family
eprovst Mar 24, 2023
9c489fd
Spend less time checking for NaN in xGEBAL
eprovst Mar 27, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
312 changes: 167 additions & 145 deletions SRC/cgebal.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,7 @@
*> \verbatim
*> ILO is INTEGER
*> \endverbatim
*>
*> \param[out] IHI
*> \verbatim
*> IHI is INTEGER
Expand Down Expand Up @@ -154,6 +155,9 @@
*>
*> Modified by Tzu-Yi Chen, Computer Science Division, University of
*> California at Berkeley, USA
*>
*> Refactored by Evert Provoost, Department of Computer Science,
*> KU Leuven, Belgium
*> \endverbatim
*>
* =====================================================================
Expand Down Expand Up @@ -183,8 +187,8 @@ SUBROUTINE CGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
PARAMETER ( FACTOR = 0.95E+0 )
* ..
* .. Local Scalars ..
LOGICAL NOCONV
INTEGER I, ICA, IEXC, IRA, J, K, L, M
LOGICAL NOCONV, CANSWAP
INTEGER I, ICA, IRA, J, K, L
REAL C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
$ SFMIN2
* ..
Expand All @@ -195,10 +199,10 @@ SUBROUTINE CGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
EXTERNAL SISNAN, LSAME, ICAMAX, SLAMCH, SCNRM2
* ..
* .. External Subroutines ..
EXTERNAL CSSCAL, CSWAP, XERBLA
EXTERNAL XERBLA, CSSCAL, CSWAP
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, AIMAG, MAX, MIN, REAL
INTRINSIC ABS, REAL, AIMAG, MAX, MIN
*
* Test the input parameters
*
Expand All @@ -216,176 +220,194 @@ SUBROUTINE CGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
RETURN
END IF
*
K = 1
L = N
* Quick returns.
*
IF( N.EQ.0 )
$ GO TO 210
IF( N.EQ.0 ) THEN
ILO = 1
IHI = 0
RETURN
END IF
*
IF( LSAME( JOB, 'N' ) ) THEN
DO 10 I = 1, N
DO I = 1, N
SCALE( I ) = ONE
10 CONTINUE
GO TO 210
END DO
ILO = 1
IHI = N
RETURN
END IF
*
IF( LSAME( JOB, 'S' ) )
$ GO TO 120
*
* Permutation to isolate eigenvalues if possible
*
GO TO 50
*
* Row and column exchange.
*
20 CONTINUE
SCALE( M ) = J
IF( J.EQ.M )
$ GO TO 30
*
CALL CSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
CALL CSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
*
30 CONTINUE
GO TO ( 40, 80 )IEXC
*
* Search for rows isolating an eigenvalue and push them down.
*
40 CONTINUE
IF( L.EQ.1 )
$ GO TO 210
L = L - 1
*
50 CONTINUE
DO 70 J = L, 1, -1
* Permutation to isolate eigenvalues if possible.
*
DO 60 I = 1, L
IF( I.EQ.J )
$ GO TO 60
IF( REAL( A( J, I ) ).NE.ZERO .OR. AIMAG( A( J, I ) ).NE.
$ ZERO )GO TO 70
60 CONTINUE
*
M = L
IEXC = 1
GO TO 20
70 CONTINUE
*
GO TO 90
K = 1
L = N
*
* Search for columns isolating an eigenvalue and push them left.
IF( .NOT.LSAME( JOB, 'S' ) ) THEN
*
80 CONTINUE
K = K + 1
* Row and column exchange.
*
90 CONTINUE
DO 110 J = K, L
NOCONV = .TRUE.
DO WHILE( NOCONV )
*
* Search for rows isolating an eigenvalue and push them down.
*
NOCONV = .FALSE.
DO I = L, 1, -1
CANSWAP = .TRUE.
DO J = 1, L
IF( I.NE.J .AND. ( REAL( A( I, J ) ).NE.ZERO .OR.
$ AIMAG( A( I, J ) ).NE.ZERO ) ) THEN
CANSWAP = .FALSE.
EXIT
END IF
END DO
*
IF( CANSWAP ) THEN
SCALE( L ) = I
IF( I.NE.L ) THEN
CALL CSWAP( L, A( 1, I ), 1, A( 1, L ), 1 )
CALL CSWAP( N-K+1, A( I, K ), LDA, A( L, K ), LDA )
END IF
NOCONV = .TRUE.
*
IF( L.EQ.1 ) THEN
ILO = 1
IHI = 1
RETURN
END IF
*
L = L - 1
END IF
END DO
*
END DO

NOCONV = .TRUE.
DO WHILE( NOCONV )
*
* Search for columns isolating an eigenvalue and push them left.
*
NOCONV = .FALSE.
DO J = K, L
CANSWAP = .TRUE.
DO I = K, L
IF( I.NE.J .AND. ( REAL( A( I, J ) ).NE.ZERO .OR.
$ AIMAG( A( I, J ) ).NE.ZERO ) ) THEN
CANSWAP = .FALSE.
EXIT
END IF
END DO
*
IF( CANSWAP ) THEN
SCALE( K ) = J
IF( J.NE.K ) THEN
CALL CSWAP( L, A( 1, J ), 1, A( 1, K ), 1 )
CALL CSWAP( N-K+1, A( J, K ), LDA, A( K, K ), LDA )
END IF
NOCONV = .TRUE.
*
K = K + 1
END IF
END DO
*
END DO
*
DO 100 I = K, L
IF( I.EQ.J )
$ GO TO 100
IF( REAL( A( I, J ) ).NE.ZERO .OR. AIMAG( A( I, J ) ).NE.
$ ZERO )GO TO 110
100 CONTINUE
END IF
*
M = K
IEXC = 2
GO TO 20
110 CONTINUE
* Initialize SCALE for non-permuted submatrix.
*
120 CONTINUE
DO 130 I = K, L
DO I = K, L
SCALE( I ) = ONE
130 CONTINUE
END DO
*
IF( LSAME( JOB, 'P' ) )
$ GO TO 210
* If we only had to permute, we are done.
*
IF( LSAME( JOB, 'P' ) ) THEN
ILO = K
IHI = L
RETURN
END IF
*
* Balance the submatrix in rows K to L.
*
* Iterative loop for norm reduction
* Iterative loop for norm reduction.
*
SFMIN1 = SLAMCH( 'S' ) / SLAMCH( 'P' )
SFMAX1 = ONE / SFMIN1
SFMIN2 = SFMIN1*SCLFAC
SFMAX2 = ONE / SFMIN2
140 CONTINUE
NOCONV = .FALSE.
*
DO 200 I = K, L
*
C = SCNRM2( L-K+1, A( K, I ), 1 )
R = SCNRM2( L-K+1, A( I , K ), LDA )
ICA = ICAMAX( L, A( 1, I ), 1 )
CA = ABS( A( ICA, I ) )
IRA = ICAMAX( N-K+1, A( I, K ), LDA )
RA = ABS( A( I, IRA+K-1 ) )
*
* Guard against zero C or R due to underflow.
*
IF( C.EQ.ZERO .OR. R.EQ.ZERO )
$ GO TO 200
G = R / SCLFAC
F = ONE
S = C + R
160 CONTINUE
IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
$ MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
IF( SISNAN( C+F+CA+R+G+RA ) ) THEN
*
* Exit if NaN to avoid infinite loop
NOCONV = .TRUE.
DO WHILE( NOCONV )
NOCONV = .FALSE.
*
INFO = -3
CALL XERBLA( 'CGEBAL', -INFO )
RETURN
END IF
F = F*SCLFAC
C = C*SCLFAC
CA = CA*SCLFAC
R = R / SCLFAC
G = G / SCLFAC
RA = RA / SCLFAC
GO TO 160
*
170 CONTINUE
G = C / SCLFAC
180 CONTINUE
IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
$ MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
F = F / SCLFAC
C = C / SCLFAC
G = G / SCLFAC
CA = CA / SCLFAC
R = R*SCLFAC
RA = RA*SCLFAC
GO TO 180
*
* Now balance.
*
190 CONTINUE
IF( ( C+R ).GE.FACTOR*S )
$ GO TO 200
IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
IF( F*SCALE( I ).LE.SFMIN1 )
$ GO TO 200
END IF
IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
IF( SCALE( I ).GE.SFMAX1 / F )
$ GO TO 200
END IF
G = ONE / F
SCALE( I ) = SCALE( I )*F
NOCONV = .TRUE.
DO I = K, L
*
CALL CSSCAL( N-K+1, G, A( I, K ), LDA )
CALL CSSCAL( L, F, A( 1, I ), 1 )
C = SCNRM2( L-K+1, A( K, I ), 1 )
R = SCNRM2( L-K+1, A( I, K ), LDA )
ICA = ICAMAX( L, A( 1, I ), 1 )
CA = ABS( A( ICA, I ) )
IRA = ICAMAX( N-K+1, A( I, K ), LDA )
RA = ABS( A( I, IRA+K-1 ) )
*
200 CONTINUE
* Guard against zero C or R due to underflow.
*
IF( C.EQ.ZERO .OR. R.EQ.ZERO ) CYCLE
*
* Exit if NaN to avoid infinite loop
*
IF( NOCONV )
$ GO TO 140
IF( SISNAN( C+CA+R+RA ) ) THEN
INFO = -3
CALL XERBLA( 'CGEBAL', -INFO )
RETURN
END IF
*
G = R / SCLFAC
F = ONE
S = C + R
*
DO WHILE( C.LT.G .AND. MAX( F, C, CA ).LT.SFMAX2 .AND.
$ MIN( R, G, RA ).GT.SFMIN2 )
F = F*SCLFAC
C = C*SCLFAC
CA = CA*SCLFAC
R = R / SCLFAC
G = G / SCLFAC
RA = RA / SCLFAC
END DO
*
G = C / SCLFAC
*
DO WHILE( G.GE.R .AND. MAX( R, RA ).LT.SFMAX2 .AND.
$ MIN( F, C, G, CA ).GT.SFMIN2 )
F = F / SCLFAC
C = C / SCLFAC
G = G / SCLFAC
CA = CA / SCLFAC
R = R*SCLFAC
RA = RA*SCLFAC
END DO
*
* Now balance.
*
IF( ( C+R ).GE.FACTOR*S ) CYCLE
IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
IF( F*SCALE( I ).LE.SFMIN1 ) CYCLE
END IF
IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
IF( SCALE( I ).GE.SFMAX1 / F ) CYCLE
END IF
G = ONE / F
SCALE( I ) = SCALE( I )*F
NOCONV = .TRUE.
*
CALL CSSCAL( N-K+1, G, A( I, K ), LDA )
CALL CSSCAL( L, F, A( 1, I ), 1 )
*
END DO
*
END DO
*
210 CONTINUE
ILO = K
IHI = L
*
Expand Down
Loading