Return to zgetc2.f CVS log

Up to [local] / rpl / lapack / lapack

En route pour la 4.0.15 !

    1:       SUBROUTINE ZGETC2( N, A, LDA, IPIV, JPIV, INFO )
    2: *
    3: *  -- LAPACK auxiliary routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2006
    7: *
    8: *     .. Scalar Arguments ..
    9:       INTEGER            INFO, LDA, N
   10: *     ..
   11: *     .. Array Arguments ..
   12:       INTEGER            IPIV( * ), JPIV( * )
   13:       COMPLEX*16         A( LDA, * )
   14: *     ..
   15: *
   16: *  Purpose
   17: *  =======
   18: *
   19: *  ZGETC2 computes an LU factorization, using complete pivoting, of the
   20: *  n-by-n matrix A. The factorization has the form A = P * L * U * Q,
   21: *  where P and Q are permutation matrices, L is lower triangular with
   22: *  unit diagonal elements and U is upper triangular.
   23: *
   24: *  This is a level 1 BLAS version of the algorithm.
   25: *
   26: *  Arguments
   27: *  =========
   28: *
   29: *  N       (input) INTEGER
   30: *          The order of the matrix A. N >= 0.
   31: *
   32: *  A       (input/output) COMPLEX*16 array, dimension (LDA, N)
   33: *          On entry, the n-by-n matrix to be factored.
   34: *          On exit, the factors L and U from the factorization
   35: *          A = P*L*U*Q; the unit diagonal elements of L are not stored.
   36: *          If U(k, k) appears to be less than SMIN, U(k, k) is given the
   37: *          value of SMIN, giving a nonsingular perturbed system.
   38: *
   39: *  LDA     (input) INTEGER
   40: *          The leading dimension of the array A.  LDA >= max(1, N).
   41: *
   42: *  IPIV    (output) INTEGER array, dimension (N).
   43: *          The pivot indices; for 1 <= i <= N, row i of the
   44: *          matrix has been interchanged with row IPIV(i).
   45: *
   46: *  JPIV    (output) INTEGER array, dimension (N).
   47: *          The pivot indices; for 1 <= j <= N, column j of the
   48: *          matrix has been interchanged with column JPIV(j).
   49: *
   50: *  INFO    (output) INTEGER
   51: *           = 0: successful exit
   52: *           > 0: if INFO = k, U(k, k) is likely to produce overflow if
   53: *                one tries to solve for x in Ax = b. So U is perturbed
   54: *                to avoid the overflow.
   55: *
   56: *  Further Details
   57: *  ===============
   58: *
   59: *  Based on contributions by
   60: *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
   61: *     Umea University, S-901 87 Umea, Sweden.
   62: *
   63: *  =====================================================================
   64: *
   65: *     .. Parameters ..
   66:       DOUBLE PRECISION   ZERO, ONE
   67:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
   68: *     ..
   69: *     .. Local Scalars ..
   70:       INTEGER            I, IP, IPV, J, JP, JPV
   71:       DOUBLE PRECISION   BIGNUM, EPS, SMIN, SMLNUM, XMAX
   72: *     ..
   73: *     .. External Subroutines ..
   74:       EXTERNAL           ZGERU, ZSWAP
   75: *     ..
   76: *     .. External Functions ..
   77:       DOUBLE PRECISION   DLAMCH
   78:       EXTERNAL           DLAMCH
   79: *     ..
   80: *     .. Intrinsic Functions ..
   81:       INTRINSIC          ABS, DCMPLX, MAX
   82: *     ..
   83: *     .. Executable Statements ..
   84: *
   85: *     Set constants to control overflow
   86: *
   87:       INFO = 0
   88:       EPS = DLAMCH( 'P' )
   89:       SMLNUM = DLAMCH( 'S' ) / EPS
   90:       BIGNUM = ONE / SMLNUM
   91:       CALL DLABAD( SMLNUM, BIGNUM )
   92: *
   93: *     Factorize A using complete pivoting.
   94: *     Set pivots less than SMIN to SMIN
   95: *
   96:       DO 40 I = 1, N - 1
   97: *
   98: *        Find max element in matrix A
   99: *
  100:          XMAX = ZERO
  101:          DO 20 IP = I, N
  102:             DO 10 JP = I, N
  103:                IF( ABS( A( IP, JP ) ).GE.XMAX ) THEN
  104:                   XMAX = ABS( A( IP, JP ) )
  105:                   IPV = IP
  106:                   JPV = JP
  107:                END IF
  108:    10       CONTINUE
  109:    20    CONTINUE
  110:          IF( I.EQ.1 )
  111:      $      SMIN = MAX( EPS*XMAX, SMLNUM )
  112: *
  113: *        Swap rows
  114: *
  115:          IF( IPV.NE.I )
  116:      $      CALL ZSWAP( N, A( IPV, 1 ), LDA, A( I, 1 ), LDA )
  117:          IPIV( I ) = IPV
  118: *
  119: *        Swap columns
  120: *
  121:          IF( JPV.NE.I )
  122:      $      CALL ZSWAP( N, A( 1, JPV ), 1, A( 1, I ), 1 )
  123:          JPIV( I ) = JPV
  124: *
  125: *        Check for singularity
  126: *
  127:          IF( ABS( A( I, I ) ).LT.SMIN ) THEN
  128:             INFO = I
  129:             A( I, I ) = DCMPLX( SMIN, ZERO )
  130:          END IF
  131:          DO 30 J = I + 1, N
  132:             A( J, I ) = A( J, I ) / A( I, I )
  133:    30    CONTINUE
  134:          CALL ZGERU( N-I, N-I, -DCMPLX( ONE ), A( I+1, I ), 1,
  135:      $               A( I, I+1 ), LDA, A( I+1, I+1 ), LDA )
  136:    40 CONTINUE
  137: *
  138:       IF( ABS( A( N, N ) ).LT.SMIN ) THEN
  139:          INFO = N
  140:          A( N, N ) = DCMPLX( SMIN, ZERO )
  141:       END IF
  142:       RETURN
  143: *
  144: *     End of ZGETC2
  145: *
  146:       END