 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.2.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * June 2010 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER JOB 00010 INTEGER IHI, ILO, INFO, LDA, N 00011 * .. 00012 * .. Array Arguments .. 00013 REAL SCALE( * ) 00014 COMPLEX A( LDA, * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * CGEBAL balances a general complex matrix A. This involves, first, 00021 * permuting A by a similarity transformation to isolate eigenvalues 00022 * in the first 1 to ILO-1 and last IHI+1 to N elements on the 00023 * diagonal; and second, applying a diagonal similarity transformation 00024 * to rows and columns ILO to IHI to make the rows and columns as 00025 * close in norm as possible. Both steps are optional. 00026 * 00027 * Balancing may reduce the 1-norm of the matrix, and improve the 00028 * accuracy of the computed eigenvalues and/or eigenvectors. 00029 * 00030 * Arguments 00031 * ========= 00032 * 00033 * JOB (input) CHARACTER*1 00034 * Specifies the operations to be performed on A: 00035 * = 'N': none: simply set ILO = 1, IHI = N, SCALE(I) = 1.0 00036 * for i = 1,...,N; 00037 * = 'P': permute only; 00038 * = 'S': scale only; 00039 * = 'B': both permute and scale. 00040 * 00041 * N (input) INTEGER 00042 * The order of the matrix A. N >= 0. 00043 * 00044 * A (input/output) COMPLEX array, dimension (LDA,N) 00045 * On entry, the input matrix A. 00046 * On exit, A is overwritten by the balanced matrix. 00047 * If JOB = 'N', A is not referenced. 00048 * See Further Details. 00049 * 00050 * LDA (input) INTEGER 00051 * The leading dimension of the array A. LDA >= max(1,N). 00052 * 00053 * ILO (output) INTEGER 00054 * IHI (output) INTEGER 00055 * ILO and IHI are set to integers such that on exit 00056 * A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N. 00057 * If JOB = 'N' or 'S', ILO = 1 and IHI = N. 00058 * 00059 * SCALE (output) REAL array, dimension (N) 00060 * Details of the permutations and scaling factors applied to 00061 * A. If P(j) is the index of the row and column interchanged 00062 * with row and column j and D(j) is the scaling factor 00063 * applied to row and column j, then 00064 * SCALE(j) = P(j) for j = 1,...,ILO-1 00065 * = D(j) for j = ILO,...,IHI 00066 * = P(j) for j = IHI+1,...,N. 00067 * The order in which the interchanges are made is N to IHI+1, 00068 * then 1 to ILO-1. 00069 * 00070 * INFO (output) INTEGER 00071 * = 0: successful exit. 00072 * < 0: if INFO = -i, the i-th argument had an illegal value. 00073 * 00074 * Further Details 00075 * =============== 00076 * 00077 * The permutations consist of row and column interchanges which put 00078 * the matrix in the form 00079 * 00080 * ( T1 X Y ) 00081 * P A P = ( 0 B Z ) 00082 * ( 0 0 T2 ) 00083 * 00084 * where T1 and T2 are upper triangular matrices whose eigenvalues lie 00085 * along the diagonal. The column indices ILO and IHI mark the starting 00086 * and ending columns of the submatrix B. Balancing consists of applying 00087 * a diagonal similarity transformation inv(D) * B * D to make the 00088 * 1-norms of each row of B and its corresponding column nearly equal. 00089 * The output matrix is 00090 * 00091 * ( T1 X*D Y ) 00092 * ( 0 inv(D)*B*D inv(D)*Z ). 00093 * ( 0 0 T2 ) 00094 * 00095 * Information about the permutations P and the diagonal matrix D is 00096 * returned in the vector SCALE. 00097 * 00098 * This subroutine is based on the EISPACK routine CBAL. 00099 * 00100 * Modified by Tzu-Yi Chen, Computer Science Division, University of 00101 * California at Berkeley, USA 00102 * 00103 * ===================================================================== 00104 * 00105 * .. Parameters .. 00106 REAL ZERO, ONE 00107 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00108 REAL SCLFAC 00109 PARAMETER ( SCLFAC = 2.0E+0 ) 00110 REAL FACTOR 00111 PARAMETER ( FACTOR = 0.95E+0 ) 00112 * .. 00113 * .. Local Scalars .. 00114 LOGICAL NOCONV 00115 INTEGER I, ICA, IEXC, IRA, J, K, L, M 00116 REAL C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1, 00117 $ SFMIN2 00118 COMPLEX CDUM 00119 * .. 00120 * .. External Functions .. 00121 LOGICAL SISNAN, LSAME 00122 INTEGER ICAMAX 00123 REAL SLAMCH 00124 EXTERNAL SISNAN, LSAME, ICAMAX, SLAMCH 00125 * .. 00126 * .. External Subroutines .. 00127 EXTERNAL CSSCAL, CSWAP, XERBLA 00128 * .. 00129 * .. Intrinsic Functions .. 00130 INTRINSIC ABS, AIMAG, MAX, MIN, REAL 00131 * .. 00132 * .. Statement Functions .. 00133 REAL CABS1 00134 * .. 00135 * .. Statement Function definitions .. 00136 CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) ) 00137 * .. 00138 * .. Executable Statements .. 00139 * 00140 * Test the input parameters 00141 * 00142 INFO = 0 00143 IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND. 00144 $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN 00145 INFO = -1 00146 ELSE IF( N.LT.0 ) THEN 00147 INFO = -2 00148 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00149 INFO = -4 00150 END IF 00151 IF( INFO.NE.0 ) THEN 00152 CALL XERBLA( 'CGEBAL', -INFO ) 00153 RETURN 00154 END IF 00155 * 00156 K = 1 00157 L = N 00158 * 00159 IF( N.EQ.0 ) 00160 $ GO TO 210 00161 * 00162 IF( LSAME( JOB, 'N' ) ) THEN 00163 DO 10 I = 1, N 00164 SCALE( I ) = ONE 00165 10 CONTINUE 00166 GO TO 210 00167 END IF 00168 * 00169 IF( LSAME( JOB, 'S' ) ) 00170 $ GO TO 120 00171 * 00172 * Permutation to isolate eigenvalues if possible 00173 * 00174 GO TO 50 00175 * 00176 * Row and column exchange. 00177 * 00178 20 CONTINUE 00179 SCALE( M ) = J 00180 IF( J.EQ.M ) 00181 $ GO TO 30 00182 * 00183 CALL CSWAP( L, A( 1, J ), 1, A( 1, M ), 1 ) 00184 CALL CSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA ) 00185 * 00186 30 CONTINUE 00187 GO TO ( 40, 80 )IEXC 00188 * 00189 * Search for rows isolating an eigenvalue and push them down. 00190 * 00191 40 CONTINUE 00192 IF( L.EQ.1 ) 00193 $ GO TO 210 00194 L = L - 1 00195 * 00196 50 CONTINUE 00197 DO 70 J = L, 1, -1 00198 * 00199 DO 60 I = 1, L 00200 IF( I.EQ.J ) 00201 $ GO TO 60 00202 IF( REAL( A( J, I ) ).NE.ZERO .OR. AIMAG( A( J, I ) ).NE. 00203 $ ZERO )GO TO 70 00204 60 CONTINUE 00205 * 00206 M = L 00207 IEXC = 1 00208 GO TO 20 00209 70 CONTINUE 00210 * 00211 GO TO 90 00212 * 00213 * Search for columns isolating an eigenvalue and push them left. 00214 * 00215 80 CONTINUE 00216 K = K + 1 00217 * 00218 90 CONTINUE 00219 DO 110 J = K, L 00220 * 00221 DO 100 I = K, L 00222 IF( I.EQ.J ) 00223 $ GO TO 100 00224 IF( REAL( A( I, J ) ).NE.ZERO .OR. AIMAG( A( I, J ) ).NE. 00225 $ ZERO )GO TO 110 00226 100 CONTINUE 00227 * 00228 M = K 00229 IEXC = 2 00230 GO TO 20 00231 110 CONTINUE 00232 * 00233 120 CONTINUE 00234 DO 130 I = K, L 00235 SCALE( I ) = ONE 00236 130 CONTINUE 00237 * 00238 IF( LSAME( JOB, 'P' ) ) 00239 $ GO TO 210 00240 * 00241 * Balance the submatrix in rows K to L. 00242 * 00243 * Iterative loop for norm reduction 00244 * 00245 SFMIN1 = SLAMCH( 'S' ) / SLAMCH( 'P' ) 00246 SFMAX1 = ONE / SFMIN1 00247 SFMIN2 = SFMIN1*SCLFAC 00248 SFMAX2 = ONE / SFMIN2 00249 140 CONTINUE 00250 NOCONV = .FALSE. 00251 * 00252 DO 200 I = K, L 00253 C = ZERO 00254 R = ZERO 00255 * 00256 DO 150 J = K, L 00257 IF( J.EQ.I ) 00258 $ GO TO 150 00259 C = C + CABS1( A( J, I ) ) 00260 R = R + CABS1( A( I, J ) ) 00261 150 CONTINUE 00262 ICA = ICAMAX( L, A( 1, I ), 1 ) 00263 CA = ABS( A( ICA, I ) ) 00264 IRA = ICAMAX( N-K+1, A( I, K ), LDA ) 00265 RA = ABS( A( I, IRA+K-1 ) ) 00266 * 00267 * Guard against zero C or R due to underflow. 00268 * 00269 IF( C.EQ.ZERO .OR. R.EQ.ZERO ) 00270 $ GO TO 200 00271 G = R / SCLFAC 00272 F = ONE 00273 S = C + R 00274 160 CONTINUE 00275 IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR. 00276 $ MIN( R, G, RA ).LE.SFMIN2 )GO TO 170 00277 IF( SISNAN( C+F+CA+R+G+RA ) ) THEN 00278 * 00279 * Exit if NaN to avoid infinite loop 00280 * 00281 INFO = -3 00282 CALL XERBLA( 'CGEBAL', -INFO ) 00283 RETURN 00284 END IF 00285 F = F*SCLFAC 00286 C = C*SCLFAC 00287 CA = CA*SCLFAC 00288 R = R / SCLFAC 00289 G = G / SCLFAC 00290 RA = RA / SCLFAC 00291 GO TO 160 00292 * 00293 170 CONTINUE 00294 G = C / SCLFAC 00295 180 CONTINUE 00296 IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR. 00297 $ MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190 00298 F = F / SCLFAC 00299 C = C / SCLFAC 00300 G = G / SCLFAC 00301 CA = CA / SCLFAC 00302 R = R*SCLFAC 00303 RA = RA*SCLFAC 00304 GO TO 180 00305 * 00306 * Now balance. 00307 * 00308 190 CONTINUE 00309 IF( ( C+R ).GE.FACTOR*S ) 00310 $ GO TO 200 00311 IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN 00312 IF( F*SCALE( I ).LE.SFMIN1 ) 00313 $ GO TO 200 00314 END IF 00315 IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN 00316 IF( SCALE( I ).GE.SFMAX1 / F ) 00317 $ GO TO 200 00318 END IF 00319 G = ONE / F 00320 SCALE( I ) = SCALE( I )*F 00321 NOCONV = .TRUE. 00322 * 00323 CALL CSSCAL( N-K+1, G, A( I, K ), LDA ) 00324 CALL CSSCAL( L, F, A( 1, I ), 1 ) 00325 * 00326 200 CONTINUE 00327 * 00328 IF( NOCONV ) 00329 $ GO TO 140 00330 * 00331 210 CONTINUE 00332 ILO = K 00333 IHI = L 00334 * 00335 RETURN 00336 * 00337 * End of CGEBAL 00338 * 00339 END
1.7.3 
