Return to zgbequ.f CVS log

Up to [local] / rpl / lapack / lapack

Commit initial.

    1:       SUBROUTINE ZGBEQU( M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND,
    2:      $                   AMAX, INFO )
    3: *
    4: *  -- LAPACK routine (version 3.2) --
    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    7: *     November 2006
    8: *
    9: *     .. Scalar Arguments ..
   10:       INTEGER            INFO, KL, KU, LDAB, M, N
   11:       DOUBLE PRECISION   AMAX, COLCND, ROWCND
   12: *     ..
   13: *     .. Array Arguments ..
   14:       DOUBLE PRECISION   C( * ), R( * )
   15:       COMPLEX*16         AB( LDAB, * )
   16: *     ..
   17: *
   18: *  Purpose
   19: *  =======
   20: *
   21: *  ZGBEQU computes row and column scalings intended to equilibrate an
   22: *  M-by-N band matrix A and reduce its condition number.  R returns the
   23: *  row scale factors and C the column scale factors, chosen to try to
   24: *  make the largest element in each row and column of the matrix B with
   25: *  elements B(i,j)=R(i)*A(i,j)*C(j) have absolute value 1.
   26: *
   27: *  R(i) and C(j) are restricted to be between SMLNUM = smallest safe
   28: *  number and BIGNUM = largest safe number.  Use of these scaling
   29: *  factors is not guaranteed to reduce the condition number of A but
   30: *  works well in practice.
   31: *
   32: *  Arguments
   33: *  =========
   34: *
   35: *  M       (input) INTEGER
   36: *          The number of rows of the matrix A.  M >= 0.
   37: *
   38: *  N       (input) INTEGER
   39: *          The number of columns of the matrix A.  N >= 0.
   40: *
   41: *  KL      (input) INTEGER
   42: *          The number of subdiagonals within the band of A.  KL >= 0.
   43: *
   44: *  KU      (input) INTEGER
   45: *          The number of superdiagonals within the band of A.  KU >= 0.
   46: *
   47: *  AB      (input) COMPLEX*16 array, dimension (LDAB,N)
   48: *          The band matrix A, stored in rows 1 to KL+KU+1.  The j-th
   49: *          column of A is stored in the j-th column of the array AB as
   50: *          follows:
   51: *          AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
   52: *
   53: *  LDAB    (input) INTEGER
   54: *          The leading dimension of the array AB.  LDAB >= KL+KU+1.
   55: *
   56: *  R       (output) DOUBLE PRECISION array, dimension (M)
   57: *          If INFO = 0, or INFO > M, R contains the row scale factors
   58: *          for A.
   59: *
   60: *  C       (output) DOUBLE PRECISION array, dimension (N)
   61: *          If INFO = 0, C contains the column scale factors for A.
   62: *
   63: *  ROWCND  (output) DOUBLE PRECISION
   64: *          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
   65: *          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
   66: *          AMAX is neither too large nor too small, it is not worth
   67: *          scaling by R.
   68: *
   69: *  COLCND  (output) DOUBLE PRECISION
   70: *          If INFO = 0, COLCND contains the ratio of the smallest
   71: *          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
   72: *          worth scaling by C.
   73: *
   74: *  AMAX    (output) DOUBLE PRECISION
   75: *          Absolute value of largest matrix element.  If AMAX is very
   76: *          close to overflow or very close to underflow, the matrix
   77: *          should be scaled.
   78: *
   79: *  INFO    (output) INTEGER
   80: *          = 0:  successful exit
   81: *          < 0:  if INFO = -i, the i-th argument had an illegal value
   82: *          > 0:  if INFO = i, and i is
   83: *                <= M:  the i-th row of A is exactly zero
   84: *                >  M:  the (i-M)-th column of A is exactly zero
   85: *
   86: *  =====================================================================
   87: *
   88: *     .. Parameters ..
   89:       DOUBLE PRECISION   ONE, ZERO
   90:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
   91: *     ..
   92: *     .. Local Scalars ..
   93:       INTEGER            I, J, KD
   94:       DOUBLE PRECISION   BIGNUM, RCMAX, RCMIN, SMLNUM
   95:       COMPLEX*16         ZDUM
   96: *     ..
   97: *     .. External Functions ..
   98:       DOUBLE PRECISION   DLAMCH
   99:       EXTERNAL           DLAMCH
  100: *     ..
  101: *     .. External Subroutines ..
  102:       EXTERNAL           XERBLA
  103: *     ..
  104: *     .. Intrinsic Functions ..
  105:       INTRINSIC          ABS, DBLE, DIMAG, MAX, MIN
  106: *     ..
  107: *     .. Statement Functions ..
  108:       DOUBLE PRECISION   CABS1
  109: *     ..
  110: *     .. Statement Function definitions ..
  111:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
  112: *     ..
  113: *     .. Executable Statements ..
  114: *
  115: *     Test the input parameters
  116: *
  117:       INFO = 0
  118:       IF( M.LT.0 ) THEN
  119:          INFO = -1
  120:       ELSE IF( N.LT.0 ) THEN
  121:          INFO = -2
  122:       ELSE IF( KL.LT.0 ) THEN
  123:          INFO = -3
  124:       ELSE IF( KU.LT.0 ) THEN
  125:          INFO = -4
  126:       ELSE IF( LDAB.LT.KL+KU+1 ) THEN
  127:          INFO = -6
  128:       END IF
  129:       IF( INFO.NE.0 ) THEN
  130:          CALL XERBLA( 'ZGBEQU', -INFO )
  131:          RETURN
  132:       END IF
  133: *
  134: *     Quick return if possible
  135: *
  136:       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
  137:          ROWCND = ONE
  138:          COLCND = ONE
  139:          AMAX = ZERO
  140:          RETURN
  141:       END IF
  142: *
  143: *     Get machine constants.
  144: *
  145:       SMLNUM = DLAMCH( 'S' )
  146:       BIGNUM = ONE / SMLNUM
  147: *
  148: *     Compute row scale factors.
  149: *
  150:       DO 10 I = 1, M
  151:          R( I ) = ZERO
  152:    10 CONTINUE
  153: *
  154: *     Find the maximum element in each row.
  155: *
  156:       KD = KU + 1
  157:       DO 30 J = 1, N
  158:          DO 20 I = MAX( J-KU, 1 ), MIN( J+KL, M )
  159:             R( I ) = MAX( R( I ), CABS1( AB( KD+I-J, J ) ) )
  160:    20    CONTINUE
  161:    30 CONTINUE
  162: *
  163: *     Find the maximum and minimum scale factors.
  164: *
  165:       RCMIN = BIGNUM
  166:       RCMAX = ZERO
  167:       DO 40 I = 1, M
  168:          RCMAX = MAX( RCMAX, R( I ) )
  169:          RCMIN = MIN( RCMIN, R( I ) )
  170:    40 CONTINUE
  171:       AMAX = RCMAX
  172: *
  173:       IF( RCMIN.EQ.ZERO ) THEN
  174: *
  175: *        Find the first zero scale factor and return an error code.
  176: *
  177:          DO 50 I = 1, M
  178:             IF( R( I ).EQ.ZERO ) THEN
  179:                INFO = I
  180:                RETURN
  181:             END IF
  182:    50    CONTINUE
  183:       ELSE
  184: *
  185: *        Invert the scale factors.
  186: *
  187:          DO 60 I = 1, M
  188:             R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM )
  189:    60    CONTINUE
  190: *
  191: *        Compute ROWCND = min(R(I)) / max(R(I))
  192: *
  193:          ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
  194:       END IF
  195: *
  196: *     Compute column scale factors
  197: *
  198:       DO 70 J = 1, N
  199:          C( J ) = ZERO
  200:    70 CONTINUE
  201: *
  202: *     Find the maximum element in each column,
  203: *     assuming the row scaling computed above.
  204: *
  205:       KD = KU + 1
  206:       DO 90 J = 1, N
  207:          DO 80 I = MAX( J-KU, 1 ), MIN( J+KL, M )
  208:             C( J ) = MAX( C( J ), CABS1( AB( KD+I-J, J ) )*R( I ) )
  209:    80    CONTINUE
  210:    90 CONTINUE
  211: *
  212: *     Find the maximum and minimum scale factors.
  213: *
  214:       RCMIN = BIGNUM
  215:       RCMAX = ZERO
  216:       DO 100 J = 1, N
  217:          RCMIN = MIN( RCMIN, C( J ) )
  218:          RCMAX = MAX( RCMAX, C( J ) )
  219:   100 CONTINUE
  220: *
  221:       IF( RCMIN.EQ.ZERO ) THEN
  222: *
  223: *        Find the first zero scale factor and return an error code.
  224: *
  225:          DO 110 J = 1, N
  226:             IF( C( J ).EQ.ZERO ) THEN
  227:                INFO = M + J
  228:                RETURN
  229:             END IF
  230:   110    CONTINUE
  231:       ELSE
  232: *
  233: *        Invert the scale factors.
  234: *
  235:          DO 120 J = 1, N
  236:             C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM )
  237:   120    CONTINUE
  238: *
  239: *        Compute COLCND = min(C(J)) / max(C(J))
  240: *
  241:          COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
  242:       END IF
  243: *
  244:       RETURN
  245: *
  246: *     End of ZGBEQU
  247: *
  248:       END