![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DGGBAL 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DGGBAL + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggbal.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggbal.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggbal.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DGGBAL( 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 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), LSCALE( * ), 00030 * $ RSCALE( * ), WORK( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> DGGBAL balances a pair of general real 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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) 00119 *> is the scaling factor 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 DOUBLE PRECISION 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) 00133 *> is the scaling factor applied to column j, then 00134 *> LSCALE(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 DOUBLE PRECISION 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 doubleGBcomputational 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 DGGBAL( 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 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), LSCALE( * ), 00191 $ RSCALE( * ), WORK( * ) 00192 * .. 00193 * 00194 * ===================================================================== 00195 * 00196 * .. Parameters .. 00197 DOUBLE PRECISION ZERO, HALF, ONE 00198 PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 ) 00199 DOUBLE PRECISION THREE, SCLFAC 00200 PARAMETER ( THREE = 3.0D+0, SCLFAC = 1.0D+1 ) 00201 * .. 00202 * .. Local Scalars .. 00203 INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1, 00204 $ K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN, 00205 $ M, NR, NRP2 00206 DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2, 00207 $ COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX, 00208 $ SFMIN, SUM, T, TA, TB, TC 00209 * .. 00210 * .. External Functions .. 00211 LOGICAL LSAME 00212 INTEGER IDAMAX 00213 DOUBLE PRECISION DDOT, DLAMCH 00214 EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH 00215 * .. 00216 * .. External Subroutines .. 00217 EXTERNAL DAXPY, DSCAL, DSWAP, XERBLA 00218 * .. 00219 * .. Intrinsic Functions .. 00220 INTRINSIC ABS, DBLE, INT, LOG10, MAX, MIN, SIGN 00221 * .. 00222 * .. Executable Statements .. 00223 * 00224 * Test the input parameters 00225 * 00226 INFO = 0 00227 IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND. 00228 $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN 00229 INFO = -1 00230 ELSE IF( N.LT.0 ) THEN 00231 INFO = -2 00232 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00233 INFO = -4 00234 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00235 INFO = -6 00236 END IF 00237 IF( INFO.NE.0 ) THEN 00238 CALL XERBLA( 'DGGBAL', -INFO ) 00239 RETURN 00240 END IF 00241 * 00242 * Quick return if possible 00243 * 00244 IF( N.EQ.0 ) THEN 00245 ILO = 1 00246 IHI = N 00247 RETURN 00248 END IF 00249 * 00250 IF( N.EQ.1 ) THEN 00251 ILO = 1 00252 IHI = N 00253 LSCALE( 1 ) = ONE 00254 RSCALE( 1 ) = ONE 00255 RETURN 00256 END IF 00257 * 00258 IF( LSAME( JOB, 'N' ) ) THEN 00259 ILO = 1 00260 IHI = N 00261 DO 10 I = 1, N 00262 LSCALE( I ) = ONE 00263 RSCALE( I ) = ONE 00264 10 CONTINUE 00265 RETURN 00266 END IF 00267 * 00268 K = 1 00269 L = N 00270 IF( LSAME( JOB, 'S' ) ) 00271 $ GO TO 190 00272 * 00273 GO TO 30 00274 * 00275 * Permute the matrices A and B to isolate the eigenvalues. 00276 * 00277 * Find row with one nonzero in columns 1 through L 00278 * 00279 20 CONTINUE 00280 L = LM1 00281 IF( L.NE.1 ) 00282 $ GO TO 30 00283 * 00284 RSCALE( 1 ) = ONE 00285 LSCALE( 1 ) = ONE 00286 GO TO 190 00287 * 00288 30 CONTINUE 00289 LM1 = L - 1 00290 DO 80 I = L, 1, -1 00291 DO 40 J = 1, LM1 00292 JP1 = J + 1 00293 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) 00294 $ GO TO 50 00295 40 CONTINUE 00296 J = L 00297 GO TO 70 00298 * 00299 50 CONTINUE 00300 DO 60 J = JP1, L 00301 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) 00302 $ GO TO 80 00303 60 CONTINUE 00304 J = JP1 - 1 00305 * 00306 70 CONTINUE 00307 M = L 00308 IFLOW = 1 00309 GO TO 160 00310 80 CONTINUE 00311 GO TO 100 00312 * 00313 * Find column with one nonzero in rows K through N 00314 * 00315 90 CONTINUE 00316 K = K + 1 00317 * 00318 100 CONTINUE 00319 DO 150 J = K, L 00320 DO 110 I = K, LM1 00321 IP1 = I + 1 00322 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) 00323 $ GO TO 120 00324 110 CONTINUE 00325 I = L 00326 GO TO 140 00327 120 CONTINUE 00328 DO 130 I = IP1, L 00329 IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO ) 00330 $ GO TO 150 00331 130 CONTINUE 00332 I = IP1 - 1 00333 140 CONTINUE 00334 M = K 00335 IFLOW = 2 00336 GO TO 160 00337 150 CONTINUE 00338 GO TO 190 00339 * 00340 * Permute rows M and I 00341 * 00342 160 CONTINUE 00343 LSCALE( M ) = I 00344 IF( I.EQ.M ) 00345 $ GO TO 170 00346 CALL DSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA ) 00347 CALL DSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB ) 00348 * 00349 * Permute columns M and J 00350 * 00351 170 CONTINUE 00352 RSCALE( M ) = J 00353 IF( J.EQ.M ) 00354 $ GO TO 180 00355 CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 ) 00356 CALL DSWAP( L, B( 1, J ), 1, B( 1, M ), 1 ) 00357 * 00358 180 CONTINUE 00359 GO TO ( 20, 90 )IFLOW 00360 * 00361 190 CONTINUE 00362 ILO = K 00363 IHI = L 00364 * 00365 IF( LSAME( JOB, 'P' ) ) THEN 00366 DO 195 I = ILO, IHI 00367 LSCALE( I ) = ONE 00368 RSCALE( I ) = ONE 00369 195 CONTINUE 00370 RETURN 00371 END IF 00372 * 00373 IF( ILO.EQ.IHI ) 00374 $ RETURN 00375 * 00376 * Balance the submatrix in rows ILO to IHI. 00377 * 00378 NR = IHI - ILO + 1 00379 DO 200 I = ILO, IHI 00380 RSCALE( I ) = ZERO 00381 LSCALE( I ) = ZERO 00382 * 00383 WORK( I ) = ZERO 00384 WORK( I+N ) = ZERO 00385 WORK( I+2*N ) = ZERO 00386 WORK( I+3*N ) = ZERO 00387 WORK( I+4*N ) = ZERO 00388 WORK( I+5*N ) = ZERO 00389 200 CONTINUE 00390 * 00391 * Compute right side vector in resulting linear equations 00392 * 00393 BASL = LOG10( SCLFAC ) 00394 DO 240 I = ILO, IHI 00395 DO 230 J = ILO, IHI 00396 TB = B( I, J ) 00397 TA = A( I, J ) 00398 IF( TA.EQ.ZERO ) 00399 $ GO TO 210 00400 TA = LOG10( ABS( TA ) ) / BASL 00401 210 CONTINUE 00402 IF( TB.EQ.ZERO ) 00403 $ GO TO 220 00404 TB = LOG10( ABS( TB ) ) / BASL 00405 220 CONTINUE 00406 WORK( I+4*N ) = WORK( I+4*N ) - TA - TB 00407 WORK( J+5*N ) = WORK( J+5*N ) - TA - TB 00408 230 CONTINUE 00409 240 CONTINUE 00410 * 00411 COEF = ONE / DBLE( 2*NR ) 00412 COEF2 = COEF*COEF 00413 COEF5 = HALF*COEF2 00414 NRP2 = NR + 2 00415 BETA = ZERO 00416 IT = 1 00417 * 00418 * Start generalized conjugate gradient iteration 00419 * 00420 250 CONTINUE 00421 * 00422 GAMMA = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) + 00423 $ DDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 ) 00424 * 00425 EW = ZERO 00426 EWC = ZERO 00427 DO 260 I = ILO, IHI 00428 EW = EW + WORK( I+4*N ) 00429 EWC = EWC + WORK( I+5*N ) 00430 260 CONTINUE 00431 * 00432 GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2 00433 IF( GAMMA.EQ.ZERO ) 00434 $ GO TO 350 00435 IF( IT.NE.1 ) 00436 $ BETA = GAMMA / PGAMMA 00437 T = COEF5*( EWC-THREE*EW ) 00438 TC = COEF5*( EW-THREE*EWC ) 00439 * 00440 CALL DSCAL( NR, BETA, WORK( ILO ), 1 ) 00441 CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 ) 00442 * 00443 CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 ) 00444 CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 ) 00445 * 00446 DO 270 I = ILO, IHI 00447 WORK( I ) = WORK( I ) + TC 00448 WORK( I+N ) = WORK( I+N ) + T 00449 270 CONTINUE 00450 * 00451 * Apply matrix to vector 00452 * 00453 DO 300 I = ILO, IHI 00454 KOUNT = 0 00455 SUM = ZERO 00456 DO 290 J = ILO, IHI 00457 IF( A( I, J ).EQ.ZERO ) 00458 $ GO TO 280 00459 KOUNT = KOUNT + 1 00460 SUM = SUM + WORK( J ) 00461 280 CONTINUE 00462 IF( B( I, J ).EQ.ZERO ) 00463 $ GO TO 290 00464 KOUNT = KOUNT + 1 00465 SUM = SUM + WORK( J ) 00466 290 CONTINUE 00467 WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM 00468 300 CONTINUE 00469 * 00470 DO 330 J = ILO, IHI 00471 KOUNT = 0 00472 SUM = ZERO 00473 DO 320 I = ILO, IHI 00474 IF( A( I, J ).EQ.ZERO ) 00475 $ GO TO 310 00476 KOUNT = KOUNT + 1 00477 SUM = SUM + WORK( I+N ) 00478 310 CONTINUE 00479 IF( B( I, J ).EQ.ZERO ) 00480 $ GO TO 320 00481 KOUNT = KOUNT + 1 00482 SUM = SUM + WORK( I+N ) 00483 320 CONTINUE 00484 WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM 00485 330 CONTINUE 00486 * 00487 SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) + 00488 $ DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 ) 00489 ALPHA = GAMMA / SUM 00490 * 00491 * Determine correction to current iteration 00492 * 00493 CMAX = ZERO 00494 DO 340 I = ILO, IHI 00495 COR = ALPHA*WORK( I+N ) 00496 IF( ABS( COR ).GT.CMAX ) 00497 $ CMAX = ABS( COR ) 00498 LSCALE( I ) = LSCALE( I ) + COR 00499 COR = ALPHA*WORK( I ) 00500 IF( ABS( COR ).GT.CMAX ) 00501 $ CMAX = ABS( COR ) 00502 RSCALE( I ) = RSCALE( I ) + COR 00503 340 CONTINUE 00504 IF( CMAX.LT.HALF ) 00505 $ GO TO 350 00506 * 00507 CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 ) 00508 CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 ) 00509 * 00510 PGAMMA = GAMMA 00511 IT = IT + 1 00512 IF( IT.LE.NRP2 ) 00513 $ GO TO 250 00514 * 00515 * End generalized conjugate gradient iteration 00516 * 00517 350 CONTINUE 00518 SFMIN = DLAMCH( 'S' ) 00519 SFMAX = ONE / SFMIN 00520 LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE ) 00521 LSFMAX = INT( LOG10( SFMAX ) / BASL ) 00522 DO 360 I = ILO, IHI 00523 IRAB = IDAMAX( N-ILO+1, A( I, ILO ), LDA ) 00524 RAB = ABS( A( I, IRAB+ILO-1 ) ) 00525 IRAB = IDAMAX( N-ILO+1, B( I, ILO ), LDB ) 00526 RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) ) 00527 LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE ) 00528 IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) ) 00529 IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB ) 00530 LSCALE( I ) = SCLFAC**IR 00531 ICAB = IDAMAX( IHI, A( 1, I ), 1 ) 00532 CAB = ABS( A( ICAB, I ) ) 00533 ICAB = IDAMAX( IHI, B( 1, I ), 1 ) 00534 CAB = MAX( CAB, ABS( B( ICAB, I ) ) ) 00535 LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE ) 00536 JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) ) 00537 JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB ) 00538 RSCALE( I ) = SCLFAC**JC 00539 360 CONTINUE 00540 * 00541 * Row scaling of matrices A and B 00542 * 00543 DO 370 I = ILO, IHI 00544 CALL DSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA ) 00545 CALL DSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB ) 00546 370 CONTINUE 00547 * 00548 * Column scaling of matrices A and B 00549 * 00550 DO 380 J = ILO, IHI 00551 CALL DSCAL( IHI, RSCALE( J ), A( 1, J ), 1 ) 00552 CALL DSCAL( IHI, RSCALE( J ), B( 1, J ), 1 ) 00553 380 CONTINUE 00554 * 00555 RETURN 00556 * 00557 * End of DGGBAL 00558 * 00559 END