Return to zgetc2.f CVS log

Up to [local] / rpl / lapack / lapack

Première mise à jour de lapack et blas.

    1: *> \brief \b ZGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix.
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download ZGETC2 + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgetc2.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgetc2.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgetc2.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZGETC2( N, A, LDA, IPIV, JPIV, INFO )
   22: *
   23: *       .. Scalar Arguments ..
   24: *       INTEGER            INFO, LDA, N
   25: *       ..
   26: *       .. Array Arguments ..
   27: *       INTEGER            IPIV( * ), JPIV( * )
   28: *       COMPLEX*16         A( LDA, * )
   29: *       ..
   30: *
   31: *
   32: *> \par Purpose:
   33: *  =============
   34: *>
   35: *> \verbatim
   36: *>
   37: *> ZGETC2 computes an LU factorization, using complete pivoting, of the
   38: *> n-by-n matrix A. The factorization has the form A = P * L * U * Q,
   39: *> where P and Q are permutation matrices, L is lower triangular with
   40: *> unit diagonal elements and U is upper triangular.
   41: *>
   42: *> This is a level 1 BLAS version of the algorithm.
   43: *> \endverbatim
   44: *
   45: *  Arguments:
   46: *  ==========
   47: *
   48: *> \param[in] N
   49: *> \verbatim
   50: *>          N is INTEGER
   51: *>          The order of the matrix A. N >= 0.
   52: *> \endverbatim
   53: *>
   54: *> \param[in,out] A
   55: *> \verbatim
   56: *>          A is COMPLEX*16 array, dimension (LDA, N)
   57: *>          On entry, the n-by-n matrix to be factored.
   58: *>          On exit, the factors L and U from the factorization
   59: *>          A = P*L*U*Q; the unit diagonal elements of L are not stored.
   60: *>          If U(k, k) appears to be less than SMIN, U(k, k) is given the
   61: *>          value of SMIN, giving a nonsingular perturbed system.
   62: *> \endverbatim
   63: *>
   64: *> \param[in] LDA
   65: *> \verbatim
   66: *>          LDA is INTEGER
   67: *>          The leading dimension of the array A.  LDA >= max(1, N).
   68: *> \endverbatim
   69: *>
   70: *> \param[out] IPIV
   71: *> \verbatim
   72: *>          IPIV is INTEGER array, dimension (N).
   73: *>          The pivot indices; for 1 <= i <= N, row i of the
   74: *>          matrix has been interchanged with row IPIV(i).
   75: *> \endverbatim
   76: *>
   77: *> \param[out] JPIV
   78: *> \verbatim
   79: *>          JPIV is INTEGER array, dimension (N).
   80: *>          The pivot indices; for 1 <= j <= N, column j of the
   81: *>          matrix has been interchanged with column JPIV(j).
   82: *> \endverbatim
   83: *>
   84: *> \param[out] INFO
   85: *> \verbatim
   86: *>          INFO is INTEGER
   87: *>           = 0: successful exit
   88: *>           > 0: if INFO = k, U(k, k) is likely to produce overflow if
   89: *>                one tries to solve for x in Ax = b. So U is perturbed
   90: *>                to avoid the overflow.
   91: *> \endverbatim
   92: *
   93: *  Authors:
   94: *  ========
   95: *
   96: *> \author Univ. of Tennessee
   97: *> \author Univ. of California Berkeley
   98: *> \author Univ. of Colorado Denver
   99: *> \author NAG Ltd.
  100: *
  101: *> \ingroup complex16GEauxiliary
  102: *
  103: *> \par Contributors:
  104: *  ==================
  105: *>
  106: *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
  107: *>     Umea University, S-901 87 Umea, Sweden.
  108: *
  109: *  =====================================================================
  110:       SUBROUTINE ZGETC2( N, A, LDA, IPIV, JPIV, INFO )
  111: *
  112: *  -- LAPACK auxiliary routine --
  113: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  114: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  115: *
  116: *     .. Scalar Arguments ..
  117:       INTEGER            INFO, LDA, N
  118: *     ..
  119: *     .. Array Arguments ..
  120:       INTEGER            IPIV( * ), JPIV( * )
  121:       COMPLEX*16         A( LDA, * )
  122: *     ..
  123: *
  124: *  =====================================================================
  125: *
  126: *     .. Parameters ..
  127:       DOUBLE PRECISION   ZERO, ONE
  128:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  129: *     ..
  130: *     .. Local Scalars ..
  131:       INTEGER            I, IP, IPV, J, JP, JPV
  132:       DOUBLE PRECISION   BIGNUM, EPS, SMIN, SMLNUM, XMAX
  133: *     ..
  134: *     .. External Subroutines ..
  135:       EXTERNAL           ZGERU, ZSWAP, DLABAD
  136: *     ..
  137: *     .. External Functions ..
  138:       DOUBLE PRECISION   DLAMCH
  139:       EXTERNAL           DLAMCH
  140: *     ..
  141: *     .. Intrinsic Functions ..
  142:       INTRINSIC          ABS, DCMPLX, MAX
  143: *     ..
  144: *     .. Executable Statements ..
  145: *
  146:       INFO = 0
  147: *
  148: *     Quick return if possible
  149: *
  150:       IF( N.EQ.0 )
  151:      $   RETURN
  152: *
  153: *     Set constants to control overflow
  154: *
  155:       EPS = DLAMCH( 'P' )
  156:       SMLNUM = DLAMCH( 'S' ) / EPS
  157:       BIGNUM = ONE / SMLNUM
  158:       CALL DLABAD( SMLNUM, BIGNUM )
  159: *
  160: *     Handle the case N=1 by itself
  161: *
  162:       IF( N.EQ.1 ) THEN
  163:          IPIV( 1 ) = 1
  164:          JPIV( 1 ) = 1
  165:          IF( ABS( A( 1, 1 ) ).LT.SMLNUM ) THEN
  166:             INFO = 1
  167:             A( 1, 1 ) = DCMPLX( SMLNUM, ZERO )
  168:          END IF
  169:          RETURN
  170:       END IF
  171: *
  172: *     Factorize A using complete pivoting.
  173: *     Set pivots less than SMIN to SMIN
  174: *
  175:       DO 40 I = 1, N - 1
  176: *
  177: *        Find max element in matrix A
  178: *
  179:          XMAX = ZERO
  180:          DO 20 IP = I, N
  181:             DO 10 JP = I, N
  182:                IF( ABS( A( IP, JP ) ).GE.XMAX ) THEN
  183:                   XMAX = ABS( A( IP, JP ) )
  184:                   IPV = IP
  185:                   JPV = JP
  186:                END IF
  187:    10       CONTINUE
  188:    20    CONTINUE
  189:          IF( I.EQ.1 )
  190:      $      SMIN = MAX( EPS*XMAX, SMLNUM )
  191: *
  192: *        Swap rows
  193: *
  194:          IF( IPV.NE.I )
  195:      $      CALL ZSWAP( N, A( IPV, 1 ), LDA, A( I, 1 ), LDA )
  196:          IPIV( I ) = IPV
  197: *
  198: *        Swap columns
  199: *
  200:          IF( JPV.NE.I )
  201:      $      CALL ZSWAP( N, A( 1, JPV ), 1, A( 1, I ), 1 )
  202:          JPIV( I ) = JPV
  203: *
  204: *        Check for singularity
  205: *
  206:          IF( ABS( A( I, I ) ).LT.SMIN ) THEN
  207:             INFO = I
  208:             A( I, I ) = DCMPLX( SMIN, ZERO )
  209:          END IF
  210:          DO 30 J = I + 1, N
  211:             A( J, I ) = A( J, I ) / A( I, I )
  212:    30    CONTINUE
  213:          CALL ZGERU( N-I, N-I, -DCMPLX( ONE ), A( I+1, I ), 1,
  214:      $               A( I, I+1 ), LDA, A( I+1, I+1 ), LDA )
  215:    40 CONTINUE
  216: *
  217:       IF( ABS( A( N, N ) ).LT.SMIN ) THEN
  218:          INFO = N
  219:          A( N, N ) = DCMPLX( SMIN, ZERO )
  220:       END IF
  221: *
  222: *     Set last pivots to N
  223: *
  224:       IPIV( N ) = N
  225:       JPIV( N ) = N
  226: *
  227:       RETURN
  228: *
  229: *     End of ZGETC2
  230: *
  231:       END