LAPACK: SRC/zggbal.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZGGBAL( 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       DOUBLE PRECISION   LSCALE( * ), RSCALE( * ), WORK( * )
00015       COMPLEX*16         A( LDA, * ), B( LDB, * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  ZGGBAL 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*16 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*16 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) 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) 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) 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       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       COMPLEX*16         CZERO
00113       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+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       DOUBLE PRECISION   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*16         CDUM
00123 *     ..
00124 *     .. External Functions ..
00125       LOGICAL            LSAME
00126       INTEGER            IZAMAX
00127       DOUBLE PRECISION   DDOT, DLAMCH
00128       EXTERNAL           LSAME, IZAMAX, DDOT, DLAMCH
00129 *     ..
00130 *     .. External Subroutines ..
00131       EXTERNAL           DAXPY, DSCAL, XERBLA, ZDSCAL, ZSWAP
00132 *     ..
00133 *     .. Intrinsic Functions ..
00134       INTRINSIC          ABS, DBLE, DIMAG, INT, LOG10, MAX, MIN, SIGN
00135 *     ..
00136 *     .. Statement Functions ..
00137       DOUBLE PRECISION   CABS1
00138 *     ..
00139 *     .. Statement Function definitions ..
00140       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( 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( 'ZGGBAL', -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 ) = 1
00205       LSCALE( 1 ) = 1
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 ZSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA )
00267       CALL ZSWAP( 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 ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
00276       CALL ZSWAP( 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 / DBLE( 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 = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) +
00347      $        DDOT( 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 DSCAL( NR, BETA, WORK( ILO ), 1 )
00365       CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 )
00366 *
00367       CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 )
00368       CALL DAXPY( 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 ) = DBLE( 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 ) = DBLE( KOUNT )*WORK( J ) + SUM
00409   330 CONTINUE
00410 *
00411       SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) +
00412      $      DDOT( 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 DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 )
00432       CALL DAXPY( 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 = DLAMCH( '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 = IZAMAX( N-ILO+1, A( I, ILO ), LDA )
00448          RAB = ABS( A( I, IRAB+ILO-1 ) )
00449          IRAB = IZAMAX( 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 = IZAMAX( IHI, A( 1, I ), 1 )
00456          CAB = ABS( A( ICAB, I ) )
00457          ICAB = IZAMAX( 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 ZDSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA )
00469          CALL ZDSCAL( 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 ZDSCAL( IHI, RSCALE( J ), A( 1, J ), 1 )
00476          CALL ZDSCAL( IHI, RSCALE( J ), B( 1, J ), 1 )
00477   380 CONTINUE
00478 *
00479       RETURN
00480 *
00481 *     End of ZGGBAL
00482 *
00483       END