Return to zgetrf2.f CVS log

Up to [local] / rpl / lapack / lapack

Première mise à jour de lapack et blas.

    1: *> \brief \b ZGETRF2
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *  Definition:
    9: *  ===========
   10: *
   11: *       RECURSIVE SUBROUTINE ZGETRF2( M, N, A, LDA, IPIV, INFO )
   12: *
   13: *       .. Scalar Arguments ..
   14: *       INTEGER            INFO, LDA, M, N
   15: *       ..
   16: *       .. Array Arguments ..
   17: *       INTEGER            IPIV( * )
   18: *       COMPLEX*16         A( LDA, * )
   19: *       ..
   20: *
   21: *
   22: *> \par Purpose:
   23: *  =============
   24: *>
   25: *> \verbatim
   26: *>
   27: *> ZGETRF2 computes an LU factorization of a general M-by-N matrix A
   28: *> using partial pivoting with row interchanges.
   29: *>
   30: *> The factorization has the form
   31: *>    A = P * L * U
   32: *> where P is a permutation matrix, L is lower triangular with unit
   33: *> diagonal elements (lower trapezoidal if m > n), and U is upper
   34: *> triangular (upper trapezoidal if m < n).
   35: *>
   36: *> This is the recursive version of the algorithm. It divides
   37: *> the matrix into four submatrices:
   38: *>
   39: *>        [  A11 | A12  ]  where A11 is n1 by n1 and A22 is n2 by n2
   40: *>    A = [ -----|----- ]  with n1 = min(m,n)/2
   41: *>        [  A21 | A22  ]       n2 = n-n1
   42: *>
   43: *>                                       [ A11 ]
   44: *> The subroutine calls itself to factor [ --- ],
   45: *>                                       [ A12 ]
   46: *>                 [ A12 ]
   47: *> do the swaps on [ --- ], solve A12, update A22,
   48: *>                 [ A22 ]
   49: *>
   50: *> then calls itself to factor A22 and do the swaps on A21.
   51: *>
   52: *> \endverbatim
   53: *
   54: *  Arguments:
   55: *  ==========
   56: *
   57: *> \param[in] M
   58: *> \verbatim
   59: *>          M is INTEGER
   60: *>          The number of rows of the matrix A.  M >= 0.
   61: *> \endverbatim
   62: *>
   63: *> \param[in] N
   64: *> \verbatim
   65: *>          N is INTEGER
   66: *>          The number of columns of the matrix A.  N >= 0.
   67: *> \endverbatim
   68: *>
   69: *> \param[in,out] A
   70: *> \verbatim
   71: *>          A is COMPLEX*16 array, dimension (LDA,N)
   72: *>          On entry, the M-by-N matrix to be factored.
   73: *>          On exit, the factors L and U from the factorization
   74: *>          A = P*L*U; the unit diagonal elements of L are not stored.
   75: *> \endverbatim
   76: *>
   77: *> \param[in] LDA
   78: *> \verbatim
   79: *>          LDA is INTEGER
   80: *>          The leading dimension of the array A.  LDA >= max(1,M).
   81: *> \endverbatim
   82: *>
   83: *> \param[out] IPIV
   84: *> \verbatim
   85: *>          IPIV is INTEGER array, dimension (min(M,N))
   86: *>          The pivot indices; for 1 <= i <= min(M,N), row i of the
   87: *>          matrix was interchanged with row IPIV(i).
   88: *> \endverbatim
   89: *>
   90: *> \param[out] INFO
   91: *> \verbatim
   92: *>          INFO is INTEGER
   93: *>          = 0:  successful exit
   94: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
   95: *>          > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
   96: *>                has been completed, but the factor U is exactly
   97: *>                singular, and division by zero will occur if it is used
   98: *>                to solve a system of equations.
   99: *> \endverbatim
  100: *
  101: *  Authors:
  102: *  ========
  103: *
  104: *> \author Univ. of Tennessee
  105: *> \author Univ. of California Berkeley
  106: *> \author Univ. of Colorado Denver
  107: *> \author NAG Ltd.
  108: *
  109: *> \ingroup complex16GEcomputational
  110: *
  111: *  =====================================================================
  112:       RECURSIVE SUBROUTINE ZGETRF2( M, N, A, LDA, IPIV, INFO )
  113: *
  114: *  -- LAPACK computational routine --
  115: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  116: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  117: *
  118: *     .. Scalar Arguments ..
  119:       INTEGER            INFO, LDA, M, N
  120: *     ..
  121: *     .. Array Arguments ..
  122:       INTEGER            IPIV( * )
  123:       COMPLEX*16         A( LDA, * )
  124: *     ..
  125: *
  126: *  =====================================================================
  127: *
  128: *     .. Parameters ..
  129:       COMPLEX*16         ONE, ZERO
  130:       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
  131:      $                     ZERO = ( 0.0D+0, 0.0D+0 ) )
  132: *     ..
  133: *     .. Local Scalars ..
  134:       DOUBLE PRECISION   SFMIN
  135:       COMPLEX*16         TEMP
  136:       INTEGER            I, IINFO, N1, N2
  137: *     ..
  138: *     .. External Functions ..
  139:       DOUBLE PRECISION   DLAMCH
  140:       INTEGER            IZAMAX
  141:       EXTERNAL           DLAMCH, IZAMAX
  142: *     ..
  143: *     .. External Subroutines ..
  144:       EXTERNAL           ZGEMM, ZSCAL, ZLASWP, ZTRSM, XERBLA
  145: *     ..
  146: *     .. Intrinsic Functions ..
  147:       INTRINSIC          MAX, MIN
  148: *     ..
  149: *     .. Executable Statements ..
  150: *
  151: *     Test the input parameters
  152: *
  153:       INFO = 0
  154:       IF( M.LT.0 ) THEN
  155:          INFO = -1
  156:       ELSE IF( N.LT.0 ) THEN
  157:          INFO = -2
  158:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  159:          INFO = -4
  160:       END IF
  161:       IF( INFO.NE.0 ) THEN
  162:          CALL XERBLA( 'ZGETRF2', -INFO )
  163:          RETURN
  164:       END IF
  165: *
  166: *     Quick return if possible
  167: *
  168:       IF( M.EQ.0 .OR. N.EQ.0 )
  169:      $   RETURN
  170: 
  171:       IF ( M.EQ.1 ) THEN
  172: *
  173: *        Use unblocked code for one row case
  174: *        Just need to handle IPIV and INFO
  175: *
  176:          IPIV( 1 ) = 1
  177:          IF ( A(1,1).EQ.ZERO )
  178:      $      INFO = 1
  179: *
  180:       ELSE IF( N.EQ.1 ) THEN
  181: *
  182: *        Use unblocked code for one column case
  183: *
  184: *
  185: *        Compute machine safe minimum
  186: *
  187:          SFMIN = DLAMCH('S')
  188: *
  189: *        Find pivot and test for singularity
  190: *
  191:          I = IZAMAX( M, A( 1, 1 ), 1 )
  192:          IPIV( 1 ) = I
  193:          IF( A( I, 1 ).NE.ZERO ) THEN
  194: *
  195: *           Apply the interchange
  196: *
  197:             IF( I.NE.1 ) THEN
  198:                TEMP = A( 1, 1 )
  199:                A( 1, 1 ) = A( I, 1 )
  200:                A( I, 1 ) = TEMP
  201:             END IF
  202: *
  203: *           Compute elements 2:M of the column
  204: *
  205:             IF( ABS(A( 1, 1 )) .GE. SFMIN ) THEN
  206:                CALL ZSCAL( M-1, ONE / A( 1, 1 ), A( 2, 1 ), 1 )
  207:             ELSE
  208:                DO 10 I = 1, M-1
  209:                   A( 1+I, 1 ) = A( 1+I, 1 ) / A( 1, 1 )
  210:    10          CONTINUE
  211:             END IF
  212: *
  213:          ELSE
  214:             INFO = 1
  215:          END IF
  216: 
  217:       ELSE
  218: *
  219: *        Use recursive code
  220: *
  221:          N1 = MIN( M, N ) / 2
  222:          N2 = N-N1
  223: *
  224: *               [ A11 ]
  225: *        Factor [ --- ]
  226: *               [ A21 ]
  227: *
  228:          CALL ZGETRF2( M, N1, A, LDA, IPIV, IINFO )
  229: 
  230:          IF ( INFO.EQ.0 .AND. IINFO.GT.0 )
  231:      $      INFO = IINFO
  232: *
  233: *                              [ A12 ]
  234: *        Apply interchanges to [ --- ]
  235: *                              [ A22 ]
  236: *
  237:          CALL ZLASWP( N2, A( 1, N1+1 ), LDA, 1, N1, IPIV, 1 )
  238: *
  239: *        Solve A12
  240: *
  241:          CALL ZTRSM( 'L', 'L', 'N', 'U', N1, N2, ONE, A, LDA,
  242:      $               A( 1, N1+1 ), LDA )
  243: *
  244: *        Update A22
  245: *
  246:          CALL ZGEMM( 'N', 'N', M-N1, N2, N1, -ONE, A( N1+1, 1 ), LDA,
  247:      $               A( 1, N1+1 ), LDA, ONE, A( N1+1, N1+1 ), LDA )
  248: *
  249: *        Factor A22
  250: *
  251:          CALL ZGETRF2( M-N1, N2, A( N1+1, N1+1 ), LDA, IPIV( N1+1 ),
  252:      $                 IINFO )
  253: *
  254: *        Adjust INFO and the pivot indices
  255: *
  256:          IF ( INFO.EQ.0 .AND. IINFO.GT.0 )
  257:      $      INFO = IINFO + N1
  258:          DO 20 I = N1+1, MIN( M, N )
  259:             IPIV( I ) = IPIV( I ) + N1
  260:    20    CONTINUE
  261: *
  262: *        Apply interchanges to A21
  263: *
  264:          CALL ZLASWP( N1, A( 1, 1 ), LDA, N1+1, MIN( M, N), IPIV, 1 )
  265: *
  266:       END IF
  267:       RETURN
  268: *
  269: *     End of ZGETRF2
  270: *
  271:       END