 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, 00002 $ RSCALE, WORK, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER JOB 00011 INTEGER IHI, ILO, INFO, LDA, LDB, N 00012 * .. 00013 * .. Array Arguments .. 00014 REAL LSCALE( * ), RSCALE( * ), WORK( * ) 00015 COMPLEX A( LDA, * ), B( LDB, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * CGGBAL balances a pair of general complex matrices (A,B). This 00022 * involves, first, permuting A and B by similarity transformations to 00023 * isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N 00024 * elements on the diagonal; and second, applying a diagonal similarity 00025 * transformation to rows and columns ILO to IHI to make the rows 00026 * and columns as close in norm as possible. Both steps are optional. 00027 * 00028 * Balancing may reduce the 1-norm of the matrices, and improve the 00029 * accuracy of the computed eigenvalues and/or eigenvectors in the 00030 * generalized eigenvalue problem A*x = lambda*B*x. 00031 * 00032 * Arguments 00033 * ========= 00034 * 00035 * JOB (input) CHARACTER*1 00036 * Specifies the operations to be performed on A and B: 00037 * = 'N': none: simply set ILO = 1, IHI = N, LSCALE(I) = 1.0 00038 * and RSCALE(I) = 1.0 for i=1,...,N; 00039 * = 'P': permute only; 00040 * = 'S': scale only; 00041 * = 'B': both permute and scale. 00042 * 00043 * N (input) INTEGER 00044 * The order of the matrices A and B. N >= 0. 00045 * 00046 * A (input/output) COMPLEX array, dimension (LDA,N) 00047 * On entry, the input matrix A. 00048 * On exit, A is overwritten by the balanced matrix. 00049 * If JOB = 'N', A is not referenced. 00050 * 00051 * LDA (input) INTEGER 00052 * The leading dimension of the array A. LDA >= max(1,N). 00053 * 00054 * B (input/output) COMPLEX array, dimension (LDB,N) 00055 * On entry, the input matrix B. 00056 * On exit, B is overwritten by the balanced matrix. 00057 * If JOB = 'N', B is not referenced. 00058 * 00059 * LDB (input) INTEGER 00060 * The leading dimension of the array B. LDB >= max(1,N). 00061 * 00062 * ILO (output) INTEGER 00063 * IHI (output) INTEGER 00064 * ILO and IHI are set to integers such that on exit 00065 * A(i,j) = 0 and B(i,j) = 0 if i > j and 00066 * j = 1,...,ILO-1 or i = IHI+1,...,N. 00067 * If JOB = 'N' or 'S', ILO = 1 and IHI = N. 00068 * 00069 * LSCALE (output) REAL array, dimension (N) 00070 * Details of the permutations and scaling factors applied 00071 * to the left side of A and B. If P(j) is the index of the 00072 * row interchanged with row j, and D(j) is the scaling factor 00073 * applied to row j, then 00074 * LSCALE(j) = P(j) for J = 1,...,ILO-1 00075 * = D(j) for J = ILO,...,IHI 00076 * = P(j) for J = IHI+1,...,N. 00077 * The order in which the interchanges are made is N to IHI+1, 00078 * then 1 to ILO-1. 00079 * 00080 * RSCALE (output) REAL array, dimension (N) 00081 * Details of the permutations and scaling factors applied 00082 * to the right side of A and B. If P(j) is the index of the 00083 * column interchanged with column j, and D(j) is the scaling 00084 * factor applied to column j, then 00085 * RSCALE(j) = P(j) for J = 1,...,ILO-1 00086 * = D(j) for J = ILO,...,IHI 00087 * = P(j) for J = IHI+1,...,N. 00088 * The order in which the interchanges are made is N to IHI+1, 00089 * then 1 to ILO-1. 00090 * 00091 * WORK (workspace) REAL array, dimension (lwork) 00092 * lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and 00093 * at least 1 when JOB = 'N' or 'P'. 00094 * 00095 * INFO (output) INTEGER 00096 * = 0: successful exit 00097 * < 0: if INFO = -i, the i-th argument had an illegal value. 00098 * 00099 * Further Details 00100 * =============== 00101 * 00102 * See R.C. WARD, Balancing the generalized eigenvalue problem, 00103 * SIAM J. Sci. Stat. Comp. 2 (1981), 141-152. 00104 * 00105 * ===================================================================== 00106 * 00107 * .. Parameters .. 00108 REAL ZERO, HALF, ONE 00109 PARAMETER ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0 ) 00110 REAL THREE, SCLFAC 00111 PARAMETER ( THREE = 3.0E+0, SCLFAC = 1.0E+1 ) 00112 COMPLEX CZERO 00113 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) ) 00114 * .. 00115 * .. Local Scalars .. 00116 INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1, 00117 $ K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN, 00118 $ M, NR, NRP2 00119 REAL ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2, 00120 $ COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX, 00121 $ SFMIN, SUM, T, TA, TB, TC 00122 COMPLEX CDUM 00123 * .. 00124 * .. External Functions .. 00125 LOGICAL LSAME 00126 INTEGER ICAMAX 00127 REAL SDOT, SLAMCH 00128 EXTERNAL LSAME, ICAMAX, SDOT, SLAMCH 00129 * .. 00130 * .. External Subroutines .. 00131 EXTERNAL CSSCAL, CSWAP, SAXPY, SSCAL, XERBLA 00132 * .. 00133 * .. Intrinsic Functions .. 00134 INTRINSIC ABS, AIMAG, INT, LOG10, MAX, MIN, REAL, SIGN 00135 * .. 00136 * .. Statement Functions .. 00137 REAL CABS1 00138 * .. 00139 * .. Statement Function definitions .. 00140 CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) ) 00141 * .. 00142 * .. Executable Statements .. 00143 * 00144 * Test the input parameters 00145 * 00146 INFO = 0 00147 IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND. 00148 $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN 00149 INFO = -1 00150 ELSE IF( N.LT.0 ) THEN 00151 INFO = -2 00152 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00153 INFO = -4 00154 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00155 INFO = -6 00156 END IF 00157 IF( INFO.NE.0 ) THEN 00158 CALL XERBLA( 'CGGBAL', -INFO ) 00159 RETURN 00160 END IF 00161 * 00162 * Quick return if possible 00163 * 00164 IF( N.EQ.0 ) THEN 00165 ILO = 1 00166 IHI = N 00167 RETURN 00168 END IF 00169 * 00170 IF( N.EQ.1 ) THEN 00171 ILO = 1 00172 IHI = N 00173 LSCALE( 1 ) = ONE 00174 RSCALE( 1 ) = ONE 00175 RETURN 00176 END IF 00177 * 00178 IF( LSAME( JOB, 'N' ) ) THEN 00179 ILO = 1 00180 IHI = N 00181 DO 10 I = 1, N 00182 LSCALE( I ) = ONE 00183 RSCALE( I ) = ONE 00184 10 CONTINUE 00185 RETURN 00186 END IF 00187 * 00188 K = 1 00189 L = N 00190 IF( LSAME( JOB, 'S' ) ) 00191 $ GO TO 190 00192 * 00193 GO TO 30 00194 * 00195 * Permute the matrices A and B to isolate the eigenvalues. 00196 * 00197 * Find row with one nonzero in columns 1 through L 00198 * 00199 20 CONTINUE 00200 L = LM1 00201 IF( L.NE.1 ) 00202 $ GO TO 30 00203 * 00204 RSCALE( 1 ) = ONE 00205 LSCALE( 1 ) = ONE 00206 GO TO 190 00207 * 00208 30 CONTINUE 00209 LM1 = L - 1 00210 DO 80 I = L, 1, -1 00211 DO 40 J = 1, LM1 00212 JP1 = J + 1 00213 IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO ) 00214 $ GO TO 50 00215 40 CONTINUE 00216 J = L 00217 GO TO 70 00218 * 00219 50 CONTINUE 00220 DO 60 J = JP1, L 00221 IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO ) 00222 $ GO TO 80 00223 60 CONTINUE 00224 J = JP1 - 1 00225 * 00226 70 CONTINUE 00227 M = L 00228 IFLOW = 1 00229 GO TO 160 00230 80 CONTINUE 00231 GO TO 100 00232 * 00233 * Find column with one nonzero in rows K through N 00234 * 00235 90 CONTINUE 00236 K = K + 1 00237 * 00238 100 CONTINUE 00239 DO 150 J = K, L 00240 DO 110 I = K, LM1 00241 IP1 = I + 1 00242 IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO ) 00243 $ GO TO 120 00244 110 CONTINUE 00245 I = L 00246 GO TO 140 00247 120 CONTINUE 00248 DO 130 I = IP1, L 00249 IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO ) 00250 $ GO TO 150 00251 130 CONTINUE 00252 I = IP1 - 1 00253 140 CONTINUE 00254 M = K 00255 IFLOW = 2 00256 GO TO 160 00257 150 CONTINUE 00258 GO TO 190 00259 * 00260 * Permute rows M and I 00261 * 00262 160 CONTINUE 00263 LSCALE( M ) = I 00264 IF( I.EQ.M ) 00265 $ GO TO 170 00266 CALL CSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA ) 00267 CALL CSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB ) 00268 * 00269 * Permute columns M and J 00270 * 00271 170 CONTINUE 00272 RSCALE( M ) = J 00273 IF( J.EQ.M ) 00274 $ GO TO 180 00275 CALL CSWAP( L, A( 1, J ), 1, A( 1, M ), 1 ) 00276 CALL CSWAP( L, B( 1, J ), 1, B( 1, M ), 1 ) 00277 * 00278 180 CONTINUE 00279 GO TO ( 20, 90 )IFLOW 00280 * 00281 190 CONTINUE 00282 ILO = K 00283 IHI = L 00284 * 00285 IF( LSAME( JOB, 'P' ) ) THEN 00286 DO 195 I = ILO, IHI 00287 LSCALE( I ) = ONE 00288 RSCALE( I ) = ONE 00289 195 CONTINUE 00290 RETURN 00291 END IF 00292 * 00293 IF( ILO.EQ.IHI ) 00294 $ RETURN 00295 * 00296 * Balance the submatrix in rows ILO to IHI. 00297 * 00298 NR = IHI - ILO + 1 00299 DO 200 I = ILO, IHI 00300 RSCALE( I ) = ZERO 00301 LSCALE( I ) = ZERO 00302 * 00303 WORK( I ) = ZERO 00304 WORK( I+N ) = ZERO 00305 WORK( I+2*N ) = ZERO 00306 WORK( I+3*N ) = ZERO 00307 WORK( I+4*N ) = ZERO 00308 WORK( I+5*N ) = ZERO 00309 200 CONTINUE 00310 * 00311 * Compute right side vector in resulting linear equations 00312 * 00313 BASL = LOG10( SCLFAC ) 00314 DO 240 I = ILO, IHI 00315 DO 230 J = ILO, IHI 00316 IF( A( I, J ).EQ.CZERO ) THEN 00317 TA = ZERO 00318 GO TO 210 00319 END IF 00320 TA = LOG10( CABS1( A( I, J ) ) ) / BASL 00321 * 00322 210 CONTINUE 00323 IF( B( I, J ).EQ.CZERO ) THEN 00324 TB = ZERO 00325 GO TO 220 00326 END IF 00327 TB = LOG10( CABS1( B( I, J ) ) ) / BASL 00328 * 00329 220 CONTINUE 00330 WORK( I+4*N ) = WORK( I+4*N ) - TA - TB 00331 WORK( J+5*N ) = WORK( J+5*N ) - TA - TB 00332 230 CONTINUE 00333 240 CONTINUE 00334 * 00335 COEF = ONE / REAL( 2*NR ) 00336 COEF2 = COEF*COEF 00337 COEF5 = HALF*COEF2 00338 NRP2 = NR + 2 00339 BETA = ZERO 00340 IT = 1 00341 * 00342 * Start generalized conjugate gradient iteration 00343 * 00344 250 CONTINUE 00345 * 00346 GAMMA = SDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) + 00347 $ SDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 ) 00348 * 00349 EW = ZERO 00350 EWC = ZERO 00351 DO 260 I = ILO, IHI 00352 EW = EW + WORK( I+4*N ) 00353 EWC = EWC + WORK( I+5*N ) 00354 260 CONTINUE 00355 * 00356 GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2 00357 IF( GAMMA.EQ.ZERO ) 00358 $ GO TO 350 00359 IF( IT.NE.1 ) 00360 $ BETA = GAMMA / PGAMMA 00361 T = COEF5*( EWC-THREE*EW ) 00362 TC = COEF5*( EW-THREE*EWC ) 00363 * 00364 CALL SSCAL( NR, BETA, WORK( ILO ), 1 ) 00365 CALL SSCAL( NR, BETA, WORK( ILO+N ), 1 ) 00366 * 00367 CALL SAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 ) 00368 CALL SAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 ) 00369 * 00370 DO 270 I = ILO, IHI 00371 WORK( I ) = WORK( I ) + TC 00372 WORK( I+N ) = WORK( I+N ) + T 00373 270 CONTINUE 00374 * 00375 * Apply matrix to vector 00376 * 00377 DO 300 I = ILO, IHI 00378 KOUNT = 0 00379 SUM = ZERO 00380 DO 290 J = ILO, IHI 00381 IF( A( I, J ).EQ.CZERO ) 00382 $ GO TO 280 00383 KOUNT = KOUNT + 1 00384 SUM = SUM + WORK( J ) 00385 280 CONTINUE 00386 IF( B( I, J ).EQ.CZERO ) 00387 $ GO TO 290 00388 KOUNT = KOUNT + 1 00389 SUM = SUM + WORK( J ) 00390 290 CONTINUE 00391 WORK( I+2*N ) = REAL( KOUNT )*WORK( I+N ) + SUM 00392 300 CONTINUE 00393 * 00394 DO 330 J = ILO, IHI 00395 KOUNT = 0 00396 SUM = ZERO 00397 DO 320 I = ILO, IHI 00398 IF( A( I, J ).EQ.CZERO ) 00399 $ GO TO 310 00400 KOUNT = KOUNT + 1 00401 SUM = SUM + WORK( I+N ) 00402 310 CONTINUE 00403 IF( B( I, J ).EQ.CZERO ) 00404 $ GO TO 320 00405 KOUNT = KOUNT + 1 00406 SUM = SUM + WORK( I+N ) 00407 320 CONTINUE 00408 WORK( J+3*N ) = REAL( KOUNT )*WORK( J ) + SUM 00409 330 CONTINUE 00410 * 00411 SUM = SDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) + 00412 $ SDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 ) 00413 ALPHA = GAMMA / SUM 00414 * 00415 * Determine correction to current iteration 00416 * 00417 CMAX = ZERO 00418 DO 340 I = ILO, IHI 00419 COR = ALPHA*WORK( I+N ) 00420 IF( ABS( COR ).GT.CMAX ) 00421 $ CMAX = ABS( COR ) 00422 LSCALE( I ) = LSCALE( I ) + COR 00423 COR = ALPHA*WORK( I ) 00424 IF( ABS( COR ).GT.CMAX ) 00425 $ CMAX = ABS( COR ) 00426 RSCALE( I ) = RSCALE( I ) + COR 00427 340 CONTINUE 00428 IF( CMAX.LT.HALF ) 00429 $ GO TO 350 00430 * 00431 CALL SAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 ) 00432 CALL SAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 ) 00433 * 00434 PGAMMA = GAMMA 00435 IT = IT + 1 00436 IF( IT.LE.NRP2 ) 00437 $ GO TO 250 00438 * 00439 * End generalized conjugate gradient iteration 00440 * 00441 350 CONTINUE 00442 SFMIN = SLAMCH( 'S' ) 00443 SFMAX = ONE / SFMIN 00444 LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE ) 00445 LSFMAX = INT( LOG10( SFMAX ) / BASL ) 00446 DO 360 I = ILO, IHI 00447 IRAB = ICAMAX( N-ILO+1, A( I, ILO ), LDA ) 00448 RAB = ABS( A( I, IRAB+ILO-1 ) ) 00449 IRAB = ICAMAX( N-ILO+1, B( I, ILO ), LDB ) 00450 RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) ) 00451 LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE ) 00452 IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) ) 00453 IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB ) 00454 LSCALE( I ) = SCLFAC**IR 00455 ICAB = ICAMAX( IHI, A( 1, I ), 1 ) 00456 CAB = ABS( A( ICAB, I ) ) 00457 ICAB = ICAMAX( IHI, B( 1, I ), 1 ) 00458 CAB = MAX( CAB, ABS( B( ICAB, I ) ) ) 00459 LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE ) 00460 JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) ) 00461 JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB ) 00462 RSCALE( I ) = SCLFAC**JC 00463 360 CONTINUE 00464 * 00465 * Row scaling of matrices A and B 00466 * 00467 DO 370 I = ILO, IHI 00468 CALL CSSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA ) 00469 CALL CSSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB ) 00470 370 CONTINUE 00471 * 00472 * Column scaling of matrices A and B 00473 * 00474 DO 380 J = ILO, IHI 00475 CALL CSSCAL( IHI, RSCALE( J ), A( 1, J ), 1 ) 00476 CALL CSSCAL( IHI, RSCALE( J ), B( 1, J ), 1 ) 00477 380 CONTINUE 00478 * 00479 RETURN 00480 * 00481 * End of CGGBAL 00482 * 00483 END
1.7.3 
