LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zgetri.f
Go to the documentation of this file.
00001 *> \brief \b ZGETRI
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZGETRI + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgetri.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgetri.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgetri.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            INFO, LDA, LWORK, N
00025 *       ..
00026 *       .. Array Arguments ..
00027 *       INTEGER            IPIV( * )
00028 *       COMPLEX*16         A( LDA, * ), WORK( * )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> ZGETRI computes the inverse of a matrix using the LU factorization
00038 *> computed by ZGETRF.
00039 *>
00040 *> This method inverts U and then computes inv(A) by solving the system
00041 *> inv(A)*L = inv(U) for inv(A).
00042 *> \endverbatim
00043 *
00044 *  Arguments:
00045 *  ==========
00046 *
00047 *> \param[in] N
00048 *> \verbatim
00049 *>          N is INTEGER
00050 *>          The order of the matrix A.  N >= 0.
00051 *> \endverbatim
00052 *>
00053 *> \param[in,out] A
00054 *> \verbatim
00055 *>          A is COMPLEX*16 array, dimension (LDA,N)
00056 *>          On entry, the factors L and U from the factorization
00057 *>          A = P*L*U as computed by ZGETRF.
00058 *>          On exit, if INFO = 0, the inverse of the original matrix A.
00059 *> \endverbatim
00060 *>
00061 *> \param[in] LDA
00062 *> \verbatim
00063 *>          LDA is INTEGER
00064 *>          The leading dimension of the array A.  LDA >= max(1,N).
00065 *> \endverbatim
00066 *>
00067 *> \param[in] IPIV
00068 *> \verbatim
00069 *>          IPIV is INTEGER array, dimension (N)
00070 *>          The pivot indices from ZGETRF; for 1<=i<=N, row i of the
00071 *>          matrix was interchanged with row IPIV(i).
00072 *> \endverbatim
00073 *>
00074 *> \param[out] WORK
00075 *> \verbatim
00076 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
00077 *>          On exit, if INFO=0, then WORK(1) returns the optimal LWORK.
00078 *> \endverbatim
00079 *>
00080 *> \param[in] LWORK
00081 *> \verbatim
00082 *>          LWORK is INTEGER
00083 *>          The dimension of the array WORK.  LWORK >= max(1,N).
00084 *>          For optimal performance LWORK >= N*NB, where NB is
00085 *>          the optimal blocksize returned by ILAENV.
00086 *>
00087 *>          If LWORK = -1, then a workspace query is assumed; the routine
00088 *>          only calculates the optimal size of the WORK array, returns
00089 *>          this value as the first entry of the WORK array, and no error
00090 *>          message related to LWORK is issued by XERBLA.
00091 *> \endverbatim
00092 *>
00093 *> \param[out] INFO
00094 *> \verbatim
00095 *>          INFO is INTEGER
00096 *>          = 0:  successful exit
00097 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00098 *>          > 0:  if INFO = i, U(i,i) is exactly zero; the matrix is
00099 *>                singular and its inverse could not be computed.
00100 *> \endverbatim
00101 *
00102 *  Authors:
00103 *  ========
00104 *
00105 *> \author Univ. of Tennessee 
00106 *> \author Univ. of California Berkeley 
00107 *> \author Univ. of Colorado Denver 
00108 *> \author NAG Ltd. 
00109 *
00110 *> \date November 2011
00111 *
00112 *> \ingroup complex16GEcomputational
00113 *
00114 *  =====================================================================
00115       SUBROUTINE ZGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
00116 *
00117 *  -- LAPACK computational routine (version 3.4.0) --
00118 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00119 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00120 *     November 2011
00121 *
00122 *     .. Scalar Arguments ..
00123       INTEGER            INFO, LDA, LWORK, N
00124 *     ..
00125 *     .. Array Arguments ..
00126       INTEGER            IPIV( * )
00127       COMPLEX*16         A( LDA, * ), WORK( * )
00128 *     ..
00129 *
00130 *  =====================================================================
00131 *
00132 *     .. Parameters ..
00133       COMPLEX*16         ZERO, ONE
00134       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
00135      $                   ONE = ( 1.0D+0, 0.0D+0 ) )
00136 *     ..
00137 *     .. Local Scalars ..
00138       LOGICAL            LQUERY
00139       INTEGER            I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB,
00140      $                   NBMIN, NN
00141 *     ..
00142 *     .. External Functions ..
00143       INTEGER            ILAENV
00144       EXTERNAL           ILAENV
00145 *     ..
00146 *     .. External Subroutines ..
00147       EXTERNAL           XERBLA, ZGEMM, ZGEMV, ZSWAP, ZTRSM, ZTRTRI
00148 *     ..
00149 *     .. Intrinsic Functions ..
00150       INTRINSIC          MAX, MIN
00151 *     ..
00152 *     .. Executable Statements ..
00153 *
00154 *     Test the input parameters.
00155 *
00156       INFO = 0
00157       NB = ILAENV( 1, 'ZGETRI', ' ', N, -1, -1, -1 )
00158       LWKOPT = N*NB
00159       WORK( 1 ) = LWKOPT
00160       LQUERY = ( LWORK.EQ.-1 )
00161       IF( N.LT.0 ) THEN
00162          INFO = -1
00163       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00164          INFO = -3
00165       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
00166          INFO = -6
00167       END IF
00168       IF( INFO.NE.0 ) THEN
00169          CALL XERBLA( 'ZGETRI', -INFO )
00170          RETURN
00171       ELSE IF( LQUERY ) THEN
00172          RETURN
00173       END IF
00174 *
00175 *     Quick return if possible
00176 *
00177       IF( N.EQ.0 )
00178      $   RETURN
00179 *
00180 *     Form inv(U).  If INFO > 0 from ZTRTRI, then U is singular,
00181 *     and the inverse is not computed.
00182 *
00183       CALL ZTRTRI( 'Upper', 'Non-unit', N, A, LDA, INFO )
00184       IF( INFO.GT.0 )
00185      $   RETURN
00186 *
00187       NBMIN = 2
00188       LDWORK = N
00189       IF( NB.GT.1 .AND. NB.LT.N ) THEN
00190          IWS = MAX( LDWORK*NB, 1 )
00191          IF( LWORK.LT.IWS ) THEN
00192             NB = LWORK / LDWORK
00193             NBMIN = MAX( 2, ILAENV( 2, 'ZGETRI', ' ', N, -1, -1, -1 ) )
00194          END IF
00195       ELSE
00196          IWS = N
00197       END IF
00198 *
00199 *     Solve the equation inv(A)*L = inv(U) for inv(A).
00200 *
00201       IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN
00202 *
00203 *        Use unblocked code.
00204 *
00205          DO 20 J = N, 1, -1
00206 *
00207 *           Copy current column of L to WORK and replace with zeros.
00208 *
00209             DO 10 I = J + 1, N
00210                WORK( I ) = A( I, J )
00211                A( I, J ) = ZERO
00212    10       CONTINUE
00213 *
00214 *           Compute current column of inv(A).
00215 *
00216             IF( J.LT.N )
00217      $         CALL ZGEMV( 'No transpose', N, N-J, -ONE, A( 1, J+1 ),
00218      $                     LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 )
00219    20    CONTINUE
00220       ELSE
00221 *
00222 *        Use blocked code.
00223 *
00224          NN = ( ( N-1 ) / NB )*NB + 1
00225          DO 50 J = NN, 1, -NB
00226             JB = MIN( NB, N-J+1 )
00227 *
00228 *           Copy current block column of L to WORK and replace with
00229 *           zeros.
00230 *
00231             DO 40 JJ = J, J + JB - 1
00232                DO 30 I = JJ + 1, N
00233                   WORK( I+( JJ-J )*LDWORK ) = A( I, JJ )
00234                   A( I, JJ ) = ZERO
00235    30          CONTINUE
00236    40       CONTINUE
00237 *
00238 *           Compute current block column of inv(A).
00239 *
00240             IF( J+JB.LE.N )
00241      $         CALL ZGEMM( 'No transpose', 'No transpose', N, JB,
00242      $                     N-J-JB+1, -ONE, A( 1, J+JB ), LDA,
00243      $                     WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA )
00244             CALL ZTRSM( 'Right', 'Lower', 'No transpose', 'Unit', N, JB,
00245      $                  ONE, WORK( J ), LDWORK, A( 1, J ), LDA )
00246    50    CONTINUE
00247       END IF
00248 *
00249 *     Apply column interchanges.
00250 *
00251       DO 60 J = N - 1, 1, -1
00252          JP = IPIV( J )
00253          IF( JP.NE.J )
00254      $      CALL ZSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 )
00255    60 CONTINUE
00256 *
00257       WORK( 1 ) = IWS
00258       RETURN
00259 *
00260 *     End of ZGETRI
00261 *
00262       END
 All Files Functions