Return to dgetri.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack.

    1: *> \brief \b DGETRI
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download DGETRI + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgetri.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgetri.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgetri.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
   22: *
   23: *       .. Scalar Arguments ..
   24: *       INTEGER            INFO, LDA, LWORK, N
   25: *       ..
   26: *       .. Array Arguments ..
   27: *       INTEGER            IPIV( * )
   28: *       DOUBLE PRECISION   A( LDA, * ), WORK( * )
   29: *       ..
   30: *
   31: *
   32: *> \par Purpose:
   33: *  =============
   34: *>
   35: *> \verbatim
   36: *>
   37: *> DGETRI computes the inverse of a matrix using the LU factorization
   38: *> computed by DGETRF.
   39: *>
   40: *> This method inverts U and then computes inv(A) by solving the system
   41: *> inv(A)*L = inv(U) for inv(A).
   42: *> \endverbatim
   43: *
   44: *  Arguments:
   45: *  ==========
   46: *
   47: *> \param[in] N
   48: *> \verbatim
   49: *>          N is INTEGER
   50: *>          The order of the matrix A.  N >= 0.
   51: *> \endverbatim
   52: *>
   53: *> \param[in,out] A
   54: *> \verbatim
   55: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
   56: *>          On entry, the factors L and U from the factorization
   57: *>          A = P*L*U as computed by DGETRF.
   58: *>          On exit, if INFO = 0, the inverse of the original matrix A.
   59: *> \endverbatim
   60: *>
   61: *> \param[in] LDA
   62: *> \verbatim
   63: *>          LDA is INTEGER
   64: *>          The leading dimension of the array A.  LDA >= max(1,N).
   65: *> \endverbatim
   66: *>
   67: *> \param[in] IPIV
   68: *> \verbatim
   69: *>          IPIV is INTEGER array, dimension (N)
   70: *>          The pivot indices from DGETRF; for 1<=i<=N, row i of the
   71: *>          matrix was interchanged with row IPIV(i).
   72: *> \endverbatim
   73: *>
   74: *> \param[out] WORK
   75: *> \verbatim
   76: *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
   77: *>          On exit, if INFO=0, then WORK(1) returns the optimal LWORK.
   78: *> \endverbatim
   79: *>
   80: *> \param[in] LWORK
   81: *> \verbatim
   82: *>          LWORK is INTEGER
   83: *>          The dimension of the array WORK.  LWORK >= max(1,N).
   84: *>          For optimal performance LWORK >= N*NB, where NB is
   85: *>          the optimal blocksize returned by ILAENV.
   86: *>
   87: *>          If LWORK = -1, then a workspace query is assumed; the routine
   88: *>          only calculates the optimal size of the WORK array, returns
   89: *>          this value as the first entry of the WORK array, and no error
   90: *>          message related to LWORK is issued by XERBLA.
   91: *> \endverbatim
   92: *>
   93: *> \param[out] INFO
   94: *> \verbatim
   95: *>          INFO is INTEGER
   96: *>          = 0:  successful exit
   97: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
   98: *>          > 0:  if INFO = i, U(i,i) is exactly zero; the matrix is
   99: *>                singular and its inverse could not be computed.
  100: *> \endverbatim
  101: *
  102: *  Authors:
  103: *  ========
  104: *
  105: *> \author Univ. of Tennessee
  106: *> \author Univ. of California Berkeley
  107: *> \author Univ. of Colorado Denver
  108: *> \author NAG Ltd.
  109: *
  110: *> \date December 2016
  111: *
  112: *> \ingroup doubleGEcomputational
  113: *
  114: *  =====================================================================
  115:       SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
  116: *
  117: *  -- LAPACK computational routine (version 3.7.0) --
  118: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  119: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  120: *     December 2016
  121: *
  122: *     .. Scalar Arguments ..
  123:       INTEGER            INFO, LDA, LWORK, N
  124: *     ..
  125: *     .. Array Arguments ..
  126:       INTEGER            IPIV( * )
  127:       DOUBLE PRECISION   A( LDA, * ), WORK( * )
  128: *     ..
  129: *
  130: *  =====================================================================
  131: *
  132: *     .. Parameters ..
  133:       DOUBLE PRECISION   ZERO, ONE
  134:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  135: *     ..
  136: *     .. Local Scalars ..
  137:       LOGICAL            LQUERY
  138:       INTEGER            I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB,
  139:      $                   NBMIN, NN
  140: *     ..
  141: *     .. External Functions ..
  142:       INTEGER            ILAENV
  143:       EXTERNAL           ILAENV
  144: *     ..
  145: *     .. External Subroutines ..
  146:       EXTERNAL           DGEMM, DGEMV, DSWAP, DTRSM, DTRTRI, XERBLA
  147: *     ..
  148: *     .. Intrinsic Functions ..
  149:       INTRINSIC          MAX, MIN
  150: *     ..
  151: *     .. Executable Statements ..
  152: *
  153: *     Test the input parameters.
  154: *
  155:       INFO = 0
  156:       NB = ILAENV( 1, 'DGETRI', ' ', N, -1, -1, -1 )
  157:       LWKOPT = N*NB
  158:       WORK( 1 ) = LWKOPT
  159:       LQUERY = ( LWORK.EQ.-1 )
  160:       IF( N.LT.0 ) THEN
  161:          INFO = -1
  162:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  163:          INFO = -3
  164:       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
  165:          INFO = -6
  166:       END IF
  167:       IF( INFO.NE.0 ) THEN
  168:          CALL XERBLA( 'DGETRI', -INFO )
  169:          RETURN
  170:       ELSE IF( LQUERY ) THEN
  171:          RETURN
  172:       END IF
  173: *
  174: *     Quick return if possible
  175: *
  176:       IF( N.EQ.0 )
  177:      $   RETURN
  178: *
  179: *     Form inv(U).  If INFO > 0 from DTRTRI, then U is singular,
  180: *     and the inverse is not computed.
  181: *
  182:       CALL DTRTRI( 'Upper', 'Non-unit', N, A, LDA, INFO )
  183:       IF( INFO.GT.0 )
  184:      $   RETURN
  185: *
  186:       NBMIN = 2
  187:       LDWORK = N
  188:       IF( NB.GT.1 .AND. NB.LT.N ) THEN
  189:          IWS = MAX( LDWORK*NB, 1 )
  190:          IF( LWORK.LT.IWS ) THEN
  191:             NB = LWORK / LDWORK
  192:             NBMIN = MAX( 2, ILAENV( 2, 'DGETRI', ' ', N, -1, -1, -1 ) )
  193:          END IF
  194:       ELSE
  195:          IWS = N
  196:       END IF
  197: *
  198: *     Solve the equation inv(A)*L = inv(U) for inv(A).
  199: *
  200:       IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN
  201: *
  202: *        Use unblocked code.
  203: *
  204:          DO 20 J = N, 1, -1
  205: *
  206: *           Copy current column of L to WORK and replace with zeros.
  207: *
  208:             DO 10 I = J + 1, N
  209:                WORK( I ) = A( I, J )
  210:                A( I, J ) = ZERO
  211:    10       CONTINUE
  212: *
  213: *           Compute current column of inv(A).
  214: *
  215:             IF( J.LT.N )
  216:      $         CALL DGEMV( 'No transpose', N, N-J, -ONE, A( 1, J+1 ),
  217:      $                     LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 )
  218:    20    CONTINUE
  219:       ELSE
  220: *
  221: *        Use blocked code.
  222: *
  223:          NN = ( ( N-1 ) / NB )*NB + 1
  224:          DO 50 J = NN, 1, -NB
  225:             JB = MIN( NB, N-J+1 )
  226: *
  227: *           Copy current block column of L to WORK and replace with
  228: *           zeros.
  229: *
  230:             DO 40 JJ = J, J + JB - 1
  231:                DO 30 I = JJ + 1, N
  232:                   WORK( I+( JJ-J )*LDWORK ) = A( I, JJ )
  233:                   A( I, JJ ) = ZERO
  234:    30          CONTINUE
  235:    40       CONTINUE
  236: *
  237: *           Compute current block column of inv(A).
  238: *
  239:             IF( J+JB.LE.N )
  240:      $         CALL DGEMM( 'No transpose', 'No transpose', N, JB,
  241:      $                     N-J-JB+1, -ONE, A( 1, J+JB ), LDA,
  242:      $                     WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA )
  243:             CALL DTRSM( 'Right', 'Lower', 'No transpose', 'Unit', N, JB,
  244:      $                  ONE, WORK( J ), LDWORK, A( 1, J ), LDA )
  245:    50    CONTINUE
  246:       END IF
  247: *
  248: *     Apply column interchanges.
  249: *
  250:       DO 60 J = N - 1, 1, -1
  251:          JP = IPIV( J )
  252:          IF( JP.NE.J )
  253:      $      CALL DSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 )
  254:    60 CONTINUE
  255: *
  256:       WORK( 1 ) = IWS
  257:       RETURN
  258: *
  259: *     End of DGETRI
  260: *
  261:       END