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