SRC/zggbal.f

   1 *> \brief \b ZGGBAL
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZGGBAL + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggbal.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggbal.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggbal.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
  22 *                          RSCALE, WORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          JOB
  26 *       INTEGER            IHI, ILO, INFO, LDA, LDB, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   LSCALE( * ), RSCALE( * ), WORK( * )
  30 *       COMPLEX*16         A( LDA, * ), B( LDB, * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> ZGGBAL balances a pair of general complex matrices (A,B).  This
  40 *> involves, first, permuting A and B by similarity transformations to
  41 *> isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
  42 *> elements on the diagonal; and second, applying a diagonal similarity
  43 *> transformation to rows and columns ILO to IHI to make the rows
  44 *> and columns as close in norm as possible. Both steps are optional.
  45 *>
  46 *> Balancing may reduce the 1-norm of the matrices, and improve the
  47 *> accuracy of the computed eigenvalues and/or eigenvectors in the
  48 *> generalized eigenvalue problem A*x = lambda*B*x.
  49 *> \endverbatim
  50 *
  51 *  Arguments:
  52 *  ==========
  53 *
  54 *> \param[in] JOB
  55 *> \verbatim
  56 *>          JOB is CHARACTER*1
  57 *>          Specifies the operations to be performed on A and B:
  58 *>          = 'N':  none:  simply set ILO = 1, IHI = N, LSCALE(I) = 1.0
  59 *>                  and RSCALE(I) = 1.0 for i=1,...,N;
  60 *>          = 'P':  permute only;
  61 *>          = 'S':  scale only;
  62 *>          = 'B':  both permute and scale.
  63 *> \endverbatim
  64 *>
  65 *> \param[in] N
  66 *> \verbatim
  67 *>          N is INTEGER
  68 *>          The order of the matrices A and B.  N >= 0.
  69 *> \endverbatim
  70 *>
  71 *> \param[in,out] A
  72 *> \verbatim
  73 *>          A is COMPLEX*16 array, dimension (LDA,N)
  74 *>          On entry, the input matrix A.
  75 *>          On exit, A is overwritten by the balanced matrix.
  76 *>          If JOB = 'N', A is not referenced.
  77 *> \endverbatim
  78 *>
  79 *> \param[in] LDA
  80 *> \verbatim
  81 *>          LDA is INTEGER
  82 *>          The leading dimension of the array A. LDA >= max(1,N).
  83 *> \endverbatim
  84 *>
  85 *> \param[in,out] B
  86 *> \verbatim
  87 *>          B is COMPLEX*16 array, dimension (LDB,N)
  88 *>          On entry, the input matrix B.
  89 *>          On exit, B is overwritten by the balanced matrix.
  90 *>          If JOB = 'N', B is not referenced.
  91 *> \endverbatim
  92 *>
  93 *> \param[in] LDB
  94 *> \verbatim
  95 *>          LDB is INTEGER
  96 *>          The leading dimension of the array B. LDB >= max(1,N).
  97 *> \endverbatim
  98 *>
  99 *> \param[out] ILO
 100 *> \verbatim
 101 *>          ILO is INTEGER
 102 *> \endverbatim
 103 *>
 104 *> \param[out] IHI
 105 *> \verbatim
 106 *>          IHI is INTEGER
 107 *>          ILO and IHI are set to integers such that on exit
 108 *>          A(i,j) = 0 and B(i,j) = 0 if i > j and
 109 *>          j = 1,...,ILO-1 or i = IHI+1,...,N.
 110 *>          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
 111 *> \endverbatim
 112 *>
 113 *> \param[out] LSCALE
 114 *> \verbatim
 115 *>          LSCALE is DOUBLE PRECISION array, dimension (N)
 116 *>          Details of the permutations and scaling factors applied
 117 *>          to the left side of A and B.  If P(j) is the index of the
 118 *>          row interchanged with row j, and D(j) is the scaling factor
 119 *>          applied to row j, then
 120 *>            LSCALE(j) = P(j)    for J = 1,...,ILO-1
 121 *>                      = D(j)    for J = ILO,...,IHI
 122 *>                      = P(j)    for J = IHI+1,...,N.
 123 *>          The order in which the interchanges are made is N to IHI+1,
 124 *>          then 1 to ILO-1.
 125 *> \endverbatim
 126 *>
 127 *> \param[out] RSCALE
 128 *> \verbatim
 129 *>          RSCALE is DOUBLE PRECISION array, dimension (N)
 130 *>          Details of the permutations and scaling factors applied
 131 *>          to the right side of A and B.  If P(j) is the index of the
 132 *>          column interchanged with column j, and D(j) is the scaling
 133 *>          factor applied to column j, then
 134 *>            RSCALE(j) = P(j)    for J = 1,...,ILO-1
 135 *>                      = D(j)    for J = ILO,...,IHI
 136 *>                      = P(j)    for J = IHI+1,...,N.
 137 *>          The order in which the interchanges are made is N to IHI+1,
 138 *>          then 1 to ILO-1.
 139 *> \endverbatim
 140 *>
 141 *> \param[out] WORK
 142 *> \verbatim
 143 *>          WORK is DOUBLE PRECISION array, dimension (lwork)
 144 *>          lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and
 145 *>          at least 1 when JOB = 'N' or 'P'.
 146 *> \endverbatim
 147 *>
 148 *> \param[out] INFO
 149 *> \verbatim
 150 *>          INFO is INTEGER
 151 *>          = 0:  successful exit
 152 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 153 *> \endverbatim
 154 *
 155 *  Authors:
 156 *  ========
 157 *
 158 *> \author Univ. of Tennessee
 159 *> \author Univ. of California Berkeley
 160 *> \author Univ. of Colorado Denver
 161 *> \author NAG Ltd.
 162 *
 163 *> \date June 2016
 164 *
 165 *> \ingroup complex16GBcomputational
 166 *
 167 *> \par Further Details:
 168 *  =====================
 169 *>
 170 *> \verbatim
 171 *>
 172 *>  See R.C. WARD, Balancing the generalized eigenvalue problem,
 173 *>                 SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
 174 *> \endverbatim
 175 *>
 176 *  =====================================================================
 177       SUBROUTINE ZGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
 178      $                   RSCALE, WORK, INFO )
 179 *
 180 *  -- LAPACK computational routine (version 3.6.1) --
 181 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 182 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 183 *     June 2016
 184 *
 185 *     .. Scalar Arguments ..
 186       CHARACTER          JOB
 187       INTEGER            IHI, ILO, INFO, LDA, LDB, N
 188 *     ..
 189 *     .. Array Arguments ..
 190       DOUBLE PRECISION   LSCALE( * ), RSCALE( * ), WORK( * )
 191       COMPLEX*16         A( LDA, * ), B( LDB, * )
 192 *     ..
 193 *
 194 *  =====================================================================
 195 *
 196 *     .. Parameters ..
 197       DOUBLE PRECISION   ZERO, HALF, ONE
 198       PARAMETER          ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 )
 199       DOUBLE PRECISION   THREE, SCLFAC
 200       PARAMETER          ( THREE = 3.0D+0, SCLFAC = 1.0D+1 )
 201       COMPLEX*16         CZERO
 202       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ) )
 203 *     ..
 204 *     .. Local Scalars ..
 205       INTEGER            I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
 206      $                   K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN,
 207      $                   M, NR, NRP2
 208       DOUBLE PRECISION   ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
 209      $                   COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX,
 210      $                   SFMIN, SUM, T, TA, TB, TC
 211       COMPLEX*16         CDUM
 212 *     ..
 213 *     .. External Functions ..
 214       LOGICAL            LSAME
 215       INTEGER            IZAMAX
 216       DOUBLE PRECISION   DDOT, DLAMCH
 217       EXTERNAL           LSAME, IZAMAX, DDOT, DLAMCH
 218 *     ..
 219 *     .. External Subroutines ..
 220       EXTERNAL           DAXPY, DSCAL, XERBLA, ZDSCAL, ZSWAP
 221 *     ..
 222 *     .. Intrinsic Functions ..
 223       INTRINSIC          ABS, DBLE, DIMAG, INT, LOG10, MAX, MIN, SIGN
 224 *     ..
 225 *     .. Statement Functions ..
 226       DOUBLE PRECISION   CABS1
 227 *     ..
 228 *     .. Statement Function definitions ..
 229       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
 230 *     ..
 231 *     .. Executable Statements ..
 232 *
 233 *     Test the input parameters
 234 *
 235       INFO = 0
 236       IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
 237      $    .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
 238          INFO = -1
 239       ELSE IF( N.LT.0 ) THEN
 240          INFO = -2
 241       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 242          INFO = -4
 243       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 244          INFO = -6
 245       END IF
 246       IF( INFO.NE.0 ) THEN
 247          CALL XERBLA( 'ZGGBAL', -INFO )
 248          RETURN
 249       END IF
 250 *
 251 *     Quick return if possible
 252 *
 253       IF( N.EQ.0 ) THEN
 254          ILO = 1
 255          IHI = N
 256          RETURN
 257       END IF
 258 *
 259       IF( N.EQ.1 ) THEN
 260          ILO = 1
 261          IHI = N
 262          LSCALE( 1 ) = ONE
 263          RSCALE( 1 ) = ONE
 264          RETURN
 265       END IF
 266 *
 267       IF( LSAME( JOB, 'N' ) ) THEN
 268          ILO = 1
 269          IHI = N
 270          DO 10 I = 1, N
 271             LSCALE( I ) = ONE
 272             RSCALE( I ) = ONE
 273    10    CONTINUE
 274          RETURN
 275       END IF
 276 *
 277       K = 1
 278       L = N
 279       IF( LSAME( JOB, 'S' ) )
 280      $   GO TO 190
 281 *
 282       GO TO 30
 283 *
 284 *     Permute the matrices A and B to isolate the eigenvalues.
 285 *
 286 *     Find row with one nonzero in columns 1 through L
 287 *
 288    20 CONTINUE
 289       L = LM1
 290       IF( L.NE.1 )
 291      $   GO TO 30
 292 *
 293       RSCALE( 1 ) = 1
 294       LSCALE( 1 ) = 1
 295       GO TO 190
 296 *
 297    30 CONTINUE
 298       LM1 = L - 1
 299       DO 80 I = L, 1, -1
 300          DO 40 J = 1, LM1
 301             JP1 = J + 1
 302             IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
 303      $         GO TO 50
 304    40    CONTINUE
 305          J = L
 306          GO TO 70
 307 *
 308    50    CONTINUE
 309          DO 60 J = JP1, L
 310             IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
 311      $         GO TO 80
 312    60    CONTINUE
 313          J = JP1 - 1
 314 *
 315    70    CONTINUE
 316          M = L
 317          IFLOW = 1
 318          GO TO 160
 319    80 CONTINUE
 320       GO TO 100
 321 *
 322 *     Find column with one nonzero in rows K through N
 323 *
 324    90 CONTINUE
 325       K = K + 1
 326 *
 327   100 CONTINUE
 328       DO 150 J = K, L
 329          DO 110 I = K, LM1
 330             IP1 = I + 1
 331             IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
 332      $         GO TO 120
 333   110    CONTINUE
 334          I = L
 335          GO TO 140
 336   120    CONTINUE
 337          DO 130 I = IP1, L
 338             IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
 339      $         GO TO 150
 340   130    CONTINUE
 341          I = IP1 - 1
 342   140    CONTINUE
 343          M = K
 344          IFLOW = 2
 345          GO TO 160
 346   150 CONTINUE
 347       GO TO 190
 348 *
 349 *     Permute rows M and I
 350 *
 351   160 CONTINUE
 352       LSCALE( M ) = I
 353       IF( I.EQ.M )
 354      $   GO TO 170
 355       CALL ZSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA )
 356       CALL ZSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB )
 357 *
 358 *     Permute columns M and J
 359 *
 360   170 CONTINUE
 361       RSCALE( M ) = J
 362       IF( J.EQ.M )
 363      $   GO TO 180
 364       CALL ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
 365       CALL ZSWAP( L, B( 1, J ), 1, B( 1, M ), 1 )
 366 *
 367   180 CONTINUE
 368       GO TO ( 20, 90 )IFLOW
 369 *
 370   190 CONTINUE
 371       ILO = K
 372       IHI = L
 373 *
 374       IF( LSAME( JOB, 'P' ) ) THEN
 375          DO 195 I = ILO, IHI
 376             LSCALE( I ) = ONE
 377             RSCALE( I ) = ONE
 378   195    CONTINUE
 379          RETURN
 380       END IF
 381 *
 382       IF( ILO.EQ.IHI )
 383      $   RETURN
 384 *
 385 *     Balance the submatrix in rows ILO to IHI.
 386 *
 387       NR = IHI - ILO + 1
 388       DO 200 I = ILO, IHI
 389          RSCALE( I ) = ZERO
 390          LSCALE( I ) = ZERO
 391 *
 392          WORK( I ) = ZERO
 393          WORK( I+N ) = ZERO
 394          WORK( I+2*N ) = ZERO
 395          WORK( I+3*N ) = ZERO
 396          WORK( I+4*N ) = ZERO
 397          WORK( I+5*N ) = ZERO
 398   200 CONTINUE
 399 *
 400 *     Compute right side vector in resulting linear equations
 401 *
 402       BASL = LOG10( SCLFAC )
 403       DO 240 I = ILO, IHI
 404          DO 230 J = ILO, IHI
 405             IF( A( I, J ).EQ.CZERO ) THEN
 406                TA = ZERO
 407                GO TO 210
 408             END IF
 409             TA = LOG10( CABS1( A( I, J ) ) ) / BASL
 410 *
 411   210       CONTINUE
 412             IF( B( I, J ).EQ.CZERO ) THEN
 413                TB = ZERO
 414                GO TO 220
 415             END IF
 416             TB = LOG10( CABS1( B( I, J ) ) ) / BASL
 417 *
 418   220       CONTINUE
 419             WORK( I+4*N ) = WORK( I+4*N ) - TA - TB
 420             WORK( J+5*N ) = WORK( J+5*N ) - TA - TB
 421   230    CONTINUE
 422   240 CONTINUE
 423 *
 424       COEF = ONE / DBLE( 2*NR )
 425       COEF2 = COEF*COEF
 426       COEF5 = HALF*COEF2
 427       NRP2 = NR + 2
 428       BETA = ZERO
 429       IT = 1
 430 *
 431 *     Start generalized conjugate gradient iteration
 432 *
 433   250 CONTINUE
 434 *
 435       GAMMA = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) +
 436      $        DDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 )
 437 *
 438       EW = ZERO
 439       EWC = ZERO
 440       DO 260 I = ILO, IHI
 441          EW = EW + WORK( I+4*N )
 442          EWC = EWC + WORK( I+5*N )
 443   260 CONTINUE
 444 *
 445       GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2
 446       IF( GAMMA.EQ.ZERO )
 447      $   GO TO 350
 448       IF( IT.NE.1 )
 449      $   BETA = GAMMA / PGAMMA
 450       T = COEF5*( EWC-THREE*EW )
 451       TC = COEF5*( EW-THREE*EWC )
 452 *
 453       CALL DSCAL( NR, BETA, WORK( ILO ), 1 )
 454       CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 )
 455 *
 456       CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 )
 457       CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 )
 458 *
 459       DO 270 I = ILO, IHI
 460          WORK( I ) = WORK( I ) + TC
 461          WORK( I+N ) = WORK( I+N ) + T
 462   270 CONTINUE
 463 *
 464 *     Apply matrix to vector
 465 *
 466       DO 300 I = ILO, IHI
 467          KOUNT = 0
 468          SUM = ZERO
 469          DO 290 J = ILO, IHI
 470             IF( A( I, J ).EQ.CZERO )
 471      $         GO TO 280
 472             KOUNT = KOUNT + 1
 473             SUM = SUM + WORK( J )
 474   280       CONTINUE
 475             IF( B( I, J ).EQ.CZERO )
 476      $         GO TO 290
 477             KOUNT = KOUNT + 1
 478             SUM = SUM + WORK( J )
 479   290    CONTINUE
 480          WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM
 481   300 CONTINUE
 482 *
 483       DO 330 J = ILO, IHI
 484          KOUNT = 0
 485          SUM = ZERO
 486          DO 320 I = ILO, IHI
 487             IF( A( I, J ).EQ.CZERO )
 488      $         GO TO 310
 489             KOUNT = KOUNT + 1
 490             SUM = SUM + WORK( I+N )
 491   310       CONTINUE
 492             IF( B( I, J ).EQ.CZERO )
 493      $         GO TO 320
 494             KOUNT = KOUNT + 1
 495             SUM = SUM + WORK( I+N )
 496   320    CONTINUE
 497          WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM
 498   330 CONTINUE
 499 *
 500       SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) +
 501      $      DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 )
 502       ALPHA = GAMMA / SUM
 503 *
 504 *     Determine correction to current iteration
 505 *
 506       CMAX = ZERO
 507       DO 340 I = ILO, IHI
 508          COR = ALPHA*WORK( I+N )
 509          IF( ABS( COR ).GT.CMAX )
 510      $      CMAX = ABS( COR )
 511          LSCALE( I ) = LSCALE( I ) + COR
 512          COR = ALPHA*WORK( I )
 513          IF( ABS( COR ).GT.CMAX )
 514      $      CMAX = ABS( COR )
 515          RSCALE( I ) = RSCALE( I ) + COR
 516   340 CONTINUE
 517       IF( CMAX.LT.HALF )
 518      $   GO TO 350
 519 *
 520       CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 )
 521       CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 )
 522 *
 523       PGAMMA = GAMMA
 524       IT = IT + 1
 525       IF( IT.LE.NRP2 )
 526      $   GO TO 250
 527 *
 528 *     End generalized conjugate gradient iteration
 529 *
 530   350 CONTINUE
 531       SFMIN = DLAMCH( 'S' )
 532       SFMAX = ONE / SFMIN
 533       LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE )
 534       LSFMAX = INT( LOG10( SFMAX ) / BASL )
 535       DO 360 I = ILO, IHI
 536          IRAB = IZAMAX( N-ILO+1, A( I, ILO ), LDA )
 537          RAB = ABS( A( I, IRAB+ILO-1 ) )
 538          IRAB = IZAMAX( N-ILO+1, B( I, ILO ), LDB )
 539          RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) )
 540          LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE )
 541          IR = INT(LSCALE( I ) + SIGN( HALF, LSCALE( I ) ))
 542          IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB )
 543          LSCALE( I ) = SCLFAC**IR
 544          ICAB = IZAMAX( IHI, A( 1, I ), 1 )
 545          CAB = ABS( A( ICAB, I ) )
 546          ICAB = IZAMAX( IHI, B( 1, I ), 1 )
 547          CAB = MAX( CAB, ABS( B( ICAB, I ) ) )
 548          LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE )
 549          JC = INT(RSCALE( I ) + SIGN( HALF, RSCALE( I ) ))
 550          JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB )
 551          RSCALE( I ) = SCLFAC**JC
 552   360 CONTINUE
 553 *
 554 *     Row scaling of matrices A and B
 555 *
 556       DO 370 I = ILO, IHI
 557          CALL ZDSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA )
 558          CALL ZDSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB )
 559   370 CONTINUE
 560 *
 561 *     Column scaling of matrices A and B
 562 *
 563       DO 380 J = ILO, IHI
 564          CALL ZDSCAL( IHI, RSCALE( J ), A( 1, J ), 1 )
 565          CALL ZDSCAL( IHI, RSCALE( J ), B( 1, J ), 1 )
 566   380 CONTINUE
 567 *
 568       RETURN
 569 *
 570 *     End of ZGGBAL
 571 *
 572       END