Refactor ?GEBAL (Reference-LAPACK PR 808)
This commit is contained in:
Martin Kroeker
GitHub
parent
5f1fb27c40
commit
d175b8f56f
|
|
@ -85,6 +85,7 @@
|
||||||
*> \verbatim
|
*> \verbatim
|
||||||
*> ILO is INTEGER
|
*> ILO is INTEGER
|
||||||
*> \endverbatim
|
*> \endverbatim
|
||||||
|
*>
|
||||||
*> \param[out] IHI
|
*> \param[out] IHI
|
||||||
*> \verbatim
|
*> \verbatim
|
||||||
*> IHI is INTEGER
|
*> IHI is INTEGER
|
||||||
|
|
@ -154,6 +155,9 @@
|
||||||
*>
|
*>
|
||||||
*> Modified by Tzu-Yi Chen, Computer Science Division, University of
|
*> Modified by Tzu-Yi Chen, Computer Science Division, University of
|
||||||
*> California at Berkeley, USA
|
*> California at Berkeley, USA
|
||||||
|
*>
|
||||||
|
*> Refactored by Evert Provoost, Department of Computer Science,
|
||||||
|
*> KU Leuven, Belgium
|
||||||
*> \endverbatim
|
*> \endverbatim
|
||||||
*>
|
*>
|
||||||
* =====================================================================
|
* =====================================================================
|
||||||
|
|
@ -183,8 +187,8 @@
|
||||||
PARAMETER ( FACTOR = 0.95E+0 )
|
PARAMETER ( FACTOR = 0.95E+0 )
|
||||||
* ..
|
* ..
|
||||||
* .. Local Scalars ..
|
* .. Local Scalars ..
|
||||||
LOGICAL NOCONV
|
LOGICAL NOCONV, CANSWAP
|
||||||
INTEGER I, ICA, IEXC, IRA, J, K, L, M
|
INTEGER I, ICA, IRA, J, K, L
|
||||||
REAL C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
|
REAL C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
|
||||||
$ SFMIN2
|
$ SFMIN2
|
||||||
* ..
|
* ..
|
||||||
|
|
@ -195,10 +199,10 @@
|
||||||
EXTERNAL SISNAN, LSAME, ICAMAX, SLAMCH, SCNRM2
|
EXTERNAL SISNAN, LSAME, ICAMAX, SLAMCH, SCNRM2
|
||||||
* ..
|
* ..
|
||||||
* .. External Subroutines ..
|
* .. External Subroutines ..
|
||||||
EXTERNAL CSSCAL, CSWAP, XERBLA
|
EXTERNAL XERBLA, CSSCAL, CSWAP
|
||||||
* ..
|
* ..
|
||||||
* .. Intrinsic Functions ..
|
* .. Intrinsic Functions ..
|
||||||
INTRINSIC ABS, AIMAG, MAX, MIN, REAL
|
INTRINSIC ABS, REAL, AIMAG, MAX, MIN
|
||||||
*
|
*
|
||||||
* Test the input parameters
|
* Test the input parameters
|
||||||
*
|
*
|
||||||
|
|
@ -215,177 +219,195 @@
|
||||||
CALL XERBLA( 'CGEBAL', -INFO )
|
CALL XERBLA( 'CGEBAL', -INFO )
|
||||||
RETURN
|
RETURN
|
||||||
END IF
|
END IF
|
||||||
|
*
|
||||||
|
* Quick returns.
|
||||||
|
*
|
||||||
|
IF( N.EQ.0 ) THEN
|
||||||
|
ILO = 1
|
||||||
|
IHI = 0
|
||||||
|
RETURN
|
||||||
|
END IF
|
||||||
|
*
|
||||||
|
IF( LSAME( JOB, 'N' ) ) THEN
|
||||||
|
DO I = 1, N
|
||||||
|
SCALE( I ) = ONE
|
||||||
|
END DO
|
||||||
|
ILO = 1
|
||||||
|
IHI = N
|
||||||
|
RETURN
|
||||||
|
END IF
|
||||||
|
*
|
||||||
|
* Permutation to isolate eigenvalues if possible.
|
||||||
*
|
*
|
||||||
K = 1
|
K = 1
|
||||||
L = N
|
L = N
|
||||||
*
|
*
|
||||||
IF( N.EQ.0 )
|
IF( .NOT.LSAME( JOB, 'S' ) ) THEN
|
||||||
$ GO TO 210
|
*
|
||||||
|
* Row and column exchange.
|
||||||
|
*
|
||||||
|
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
|
||||||
*
|
*
|
||||||
IF( LSAME( JOB, 'N' ) ) THEN
|
|
||||||
DO 10 I = 1, N
|
|
||||||
SCALE( I ) = ONE
|
|
||||||
10 CONTINUE
|
|
||||||
GO TO 210
|
|
||||||
END IF
|
END IF
|
||||||
*
|
*
|
||||||
IF( LSAME( JOB, 'S' ) )
|
* Initialize SCALE for non-permuted submatrix.
|
||||||
$ GO TO 120
|
|
||||||
*
|
*
|
||||||
* Permutation to isolate eigenvalues if possible
|
DO I = K, L
|
||||||
*
|
|
||||||
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
|
|
||||||
*
|
|
||||||
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
|
|
||||||
*
|
|
||||||
* Search for columns isolating an eigenvalue and push them left.
|
|
||||||
*
|
|
||||||
80 CONTINUE
|
|
||||||
K = K + 1
|
|
||||||
*
|
|
||||||
90 CONTINUE
|
|
||||||
DO 110 J = K, L
|
|
||||||
*
|
|
||||||
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
|
|
||||||
*
|
|
||||||
M = K
|
|
||||||
IEXC = 2
|
|
||||||
GO TO 20
|
|
||||||
110 CONTINUE
|
|
||||||
*
|
|
||||||
120 CONTINUE
|
|
||||||
DO 130 I = K, L
|
|
||||||
SCALE( I ) = ONE
|
SCALE( I ) = ONE
|
||||||
130 CONTINUE
|
END DO
|
||||||
*
|
*
|
||||||
IF( LSAME( JOB, 'P' ) )
|
* If we only had to permute, we are done.
|
||||||
$ GO TO 210
|
*
|
||||||
|
IF( LSAME( JOB, 'P' ) ) THEN
|
||||||
|
ILO = K
|
||||||
|
IHI = L
|
||||||
|
RETURN
|
||||||
|
END IF
|
||||||
*
|
*
|
||||||
* Balance the submatrix in rows K to L.
|
* Balance the submatrix in rows K to L.
|
||||||
*
|
*
|
||||||
* Iterative loop for norm reduction
|
* Iterative loop for norm reduction.
|
||||||
*
|
*
|
||||||
SFMIN1 = SLAMCH( 'S' ) / SLAMCH( 'P' )
|
SFMIN1 = SLAMCH( 'S' ) / SLAMCH( 'P' )
|
||||||
SFMAX1 = ONE / SFMIN1
|
SFMAX1 = ONE / SFMIN1
|
||||||
SFMIN2 = SFMIN1*SCLFAC
|
SFMIN2 = SFMIN1*SCLFAC
|
||||||
SFMAX2 = ONE / SFMIN2
|
SFMAX2 = ONE / SFMIN2
|
||||||
140 CONTINUE
|
|
||||||
NOCONV = .FALSE.
|
|
||||||
*
|
*
|
||||||
DO 200 I = K, L
|
NOCONV = .TRUE.
|
||||||
|
DO WHILE( NOCONV )
|
||||||
|
NOCONV = .FALSE.
|
||||||
*
|
*
|
||||||
C = SCNRM2( L-K+1, A( K, I ), 1 )
|
DO I = K, L
|
||||||
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.
|
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 ) )
|
||||||
*
|
*
|
||||||
IF( C.EQ.ZERO .OR. R.EQ.ZERO )
|
* Guard against zero C or R due to underflow.
|
||||||
$ GO TO 200
|
*
|
||||||
G = R / SCLFAC
|
IF( C.EQ.ZERO .OR. R.EQ.ZERO ) CYCLE
|
||||||
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
|
* Exit if NaN to avoid infinite loop
|
||||||
*
|
*
|
||||||
INFO = -3
|
IF( SISNAN( C+CA+R+RA ) ) THEN
|
||||||
CALL XERBLA( 'CGEBAL', -INFO )
|
INFO = -3
|
||||||
RETURN
|
CALL XERBLA( 'CGEBAL', -INFO )
|
||||||
END IF
|
RETURN
|
||||||
F = F*SCLFAC
|
END IF
|
||||||
C = C*SCLFAC
|
|
||||||
CA = CA*SCLFAC
|
|
||||||
R = R / SCLFAC
|
|
||||||
G = G / SCLFAC
|
|
||||||
RA = RA / SCLFAC
|
|
||||||
GO TO 160
|
|
||||||
*
|
*
|
||||||
170 CONTINUE
|
G = R / SCLFAC
|
||||||
G = C / SCLFAC
|
F = ONE
|
||||||
180 CONTINUE
|
S = C + R
|
||||||
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.
|
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
|
||||||
*
|
*
|
||||||
190 CONTINUE
|
G = C / SCLFAC
|
||||||
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.
|
|
||||||
*
|
*
|
||||||
CALL CSSCAL( N-K+1, G, A( I, K ), LDA )
|
DO WHILE( G.GE.R .AND. MAX( R, RA ).LT.SFMAX2 .AND.
|
||||||
CALL CSSCAL( L, F, A( 1, I ), 1 )
|
$ 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
|
||||||
*
|
*
|
||||||
200 CONTINUE
|
* Now balance.
|
||||||
*
|
*
|
||||||
IF( NOCONV )
|
IF( ( C+R ).GE.FACTOR*S ) CYCLE
|
||||||
$ GO TO 140
|
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
|
ILO = K
|
||||||
IHI = L
|
IHI = L
|
||||||
*
|
*
|
||||||
|
|
|
||||||
|
|
@ -153,6 +153,9 @@
|
||||||
*>
|
*>
|
||||||
*> Modified by Tzu-Yi Chen, Computer Science Division, University of
|
*> Modified by Tzu-Yi Chen, Computer Science Division, University of
|
||||||
*> California at Berkeley, USA
|
*> California at Berkeley, USA
|
||||||
|
*>
|
||||||
|
*> Refactored by Evert Provoost, Department of Computer Science,
|
||||||
|
*> KU Leuven, Belgium
|
||||||
*> \endverbatim
|
*> \endverbatim
|
||||||
*>
|
*>
|
||||||
* =====================================================================
|
* =====================================================================
|
||||||
|
|
@ -181,8 +184,8 @@
|
||||||
PARAMETER ( FACTOR = 0.95D+0 )
|
PARAMETER ( FACTOR = 0.95D+0 )
|
||||||
* ..
|
* ..
|
||||||
* .. Local Scalars ..
|
* .. Local Scalars ..
|
||||||
LOGICAL NOCONV
|
LOGICAL NOCONV, CANSWAP
|
||||||
INTEGER I, ICA, IEXC, IRA, J, K, L, M
|
INTEGER I, ICA, IRA, J, K, L
|
||||||
DOUBLE PRECISION C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
|
DOUBLE PRECISION C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
|
||||||
$ SFMIN2
|
$ SFMIN2
|
||||||
* ..
|
* ..
|
||||||
|
|
@ -213,178 +216,193 @@
|
||||||
CALL XERBLA( 'DGEBAL', -INFO )
|
CALL XERBLA( 'DGEBAL', -INFO )
|
||||||
RETURN
|
RETURN
|
||||||
END IF
|
END IF
|
||||||
|
*
|
||||||
|
* Quick returns.
|
||||||
|
*
|
||||||
|
IF( N.EQ.0 ) THEN
|
||||||
|
ILO = 1
|
||||||
|
IHI = 0
|
||||||
|
RETURN
|
||||||
|
END IF
|
||||||
|
*
|
||||||
|
IF( LSAME( JOB, 'N' ) ) THEN
|
||||||
|
DO I = 1, N
|
||||||
|
SCALE( I ) = ONE
|
||||||
|
END DO
|
||||||
|
ILO = 1
|
||||||
|
IHI = N
|
||||||
|
RETURN
|
||||||
|
END IF
|
||||||
|
*
|
||||||
|
* Permutation to isolate eigenvalues if possible.
|
||||||
*
|
*
|
||||||
K = 1
|
K = 1
|
||||||
L = N
|
L = N
|
||||||
*
|
*
|
||||||
IF( N.EQ.0 )
|
IF( .NOT.LSAME( JOB, 'S' ) ) THEN
|
||||||
$ GO TO 210
|
*
|
||||||
|
* Row and column exchange.
|
||||||
|
*
|
||||||
|
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. 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 DSWAP( L, A( 1, I ), 1, A( 1, L ), 1 )
|
||||||
|
CALL DSWAP( 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. 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 DSWAP( L, A( 1, J ), 1, A( 1, K ), 1 )
|
||||||
|
CALL DSWAP( N-K+1, A( J, K ), LDA, A( K, K ), LDA )
|
||||||
|
END IF
|
||||||
|
NOCONV = .TRUE.
|
||||||
|
*
|
||||||
|
K = K + 1
|
||||||
|
END IF
|
||||||
|
END DO
|
||||||
|
*
|
||||||
|
END DO
|
||||||
*
|
*
|
||||||
IF( LSAME( JOB, 'N' ) ) THEN
|
|
||||||
DO 10 I = 1, N
|
|
||||||
SCALE( I ) = ONE
|
|
||||||
10 CONTINUE
|
|
||||||
GO TO 210
|
|
||||||
END IF
|
END IF
|
||||||
*
|
*
|
||||||
IF( LSAME( JOB, 'S' ) )
|
* Initialize SCALE for non-permuted submatrix.
|
||||||
$ GO TO 120
|
|
||||||
*
|
*
|
||||||
* Permutation to isolate eigenvalues if possible
|
DO I = K, L
|
||||||
*
|
|
||||||
GO TO 50
|
|
||||||
*
|
|
||||||
* Row and column exchange.
|
|
||||||
*
|
|
||||||
20 CONTINUE
|
|
||||||
SCALE( M ) = J
|
|
||||||
IF( J.EQ.M )
|
|
||||||
$ GO TO 30
|
|
||||||
*
|
|
||||||
CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
|
|
||||||
CALL DSWAP( 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
|
|
||||||
*
|
|
||||||
DO 60 I = 1, L
|
|
||||||
IF( I.EQ.J )
|
|
||||||
$ GO TO 60
|
|
||||||
IF( A( J, I ).NE.ZERO )
|
|
||||||
$ GO TO 70
|
|
||||||
60 CONTINUE
|
|
||||||
*
|
|
||||||
M = L
|
|
||||||
IEXC = 1
|
|
||||||
GO TO 20
|
|
||||||
70 CONTINUE
|
|
||||||
*
|
|
||||||
GO TO 90
|
|
||||||
*
|
|
||||||
* Search for columns isolating an eigenvalue and push them left.
|
|
||||||
*
|
|
||||||
80 CONTINUE
|
|
||||||
K = K + 1
|
|
||||||
*
|
|
||||||
90 CONTINUE
|
|
||||||
DO 110 J = K, L
|
|
||||||
*
|
|
||||||
DO 100 I = K, L
|
|
||||||
IF( I.EQ.J )
|
|
||||||
$ GO TO 100
|
|
||||||
IF( A( I, J ).NE.ZERO )
|
|
||||||
$ GO TO 110
|
|
||||||
100 CONTINUE
|
|
||||||
*
|
|
||||||
M = K
|
|
||||||
IEXC = 2
|
|
||||||
GO TO 20
|
|
||||||
110 CONTINUE
|
|
||||||
*
|
|
||||||
120 CONTINUE
|
|
||||||
DO 130 I = K, L
|
|
||||||
SCALE( I ) = ONE
|
SCALE( I ) = ONE
|
||||||
130 CONTINUE
|
END DO
|
||||||
*
|
*
|
||||||
IF( LSAME( JOB, 'P' ) )
|
* If we only had to permute, we are done.
|
||||||
$ GO TO 210
|
*
|
||||||
|
IF( LSAME( JOB, 'P' ) ) THEN
|
||||||
|
ILO = K
|
||||||
|
IHI = L
|
||||||
|
RETURN
|
||||||
|
END IF
|
||||||
*
|
*
|
||||||
* Balance the submatrix in rows K to L.
|
* Balance the submatrix in rows K to L.
|
||||||
*
|
*
|
||||||
* Iterative loop for norm reduction
|
* Iterative loop for norm reduction.
|
||||||
*
|
*
|
||||||
SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
|
SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
|
||||||
SFMAX1 = ONE / SFMIN1
|
SFMAX1 = ONE / SFMIN1
|
||||||
SFMIN2 = SFMIN1*SCLFAC
|
SFMIN2 = SFMIN1*SCLFAC
|
||||||
SFMAX2 = ONE / SFMIN2
|
SFMAX2 = ONE / SFMIN2
|
||||||
*
|
*
|
||||||
140 CONTINUE
|
NOCONV = .TRUE.
|
||||||
NOCONV = .FALSE.
|
DO WHILE( NOCONV )
|
||||||
|
NOCONV = .FALSE.
|
||||||
*
|
*
|
||||||
DO 200 I = K, L
|
DO I = K, L
|
||||||
*
|
*
|
||||||
C = DNRM2( L-K+1, A( K, I ), 1 )
|
C = DNRM2( L-K+1, A( K, I ), 1 )
|
||||||
R = DNRM2( L-K+1, A( I, K ), LDA )
|
R = DNRM2( L-K+1, A( I, K ), LDA )
|
||||||
ICA = IDAMAX( L, A( 1, I ), 1 )
|
ICA = IDAMAX( L, A( 1, I ), 1 )
|
||||||
CA = ABS( A( ICA, I ) )
|
CA = ABS( A( ICA, I ) )
|
||||||
IRA = IDAMAX( N-K+1, A( I, K ), LDA )
|
IRA = IDAMAX( N-K+1, A( I, K ), LDA )
|
||||||
RA = ABS( A( I, IRA+K-1 ) )
|
RA = ABS( A( I, IRA+K-1 ) )
|
||||||
*
|
*
|
||||||
* Guard against zero C or R due to underflow.
|
* Guard against zero C or R due to underflow.
|
||||||
*
|
*
|
||||||
IF( C.EQ.ZERO .OR. R.EQ.ZERO )
|
IF( C.EQ.ZERO .OR. R.EQ.ZERO ) CYCLE
|
||||||
$ 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( DISNAN( C+F+CA+R+G+RA ) ) THEN
|
|
||||||
*
|
*
|
||||||
* Exit if NaN to avoid infinite loop
|
* Exit if NaN to avoid infinite loop
|
||||||
*
|
*
|
||||||
INFO = -3
|
IF( DISNAN( C+CA+R+RA ) ) THEN
|
||||||
CALL XERBLA( 'DGEBAL', -INFO )
|
INFO = -3
|
||||||
RETURN
|
CALL XERBLA( 'DGEBAL', -INFO )
|
||||||
END IF
|
RETURN
|
||||||
F = F*SCLFAC
|
END IF
|
||||||
C = C*SCLFAC
|
|
||||||
CA = CA*SCLFAC
|
|
||||||
R = R / SCLFAC
|
|
||||||
G = G / SCLFAC
|
|
||||||
RA = RA / SCLFAC
|
|
||||||
GO TO 160
|
|
||||||
*
|
*
|
||||||
170 CONTINUE
|
G = R / SCLFAC
|
||||||
G = C / SCLFAC
|
F = ONE
|
||||||
180 CONTINUE
|
S = C + R
|
||||||
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.
|
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
|
||||||
*
|
*
|
||||||
190 CONTINUE
|
G = C / SCLFAC
|
||||||
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.
|
|
||||||
*
|
*
|
||||||
CALL DSCAL( N-K+1, G, A( I, K ), LDA )
|
DO WHILE( G.GE.R .AND. MAX( R, RA ).LT.SFMAX2 .AND.
|
||||||
CALL DSCAL( L, F, A( 1, I ), 1 )
|
$ 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
|
||||||
*
|
*
|
||||||
200 CONTINUE
|
* Now balance.
|
||||||
*
|
*
|
||||||
IF( NOCONV )
|
IF( ( C+R ).GE.FACTOR*S ) CYCLE
|
||||||
$ GO TO 140
|
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 DSCAL( N-K+1, G, A( I, K ), LDA )
|
||||||
|
CALL DSCAL( L, F, A( 1, I ), 1 )
|
||||||
|
*
|
||||||
|
END DO
|
||||||
|
*
|
||||||
|
END DO
|
||||||
*
|
*
|
||||||
210 CONTINUE
|
|
||||||
ILO = K
|
ILO = K
|
||||||
IHI = L
|
IHI = L
|
||||||
*
|
*
|
||||||
|
|
|
||||||
|
|
@ -153,6 +153,9 @@
|
||||||
*>
|
*>
|
||||||
*> Modified by Tzu-Yi Chen, Computer Science Division, University of
|
*> Modified by Tzu-Yi Chen, Computer Science Division, University of
|
||||||
*> California at Berkeley, USA
|
*> California at Berkeley, USA
|
||||||
|
*>
|
||||||
|
*> Refactored by Evert Provoost, Department of Computer Science,
|
||||||
|
*> KU Leuven, Belgium
|
||||||
*> \endverbatim
|
*> \endverbatim
|
||||||
*>
|
*>
|
||||||
* =====================================================================
|
* =====================================================================
|
||||||
|
|
@ -181,8 +184,8 @@
|
||||||
PARAMETER ( FACTOR = 0.95E+0 )
|
PARAMETER ( FACTOR = 0.95E+0 )
|
||||||
* ..
|
* ..
|
||||||
* .. Local Scalars ..
|
* .. Local Scalars ..
|
||||||
LOGICAL NOCONV
|
LOGICAL NOCONV, CANSWAP
|
||||||
INTEGER I, ICA, IEXC, IRA, J, K, L, M
|
INTEGER I, ICA, IRA, J, K, L
|
||||||
REAL C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
|
REAL C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
|
||||||
$ SFMIN2
|
$ SFMIN2
|
||||||
* ..
|
* ..
|
||||||
|
|
@ -197,7 +200,7 @@
|
||||||
* ..
|
* ..
|
||||||
* .. Intrinsic Functions ..
|
* .. Intrinsic Functions ..
|
||||||
INTRINSIC ABS, MAX, MIN
|
INTRINSIC ABS, MAX, MIN
|
||||||
*
|
* ..
|
||||||
* Test the input parameters
|
* Test the input parameters
|
||||||
*
|
*
|
||||||
INFO = 0
|
INFO = 0
|
||||||
|
|
@ -213,177 +216,193 @@
|
||||||
CALL XERBLA( 'SGEBAL', -INFO )
|
CALL XERBLA( 'SGEBAL', -INFO )
|
||||||
RETURN
|
RETURN
|
||||||
END IF
|
END IF
|
||||||
|
*
|
||||||
|
* Quick returns.
|
||||||
|
*
|
||||||
|
IF( N.EQ.0 ) THEN
|
||||||
|
ILO = 1
|
||||||
|
IHI = 0
|
||||||
|
RETURN
|
||||||
|
END IF
|
||||||
|
*
|
||||||
|
IF( LSAME( JOB, 'N' ) ) THEN
|
||||||
|
DO I = 1, N
|
||||||
|
SCALE( I ) = ONE
|
||||||
|
END DO
|
||||||
|
ILO = 1
|
||||||
|
IHI = N
|
||||||
|
RETURN
|
||||||
|
END IF
|
||||||
|
*
|
||||||
|
* Permutation to isolate eigenvalues if possible.
|
||||||
*
|
*
|
||||||
K = 1
|
K = 1
|
||||||
L = N
|
L = N
|
||||||
*
|
*
|
||||||
IF( N.EQ.0 )
|
IF( .NOT.LSAME( JOB, 'S' ) ) THEN
|
||||||
$ GO TO 210
|
*
|
||||||
|
* Row and column exchange.
|
||||||
|
*
|
||||||
|
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. 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 SSWAP( L, A( 1, I ), 1, A( 1, L ), 1 )
|
||||||
|
CALL SSWAP( 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. 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 SSWAP( L, A( 1, J ), 1, A( 1, K ), 1 )
|
||||||
|
CALL SSWAP( N-K+1, A( J, K ), LDA, A( K, K ), LDA )
|
||||||
|
END IF
|
||||||
|
NOCONV = .TRUE.
|
||||||
|
*
|
||||||
|
K = K + 1
|
||||||
|
END IF
|
||||||
|
END DO
|
||||||
|
*
|
||||||
|
END DO
|
||||||
*
|
*
|
||||||
IF( LSAME( JOB, 'N' ) ) THEN
|
|
||||||
DO 10 I = 1, N
|
|
||||||
SCALE( I ) = ONE
|
|
||||||
10 CONTINUE
|
|
||||||
GO TO 210
|
|
||||||
END IF
|
END IF
|
||||||
*
|
*
|
||||||
IF( LSAME( JOB, 'S' ) )
|
* Initialize SCALE for non-permuted submatrix.
|
||||||
$ GO TO 120
|
|
||||||
*
|
*
|
||||||
* Permutation to isolate eigenvalues if possible
|
DO I = K, L
|
||||||
*
|
|
||||||
GO TO 50
|
|
||||||
*
|
|
||||||
* Row and column exchange.
|
|
||||||
*
|
|
||||||
20 CONTINUE
|
|
||||||
SCALE( M ) = J
|
|
||||||
IF( J.EQ.M )
|
|
||||||
$ GO TO 30
|
|
||||||
*
|
|
||||||
CALL SSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
|
|
||||||
CALL SSWAP( 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
|
|
||||||
*
|
|
||||||
DO 60 I = 1, L
|
|
||||||
IF( I.EQ.J )
|
|
||||||
$ GO TO 60
|
|
||||||
IF( A( J, I ).NE.ZERO )
|
|
||||||
$ GO TO 70
|
|
||||||
60 CONTINUE
|
|
||||||
*
|
|
||||||
M = L
|
|
||||||
IEXC = 1
|
|
||||||
GO TO 20
|
|
||||||
70 CONTINUE
|
|
||||||
*
|
|
||||||
GO TO 90
|
|
||||||
*
|
|
||||||
* Search for columns isolating an eigenvalue and push them left.
|
|
||||||
*
|
|
||||||
80 CONTINUE
|
|
||||||
K = K + 1
|
|
||||||
*
|
|
||||||
90 CONTINUE
|
|
||||||
DO 110 J = K, L
|
|
||||||
*
|
|
||||||
DO 100 I = K, L
|
|
||||||
IF( I.EQ.J )
|
|
||||||
$ GO TO 100
|
|
||||||
IF( A( I, J ).NE.ZERO )
|
|
||||||
$ GO TO 110
|
|
||||||
100 CONTINUE
|
|
||||||
*
|
|
||||||
M = K
|
|
||||||
IEXC = 2
|
|
||||||
GO TO 20
|
|
||||||
110 CONTINUE
|
|
||||||
*
|
|
||||||
120 CONTINUE
|
|
||||||
DO 130 I = K, L
|
|
||||||
SCALE( I ) = ONE
|
SCALE( I ) = ONE
|
||||||
130 CONTINUE
|
END DO
|
||||||
*
|
*
|
||||||
IF( LSAME( JOB, 'P' ) )
|
* If we only had to permute, we are done.
|
||||||
$ GO TO 210
|
*
|
||||||
|
IF( LSAME( JOB, 'P' ) ) THEN
|
||||||
|
ILO = K
|
||||||
|
IHI = L
|
||||||
|
RETURN
|
||||||
|
END IF
|
||||||
*
|
*
|
||||||
* Balance the submatrix in rows K to L.
|
* Balance the submatrix in rows K to L.
|
||||||
*
|
*
|
||||||
* Iterative loop for norm reduction
|
* Iterative loop for norm reduction.
|
||||||
*
|
*
|
||||||
SFMIN1 = SLAMCH( 'S' ) / SLAMCH( 'P' )
|
SFMIN1 = SLAMCH( 'S' ) / SLAMCH( 'P' )
|
||||||
SFMAX1 = ONE / SFMIN1
|
SFMAX1 = ONE / SFMIN1
|
||||||
SFMIN2 = SFMIN1*SCLFAC
|
SFMIN2 = SFMIN1*SCLFAC
|
||||||
SFMAX2 = ONE / SFMIN2
|
SFMAX2 = ONE / SFMIN2
|
||||||
140 CONTINUE
|
|
||||||
NOCONV = .FALSE.
|
|
||||||
*
|
*
|
||||||
DO 200 I = K, L
|
NOCONV = .TRUE.
|
||||||
|
DO WHILE( NOCONV )
|
||||||
|
NOCONV = .FALSE.
|
||||||
*
|
*
|
||||||
C = SNRM2( L-K+1, A( K, I ), 1 )
|
DO I = K, L
|
||||||
R = SNRM2( L-K+1, A( I, K ), LDA )
|
|
||||||
ICA = ISAMAX( L, A( 1, I ), 1 )
|
|
||||||
CA = ABS( A( ICA, I ) )
|
|
||||||
IRA = ISAMAX( N-K+1, A( I, K ), LDA )
|
|
||||||
RA = ABS( A( I, IRA+K-1 ) )
|
|
||||||
*
|
*
|
||||||
* Guard against zero C or R due to underflow.
|
C = SNRM2( L-K+1, A( K, I ), 1 )
|
||||||
|
R = SNRM2( L-K+1, A( I, K ), LDA )
|
||||||
|
ICA = ISAMAX( L, A( 1, I ), 1 )
|
||||||
|
CA = ABS( A( ICA, I ) )
|
||||||
|
IRA = ISAMAX( N-K+1, A( I, K ), LDA )
|
||||||
|
RA = ABS( A( I, IRA+K-1 ) )
|
||||||
*
|
*
|
||||||
IF( C.EQ.ZERO .OR. R.EQ.ZERO )
|
* Guard against zero C or R due to underflow.
|
||||||
$ 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
|
|
||||||
F = F*SCLFAC
|
|
||||||
C = C*SCLFAC
|
|
||||||
CA = CA*SCLFAC
|
|
||||||
R = R / SCLFAC
|
|
||||||
G = G / SCLFAC
|
|
||||||
RA = RA / SCLFAC
|
|
||||||
GO TO 160
|
|
||||||
*
|
*
|
||||||
170 CONTINUE
|
IF( C.EQ.ZERO .OR. R.EQ.ZERO ) CYCLE
|
||||||
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
|
|
||||||
IF( SISNAN( C+F+CA+R+G+RA ) ) THEN
|
|
||||||
*
|
*
|
||||||
* Exit if NaN to avoid infinite loop
|
* Exit if NaN to avoid infinite loop
|
||||||
*
|
*
|
||||||
INFO = -3
|
IF( SISNAN( C+CA+R+RA ) ) THEN
|
||||||
CALL XERBLA( 'SGEBAL', -INFO )
|
INFO = -3
|
||||||
RETURN
|
CALL XERBLA( 'SGEBAL', -INFO )
|
||||||
END IF
|
RETURN
|
||||||
F = F / SCLFAC
|
END IF
|
||||||
C = C / SCLFAC
|
|
||||||
G = G / SCLFAC
|
|
||||||
CA = CA / SCLFAC
|
|
||||||
R = R*SCLFAC
|
|
||||||
RA = RA*SCLFAC
|
|
||||||
GO TO 180
|
|
||||||
*
|
*
|
||||||
* Now balance.
|
G = R / SCLFAC
|
||||||
|
F = ONE
|
||||||
|
S = C + R
|
||||||
*
|
*
|
||||||
190 CONTINUE
|
DO WHILE( C.LT.G .AND. MAX( F, C, CA ).LT.SFMAX2 .AND.
|
||||||
IF( ( C+R ).GE.FACTOR*S )
|
$ MIN( R, G, RA ).GT.SFMIN2 )
|
||||||
$ GO TO 200
|
F = F*SCLFAC
|
||||||
IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
|
C = C*SCLFAC
|
||||||
IF( F*SCALE( I ).LE.SFMIN1 )
|
CA = CA*SCLFAC
|
||||||
$ GO TO 200
|
R = R / SCLFAC
|
||||||
END IF
|
G = G / SCLFAC
|
||||||
IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
|
RA = RA / SCLFAC
|
||||||
IF( SCALE( I ).GE.SFMAX1 / F )
|
END DO
|
||||||
$ GO TO 200
|
|
||||||
END IF
|
|
||||||
G = ONE / F
|
|
||||||
SCALE( I ) = SCALE( I )*F
|
|
||||||
NOCONV = .TRUE.
|
|
||||||
*
|
*
|
||||||
CALL SSCAL( N-K+1, G, A( I, K ), LDA )
|
G = C / SCLFAC
|
||||||
CALL SSCAL( L, F, A( 1, I ), 1 )
|
|
||||||
*
|
*
|
||||||
200 CONTINUE
|
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
|
||||||
*
|
*
|
||||||
IF( NOCONV )
|
* Now balance.
|
||||||
$ GO TO 140
|
*
|
||||||
|
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 SSCAL( N-K+1, G, A( I, K ), LDA )
|
||||||
|
CALL SSCAL( L, F, A( 1, I ), 1 )
|
||||||
|
*
|
||||||
|
END DO
|
||||||
|
*
|
||||||
|
END DO
|
||||||
*
|
*
|
||||||
210 CONTINUE
|
|
||||||
ILO = K
|
ILO = K
|
||||||
IHI = L
|
IHI = L
|
||||||
*
|
*
|
||||||
|
|
|
||||||
|
|
@ -89,7 +89,7 @@
|
||||||
*> \param[out] IHI
|
*> \param[out] IHI
|
||||||
*> \verbatim
|
*> \verbatim
|
||||||
*> IHI is INTEGER
|
*> IHI is INTEGER
|
||||||
*> ILO and IHI are set to INTEGER such that on exit
|
*> ILO and IHI are set to integers such that on exit
|
||||||
*> A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
|
*> A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
|
||||||
*> If JOB = 'N' or 'S', ILO = 1 and IHI = N.
|
*> If JOB = 'N' or 'S', ILO = 1 and IHI = N.
|
||||||
*> \endverbatim
|
*> \endverbatim
|
||||||
|
|
@ -155,6 +155,9 @@
|
||||||
*>
|
*>
|
||||||
*> Modified by Tzu-Yi Chen, Computer Science Division, University of
|
*> Modified by Tzu-Yi Chen, Computer Science Division, University of
|
||||||
*> California at Berkeley, USA
|
*> California at Berkeley, USA
|
||||||
|
*>
|
||||||
|
*> Refactored by Evert Provoost, Department of Computer Science,
|
||||||
|
*> KU Leuven, Belgium
|
||||||
*> \endverbatim
|
*> \endverbatim
|
||||||
*>
|
*>
|
||||||
* =====================================================================
|
* =====================================================================
|
||||||
|
|
@ -184,8 +187,8 @@
|
||||||
PARAMETER ( FACTOR = 0.95D+0 )
|
PARAMETER ( FACTOR = 0.95D+0 )
|
||||||
* ..
|
* ..
|
||||||
* .. Local Scalars ..
|
* .. Local Scalars ..
|
||||||
LOGICAL NOCONV
|
LOGICAL NOCONV, CANSWAP
|
||||||
INTEGER I, ICA, IEXC, IRA, J, K, L, M
|
INTEGER I, ICA, IRA, J, K, L
|
||||||
DOUBLE PRECISION C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
|
DOUBLE PRECISION C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
|
||||||
$ SFMIN2
|
$ SFMIN2
|
||||||
* ..
|
* ..
|
||||||
|
|
@ -216,177 +219,195 @@
|
||||||
CALL XERBLA( 'ZGEBAL', -INFO )
|
CALL XERBLA( 'ZGEBAL', -INFO )
|
||||||
RETURN
|
RETURN
|
||||||
END IF
|
END IF
|
||||||
|
*
|
||||||
|
* Quick returns.
|
||||||
|
*
|
||||||
|
IF( N.EQ.0 ) THEN
|
||||||
|
ILO = 1
|
||||||
|
IHI = 0
|
||||||
|
RETURN
|
||||||
|
END IF
|
||||||
|
*
|
||||||
|
IF( LSAME( JOB, 'N' ) ) THEN
|
||||||
|
DO I = 1, N
|
||||||
|
SCALE( I ) = ONE
|
||||||
|
END DO
|
||||||
|
ILO = 1
|
||||||
|
IHI = N
|
||||||
|
RETURN
|
||||||
|
END IF
|
||||||
|
*
|
||||||
|
* Permutation to isolate eigenvalues if possible.
|
||||||
*
|
*
|
||||||
K = 1
|
K = 1
|
||||||
L = N
|
L = N
|
||||||
*
|
*
|
||||||
IF( N.EQ.0 )
|
IF( .NOT.LSAME( JOB, 'S' ) ) THEN
|
||||||
$ GO TO 210
|
*
|
||||||
|
* Row and column exchange.
|
||||||
|
*
|
||||||
|
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. ( DBLE( A( I, J ) ).NE.ZERO .OR.
|
||||||
|
$ DIMAG( 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 ZSWAP( L, A( 1, I ), 1, A( 1, L ), 1 )
|
||||||
|
CALL ZSWAP( 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. ( DBLE( A( I, J ) ).NE.ZERO .OR.
|
||||||
|
$ DIMAG( 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 ZSWAP( L, A( 1, J ), 1, A( 1, K ), 1 )
|
||||||
|
CALL ZSWAP( N-K+1, A( J, K ), LDA, A( K, K ), LDA )
|
||||||
|
END IF
|
||||||
|
NOCONV = .TRUE.
|
||||||
|
*
|
||||||
|
K = K + 1
|
||||||
|
END IF
|
||||||
|
END DO
|
||||||
|
*
|
||||||
|
END DO
|
||||||
*
|
*
|
||||||
IF( LSAME( JOB, 'N' ) ) THEN
|
|
||||||
DO 10 I = 1, N
|
|
||||||
SCALE( I ) = ONE
|
|
||||||
10 CONTINUE
|
|
||||||
GO TO 210
|
|
||||||
END IF
|
END IF
|
||||||
*
|
*
|
||||||
IF( LSAME( JOB, 'S' ) )
|
* Initialize SCALE for non-permuted submatrix.
|
||||||
$ GO TO 120
|
|
||||||
*
|
*
|
||||||
* Permutation to isolate eigenvalues if possible
|
DO I = K, L
|
||||||
*
|
|
||||||
GO TO 50
|
|
||||||
*
|
|
||||||
* Row and column exchange.
|
|
||||||
*
|
|
||||||
20 CONTINUE
|
|
||||||
SCALE( M ) = J
|
|
||||||
IF( J.EQ.M )
|
|
||||||
$ GO TO 30
|
|
||||||
*
|
|
||||||
CALL ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
|
|
||||||
CALL ZSWAP( 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
|
|
||||||
*
|
|
||||||
DO 60 I = 1, L
|
|
||||||
IF( I.EQ.J )
|
|
||||||
$ GO TO 60
|
|
||||||
IF( DBLE( A( J, I ) ).NE.ZERO .OR. DIMAG( A( J, I ) ).NE.
|
|
||||||
$ ZERO )GO TO 70
|
|
||||||
60 CONTINUE
|
|
||||||
*
|
|
||||||
M = L
|
|
||||||
IEXC = 1
|
|
||||||
GO TO 20
|
|
||||||
70 CONTINUE
|
|
||||||
*
|
|
||||||
GO TO 90
|
|
||||||
*
|
|
||||||
* Search for columns isolating an eigenvalue and push them left.
|
|
||||||
*
|
|
||||||
80 CONTINUE
|
|
||||||
K = K + 1
|
|
||||||
*
|
|
||||||
90 CONTINUE
|
|
||||||
DO 110 J = K, L
|
|
||||||
*
|
|
||||||
DO 100 I = K, L
|
|
||||||
IF( I.EQ.J )
|
|
||||||
$ GO TO 100
|
|
||||||
IF( DBLE( A( I, J ) ).NE.ZERO .OR. DIMAG( A( I, J ) ).NE.
|
|
||||||
$ ZERO )GO TO 110
|
|
||||||
100 CONTINUE
|
|
||||||
*
|
|
||||||
M = K
|
|
||||||
IEXC = 2
|
|
||||||
GO TO 20
|
|
||||||
110 CONTINUE
|
|
||||||
*
|
|
||||||
120 CONTINUE
|
|
||||||
DO 130 I = K, L
|
|
||||||
SCALE( I ) = ONE
|
SCALE( I ) = ONE
|
||||||
130 CONTINUE
|
END DO
|
||||||
*
|
*
|
||||||
IF( LSAME( JOB, 'P' ) )
|
* If we only had to permute, we are done.
|
||||||
$ GO TO 210
|
*
|
||||||
|
IF( LSAME( JOB, 'P' ) ) THEN
|
||||||
|
ILO = K
|
||||||
|
IHI = L
|
||||||
|
RETURN
|
||||||
|
END IF
|
||||||
*
|
*
|
||||||
* Balance the submatrix in rows K to L.
|
* Balance the submatrix in rows K to L.
|
||||||
*
|
*
|
||||||
* Iterative loop for norm reduction
|
* Iterative loop for norm reduction.
|
||||||
*
|
*
|
||||||
SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
|
SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
|
||||||
SFMAX1 = ONE / SFMIN1
|
SFMAX1 = ONE / SFMIN1
|
||||||
SFMIN2 = SFMIN1*SCLFAC
|
SFMIN2 = SFMIN1*SCLFAC
|
||||||
SFMAX2 = ONE / SFMIN2
|
SFMAX2 = ONE / SFMIN2
|
||||||
140 CONTINUE
|
|
||||||
NOCONV = .FALSE.
|
|
||||||
*
|
*
|
||||||
DO 200 I = K, L
|
NOCONV = .TRUE.
|
||||||
|
DO WHILE( NOCONV )
|
||||||
|
NOCONV = .FALSE.
|
||||||
*
|
*
|
||||||
C = DZNRM2( L-K+1, A( K, I ), 1 )
|
DO I = K, L
|
||||||
R = DZNRM2( L-K+1, A( I, K ), LDA )
|
|
||||||
ICA = IZAMAX( L, A( 1, I ), 1 )
|
|
||||||
CA = ABS( A( ICA, I ) )
|
|
||||||
IRA = IZAMAX( N-K+1, A( I, K ), LDA )
|
|
||||||
RA = ABS( A( I, IRA+K-1 ) )
|
|
||||||
*
|
*
|
||||||
* Guard against zero C or R due to underflow.
|
C = DZNRM2( L-K+1, A( K, I ), 1 )
|
||||||
|
R = DZNRM2( L-K+1, A( I, K ), LDA )
|
||||||
|
ICA = IZAMAX( L, A( 1, I ), 1 )
|
||||||
|
CA = ABS( A( ICA, I ) )
|
||||||
|
IRA = IZAMAX( N-K+1, A( I, K ), LDA )
|
||||||
|
RA = ABS( A( I, IRA+K-1 ) )
|
||||||
*
|
*
|
||||||
IF( C.EQ.ZERO .OR. R.EQ.ZERO )
|
* Guard against zero C or R due to underflow.
|
||||||
$ GO TO 200
|
*
|
||||||
G = R / SCLFAC
|
IF( C.EQ.ZERO .OR. R.EQ.ZERO ) CYCLE
|
||||||
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( DISNAN( C+F+CA+R+G+RA ) ) THEN
|
|
||||||
*
|
*
|
||||||
* Exit if NaN to avoid infinite loop
|
* Exit if NaN to avoid infinite loop
|
||||||
*
|
*
|
||||||
INFO = -3
|
IF( DISNAN( C+CA+R+RA ) ) THEN
|
||||||
CALL XERBLA( 'ZGEBAL', -INFO )
|
INFO = -3
|
||||||
RETURN
|
CALL XERBLA( 'ZGEBAL', -INFO )
|
||||||
END IF
|
RETURN
|
||||||
F = F*SCLFAC
|
END IF
|
||||||
C = C*SCLFAC
|
|
||||||
CA = CA*SCLFAC
|
|
||||||
R = R / SCLFAC
|
|
||||||
G = G / SCLFAC
|
|
||||||
RA = RA / SCLFAC
|
|
||||||
GO TO 160
|
|
||||||
*
|
*
|
||||||
170 CONTINUE
|
G = R / SCLFAC
|
||||||
G = C / SCLFAC
|
F = ONE
|
||||||
180 CONTINUE
|
S = C + R
|
||||||
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.
|
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
|
||||||
*
|
*
|
||||||
190 CONTINUE
|
G = C / SCLFAC
|
||||||
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.
|
|
||||||
*
|
*
|
||||||
CALL ZDSCAL( N-K+1, G, A( I, K ), LDA )
|
DO WHILE( G.GE.R .AND. MAX( R, RA ).LT.SFMAX2 .AND.
|
||||||
CALL ZDSCAL( L, F, A( 1, I ), 1 )
|
$ 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
|
||||||
*
|
*
|
||||||
200 CONTINUE
|
* Now balance.
|
||||||
*
|
*
|
||||||
IF( NOCONV )
|
IF( ( C+R ).GE.FACTOR*S ) CYCLE
|
||||||
$ GO TO 140
|
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 ZDSCAL( N-K+1, G, A( I, K ), LDA )
|
||||||
|
CALL ZDSCAL( L, F, A( 1, I ), 1 )
|
||||||
|
*
|
||||||
|
END DO
|
||||||
|
*
|
||||||
|
END DO
|
||||||
*
|
*
|
||||||
210 CONTINUE
|
|
||||||
ILO = K
|
ILO = K
|
||||||
IHI = L
|
IHI = L
|
||||||
*
|
*
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue