Return to zgebal.f CVS log

Up to [local] / rpl / lapack / lapack

Initial revision

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