Return to dgebal.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack vers la version 3.3.0.

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