SRC/cgebal.f

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