File:  [local] / rpl / lapack / lapack / zgeequb.f
Revision 1.14: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:39:16 2023 UTC (9 months, 1 week ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

    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.  Barring over- and underflow, scaling by
   53: *> these factors introduces no additional rounding errors.  However, the
   54: *> scaled entries' magnitudes 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: *> \ingroup complex16GEcomputational
  143: *
  144: *  =====================================================================
  145:       SUBROUTINE ZGEEQUB( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
  146:      $                    INFO )
  147: *
  148: *  -- LAPACK computational routine --
  149: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  150: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  151: *
  152: *     .. Scalar Arguments ..
  153:       INTEGER            INFO, LDA, M, N
  154:       DOUBLE PRECISION   AMAX, COLCND, ROWCND
  155: *     ..
  156: *     .. Array Arguments ..
  157:       DOUBLE PRECISION   C( * ), R( * )
  158:       COMPLEX*16         A( LDA, * )
  159: *     ..
  160: *
  161: *  =====================================================================
  162: *
  163: *     .. Parameters ..
  164:       DOUBLE PRECISION   ONE, ZERO
  165:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  166: *     ..
  167: *     .. Local Scalars ..
  168:       INTEGER            I, J
  169:       DOUBLE PRECISION   BIGNUM, RCMAX, RCMIN, SMLNUM, RADIX, LOGRDX
  170:       COMPLEX*16         ZDUM
  171: *     ..
  172: *     .. External Functions ..
  173:       DOUBLE PRECISION   DLAMCH
  174:       EXTERNAL           DLAMCH
  175: *     ..
  176: *     .. External Subroutines ..
  177:       EXTERNAL           XERBLA
  178: *     ..
  179: *     .. Intrinsic Functions ..
  180:       INTRINSIC          ABS, MAX, MIN, LOG, DBLE, DIMAG
  181: *     ..
  182: *     .. Statement Functions ..
  183:       DOUBLE PRECISION   CABS1
  184: *     ..
  185: *     .. Statement Function definitions ..
  186:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
  187: *     ..
  188: *     .. Executable Statements ..
  189: *
  190: *     Test the input parameters.
  191: *
  192:       INFO = 0
  193:       IF( M.LT.0 ) THEN
  194:          INFO = -1
  195:       ELSE IF( N.LT.0 ) THEN
  196:          INFO = -2
  197:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  198:          INFO = -4
  199:       END IF
  200:       IF( INFO.NE.0 ) THEN
  201:          CALL XERBLA( 'ZGEEQUB', -INFO )
  202:          RETURN
  203:       END IF
  204: *
  205: *     Quick return if possible.
  206: *
  207:       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
  208:          ROWCND = ONE
  209:          COLCND = ONE
  210:          AMAX = ZERO
  211:          RETURN
  212:       END IF
  213: *
  214: *     Get machine constants.  Assume SMLNUM is a power of the radix.
  215: *
  216:       SMLNUM = DLAMCH( 'S' )
  217:       BIGNUM = ONE / SMLNUM
  218:       RADIX = DLAMCH( 'B' )
  219:       LOGRDX = LOG( RADIX )
  220: *
  221: *     Compute row scale factors.
  222: *
  223:       DO 10 I = 1, M
  224:          R( I ) = ZERO
  225:    10 CONTINUE
  226: *
  227: *     Find the maximum element in each row.
  228: *
  229:       DO 30 J = 1, N
  230:          DO 20 I = 1, M
  231:             R( I ) = MAX( R( I ), CABS1( A( I, J ) ) )
  232:    20    CONTINUE
  233:    30 CONTINUE
  234:       DO I = 1, M
  235:          IF( R( I ).GT.ZERO ) THEN
  236:             R( I ) = RADIX**INT( LOG(R( I ) ) / LOGRDX )
  237:          END IF
  238:       END DO
  239: *
  240: *     Find the maximum and minimum scale factors.
  241: *
  242:       RCMIN = BIGNUM
  243:       RCMAX = ZERO
  244:       DO 40 I = 1, M
  245:          RCMAX = MAX( RCMAX, R( I ) )
  246:          RCMIN = MIN( RCMIN, R( I ) )
  247:    40 CONTINUE
  248:       AMAX = RCMAX
  249: *
  250:       IF( RCMIN.EQ.ZERO ) THEN
  251: *
  252: *        Find the first zero scale factor and return an error code.
  253: *
  254:          DO 50 I = 1, M
  255:             IF( R( I ).EQ.ZERO ) THEN
  256:                INFO = I
  257:                RETURN
  258:             END IF
  259:    50    CONTINUE
  260:       ELSE
  261: *
  262: *        Invert the scale factors.
  263: *
  264:          DO 60 I = 1, M
  265:             R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM )
  266:    60    CONTINUE
  267: *
  268: *        Compute ROWCND = min(R(I)) / max(R(I)).
  269: *
  270:          ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
  271:       END IF
  272: *
  273: *     Compute column scale factors.
  274: *
  275:       DO 70 J = 1, N
  276:          C( J ) = ZERO
  277:    70 CONTINUE
  278: *
  279: *     Find the maximum element in each column,
  280: *     assuming the row scaling computed above.
  281: *
  282:       DO 90 J = 1, N
  283:          DO 80 I = 1, M
  284:             C( J ) = MAX( C( J ), CABS1( A( I, J ) )*R( I ) )
  285:    80    CONTINUE
  286:          IF( C( J ).GT.ZERO ) THEN
  287:             C( J ) = RADIX**INT( LOG( C( J ) ) / LOGRDX )
  288:          END IF
  289:    90 CONTINUE
  290: *
  291: *     Find the maximum and minimum scale factors.
  292: *
  293:       RCMIN = BIGNUM
  294:       RCMAX = ZERO
  295:       DO 100 J = 1, N
  296:          RCMIN = MIN( RCMIN, C( J ) )
  297:          RCMAX = MAX( RCMAX, C( J ) )
  298:   100 CONTINUE
  299: *
  300:       IF( RCMIN.EQ.ZERO ) THEN
  301: *
  302: *        Find the first zero scale factor and return an error code.
  303: *
  304:          DO 110 J = 1, N
  305:             IF( C( J ).EQ.ZERO ) THEN
  306:                INFO = M + J
  307:                RETURN
  308:             END IF
  309:   110    CONTINUE
  310:       ELSE
  311: *
  312: *        Invert the scale factors.
  313: *
  314:          DO 120 J = 1, N
  315:             C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM )
  316:   120    CONTINUE
  317: *
  318: *        Compute COLCND = min(C(J)) / max(C(J)).
  319: *
  320:          COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
  321:       END IF
  322: *
  323:       RETURN
  324: *
  325: *     End of ZGEEQUB
  326: *
  327:       END

CVSweb interface <joel.bertrand@systella.fr>