LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cggbal.f
Go to the documentation of this file.
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
 All Files Functions