SRC/zgebal.f

   1 *> \brief \b ZGEBAL
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZGEBAL + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgebal.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgebal.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgebal.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          JOB
  25 *       INTEGER            IHI, ILO, INFO, LDA, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       DOUBLE PRECISION   SCALE( * )
  29 *       COMPLEX*16         A( LDA, * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> ZGEBAL balances a general complex matrix A.  This involves, first,
  39 *> permuting A by a similarity transformation to isolate eigenvalues
  40 *> in the first 1 to ILO-1 and last IHI+1 to N elements on the
  41 *> diagonal; and second, applying a diagonal similarity transformation
  42 *> to rows and columns ILO to IHI to make the rows and columns as
  43 *> close in norm as possible.  Both steps are optional.
  44 *>
  45 *> Balancing may reduce the 1-norm of the matrix, and improve the
  46 *> accuracy of the computed eigenvalues and/or eigenvectors.
  47 *> \endverbatim
  48 *
  49 *  Arguments:
  50 *  ==========
  51 *
  52 *> \param[in] JOB
  53 *> \verbatim
  54 *>          JOB is CHARACTER*1
  55 *>          Specifies the operations to be performed on A:
  56 *>          = 'N':  none:  simply set ILO = 1, IHI = N, SCALE(I) = 1.0
  57 *>                  for i = 1,...,N;
  58 *>          = 'P':  permute only;
  59 *>          = 'S':  scale only;
  60 *>          = 'B':  both permute and scale.
  61 *> \endverbatim
  62 *>
  63 *> \param[in] N
  64 *> \verbatim
  65 *>          N is INTEGER
  66 *>          The order of the matrix A.  N >= 0.
  67 *> \endverbatim
  68 *>
  69 *> \param[in,out] A
  70 *> \verbatim
  71 *>          A is COMPLEX*16 array, dimension (LDA,N)
  72 *>          On entry, the input matrix A.
  73 *>          On exit,  A is overwritten by the balanced matrix.
  74 *>          If JOB = 'N', A is not referenced.
  75 *>          See Further Details.
  76 *> \endverbatim
  77 *>
  78 *> \param[in] LDA
  79 *> \verbatim
  80 *>          LDA is INTEGER
  81 *>          The leading dimension of the array A.  LDA >= max(1,N).
  82 *> \endverbatim
  83 *>
  84 *> \param[out] ILO
  85 *> \verbatim
  86 *> \endverbatim
  87 *>
  88 *> \param[out] IHI
  89 *> \verbatim
  90 *>          ILO and IHI are set to INTEGER such that on exit
  91 *>          A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
  92 *>          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
  93 *> \endverbatim
  94 *>
  95 *> \param[out] SCALE
  96 *> \verbatim
  97 *>          SCALE is DOUBLE PRECISION array, dimension (N)
  98 *>          Details of the permutations and scaling factors applied to
  99 *>          A.  If P(j) is the index of the row and column interchanged
 100 *>          with row and column j and D(j) is the scaling factor
 101 *>          applied to row and column j, then
 102 *>          SCALE(j) = P(j)    for j = 1,...,ILO-1
 103 *>                   = D(j)    for j = ILO,...,IHI
 104 *>                   = P(j)    for j = IHI+1,...,N.
 105 *>          The order in which the interchanges are made is N to IHI+1,
 106 *>          then 1 to ILO-1.
 107 *> \endverbatim
 108 *>
 109 *> \param[out] INFO
 110 *> \verbatim
 111 *>          INFO is INTEGER
 112 *>          = 0:  successful exit.
 113 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 114 *> \endverbatim
 115 *
 116 *  Authors:
 117 *  ========
 118 *
 119 *> \author Univ. of Tennessee
 120 *> \author Univ. of California Berkeley
 121 *> \author Univ. of Colorado Denver
 122 *> \author NAG Ltd.
 123 *
 124 *> \date November 2015
 125 *
 126 *> \ingroup complex16GEcomputational
 127 *
 128 *> \par Further Details:
 129 *  =====================
 130 *>
 131 *> \verbatim
 132 *>
 133 *>  The permutations consist of row and column interchanges which put
 134 *>  the matrix in the form
 135 *>
 136 *>             ( T1   X   Y  )
 137 *>     P A P = (  0   B   Z  )
 138 *>             (  0   0   T2 )
 139 *>
 140 *>  where T1 and T2 are upper triangular matrices whose eigenvalues lie
 141 *>  along the diagonal.  The column indices ILO and IHI mark the starting
 142 *>  and ending columns of the submatrix B. Balancing consists of applying
 143 *>  a diagonal similarity transformation inv(D) * B * D to make the
 144 *>  1-norms of each row of B and its corresponding column nearly equal.
 145 *>  The output matrix is
 146 *>
 147 *>     ( T1     X*D          Y    )
 148 *>     (  0  inv(D)*B*D  inv(D)*Z ).
 149 *>     (  0      0           T2   )
 150 *>
 151 *>  Information about the permutations P and the diagonal matrix D is
 152 *>  returned in the vector SCALE.
 153 *>
 154 *>  This subroutine is based on the EISPACK routine CBAL.
 155 *>
 156 *>  Modified by Tzu-Yi Chen, Computer Science Division, University of
 157 *>    California at Berkeley, USA
 158 *> \endverbatim
 159 *>
 160 *  =====================================================================
 161       SUBROUTINE ZGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
 162 *
 163 *  -- LAPACK computational routine (version 3.6.0) --
 164 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 165 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 166 *     November 2015
 167 *
 168 *     .. Scalar Arguments ..
 169       CHARACTER          JOB
 170       INTEGER            IHI, ILO, INFO, LDA, N
 171 *     ..
 172 *     .. Array Arguments ..
 173       DOUBLE PRECISION   SCALE( * )
 174       COMPLEX*16         A( LDA, * )
 175 *     ..
 176 *
 177 *  =====================================================================
 178 *
 179 *     .. Parameters ..
 180       DOUBLE PRECISION   ZERO, ONE
 181       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 182       DOUBLE PRECISION   SCLFAC
 183       PARAMETER          ( SCLFAC = 2.0D+0 )
 184       DOUBLE PRECISION   FACTOR
 185       PARAMETER          ( FACTOR = 0.95D+0 )
 186 *     ..
 187 *     .. Local Scalars ..
 188       LOGICAL            NOCONV
 189       INTEGER            I, ICA, IEXC, IRA, J, K, L, M
 190       DOUBLE PRECISION   C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
 191      $                   SFMIN2
 192       COMPLEX*16         CDUM
 193 *     ..
 194 *     .. External Functions ..
 195       LOGICAL            DISNAN, LSAME
 196       INTEGER            IZAMAX
 197       DOUBLE PRECISION   DLAMCH, DZNRM2
 198       EXTERNAL           DISNAN, LSAME, IZAMAX, DLAMCH, DZNRM2
 199 *     ..
 200 *     .. External Subroutines ..
 201       EXTERNAL           XERBLA, ZDSCAL, ZSWAP
 202 *     ..
 203 *     .. Intrinsic Functions ..
 204       INTRINSIC          ABS, DBLE, DIMAG, MAX, MIN
 205 *
 206 *     Test the input parameters
 207 *
 208       INFO = 0
 209       IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
 210      $    .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
 211          INFO = -1
 212       ELSE IF( N.LT.0 ) THEN
 213          INFO = -2
 214       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 215          INFO = -4
 216       END IF
 217       IF( INFO.NE.0 ) THEN
 218          CALL XERBLA( 'ZGEBAL', -INFO )
 219          RETURN
 220       END IF
 221 *
 222       K = 1
 223       L = N
 224 *
 225       IF( N.EQ.0 )
 226      $   GO TO 210
 227 *
 228       IF( LSAME( JOB, 'N' ) ) THEN
 229          DO 10 I = 1, N
 230             SCALE( I ) = ONE
 231    10    CONTINUE
 232          GO TO 210
 233       END IF
 234 *
 235       IF( LSAME( JOB, 'S' ) )
 236      $   GO TO 120
 237 *
 238 *     Permutation to isolate eigenvalues if possible
 239 *
 240       GO TO 50
 241 *
 242 *     Row and column exchange.
 243 *
 244    20 CONTINUE
 245       SCALE( M ) = J
 246       IF( J.EQ.M )
 247      $   GO TO 30
 248 *
 249       CALL ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
 250       CALL ZSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
 251 *
 252    30 CONTINUE
 253       GO TO ( 40, 80 )IEXC
 254 *
 255 *     Search for rows isolating an eigenvalue and push them down.
 256 *
 257    40 CONTINUE
 258       IF( L.EQ.1 )
 259      $   GO TO 210
 260       L = L - 1
 261 *
 262    50 CONTINUE
 263       DO 70 J = L, 1, -1
 264 *
 265          DO 60 I = 1, L
 266             IF( I.EQ.J )
 267      $         GO TO 60
 268             IF( DBLE( A( J, I ) ).NE.ZERO .OR. DIMAG( A( J, I ) ).NE.
 269      $          ZERO )GO TO 70
 270    60    CONTINUE
 271 *
 272          M = L
 273          IEXC = 1
 274          GO TO 20
 275    70 CONTINUE
 276 *
 277       GO TO 90
 278 *
 279 *     Search for columns isolating an eigenvalue and push them left.
 280 *
 281    80 CONTINUE
 282       K = K + 1
 283 *
 284    90 CONTINUE
 285       DO 110 J = K, L
 286 *
 287          DO 100 I = K, L
 288             IF( I.EQ.J )
 289      $         GO TO 100
 290             IF( DBLE( A( I, J ) ).NE.ZERO .OR. DIMAG( A( I, J ) ).NE.
 291      $          ZERO )GO TO 110
 292   100    CONTINUE
 293 *
 294          M = K
 295          IEXC = 2
 296          GO TO 20
 297   110 CONTINUE
 298 *
 299   120 CONTINUE
 300       DO 130 I = K, L
 301          SCALE( I ) = ONE
 302   130 CONTINUE
 303 *
 304       IF( LSAME( JOB, 'P' ) )
 305      $   GO TO 210
 306 *
 307 *     Balance the submatrix in rows K to L.
 308 *
 309 *     Iterative loop for norm reduction
 310 *
 311       SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
 312       SFMAX1 = ONE / SFMIN1
 313       SFMIN2 = SFMIN1*SCLFAC
 314       SFMAX2 = ONE / SFMIN2
 315   140 CONTINUE
 316       NOCONV = .FALSE.
 317 *
 318       DO 200 I = K, L
 319 *
 320          C = DZNRM2( L-K+1, A( K, I ), 1 )
 321          R = DZNRM2( L-K+1, A( I, K ), LDA )
 322          ICA = IZAMAX( L, A( 1, I ), 1 )
 323          CA = ABS( A( ICA, I ) )
 324          IRA = IZAMAX( N-K+1, A( I, K ), LDA )
 325          RA = ABS( A( I, IRA+K-1 ) )
 326 *
 327 *        Guard against zero C or R due to underflow.
 328 *
 329          IF( C.EQ.ZERO .OR. R.EQ.ZERO )
 330      $      GO TO 200
 331          G = R / SCLFAC
 332          F = ONE
 333          S = C + R
 334   160    CONTINUE
 335          IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
 336      $       MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
 337             IF( DISNAN( C+F+CA+R+G+RA ) ) THEN
 338 *
 339 *           Exit if NaN to avoid infinite loop
 340 *
 341             INFO = -3
 342             CALL XERBLA( 'ZGEBAL', -INFO )
 343             RETURN
 344          END IF
 345          F = F*SCLFAC
 346          C = C*SCLFAC
 347          CA = CA*SCLFAC
 348          R = R / SCLFAC
 349          G = G / SCLFAC
 350          RA = RA / SCLFAC
 351          GO TO 160
 352 *
 353   170    CONTINUE
 354          G = C / SCLFAC
 355   180    CONTINUE
 356          IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
 357      $       MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
 358          F = F / SCLFAC
 359          C = C / SCLFAC
 360          G = G / SCLFAC
 361          CA = CA / SCLFAC
 362          R = R*SCLFAC
 363          RA = RA*SCLFAC
 364          GO TO 180
 365 *
 366 *        Now balance.
 367 *
 368   190    CONTINUE
 369          IF( ( C+R ).GE.FACTOR*S )
 370      $      GO TO 200
 371          IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
 372             IF( F*SCALE( I ).LE.SFMIN1 )
 373      $         GO TO 200
 374          END IF
 375          IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
 376             IF( SCALE( I ).GE.SFMAX1 / F )
 377      $         GO TO 200
 378          END IF
 379          G = ONE / F
 380          SCALE( I ) = SCALE( I )*F
 381          NOCONV = .TRUE.
 382 *
 383          CALL ZDSCAL( N-K+1, G, A( I, K ), LDA )
 384          CALL ZDSCAL( L, F, A( 1, I ), 1 )
 385 *
 386   200 CONTINUE
 387 *
 388       IF( NOCONV )
 389      $   GO TO 140
 390 *
 391   210 CONTINUE
 392       ILO = K
 393       IHI = L
 394 *
 395       RETURN
 396 *
 397 *     End of ZGEBAL
 398 *
 399       END