Return to dgetrf2.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de Lapack.

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