Return to dgbequb.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

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