 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, 00002 $ RSCALE, WORK, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2.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 * June 2010 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER JOB 00011 INTEGER IHI, ILO, INFO, LDA, LDB, N 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), LSCALE( * ), 00015 $ RSCALE( * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DGGBAL balances a pair of general real 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) 00073 * is the scaling factor 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) DOUBLE PRECISION 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) 00084 * is the scaling factor applied to column j, then 00085 * LSCALE(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) DOUBLE PRECISION 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 DOUBLE PRECISION ZERO, HALF, ONE 00109 PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 ) 00110 DOUBLE PRECISION THREE, SCLFAC 00111 PARAMETER ( THREE = 3.0D+0, SCLFAC = 1.0D+1 ) 00112 * .. 00113 * .. Local Scalars .. 00114 INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1, 00115 $ K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN, 00116 $ M, NR, NRP2 00117 DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2, 00118 $ COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX, 00119 $ SFMIN, SUM, T, TA, TB, TC 00120 * .. 00121 * .. External Functions .. 00122 LOGICAL LSAME 00123 INTEGER IDAMAX 00124 DOUBLE PRECISION DDOT, DLAMCH 00125 EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH 00126 * .. 00127 * .. External Subroutines .. 00128 EXTERNAL DAXPY, DSCAL, DSWAP, XERBLA 00129 * .. 00130 * .. Intrinsic Functions .. 00131 INTRINSIC ABS, DBLE, INT, LOG10, MAX, MIN, SIGN 00132 * .. 00133 * .. Executable Statements .. 00134 * 00135 * Test the input parameters 00136 * 00137 INFO = 0 00138 IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND. 00139 $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN 00140 INFO = -1 00141 ELSE IF( N.LT.0 ) THEN 00142 INFO = -2 00143 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00144 INFO = -4 00145 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00146 INFO = -6 00147 END IF 00148 IF( INFO.NE.0 ) THEN 00149 CALL XERBLA( 'DGGBAL', -INFO ) 00150 RETURN 00151 END IF 00152 * 00153 * Quick return if possible 00154 * 00155 IF( N.EQ.0 ) THEN 00156 ILO = 1 00157 IHI = N 00158 RETURN 00159 END IF 00160 * 00161 IF( N.EQ.1 ) THEN 00162 ILO = 1 00163 IHI = N 00164 LSCALE( 1 ) = ONE 00165 RSCALE( 1 ) = ONE 00166 RETURN 00167 END IF 00168 * 00169 IF( LSAME( JOB, 'N' ) ) THEN 00170 ILO = 1 00171 IHI = N 00172 DO 10 I = 1, N 00173 LSCALE( I ) = ONE 00174 RSCALE( I ) = ONE 00175 10 CONTINUE 00176 RETURN 00177 END IF 00178 * 00179 K = 1 00180 L = N 00181 IF( LSAME( JOB, 'S' ) ) 00182 $ GO TO 190 00183 * 00184 GO TO 30 00185 * 00186 * Permute the matrices A and B to isolate the eigenvalues. 00187 * 00188 * Find row with one nonzero in columns 1 through L 00189 * 00190 20 CONTINUE 00191 L = LM1 00192 IF( L.NE.1 ) 00193 $ GO TO 30 00194 * 00195 RSCALE( 1 ) = ONE 00196 LSCALE( 1 ) = ONE 00197 GO TO 190 00198 * 00199 30 CONTINUE 00200 LM1 = L - 1 00201 DO 80 I = L, 1, -1 00202 DO 40 J = 1, LM1 00203 JP1 = J + 1 00204 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) 00205 $ GO TO 50 00206 40 CONTINUE 00207 J = L 00208 GO TO 70 00209 * 00210 50 CONTINUE 00211 DO 60 J = JP1, L 00212 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) 00213 $ GO TO 80 00214 60 CONTINUE 00215 J = JP1 - 1 00216 * 00217 70 CONTINUE 00218 M = L 00219 IFLOW = 1 00220 GO TO 160 00221 80 CONTINUE 00222 GO TO 100 00223 * 00224 * Find column with one nonzero in rows K through N 00225 * 00226 90 CONTINUE 00227 K = K + 1 00228 * 00229 100 CONTINUE 00230 DO 150 J = K, L 00231 DO 110 I = K, LM1 00232 IP1 = I + 1 00233 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) 00234 $ GO TO 120 00235 110 CONTINUE 00236 I = L 00237 GO TO 140 00238 120 CONTINUE 00239 DO 130 I = IP1, L 00240 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) 00241 $ GO TO 150 00242 130 CONTINUE 00243 I = IP1 - 1 00244 140 CONTINUE 00245 M = K 00246 IFLOW = 2 00247 GO TO 160 00248 150 CONTINUE 00249 GO TO 190 00250 * 00251 * Permute rows M and I 00252 * 00253 160 CONTINUE 00254 LSCALE( M ) = I 00255 IF( I.EQ.M ) 00256 $ GO TO 170 00257 CALL DSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA ) 00258 CALL DSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB ) 00259 * 00260 * Permute columns M and J 00261 * 00262 170 CONTINUE 00263 RSCALE( M ) = J 00264 IF( J.EQ.M ) 00265 $ GO TO 180 00266 CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 ) 00267 CALL DSWAP( L, B( 1, J ), 1, B( 1, M ), 1 ) 00268 * 00269 180 CONTINUE 00270 GO TO ( 20, 90 )IFLOW 00271 * 00272 190 CONTINUE 00273 ILO = K 00274 IHI = L 00275 * 00276 IF( LSAME( JOB, 'P' ) ) THEN 00277 DO 195 I = ILO, IHI 00278 LSCALE( I ) = ONE 00279 RSCALE( I ) = ONE 00280 195 CONTINUE 00281 RETURN 00282 END IF 00283 * 00284 IF( ILO.EQ.IHI ) 00285 $ RETURN 00286 * 00287 * Balance the submatrix in rows ILO to IHI. 00288 * 00289 NR = IHI - ILO + 1 00290 DO 200 I = ILO, IHI 00291 RSCALE( I ) = ZERO 00292 LSCALE( I ) = ZERO 00293 * 00294 WORK( I ) = ZERO 00295 WORK( I+N ) = ZERO 00296 WORK( I+2*N ) = ZERO 00297 WORK( I+3*N ) = ZERO 00298 WORK( I+4*N ) = ZERO 00299 WORK( I+5*N ) = ZERO 00300 200 CONTINUE 00301 * 00302 * Compute right side vector in resulting linear equations 00303 * 00304 BASL = LOG10( SCLFAC ) 00305 DO 240 I = ILO, IHI 00306 DO 230 J = ILO, IHI 00307 TB = B( I, J ) 00308 TA = A( I, J ) 00309 IF( TA.EQ.ZERO ) 00310 $ GO TO 210 00311 TA = LOG10( ABS( TA ) ) / BASL 00312 210 CONTINUE 00313 IF( TB.EQ.ZERO ) 00314 $ GO TO 220 00315 TB = LOG10( ABS( TB ) ) / BASL 00316 220 CONTINUE 00317 WORK( I+4*N ) = WORK( I+4*N ) - TA - TB 00318 WORK( J+5*N ) = WORK( J+5*N ) - TA - TB 00319 230 CONTINUE 00320 240 CONTINUE 00321 * 00322 COEF = ONE / DBLE( 2*NR ) 00323 COEF2 = COEF*COEF 00324 COEF5 = HALF*COEF2 00325 NRP2 = NR + 2 00326 BETA = ZERO 00327 IT = 1 00328 * 00329 * Start generalized conjugate gradient iteration 00330 * 00331 250 CONTINUE 00332 * 00333 GAMMA = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) + 00334 $ DDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 ) 00335 * 00336 EW = ZERO 00337 EWC = ZERO 00338 DO 260 I = ILO, IHI 00339 EW = EW + WORK( I+4*N ) 00340 EWC = EWC + WORK( I+5*N ) 00341 260 CONTINUE 00342 * 00343 GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2 00344 IF( GAMMA.EQ.ZERO ) 00345 $ GO TO 350 00346 IF( IT.NE.1 ) 00347 $ BETA = GAMMA / PGAMMA 00348 T = COEF5*( EWC-THREE*EW ) 00349 TC = COEF5*( EW-THREE*EWC ) 00350 * 00351 CALL DSCAL( NR, BETA, WORK( ILO ), 1 ) 00352 CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 ) 00353 * 00354 CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 ) 00355 CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 ) 00356 * 00357 DO 270 I = ILO, IHI 00358 WORK( I ) = WORK( I ) + TC 00359 WORK( I+N ) = WORK( I+N ) + T 00360 270 CONTINUE 00361 * 00362 * Apply matrix to vector 00363 * 00364 DO 300 I = ILO, IHI 00365 KOUNT = 0 00366 SUM = ZERO 00367 DO 290 J = ILO, IHI 00368 IF( A( I, J ).EQ.ZERO ) 00369 $ GO TO 280 00370 KOUNT = KOUNT + 1 00371 SUM = SUM + WORK( J ) 00372 280 CONTINUE 00373 IF( B( I, J ).EQ.ZERO ) 00374 $ GO TO 290 00375 KOUNT = KOUNT + 1 00376 SUM = SUM + WORK( J ) 00377 290 CONTINUE 00378 WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM 00379 300 CONTINUE 00380 * 00381 DO 330 J = ILO, IHI 00382 KOUNT = 0 00383 SUM = ZERO 00384 DO 320 I = ILO, IHI 00385 IF( A( I, J ).EQ.ZERO ) 00386 $ GO TO 310 00387 KOUNT = KOUNT + 1 00388 SUM = SUM + WORK( I+N ) 00389 310 CONTINUE 00390 IF( B( I, J ).EQ.ZERO ) 00391 $ GO TO 320 00392 KOUNT = KOUNT + 1 00393 SUM = SUM + WORK( I+N ) 00394 320 CONTINUE 00395 WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM 00396 330 CONTINUE 00397 * 00398 SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) + 00399 $ DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 ) 00400 ALPHA = GAMMA / SUM 00401 * 00402 * Determine correction to current iteration 00403 * 00404 CMAX = ZERO 00405 DO 340 I = ILO, IHI 00406 COR = ALPHA*WORK( I+N ) 00407 IF( ABS( COR ).GT.CMAX ) 00408 $ CMAX = ABS( COR ) 00409 LSCALE( I ) = LSCALE( I ) + COR 00410 COR = ALPHA*WORK( I ) 00411 IF( ABS( COR ).GT.CMAX ) 00412 $ CMAX = ABS( COR ) 00413 RSCALE( I ) = RSCALE( I ) + COR 00414 340 CONTINUE 00415 IF( CMAX.LT.HALF ) 00416 $ GO TO 350 00417 * 00418 CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 ) 00419 CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 ) 00420 * 00421 PGAMMA = GAMMA 00422 IT = IT + 1 00423 IF( IT.LE.NRP2 ) 00424 $ GO TO 250 00425 * 00426 * End generalized conjugate gradient iteration 00427 * 00428 350 CONTINUE 00429 SFMIN = DLAMCH( 'S' ) 00430 SFMAX = ONE / SFMIN 00431 LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE ) 00432 LSFMAX = INT( LOG10( SFMAX ) / BASL ) 00433 DO 360 I = ILO, IHI 00434 IRAB = IDAMAX( N-ILO+1, A( I, ILO ), LDA ) 00435 RAB = ABS( A( I, IRAB+ILO-1 ) ) 00436 IRAB = IDAMAX( N-ILO+1, B( I, ILO ), LDB ) 00437 RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) ) 00438 LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE ) 00439 IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) ) 00440 IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB ) 00441 LSCALE( I ) = SCLFAC**IR 00442 ICAB = IDAMAX( IHI, A( 1, I ), 1 ) 00443 CAB = ABS( A( ICAB, I ) ) 00444 ICAB = IDAMAX( IHI, B( 1, I ), 1 ) 00445 CAB = MAX( CAB, ABS( B( ICAB, I ) ) ) 00446 LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE ) 00447 JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) ) 00448 JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB ) 00449 RSCALE( I ) = SCLFAC**JC 00450 360 CONTINUE 00451 * 00452 * Row scaling of matrices A and B 00453 * 00454 DO 370 I = ILO, IHI 00455 CALL DSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA ) 00456 CALL DSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB ) 00457 370 CONTINUE 00458 * 00459 * Column scaling of matrices A and B 00460 * 00461 DO 380 J = ILO, IHI 00462 CALL DSCAL( IHI, RSCALE( J ), A( 1, J ), 1 ) 00463 CALL DSCAL( IHI, RSCALE( J ), B( 1, J ), 1 ) 00464 380 CONTINUE 00465 * 00466 RETURN 00467 * 00468 * End of DGGBAL 00469 * 00470 END
1.7.3 
