Return to zgeequb.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence Lapack.

    1: *> \brief \b ZGEEQUB
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download ZGEEQUB + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgeequb.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgeequb.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgeequb.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZGEEQUB( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
   22: *                           INFO )
   23: * 
   24: *       .. Scalar Arguments ..
   25: *       INTEGER            INFO, LDA, M, N
   26: *       DOUBLE PRECISION   AMAX, COLCND, ROWCND
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       DOUBLE PRECISION   C( * ), R( * )
   30: *       COMPLEX*16         A( LDA, * )
   31: *       ..
   32: *  
   33: *
   34: *> \par Purpose:
   35: *  =============
   36: *>
   37: *> \verbatim
   38: *>
   39: *> ZGEEQUB computes row and column scalings intended to equilibrate an
   40: *> M-by-N matrix A and reduce its condition number.  R returns the row
   41: *> scale factors and C the column scale factors, chosen to try to make
   42: *> the largest element in each row and column of the matrix B with
   43: *> elements B(i,j)=R(i)*A(i,j)*C(j) have an absolute value of at most
   44: *> the radix.
   45: *>
   46: *> R(i) and C(j) are restricted to be a power of the radix between
   47: *> SMLNUM = smallest safe number and BIGNUM = largest safe number.  Use
   48: *> of these scaling factors is not guaranteed to reduce the condition
   49: *> number of A but works well in practice.
   50: *>
   51: *> This routine differs from ZGEEQU by restricting the scaling factors
   52: *> to a power of the radix.  Baring over- and underflow, scaling by
   53: *> these factors introduces no additional rounding errors.  However, the
   54: *> scaled entries' magnitured are no longer approximately 1 but lie
   55: *> between sqrt(radix) and 1/sqrt(radix).
   56: *> \endverbatim
   57: *
   58: *  Arguments:
   59: *  ==========
   60: *
   61: *> \param[in] M
   62: *> \verbatim
   63: *>          M is INTEGER
   64: *>          The number of rows of the matrix A.  M >= 0.
   65: *> \endverbatim
   66: *>
   67: *> \param[in] N
   68: *> \verbatim
   69: *>          N is INTEGER
   70: *>          The number of columns of the matrix A.  N >= 0.
   71: *> \endverbatim
   72: *>
   73: *> \param[in] A
   74: *> \verbatim
   75: *>          A is COMPLEX*16 array, dimension (LDA,N)
   76: *>          The M-by-N matrix whose equilibration factors are
   77: *>          to be computed.
   78: *> \endverbatim
   79: *>
   80: *> \param[in] LDA
   81: *> \verbatim
   82: *>          LDA is INTEGER
   83: *>          The leading dimension of the array A.  LDA >= max(1,M).
   84: *> \endverbatim
   85: *>
   86: *> \param[out] R
   87: *> \verbatim
   88: *>          R is DOUBLE PRECISION array, dimension (M)
   89: *>          If INFO = 0 or INFO > M, R contains the row scale factors
   90: *>          for A.
   91: *> \endverbatim
   92: *>
   93: *> \param[out] C
   94: *> \verbatim
   95: *>          C is DOUBLE PRECISION array, dimension (N)
   96: *>          If INFO = 0,  C contains the column scale factors for A.
   97: *> \endverbatim
   98: *>
   99: *> \param[out] ROWCND
  100: *> \verbatim
  101: *>          ROWCND is DOUBLE PRECISION
  102: *>          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
  103: *>          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
  104: *>          AMAX is neither too large nor too small, it is not worth
  105: *>          scaling by R.
  106: *> \endverbatim
  107: *>
  108: *> \param[out] COLCND
  109: *> \verbatim
  110: *>          COLCND is DOUBLE PRECISION
  111: *>          If INFO = 0, COLCND contains the ratio of the smallest
  112: *>          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
  113: *>          worth scaling by C.
  114: *> \endverbatim
  115: *>
  116: *> \param[out] AMAX
  117: *> \verbatim
  118: *>          AMAX is DOUBLE PRECISION
  119: *>          Absolute value of largest matrix element.  If AMAX is very
  120: *>          close to overflow or very close to underflow, the matrix
  121: *>          should be scaled.
  122: *> \endverbatim
  123: *>
  124: *> \param[out] INFO
  125: *> \verbatim
  126: *>          INFO is INTEGER
  127: *>          = 0:  successful exit
  128: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  129: *>          > 0:  if INFO = i,  and i is
  130: *>                <= M:  the i-th row of A is exactly zero
  131: *>                >  M:  the (i-M)-th column of A is exactly zero
  132: *> \endverbatim
  133: *
  134: *  Authors:
  135: *  ========
  136: *
  137: *> \author Univ. of Tennessee 
  138: *> \author Univ. of California Berkeley 
  139: *> \author Univ. of Colorado Denver 
  140: *> \author NAG Ltd. 
  141: *
  142: *> \date November 2011
  143: *
  144: *> \ingroup complex16GEcomputational
  145: *
  146: *  =====================================================================
  147:       SUBROUTINE ZGEEQUB( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
  148:      $                    INFO )
  149: *
  150: *  -- LAPACK computational routine (version 3.4.0) --
  151: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  152: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  153: *     November 2011
  154: *
  155: *     .. Scalar Arguments ..
  156:       INTEGER            INFO, LDA, M, N
  157:       DOUBLE PRECISION   AMAX, COLCND, ROWCND
  158: *     ..
  159: *     .. Array Arguments ..
  160:       DOUBLE PRECISION   C( * ), R( * )
  161:       COMPLEX*16         A( LDA, * )
  162: *     ..
  163: *
  164: *  =====================================================================
  165: *
  166: *     .. Parameters ..
  167:       DOUBLE PRECISION   ONE, ZERO
  168:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  169: *     ..
  170: *     .. Local Scalars ..
  171:       INTEGER            I, J
  172:       DOUBLE PRECISION   BIGNUM, RCMAX, RCMIN, SMLNUM, RADIX, LOGRDX
  173:       COMPLEX*16         ZDUM
  174: *     ..
  175: *     .. External Functions ..
  176:       DOUBLE PRECISION   DLAMCH
  177:       EXTERNAL           DLAMCH
  178: *     ..
  179: *     .. External Subroutines ..
  180:       EXTERNAL           XERBLA
  181: *     ..
  182: *     .. Intrinsic Functions ..
  183:       INTRINSIC          ABS, MAX, MIN, LOG, DBLE, DIMAG
  184: *     ..
  185: *     .. Statement Functions ..
  186:       DOUBLE PRECISION   CABS1
  187: *     ..
  188: *     .. Statement Function definitions ..
  189:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
  190: *     ..
  191: *     .. Executable Statements ..
  192: *
  193: *     Test the input parameters.
  194: *
  195:       INFO = 0
  196:       IF( M.LT.0 ) THEN
  197:          INFO = -1
  198:       ELSE IF( N.LT.0 ) THEN
  199:          INFO = -2
  200:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  201:          INFO = -4
  202:       END IF
  203:       IF( INFO.NE.0 ) THEN
  204:          CALL XERBLA( 'ZGEEQUB', -INFO )
  205:          RETURN
  206:       END IF
  207: *
  208: *     Quick return if possible.
  209: *
  210:       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
  211:          ROWCND = ONE
  212:          COLCND = ONE
  213:          AMAX = ZERO
  214:          RETURN
  215:       END IF
  216: *
  217: *     Get machine constants.  Assume SMLNUM is a power of the radix.
  218: *
  219:       SMLNUM = DLAMCH( 'S' )
  220:       BIGNUM = ONE / SMLNUM
  221:       RADIX = DLAMCH( 'B' )
  222:       LOGRDX = LOG( RADIX )
  223: *
  224: *     Compute row scale factors.
  225: *
  226:       DO 10 I = 1, M
  227:          R( I ) = ZERO
  228:    10 CONTINUE
  229: *
  230: *     Find the maximum element in each row.
  231: *
  232:       DO 30 J = 1, N
  233:          DO 20 I = 1, M
  234:             R( I ) = MAX( R( I ), CABS1( A( I, J ) ) )
  235:    20    CONTINUE
  236:    30 CONTINUE
  237:       DO I = 1, M
  238:          IF( R( I ).GT.ZERO ) THEN
  239:             R( I ) = RADIX**INT( LOG(R( I ) ) / LOGRDX )
  240:          END IF
  241:       END DO
  242: *
  243: *     Find the maximum and minimum scale factors.
  244: *
  245:       RCMIN = BIGNUM
  246:       RCMAX = ZERO
  247:       DO 40 I = 1, M
  248:          RCMAX = MAX( RCMAX, R( I ) )
  249:          RCMIN = MIN( RCMIN, R( I ) )
  250:    40 CONTINUE
  251:       AMAX = RCMAX
  252: *
  253:       IF( RCMIN.EQ.ZERO ) THEN
  254: *
  255: *        Find the first zero scale factor and return an error code.
  256: *
  257:          DO 50 I = 1, M
  258:             IF( R( I ).EQ.ZERO ) THEN
  259:                INFO = I
  260:                RETURN
  261:             END IF
  262:    50    CONTINUE
  263:       ELSE
  264: *
  265: *        Invert the scale factors.
  266: *
  267:          DO 60 I = 1, M
  268:             R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM )
  269:    60    CONTINUE
  270: *
  271: *        Compute ROWCND = min(R(I)) / max(R(I)).
  272: *
  273:          ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
  274:       END IF
  275: *
  276: *     Compute column scale factors.
  277: *
  278:       DO 70 J = 1, N
  279:          C( J ) = ZERO
  280:    70 CONTINUE
  281: *
  282: *     Find the maximum element in each column,
  283: *     assuming the row scaling computed above.
  284: *
  285:       DO 90 J = 1, N
  286:          DO 80 I = 1, M
  287:             C( J ) = MAX( C( J ), CABS1( A( I, J ) )*R( I ) )
  288:    80    CONTINUE
  289:          IF( C( J ).GT.ZERO ) THEN
  290:             C( J ) = RADIX**INT( LOG( C( J ) ) / LOGRDX )
  291:          END IF
  292:    90 CONTINUE
  293: *
  294: *     Find the maximum and minimum scale factors.
  295: *
  296:       RCMIN = BIGNUM
  297:       RCMAX = ZERO
  298:       DO 100 J = 1, N
  299:          RCMIN = MIN( RCMIN, C( J ) )
  300:          RCMAX = MAX( RCMAX, C( J ) )
  301:   100 CONTINUE
  302: *
  303:       IF( RCMIN.EQ.ZERO ) THEN
  304: *
  305: *        Find the first zero scale factor and return an error code.
  306: *
  307:          DO 110 J = 1, N
  308:             IF( C( J ).EQ.ZERO ) THEN
  309:                INFO = M + J
  310:                RETURN
  311:             END IF
  312:   110    CONTINUE
  313:       ELSE
  314: *
  315: *        Invert the scale factors.
  316: *
  317:          DO 120 J = 1, N
  318:             C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM )
  319:   120    CONTINUE
  320: *
  321: *        Compute COLCND = min(C(J)) / max(C(J)).
  322: *
  323:          COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
  324:       END IF
  325: *
  326:       RETURN
  327: *
  328: *     End of ZGEEQUB
  329: *
  330:       END