Return to zgetri.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack.

    1: *> \brief \b ZGETRI
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download ZGETRI + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgetri.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgetri.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgetri.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZGETRI( 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: *       COMPLEX*16         A( LDA, * ), WORK( * )
   29: *       ..
   30: *  
   31: *
   32: *> \par Purpose:
   33: *  =============
   34: *>
   35: *> \verbatim
   36: *>
   37: *> ZGETRI computes the inverse of a matrix using the LU factorization
   38: *> computed by ZGETRF.
   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 COMPLEX*16 array, dimension (LDA,N)
   56: *>          On entry, the factors L and U from the factorization
   57: *>          A = P*L*U as computed by ZGETRF.
   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 ZGETRF; 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 COMPLEX*16 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 November 2011
  111: *
  112: *> \ingroup complex16GEcomputational
  113: *
  114: *  =====================================================================
  115:       SUBROUTINE ZGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
  116: *
  117: *  -- LAPACK computational routine (version 3.4.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: *     November 2011
  121: *
  122: *     .. Scalar Arguments ..
  123:       INTEGER            INFO, LDA, LWORK, N
  124: *     ..
  125: *     .. Array Arguments ..
  126:       INTEGER            IPIV( * )
  127:       COMPLEX*16         A( LDA, * ), WORK( * )
  128: *     ..
  129: *
  130: *  =====================================================================
  131: *
  132: *     .. Parameters ..
  133:       COMPLEX*16         ZERO, ONE
  134:       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
  135:      $                   ONE = ( 1.0D+0, 0.0D+0 ) )
  136: *     ..
  137: *     .. Local Scalars ..
  138:       LOGICAL            LQUERY
  139:       INTEGER            I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB,
  140:      $                   NBMIN, NN
  141: *     ..
  142: *     .. External Functions ..
  143:       INTEGER            ILAENV
  144:       EXTERNAL           ILAENV
  145: *     ..
  146: *     .. External Subroutines ..
  147:       EXTERNAL           XERBLA, ZGEMM, ZGEMV, ZSWAP, ZTRSM, ZTRTRI
  148: *     ..
  149: *     .. Intrinsic Functions ..
  150:       INTRINSIC          MAX, MIN
  151: *     ..
  152: *     .. Executable Statements ..
  153: *
  154: *     Test the input parameters.
  155: *
  156:       INFO = 0
  157:       NB = ILAENV( 1, 'ZGETRI', ' ', N, -1, -1, -1 )
  158:       LWKOPT = N*NB
  159:       WORK( 1 ) = LWKOPT
  160:       LQUERY = ( LWORK.EQ.-1 )
  161:       IF( N.LT.0 ) THEN
  162:          INFO = -1
  163:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  164:          INFO = -3
  165:       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
  166:          INFO = -6
  167:       END IF
  168:       IF( INFO.NE.0 ) THEN
  169:          CALL XERBLA( 'ZGETRI', -INFO )
  170:          RETURN
  171:       ELSE IF( LQUERY ) THEN
  172:          RETURN
  173:       END IF
  174: *
  175: *     Quick return if possible
  176: *
  177:       IF( N.EQ.0 )
  178:      $   RETURN
  179: *
  180: *     Form inv(U).  If INFO > 0 from ZTRTRI, then U is singular,
  181: *     and the inverse is not computed.
  182: *
  183:       CALL ZTRTRI( 'Upper', 'Non-unit', N, A, LDA, INFO )
  184:       IF( INFO.GT.0 )
  185:      $   RETURN
  186: *
  187:       NBMIN = 2
  188:       LDWORK = N
  189:       IF( NB.GT.1 .AND. NB.LT.N ) THEN
  190:          IWS = MAX( LDWORK*NB, 1 )
  191:          IF( LWORK.LT.IWS ) THEN
  192:             NB = LWORK / LDWORK
  193:             NBMIN = MAX( 2, ILAENV( 2, 'ZGETRI', ' ', N, -1, -1, -1 ) )
  194:          END IF
  195:       ELSE
  196:          IWS = N
  197:       END IF
  198: *
  199: *     Solve the equation inv(A)*L = inv(U) for inv(A).
  200: *
  201:       IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN
  202: *
  203: *        Use unblocked code.
  204: *
  205:          DO 20 J = N, 1, -1
  206: *
  207: *           Copy current column of L to WORK and replace with zeros.
  208: *
  209:             DO 10 I = J + 1, N
  210:                WORK( I ) = A( I, J )
  211:                A( I, J ) = ZERO
  212:    10       CONTINUE
  213: *
  214: *           Compute current column of inv(A).
  215: *
  216:             IF( J.LT.N )
  217:      $         CALL ZGEMV( 'No transpose', N, N-J, -ONE, A( 1, J+1 ),
  218:      $                     LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 )
  219:    20    CONTINUE
  220:       ELSE
  221: *
  222: *        Use blocked code.
  223: *
  224:          NN = ( ( N-1 ) / NB )*NB + 1
  225:          DO 50 J = NN, 1, -NB
  226:             JB = MIN( NB, N-J+1 )
  227: *
  228: *           Copy current block column of L to WORK and replace with
  229: *           zeros.
  230: *
  231:             DO 40 JJ = J, J + JB - 1
  232:                DO 30 I = JJ + 1, N
  233:                   WORK( I+( JJ-J )*LDWORK ) = A( I, JJ )
  234:                   A( I, JJ ) = ZERO
  235:    30          CONTINUE
  236:    40       CONTINUE
  237: *
  238: *           Compute current block column of inv(A).
  239: *
  240:             IF( J+JB.LE.N )
  241:      $         CALL ZGEMM( 'No transpose', 'No transpose', N, JB,
  242:      $                     N-J-JB+1, -ONE, A( 1, J+JB ), LDA,
  243:      $                     WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA )
  244:             CALL ZTRSM( 'Right', 'Lower', 'No transpose', 'Unit', N, JB,
  245:      $                  ONE, WORK( J ), LDWORK, A( 1, J ), LDA )
  246:    50    CONTINUE
  247:       END IF
  248: *
  249: *     Apply column interchanges.
  250: *
  251:       DO 60 J = N - 1, 1, -1
  252:          JP = IPIV( J )
  253:          IF( JP.NE.J )
  254:      $      CALL ZSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 )
  255:    60 CONTINUE
  256: *
  257:       WORK( 1 ) = IWS
  258:       RETURN
  259: *
  260: *     End of ZGETRI
  261: *
  262:       END