Return to dgbequ.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

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