diff options
Diffstat (limited to 'SRC/cgebal.f')
| -rw-r--r-- | SRC/cgebal.f | 330 |
1 files changed, 330 insertions, 0 deletions
diff --git a/SRC/cgebal.f b/SRC/cgebal.f new file mode 100644 index 00000000..12394eac --- /dev/null +++ b/SRC/cgebal.f @@ -0,0 +1,330 @@ + SUBROUTINE CGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO ) +* +* -- LAPACK routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER JOB + INTEGER IHI, ILO, INFO, LDA, N +* .. +* .. Array Arguments .. + REAL SCALE( * ) + COMPLEX A( LDA, * ) +* .. +* +* Purpose +* ======= +* +* CGEBAL balances a general complex matrix A. This involves, first, +* permuting A by a similarity transformation to isolate eigenvalues +* in the first 1 to ILO-1 and last IHI+1 to N elements on the +* diagonal; and second, applying a diagonal similarity transformation +* to rows and columns ILO to IHI to make the rows and columns as +* close in norm as possible. Both steps are optional. +* +* Balancing may reduce the 1-norm of the matrix, and improve the +* accuracy of the computed eigenvalues and/or eigenvectors. +* +* Arguments +* ========= +* +* JOB (input) CHARACTER*1 +* Specifies the operations to be performed on A: +* = 'N': none: simply set ILO = 1, IHI = N, SCALE(I) = 1.0 +* for i = 1,...,N; +* = 'P': permute only; +* = 'S': scale only; +* = 'B': both permute and scale. +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input/output) COMPLEX array, dimension (LDA,N) +* On entry, the input matrix A. +* On exit, A is overwritten by the balanced matrix. +* If JOB = 'N', A is not referenced. +* See Further Details. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* ILO (output) INTEGER +* IHI (output) INTEGER +* 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. +* If JOB = 'N' or 'S', ILO = 1 and IHI = N. +* +* SCALE (output) REAL array, dimension (N) +* Details of the permutations and scaling factors applied to +* A. If P(j) is the index of the row and column interchanged +* with row and column j and D(j) is the scaling factor +* applied to row and column j, then +* SCALE(j) = P(j) for j = 1,...,ILO-1 +* = D(j) for j = ILO,...,IHI +* = P(j) for j = IHI+1,...,N. +* The order in which the interchanges are made is N to IHI+1, +* then 1 to ILO-1. +* +* INFO (output) INTEGER +* = 0: successful exit. +* < 0: if INFO = -i, the i-th argument had an illegal value. +* +* Further Details +* =============== +* +* The permutations consist of row and column interchanges which put +* the matrix in the form +* +* ( T1 X Y ) +* P A P = ( 0 B Z ) +* ( 0 0 T2 ) +* +* where T1 and T2 are upper triangular matrices whose eigenvalues lie +* along the diagonal. The column indices ILO and IHI mark the starting +* and ending columns of the submatrix B. Balancing consists of applying +* a diagonal similarity transformation inv(D) * B * D to make the +* 1-norms of each row of B and its corresponding column nearly equal. +* The output matrix is +* +* ( T1 X*D Y ) +* ( 0 inv(D)*B*D inv(D)*Z ). +* ( 0 0 T2 ) +* +* Information about the permutations P and the diagonal matrix D is +* returned in the vector SCALE. +* +* This subroutine is based on the EISPACK routine CBAL. +* +* Modified by Tzu-Yi Chen, Computer Science Division, University of +* California at Berkeley, USA +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) + REAL SCLFAC + PARAMETER ( SCLFAC = 2.0E+0 ) + REAL FACTOR + PARAMETER ( FACTOR = 0.95E+0 ) +* .. +* .. Local Scalars .. + LOGICAL NOCONV + INTEGER I, ICA, IEXC, IRA, J, K, L, M + REAL C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1, + $ SFMIN2 + COMPLEX CDUM +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ICAMAX + REAL SLAMCH + EXTERNAL LSAME, ICAMAX, SLAMCH +* .. +* .. External Subroutines .. + EXTERNAL CSSCAL, CSWAP, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, AIMAG, MAX, MIN, REAL +* .. +* .. Statement Functions .. + REAL CABS1 +* .. +* .. Statement Function definitions .. + CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) ) +* .. +* .. Executable Statements .. +* +* Test the input parameters +* + INFO = 0 + IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND. + $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CGEBAL', -INFO ) + RETURN + END IF +* + K = 1 + L = N +* + IF( N.EQ.0 ) + $ GO TO 210 +* + IF( LSAME( JOB, 'N' ) ) THEN + DO 10 I = 1, N + SCALE( I ) = ONE + 10 CONTINUE + GO TO 210 + 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 +* + 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 + 130 CONTINUE +* + IF( LSAME( JOB, 'P' ) ) + $ GO TO 210 +* +* Balance the submatrix in rows K to L. +* +* 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 = ZERO + R = ZERO +* + DO 150 J = K, L + IF( J.EQ.I ) + $ GO TO 150 + C = C + CABS1( A( J, I ) ) + R = R + CABS1( A( I, J ) ) + 150 CONTINUE + 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 + 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. +* + CALL CSSCAL( N-K+1, G, A( I, K ), LDA ) + CALL CSSCAL( L, F, A( 1, I ), 1 ) +* + 200 CONTINUE +* + IF( NOCONV ) + $ GO TO 140 +* + 210 CONTINUE + ILO = K + IHI = L +* + RETURN +* +* End of CGEBAL +* + END |