Return to zgetrf2.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack.

    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: *> \date June 2016
  110: *
  111: *> \ingroup complex16GEcomputational
  112: *
  113: *  =====================================================================
  114:       RECURSIVE SUBROUTINE ZGETRF2( M, N, A, LDA, IPIV, INFO )
  115: *
  116: *  -- LAPACK computational routine (version 3.6.1) --
  117: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  118: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  119: *     June 2016
  120: *
  121: *     .. Scalar Arguments ..
  122:       INTEGER            INFO, LDA, M, N
  123: *     ..
  124: *     .. Array Arguments ..
  125:       INTEGER            IPIV( * )
  126:       COMPLEX*16         A( LDA, * )
  127: *     ..
  128: *
  129: *  =====================================================================
  130: *
  131: *     .. Parameters ..
  132:       COMPLEX*16         ONE, ZERO
  133:       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
  134:      $                     ZERO = ( 0.0D+0, 0.0D+0 ) )
  135: *     ..
  136: *     .. Local Scalars ..
  137:       DOUBLE PRECISION   SFMIN
  138:       COMPLEX*16         TEMP
  139:       INTEGER            I, IINFO, N1, N2
  140: *     ..
  141: *     .. External Functions ..
  142:       DOUBLE PRECISION   DLAMCH
  143:       INTEGER            IZAMAX
  144:       EXTERNAL           DLAMCH, IZAMAX
  145: *     ..
  146: *     .. External Subroutines ..
  147:       EXTERNAL           ZGEMM, ZSCAL, ZLASWP, ZTRSM, XERBLA
  148: *     ..
  149: *     .. Intrinsic Functions ..
  150:       INTRINSIC          MAX, MIN
  151: *     ..
  152: *     .. Executable Statements ..
  153: *
  154: *     Test the input parameters
  155: *
  156:       INFO = 0
  157:       IF( M.LT.0 ) THEN
  158:          INFO = -1
  159:       ELSE IF( N.LT.0 ) THEN
  160:          INFO = -2
  161:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  162:          INFO = -4
  163:       END IF
  164:       IF( INFO.NE.0 ) THEN
  165:          CALL XERBLA( 'ZGETRF2', -INFO )
  166:          RETURN
  167:       END IF
  168: *
  169: *     Quick return if possible
  170: *
  171:       IF( M.EQ.0 .OR. N.EQ.0 )
  172:      $   RETURN
  173: 
  174:       IF ( M.EQ.1 ) THEN
  175: *
  176: *        Use unblocked code for one row case
  177: *        Just need to handle IPIV and INFO
  178: *
  179:          IPIV( 1 ) = 1
  180:          IF ( A(1,1).EQ.ZERO )
  181:      $      INFO = 1
  182: *
  183:       ELSE IF( N.EQ.1 ) THEN
  184: *
  185: *        Use unblocked code for one column case
  186: *
  187: *
  188: *        Compute machine safe minimum
  189: *
  190:          SFMIN = DLAMCH('S')
  191: *
  192: *        Find pivot and test for singularity
  193: *
  194:          I = IZAMAX( M, A( 1, 1 ), 1 )
  195:          IPIV( 1 ) = I
  196:          IF( A( I, 1 ).NE.ZERO ) THEN
  197: *
  198: *           Apply the interchange
  199: *
  200:             IF( I.NE.1 ) THEN
  201:                TEMP = A( 1, 1 )
  202:                A( 1, 1 ) = A( I, 1 )
  203:                A( I, 1 ) = TEMP
  204:             END IF
  205: *
  206: *           Compute elements 2:M of the column
  207: *
  208:             IF( ABS(A( 1, 1 )) .GE. SFMIN ) THEN
  209:                CALL ZSCAL( M-1, ONE / A( 1, 1 ), A( 2, 1 ), 1 )
  210:             ELSE
  211:                DO 10 I = 1, M-1
  212:                   A( 1+I, 1 ) = A( 1+I, 1 ) / A( 1, 1 )
  213:    10          CONTINUE
  214:             END IF
  215: *
  216:          ELSE
  217:             INFO = 1
  218:          END IF
  219: 
  220:       ELSE
  221: *
  222: *        Use recursive code
  223: *
  224:          N1 = MIN( M, N ) / 2
  225:          N2 = N-N1
  226: *
  227: *               [ A11 ]
  228: *        Factor [ --- ]
  229: *               [ A21 ]
  230: *
  231:          CALL ZGETRF2( M, N1, A, LDA, IPIV, IINFO )
  232: 
  233:          IF ( INFO.EQ.0 .AND. IINFO.GT.0 )
  234:      $      INFO = IINFO
  235: *
  236: *                              [ A12 ]
  237: *        Apply interchanges to [ --- ]
  238: *                              [ A22 ]
  239: *
  240:          CALL ZLASWP( N2, A( 1, N1+1 ), LDA, 1, N1, IPIV, 1 )
  241: *
  242: *        Solve A12
  243: *
  244:          CALL ZTRSM( 'L', 'L', 'N', 'U', N1, N2, ONE, A, LDA, 
  245:      $               A( 1, N1+1 ), LDA )
  246: *
  247: *        Update A22
  248: *
  249:          CALL ZGEMM( 'N', 'N', M-N1, N2, N1, -ONE, A( N1+1, 1 ), LDA, 
  250:      $               A( 1, N1+1 ), LDA, ONE, A( N1+1, N1+1 ), LDA )
  251: *
  252: *        Factor A22
  253: *
  254:          CALL ZGETRF2( M-N1, N2, A( N1+1, N1+1 ), LDA, IPIV( N1+1 ),
  255:      $                 IINFO )
  256: *
  257: *        Adjust INFO and the pivot indices
  258: *
  259:          IF ( INFO.EQ.0 .AND. IINFO.GT.0 )
  260:      $      INFO = IINFO + N1
  261:          DO 20 I = N1+1, MIN( M, N )
  262:             IPIV( I ) = IPIV( I ) + N1
  263:    20    CONTINUE
  264: *
  265: *        Apply interchanges to A21
  266: *
  267:          CALL ZLASWP( N1, A( 1, 1 ), LDA, N1+1, MIN( M, N), IPIV, 1 )
  268: *
  269:       END IF
  270:       RETURN
  271: *
  272: *     End of ZGETRF2
  273: *
  274:       END